Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/task_polfromgain.py: 93%
178 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-01 07:19 +0000
1# enable local tools
2import os
3from math import pi,floor,atan2,sin,cos,sqrt
4import pylab as mypl
6import casatools
7from casatasks import casalog
9mytb=casatools.table()
10myme=casatools.measures()
11mymd=casatools.msmetadata()
13def polfromgain(vis,tablein,caltable,paoffset,minpacov):
15 casalog.origin('polfromgain')
17 casalog.post("Deriving calibrator linear polarization from gain ratios.")
19 casalog.post("Requiring at least "+str(minpacov)+" deg of parallactic angle coverage for each antenna solution.")
20 minpacovR=minpacov*pi/180.0
22 try:
24 if ((type(vis)==str) & (os.path.exists(vis))):
25 mymd.open(vis)
26 else:
27 raise ValueError('Visibility data set not found - please verify the name')
29 nant=mymd.nantennas()
30 nfld=mymd.nfields()
31 fldnames=mymd.fieldnames()
32 nspw=mymd.nspw()
34 rempol=False
35 if ((type(tablein)==str) & (os.path.exists(tablein))):
36 if type(caltable)==str and len(caltable)>0:
38 if os.path.exists(caltable):
39 raise ValueError('Output caltable='+caltable+' exists. Choose another name or delete it.')
41 casalog.post("New caltable, "+caltable+", corrected for linear polarization, will be generated.")
42 mytb.open(tablein)
43 myout=mytb.copy(newtablename=caltable,deep=True)
44 mytb.close()
45 myout.close()
46 rempol=True
47 else:
48 casalog.post("No new caltable will be generated")
49 caltable=tablein
50 else:
51 raise ValueError('input calibration table not found - please verify the name')
54 if paoffset!=0.0:
55 casalog.post("NB: default band position angle will be offset by "+str(paoffset)+"deg.")
57 # Field coords
58 mytb.open(caltable+'/FIELD')
59 dirs=mytb.getcol('DELAY_DIR')[:,0,:]
60 mytb.close()
62 # Must retrieve nominal feed angles from MS.FEED!
63 mytb.open(vis+'/FEED')
64 nfeed=mytb.nrows()
65 fang=mytb.getcol('RECEPTOR_ANGLE')
66 fspw=mytb.getcol('SPECTRAL_WINDOW_ID')
67 fant=mytb.getcol('ANTENNA_ID')
68 rang=mypl.zeros((nant,nspw));
69 for ifeed in range(nfeed):
70 rang[fant[ifeed],fspw[ifeed]]=fang[0,ifeed]
71 mytb.close()
73 R=mypl.zeros((nspw,nfld))
74 Q=mypl.zeros((nspw,nfld))
75 U=mypl.zeros((nspw,nfld))
76 mask=mypl.zeros((nspw,nfld),dtype=bool)
78 IQUV={}
79 nomod=not rempol
80 mytb.open(caltable,nomodify=nomod)
81 uflds=mypl.unique(mytb.getcol('FIELD_ID'))
82 uspws=mypl.unique(mytb.getcol('SPECTRAL_WINDOW_ID'))
83 for ifld in uflds:
84 rah=dirs[0,ifld]*12.0/pi
85 decr=dirs[1,ifld]
86 IQUV[fldnames[ifld]]={}
87 for ispw in uspws:
89 r=mypl.zeros(nant)
90 q=mypl.zeros(nant)
91 u=mypl.zeros(nant)
92 antok=mypl.zeros(nant,dtype=bool)
93 parangbyant={} # we will remember parang outside main iant loop (for rempol on bad ants)
95 casalog.post('Fld='+fldnames[ifld]+' Spw='+str(ispw)+':')
97 for iant in range(nant):
98 qstring='FIELD_ID=='+str(ifld)+' && SPECTRAL_WINDOW_ID=='+str(ispw)+' && ANTENNA1=='+str(iant)
99 st=mytb.query(query=qstring)
100 nrows=st.nrows()
101 if nrows > 0:
103 times=st.getcol('TIME')
104 gains=st.getcol('CPARAM')
105 flags=st.getcol('FLAG')
106 flags=mypl.logical_or(flags[0,0,:],flags[1,0,:]) # 1D
108 # parang
109 parang=mypl.zeros(len(times))
111 apos=mymd.antennaposition(iant)
112 latr=myme.measure(apos,'WGS84')['m1']['value']
113 myme.doframe(apos)
114 har=mypl.zeros(nrows)
115 for itim in range(len(times)):
116 tm=myme.epoch('UTC',str(times[itim])+'s')
117 last=myme.measure(tm,'LAST')['m0']['value']
118 last-=floor(last) # days
119 last*=24.0 # hours
120 ha=last-rah # hours
121 har[itim]=ha*2.0*pi/24.0 # radians
123 parang=mypl.arctan2( (cos(latr)*mypl.sin(har)),
124 (sin(latr)*cos(decr)-cos(latr)*sin(decr)*mypl.cos(har)) )
125 # correct for cycle at +/-180.
126 # (makes values crossing +/-180 deg all positive, and thus continuous)
127 # (hmm, this will still fail at inferior circumpolar transit, but that is very rare)
128 if (latr<decr):
129 parang[parang<0.0]+=(2*pi)
131 parang+=rang[iant,ispw]
132 parang+=(paoffset*pi/180.) # manual feed pa offset
134 # save parang values for this antenna
135 parangbyant[iant]=parang
137 # Escape if insufficient samples
138 nsamp=nrows-mypl.sum(flags)
139 if (nsamp<3):
140 antok[iant]=False
141 casalog.post(' Ant='+str(iant)+' has insuffiencient sampling: nsamp='+
142 str(nsamp)+' < 3')
143 st.close()
144 continue
146 # Check parang coverage
147 dparang=abs(parang[~flags].max()-parang[~flags].min()) # rad
148 if dparang<minpacovR:
149 antok[iant]=False
150 casalog.post(' Ant='+str(iant)+' has insuffiencient parang cov: '+
151 str(round(dparang*180/pi,2))+' < '+str(minpacov)+'deg')
152 continue
156 # indep var matrix
157 A=mypl.ones((nrows,3))
158 A[:,1]=mypl.cos(2*parang)
159 A[:,2]=mypl.sin(2*parang)
160 A[flags,:]=0.0 # zero flagged rows
162 # squared gain amplitude ratio
163 amps=mypl.absolute(gains)
164 amps[amps==0.0]=1.0
165 gratio2=mypl.square(amps[0,0,:]/amps[1,0,:])
166 gratio2[flags]=0.0 # zero flagged samples
168 fit=mypl.lstsq(A,gratio2,rcond=None)
170 r2=fit[0][0]
171 if r2<0.0:
172 casalog.post(' Ant='+str(iant)+' yielded an unphysical solution; skipping.')
173 continue
175 # Reaching here, we have a nominally good solution
176 antok[iant]=True;
178 q[iant]=fit[0][1]/r2/2.0
179 u[iant]=fit[0][2]/r2/2.0
180 p=sqrt(q[iant]**2+u[iant]**2)
181 x=0.5*atan2(u[iant],q[iant])*180/pi
183 casalog.post(' Ant='+str(iant)+
184 ' (PA offset='+str(round(rang[iant,ispw]*180/pi+paoffset,2))+'deg)'+
185 ' q='+str(round(q[iant],4))+' u='+str(round(u[iant],4))+' p='+str(round(p,4))+' x='+str(round(x,3))+
186 ' Gx/Gy='+str(round(sqrt(r2),4))+
187 ' (parang cov='+str(round(dparang*180/pi,1))+'deg; nsamp='+str(nsamp))
189 if rempol:
190 if p<1.0:
191 Qpsi=q[iant]*mypl.cos(2*parang) + u[iant]*mypl.sin(2*parang)
192 gains[0,0,:]/=mypl.sqrt(1.0+Qpsi)
193 gains[1,0,:]/=mypl.sqrt(1.0-Qpsi)
194 st.putcol('CPARAM',gains)
195 else:
196 st.close()
197 raise RuntimeError('Spurious fractional polarization!')
198 st.close()
200 nantok=mypl.sum(antok)
201 if nantok==0:
202 casalog.post('Found no good polarization solutions for Fld='+fldnames[ifld]+' Spw='+str(ispw),'WARN')
203 mask[ispw,ifld]=False
204 else:
205 Q[ispw,ifld]=mypl.sum(q)/nantok
206 U[ispw,ifld]=mypl.sum(u)/nantok
207 R[ispw,ifld]=mypl.sum(r)/nantok
208 mask[ispw,ifld]=True
210 P=sqrt(Q[ispw,ifld]**2+U[ispw,ifld]**2)
211 X=0.5*atan2(U[ispw,ifld],Q[ispw,ifld])*180/pi
213 casalog.post(' Ant=<*> '+
214 ' Q='+str(round(Q[ispw,ifld],4))+
215 ' U='+str(round(U[ispw,ifld],4))+
216 ' P='+str(round(P,4))+' X='+str(round(X,3)))
218 IQUV[fldnames[ifld]]['Spw'+str(ispw)]=[1.0,Q[ispw,ifld],U[ispw,ifld],0.0]
220 # if required, remove ant-averaged polarization from non-ok antennas (if any)
221 if rempol and nantok<nant:
222 badantlist=[i for i in range(nant) if not antok[i]]
223 casalog.post(' (Correcting undersampled antennas ('+str(badantlist)+
224 ') with <*> solution.)')
225 for iant in badantlist:
226 if iant in parangbyant.keys():
227 qstring='FIELD_ID=='+str(ifld)+' && SPECTRAL_WINDOW_ID=='+str(ispw)+' && ANTENNA1=='+str(iant)
228 st=mytb.query(query=qstring)
229 if st.nrows()>0 and P<1.0:
230 gains=st.getcol('CPARAM')
231 parang=parangbyant[iant]
232 Qpsi=Q[ispw,ifld]*mypl.cos(2*parang) + U[ispw,ifld]*mypl.sin(2*parang)
233 gains[0,0,:]/=mypl.sqrt(1.0+Qpsi)
234 gains[1,0,:]/=mypl.sqrt(1.0-Qpsi)
235 st.putcol('CPARAM',gains)
236 st.close()
239 if sum(mask[:,ifld])>0:
240 casalog.post('For field='+fldnames[ifld]+' there are '+str(sum(mask[:,ifld]))+' good spws.')
241 Qm=mypl.mean(Q[mask[:,ifld],ifld])
242 Um=mypl.mean(U[mask[:,ifld],ifld])
243 IQUV[fldnames[ifld]]['SpwAve']=[1.0,Qm,Um,0.0]
244 Qe=mypl.std(Q[mask[:,ifld],ifld])
245 Ue=mypl.std(U[mask[:,ifld],ifld])
246 Pm=sqrt(Qm**2+Um**2)
247 Xm=0.5*atan2(Um,Qm)*180/pi
248 casalog.post('Spw mean: Fld='+fldnames[ifld]+' Q='+str(round(Qm,4))+' U='+str(round(Um,4))+' P='+str(round(Pm,4))+' X='+str(round(Xm,3)))
249 else:
250 casalog.post('Found no good polarization solutions for Fld='+fldnames[ifld]+' in any spw.','WARN')
252 mytb.close()
254 casalog.post("NB: Returning dictionary containing fractional Stokes results.")
255 return IQUV
257 finally:
258 mymd.close()
259 myme.done()