|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 #!/usr/bin/env python 00002 # 00003 # linfeedpolhelpers.py 00004 # 00005 # History: 00006 # v1.0 (gmoellen, 2012Oct24) == initial version 00007 # 00008 # This script defines several methods currently needed for 00009 # instrumental polarization calibration for the linear 00010 # feed basis. 00011 # To import these functions, type (at the CASA prompt): 00012 # 00013 # from recipes.linfeedpolhelpers import * 00014 # 00015 # 00016 import os 00017 00018 def qufromgain(gt,badspw=[]): 00019 00020 me.doframe(me.observatory('alma')) 00021 00022 tb.open(gt+'/ANTENNA') 00023 pos=tb.getcol('POSITION') 00024 meanpos=pl.mean(pos,1) 00025 frame=tb.getcolkeyword('POSITION','MEASINFO')['Ref'] 00026 units=tb.getcolkeyword('POSITION','QuantumUnits') 00027 mpos=me.position(frame, 00028 str(meanpos[0])+units[0], 00029 str(meanpos[1])+units[1], 00030 str(meanpos[2])+units[2]) 00031 me.doframe(mpos) 00032 00033 # _geodetic_ latitude 00034 latr=me.measure(mpos,'WGS84')['m1']['value'] 00035 00036 tb.open(gt+'/FIELD') 00037 nfld=tb.nrows() 00038 dirs=tb.getcol('DELAY_DIR')[:,0,:] 00039 tb.close() 00040 print 'Found as many as '+str(nfld)+' fields.' 00041 00042 tb.open(gt+'/SPECTRAL_WINDOW') 00043 nspw=tb.nrows() 00044 tb.close() 00045 print 'Found as many as '+str(nspw)+' spws.' 00046 00047 R=pl.zeros((nspw,nfld)) 00048 Q=pl.zeros((nspw,nfld)) 00049 U=pl.zeros((nspw,nfld)) 00050 mask=pl.ones((nspw,nfld),dtype=bool) 00051 00052 if (len(badspw)>0): 00053 mask[badspw,:]=False 00054 00055 tb.open(gt) 00056 for ifld in range(nfld): 00057 for ispw in range(nspw): 00058 if not mask[ispw,ifld]: 00059 continue 00060 st=tb.query('FIELD_ID=='+str(ifld)+' && SPECTRAL_WINDOW_ID=='+str(ispw)) 00061 nrows=st.nrows() 00062 if nrows > 0: 00063 rah=dirs[0,ifld]*12.0/pi 00064 decr=dirs[1,ifld] 00065 times=st.getcol('TIME') 00066 gains=st.getcol('CPARAM') 00067 ants=st.getcol('ANTENNA1') 00068 00069 ntimes=len(pl.unique(times)) 00070 nants=ants.max()+1 00071 00072 print nrows, ntimes, nants 00073 00074 # times 00075 time0=86400.0*floor(times[0]/86400.0) 00076 rtimes=times-time0 00077 00078 # amplitude ratio 00079 amps=pl.absolute(gains) 00080 amps[amps==0.0]=1.0 00081 ratio=amps[0,0,:]/amps[1,0,:] 00082 00083 00084 # parang 00085 parang=pl.zeros(len(times)) 00086 00087 for itim in range(len(times)): 00088 tm=me.epoch('UTC',str(times[itim])+'s') 00089 last=me.measure(tm,'LAST')['m0']['value'] 00090 last-=floor(last) # days 00091 last*=24.0 # hours 00092 ha=last-rah # hours 00093 har=ha*2.0*pi/24.0 00094 00095 parang[itim]=atan2( (cos(latr)*sin(har)), 00096 (sin(latr)*cos(decr)-cos(latr)*sin(decr)*cos(har)) ) 00097 00098 ratio.resize(nrows/nants,nants) 00099 parang.resize(nrows/nants,nants) 00100 parangd=parang*(180.0/pi) 00101 00102 A=pl.ones((nrows/nants,3)) 00103 A[:,1]=pl.cos(2*parang[:,0]) 00104 A[:,2]=pl.sin(2*parang[:,0]) 00105 00106 fit=pl.lstsq(A,pl.square(ratio)) 00107 00108 Rant=fit[0][0] 00109 Qant=fit[0][1]/fit[0][0]/2 00110 Uant=fit[0][2]/fit[0][0]/2 00111 Pant=pl.sqrt(Qant*Qant + Uant*Uant) 00112 Xant=0.5*pl.arctan2(Uant,Qant)*180/pi 00113 00114 print 'By antenna:' 00115 print ' R = ', Rant,pl.mean(Rant) 00116 print ' Q = ', Qant,pl.mean(Qant) 00117 print ' U = ', Uant,pl.mean(Uant) 00118 print ' P = ', Pant,pl.mean(Pant) 00119 print ' X = ', Xant,pl.mean(Xant) 00120 00121 pl.plot(Qant,Uant,',') 00122 00123 00124 ants0=range(nants) 00125 rsum=pl.sum(ratio[:,ants0],1) 00126 rsum/=len(ants0) 00127 00128 fit=pl.lstsq(A,pl.square(rsum)) 00129 00130 R[ispw,ifld]=fit[0][0] 00131 Q[ispw,ifld]=fit[0][1]/R[ispw,ifld]/2.0 00132 U[ispw,ifld]=fit[0][2]/R[ispw,ifld]/2.0 00133 P=sqrt(Q[ispw,ifld]**2+U[ispw,ifld]**2) 00134 X=0.5*atan2(U[ispw,ifld],Q[ispw,ifld])*180/pi 00135 00136 print 'Fld=',ifld,'Spw=',ispw,'Gx/Gy=',R[ispw,ifld],'Q=',Q[ispw,ifld],'U=',U[ispw,ifld],'P=',P,'X=',X 00137 00138 pl.plot(Q[ispw,ifld],U[ispw,ifld],'o') 00139 00140 00141 else: 00142 mask[ispw,ifld]=False 00143 00144 st.close() 00145 00146 print 'For field id = ',ifld,' there are ',sum(mask[:,ifld]),'good spws.' 00147 00148 Qm=pl.mean(Q[mask[:,ifld],ifld]) 00149 Um=pl.mean(U[mask[:,ifld],ifld]) 00150 Qe=pl.std(Q[mask[:,ifld],ifld]) 00151 Ue=pl.std(U[mask[:,ifld],ifld]) 00152 Pm=sqrt(Qm**2+Um**2) 00153 Xm=0.5*atan2(Um,Qm)*180/pi 00154 print 'Spw mean: Fld=', ifld,' Fractional: Q=',Qm,'U=',Um,'(rms=',Qe,Ue,')','P=',Pm,'X=',Xm 00155 tb.close() 00156 00157 pl.plot(Qm,Um,'*') 00158 00159 return pl.array([Qm,Um]) 00160 00161 00162 def xyamb(xy,qu,xyout=''): 00163 if xyout=='': 00164 xyout=xy 00165 if xyout!=xy: 00166 os.system('cp -r '+xy+' '+xyout) 00167 00168 tb.open(xyout,nomodify=F) 00169 00170 QU=tb.getkeyword('QU')['QU'] 00171 P=pl.sqrt(QU[0,:]**2+QU[1,:]**2) 00172 00173 nspw=P.shape[0] 00174 for ispw in range(nspw): 00175 q=QU[0,ispw] 00176 u=QU[1,ispw] 00177 if ( (abs(q)>0.0 and abs(qu[0])>0.0 and (q/qu[0])<0.0) or 00178 (abs(u)>0.0 and abs(qu[1])>0.0 and (u/qu[1])<0.0) ): 00179 print 'Fixing ambiguity in spw',ispw 00180 st=tb.query('SPECTRAL_WINDOW_ID=='+str(ispw)) 00181 c=st.getcol('CPARAM') 00182 c[0,:,:]*=-1.0 00183 st.putcol('CPARAM',c) 00184 st.close() 00185 QU[:,ispw]*=-1 00186 QUr={} 00187 QUr['QU']=QU 00188 tb.putkeyword('QU',QUr) 00189 tb.close() 00190 QUm=pl.mean(QU[:,P>0],1) 00191 QUe=pl.std(QU[:,P>0],1) 00192 print 'mean ambiguity-resolved fractional QU:',QUm,' +/- ',QUe 00193 00194 return pl.array([1.0,QUm[0],QUm[1],0.0]) 00195 00196 00197 00198 def dxy(dtab,xtab,dout): 00199 00200 os.system('cp -rf '+dtab+' '+dout) 00201 00202 # How many spws 00203 tb.open(dtab+'/SPECTRAL_WINDOW') 00204 nspw=tb.nrows() 00205 tb.close() 00206 00207 00208 for ispw in range(nspw): 00209 print 'Spw = ',ispw 00210 tb.open(xtab) 00211 x=[] 00212 st=tb.query('SPECTRAL_WINDOW_ID=='+str(ispw)) 00213 x=st.getcol('CPARAM') 00214 st.close() 00215 tb.close() 00216 #print ' x.shape = ',x.shape, len(x) 00217 00218 if (len(x)>0): 00219 00220 tb.open(dout,nomodify=F) 00221 st=tb.query('SPECTRAL_WINDOW_ID=='+str(ispw)) 00222 d=st.getcol('CPARAM') 00223 #print ' d.shape = ',d.shape, len(d), d.shape==x.shape 00224 00225 00226 # the following assumes all antennas and chans same in both tables. 00227 00228 if (len(d)>0): 00229 if (d.shape[1:3]==x.shape[1:3]): 00230 # Xinv.D.X: 00231 d[0,:,:]*=pl.conj(x[0,:,:]) 00232 d[1,:,:]*=x[0,:,:] 00233 st.putcol('CPARAM',d) 00234 else: 00235 print ' D and X shapes do not match!' 00236 else: 00237 print ' No D solutions for this spw' 00238 st.close() 00239 tb.close() 00240 else: 00241 print ' No X solutions for this spw' 00242 00243 00244 00245 def zeromeanD(dtab,dout): 00246 00247 os.system('cp -rf '+dtab+' '+dout) 00248 00249 tb.open(dout+'/SPECTRAL_WINDOW') 00250 # How many spws 00251 nspw=tb.nrows() 00252 tb.close() 00253 00254 tb.open(dout,nomodify=F) 00255 for ispw in range(nspw): 00256 print 'Spw = ',ispw 00257 00258 st=tb.query('SPECTRAL_WINDOW_ID=='+str(ispw)) 00259 00260 if (st.nrows()>0): 00261 d=st.getcol('CPARAM') 00262 fl=st.getcol('FLAG') 00263 wt=pl.ones(d.shape) 00264 wt[fl]=0.0 00265 drm=pl.average(d[0,:,:],1,wt[0,:,:]) 00266 dlm=pl.average(d[1,:,:],1,wt[1,:,:]) 00267 a=(drm-pl.conj(dlm))/2.0 00268 a=a.reshape((len(a),1)) 00269 d[0,:,:]=d[0,:,:]-a 00270 d[1,:,:]=d[1,:,:]+pl.conj(a) 00271 st.putcol('CPARAM',d) 00272 else: 00273 print ' No D in this spw' 00274 st.close() 00275 tb.close() 00276 00277 00278 00279 def nonorthogD(dtab,dout): 00280 00281 os.system('cp -rf '+dtab+' '+dout) 00282 00283 tb.open(dout+'/SPECTRAL_WINDOW') 00284 # How many spws 00285 nspw=tb.nrows() 00286 tb.close() 00287 00288 tb.open(dout,nomodify=F) 00289 for ispw in range(nspw): 00290 print 'Spw = ',ispw 00291 00292 st=tb.query('SPECTRAL_WINDOW_ID=='+str(ispw)) 00293 00294 if (st.nrows()>0): 00295 d=st.getcol('CPARAM') 00296 dorth=(d[0,:,:]-pl.conj(d[1,:,:]))/2.0 00297 d[0,:,:]-=dorth 00298 d[1,:,:]+=pl.conj(dorth) 00299 st.putcol('CPARAM',d) 00300 else: 00301 print ' No D in this spw' 00302 st.close() 00303 tb.close()
1.8.0
