#!/usr/bin/env python from numpy import * from numpy.linalg import lstsq import pyfits from pylab import * from pyraf import iraf # TRANSFORMATION PIPELINE (TP2) FOR q957g FRAMES fin = open('q957g_BDR','r') # basic data release (input) imlist = 'imlist' # list of frames (input) fimlist = open(imlist,'r') wfile = 'q957g_FDR2' # final data release (output) N = 200 # number of frames (total number of available frames = 200) Nstar = 7 # number of stars (total number of available stars = 19): X, G, F, H, E, D, R, ... ind = zeros(N) # frame number [0] (in BDR) date = zeros(N,float) # frame label [1] t = zeros(N,float) # MJD - 50000 (days) [2] fwhm = zeros(N,float) # FWHM (") of H star (from starfind)[3] ellip = zeros(N,float) # ellipticity of H star [4] snr = zeros(N,float) # signal-to-noise ratio (A) [5] qA = zeros(N,float) # flux of A (counts) [6] qB = zeros(N,float) # flux of B (counts) [7] chi2 = zeros(N,float) # chi2 (fit to A+B subframe) [8] stars = zeros((Nstar,N),float) # fluxes of stars (counts) [9-15] xs = zeros((Nstar,N),float) # x coordinates of stars (origin in the bottom left corner of the FOV) ys = zeros((Nstar,N),float) # y coordinates of stars xqA = zeros(N,float) # x-coordinate of A yqA = zeros(N,float) # y-coordinate of A xqB = zeros(N,float) # x-coordinate of B yqB = zeros(N,float) # y-coordinate of B exptime = zeros(N,float) # exposure time (to normalize the instrumental fluxes -> counts/sec) for i in range(N): filename = fimlist.readline() filename = filename[:-1] hdr = pyfits.getheader(filename+'.fits') exptime[i] = hdr['exptime'] fcoo1 = open(filename+'.coo.1') # IRAF initial coordinates of stars (input) for j in range(Nstar): s = fcoo1.readline() words = s.split() xs[j,i] = float(words[0]) ys[j,i] = float(words[1]) fcoo1.close() fcoo2 = open(filename+'.coo.2') # IRAF initial coordinates of quasar images (input) s = fcoo2.readline() words = s.split() xqA[i] = float(words[0]) yqA[i] = float(words[1]) s = fcoo2.readline() words = s.split() xqB[i] = float(words[0]) yqB[i] = float(words[1]) fcoo2.close() s = fin.readline() words = s.split() ind[i] = int(words[0]) date[i] = int(words[1]) t[i] = float(words[2]) fwhm[i] = float(words[3]) ellip[i] = float(words[4]) snr[i] = float(words[5]) qA[i] = float(words[6]) qB[i] = float(words[7]) chi2[i] = float(words[8]) for j in range(Nstar): stars[j,i] = float(words[9+j]) fimlist.close() fin.close() # SDSS stellar magnitudes (input). Note: the r-magnitude of D was derived from a r = r(VR) relationship g_cat = zeros(Nstar,float) r_cat = zeros(Nstar,float) fcat = open('SDSS19.dat','r') for j in range(Nstar): s = fcat.readline() words = s.split() g_cat[j] = float(words[4]) r_cat[j] = float(words[6]) fcat.close() # Least squares fit (PHOTOMETRIC WEIGHTS): # deviation in magnitude (delta = inst - cat) vs. zero point+colour+inhomogeneity model xs[:,:] = xs - 512 # origin of stellar coordinates in the centre of the FOV ys[:,:] = ys - 512 isel = zeros(N) isel[:] = 0 Ngood = 0 Nsima = 0 for i in range(N): if (ind[i] != 7 and snr[i] >= 100): # selection of good frames: Ngood < N isel[i] = 1 Ngood = Ngood + 1 for j in range(Nstar): if (stars[j,i] != 0): # stellar images in good frames: Nsima < Nstar * Ngood Nsima = Nsima + 1 print Ngood,Nsima delta = zeros(Nsima, float) A = zeros((Nsima,2*Ngood+5), float) k = 0 l = -1 for i in range(N): if (isel[i] == 1): l = l + 1 for j in range(Nstar): if (stars[j,i] != 0): delta[k] = (-2.5*log10(stars[j,i]/exptime[i])- g_cat[j])*sqrt(stars[j,i]) # deviation A[k,l] = sqrt(stars[j,i]) # zero point term A[k,Ngood+l] = (g_cat[j] - r_cat[j])*sqrt(stars[j,i]) # colour term A[k,2*Ngood] = xs[j,i]*sqrt(stars[j,i]) # linear inhomogeneity terms A[k,2*Ngood+1] = ys[j,i]*sqrt(stars[j,i]) A[k,2*Ngood+2] = xs[j,i]**2*sqrt(stars[j,i]) # quadratic inhomogeneity terms A[k,2*Ngood+3] = xs[j,i]*ys[j,i]*sqrt(stars[j,i]) A[k,2*Ngood+4] = ys[j,i]**2*sqrt(stars[j,i]) k = k + 1 (p, residuals, rank, s) = lstsq(A,delta) # p are the parameters that minimize the sum of (deltak - modelk)^2 (= residuals) # Output of the optimization block a0 = p[0:Ngood] C = p[Ngood:2*Ngood] I10 = p[2*Ngood] I01 = p[2*Ngood+1] I20 = p[2*Ngood+2] I11 = p[2*Ngood+3] I02 = p[2*Ngood+4] resdof = residuals/(Nsima-2*Ngood-5) print a0 print C print I10,I01,I20,I11,I02 print resdof # Average colour coefficient (checking the results: is expected to be around -0.09) colour = 0 for k in range(Ngood): colour = colour + C[k] print colour/Ngood #2D inhomogeneity pattern iraf.imdel('surface2') Nsize = 1024 im = zeros((Nsize,Nsize), float) for i in range(Nsize): y = float(i+1-Nsize/2) for j in range(Nsize): x = float(j+1-Nsize/2) im[i,j] = I10*x + I01*y + I20*x**2 + I11*x*y + I02*y**2 pyfits.writeto('surface2.fits',im) # Final stellar magnitudes mag = zeros((Nstar,Ngood), float) mag[:,:] = 0. k = -1 for i in range(N): if (isel[i] == 1): k = k + 1 for j in range(Nstar): if (stars[j,i] != 0): mag[j,k] = -2.5*log10(stars[j,i]/exptime[i]) - a0[k] - C[k]*(g_cat[j] - r_cat[j]) - I10*xs[j,i] - I01*ys[j,i] - I20*xs[j,i]**2 - I11*xs[j,i]*ys[j,i] - I02*ys[j,i]**2 for j in range(Nstar): mask = where(mag[j,:] > 0, 1, 0) magok = compress(mask, mag[j,:]) scatter = magok.std() print j,scatter # Final quasar magnitudes xqA[:] = xqA - 512 yqA[:] = yqA - 512 xqB[:] = xqB - 512 yqB[:] = yqB - 512 night = zeros(Ngood,float) time = zeros(Ngood,float) seeing = zeros(Ngood,float) sym = zeros(Ngood,float) signal = zeros(Ngood,float) chi2sys = zeros(Ngood,float) magqA = zeros(Ngood,float) magqB = zeros(Ngood,float) k = -1 for i in range(N): if (isel[i] == 1): k = k + 1 night[k] = date[i] time[k] = t[i] seeing[k] = fwhm[i] sym[k] = ellip[i] signal[k] = snr[i] chi2sys[k] = chi2[i] magqA[k] = -2.5*log10(qA[i]/exptime[i]) - a0[k] - C[k]*0.229 - I10*xqA[i] - I01*yqA[i] - I20*xqA[i]**2 - I11*xqA[i]*yqA[i] - I02*yqA[i]**2 magqB[k] = -2.5*log10(qB[i]/exptime[i]) - a0[k] - C[k]*0.165 - I10*xqB[i] - I01*yqB[i] - I20*xqB[i]**2 - I11*xqB[i]*yqB[i] - I02*yqB[i]**2 # Making the q957g_FDR2 file: (1) frame label, (2) MJD - 50000 (days), (3) FWHM ("), (4) ellipticity, (5) SNR, (6) chi2 (imfitfits), (7) zero point term, (8) colour coefficient, (9-15) magnitudes of X, G, F, H, E, D, R, and (16-17) magnitudes of A and B (quasar images). This file also includes the five inhomogeneity coefficients fw = open(wfile,'w') for k in range(0,Ngood): print >> fw, str(int(night[k])), str(time[k])[0:9], str(seeing[k])[0:4], str(sym[k])[0:5], str(int(signal[k])), str(chi2sys[k])[0:4], str(a0[k])[0:7], str(C[k])[0:7], for j in range(Nstar): print >> fw, str(mag[j,k])[0:6], print >> fw, str(magqA[k])[0:6], str(magqB[k])[0:6] print >> fw, I10, I01, I20, I11, I02 fw.close()