Commits
George Moellenbrock authored and Ville Suoranta committed 9595e18f176 Merge
10 10 | mytb=casatools.table() |
11 11 | myme=casatools.measures() |
12 12 | mymd=casatools.msmetadata() |
13 13 | except ImportError: |
14 14 | from taskinit import * |
15 15 | |
16 16 | mytb=tbtool() |
17 17 | myme=metool() |
18 18 | mymd=msmdtool() |
19 19 | |
20 - | def polfromgain(vis,tablein,caltable,paoffset): |
20 + | def polfromgain(vis,tablein,caltable,paoffset,minpacov): |
21 21 | |
22 22 | casalog.origin('polfromgain') |
23 23 | |
24 24 | casalog.post("Deriving calibrator linear polarization from gain ratios.") |
25 25 | |
26 + | casalog.post("Requiring at least "+str(minpacov)+" deg of parallactic angle coverage for each antenna solution.") |
27 + | minpacovR=minpacov*pi/180.0 |
28 + | |
26 29 | try: |
27 30 | |
28 31 | if ((type(vis)==str) & (os.path.exists(vis))): |
29 32 | mymd.open(vis) |
30 33 | else: |
31 34 | raise ValueError('Visibility data set not found - please verify the name') |
32 35 | |
33 36 | nant=mymd.nantennas() |
34 37 | nfld=mymd.nfields() |
35 38 | fldnames=mymd.fieldnames() |
87 90 | for ifld in uflds: |
88 91 | rah=dirs[0,ifld]*12.0/pi |
89 92 | decr=dirs[1,ifld] |
90 93 | IQUV[fldnames[ifld]]={} |
91 94 | for ispw in uspws: |
92 95 | |
93 96 | r=mypl.zeros(nant) |
94 97 | q=mypl.zeros(nant) |
95 98 | u=mypl.zeros(nant) |
96 99 | antok=mypl.zeros(nant,dtype=bool) |
100 + | parangbyant={} # we will remember parang outside main iant loop (for rempol on bad ants) |
101 + | |
102 + | casalog.post('Fld='+fldnames[ifld]+' Spw='+str(ispw)+':') |
97 103 | |
98 104 | for iant in range(nant): |
99 105 | qstring='FIELD_ID=='+str(ifld)+' && SPECTRAL_WINDOW_ID=='+str(ispw)+' && ANTENNA1=='+str(iant) |
100 106 | st=mytb.query(query=qstring) |
101 107 | nrows=st.nrows() |
102 108 | if nrows > 0: |
103 109 | |
104 110 | times=st.getcol('TIME') |
105 111 | gains=st.getcol('CPARAM') |
106 112 | flags=st.getcol('FLAG') |
107 113 | flags=mypl.logical_or(flags[0,0,:],flags[1,0,:]) # 1D |
108 114 | |
109 - | # Escape if insufficient data |
110 - | if (nrows-mypl.sum(flags))<3: |
111 - | antok[iant]=False |
112 - | st.close() |
113 - | continue |
114 - | |
115 - | |
116 115 | # parang |
117 116 | parang=mypl.zeros(len(times)) |
118 117 | |
119 118 | apos=mymd.antennaposition(iant) |
120 119 | latr=myme.measure(apos,'WGS84')['m1']['value'] |
121 120 | myme.doframe(apos) |
122 121 | har=mypl.zeros(nrows) |
123 122 | for itim in range(len(times)): |
124 123 | tm=myme.epoch('UTC',str(times[itim])+'s') |
125 124 | last=myme.measure(tm,'LAST')['m0']['value'] |
126 125 | last-=floor(last) # days |
127 126 | last*=24.0 # hours |
128 127 | ha=last-rah # hours |
129 128 | har[itim]=ha*2.0*pi/24.0 # radians |
130 129 | |
131 130 | parang=mypl.arctan2( (cos(latr)*mypl.sin(har)), |
132 131 | (sin(latr)*cos(decr)-cos(latr)*sin(decr)*mypl.cos(har)) ) |
132 + | # correct for cycle at +/-180. |
133 + | # (makes values crossing +/-180 deg all positive, and thus continuous) |
134 + | # (hmm, this will still fail at inferior circumpolar transit, but that is very rare) |
135 + | if (latr<decr): |
136 + | parang[parang<0.0]+=(2*pi) |
133 137 | |
134 138 | parang+=rang[iant,ispw] |
135 139 | parang+=(paoffset*pi/180.) # manual feed pa offset |
136 - | |
140 + | |
141 + | # save parang values for this antenna |
142 + | parangbyant[iant]=parang |
143 + | |
144 + | # Escape if insufficient samples |
145 + | nsamp=nrows-mypl.sum(flags) |
146 + | if (nsamp<3): |
147 + | antok[iant]=False |
148 + | casalog.post(' Ant='+str(iant)+' has insuffiencient sampling: nsamp='+ |
149 + | str(nsamp)+' < 3') |
150 + | st.close() |
151 + | continue |
152 + | |
153 + | # Check parang coverage |
154 + | dparang=abs(parang[~flags].max()-parang[~flags].min()) # rad |
155 + | if dparang<minpacovR: |
156 + | antok[iant]=False |
157 + | casalog.post(' Ant='+str(iant)+' has insuffiencient parang cov: '+ |
158 + | str(round(dparang*180/pi,2))+' < '+str(minpacov)+'deg') |
159 + | continue |
160 + | |
161 + | |
162 + | |
137 163 | # indep var matrix |
138 164 | A=mypl.ones((nrows,3)) |
139 165 | A[:,1]=mypl.cos(2*parang) |
140 166 | A[:,2]=mypl.sin(2*parang) |
141 167 | A[flags,:]=0.0 # zero flagged rows |
142 168 | |
143 169 | # squared gain amplitude ratio |
144 170 | amps=mypl.absolute(gains) |
145 171 | amps[amps==0.0]=1.0 |
146 172 | gratio2=mypl.square(amps[0,0,:]/amps[1,0,:]) |
147 173 | gratio2[flags]=0.0 # zero flagged samples |
148 174 | |
149 - | fit=mypl.lstsq(A,gratio2) |
175 + | fit=mypl.lstsq(A,gratio2,rcond=None) |
150 176 | |
151 - | r[iant]=fit[0][0] |
152 - | q[iant]=fit[0][1]/r[iant]/2.0 |
153 - | u[iant]=fit[0][2]/r[iant]/2.0 |
154 - | r[iant]=sqrt(r[iant]) |
177 + | r2=fit[0][0] |
178 + | if r2<0.0: |
179 + | casalog.post(' Ant='+str(iant)+' yielded an unphysical solution; skipping.') |
180 + | continue |
181 + | |
182 + | # Reaching here, we have a nominally good solution |
183 + | antok[iant]=True; |
184 + | |
185 + | q[iant]=fit[0][1]/r2/2.0 |
186 + | u[iant]=fit[0][2]/r2/2.0 |
155 187 | p=sqrt(q[iant]**2+u[iant]**2) |
156 188 | x=0.5*atan2(u[iant],q[iant])*180/pi |
157 - | |
158 - | antok[iant]=True; |
159 189 | |
160 - | # casalog.post('Fld='+fldnames[ifld],'Spw='+str(ispw),'Ant='+str(iant), '(PA offset='+str(rang[iant,ispw]*180/pi+paoffset)+'deg)','Gx/Gy='+str(r[iant]),'q='+str(q[iant]),'u='+str(u[iant]),'p='+str(p),'x='+str(x)) |
161 - | casalog.post('Fld='+fldnames[ifld]+' Spw='+str(ispw)+' Ant='+str(iant)+' (PA offset='+str(rang[iant,ispw]*180/pi+paoffset)+'deg)'+' q='+str(q[iant])+' u='+str(u[iant])+' p='+str(p)+' x='+str(x)+' Gx/Gy='+str(sqrt(r[iant]))) |
190 + | casalog.post(' Ant='+str(iant)+ |
191 + | ' (PA offset='+str(round(rang[iant,ispw]*180/pi+paoffset,2))+'deg)'+ |
192 + | ' q='+str(round(q[iant],4))+' u='+str(round(u[iant],4))+' p='+str(round(p,4))+' x='+str(round(x,3))+ |
193 + | ' Gx/Gy='+str(round(sqrt(r2),4))+ |
194 + | ' (parang cov='+str(round(dparang*180/pi,1))+'deg; nsamp='+str(nsamp)) |
162 195 | |
163 196 | if rempol: |
164 197 | if p<1.0: |
165 198 | Qpsi=q[iant]*mypl.cos(2*parang) + u[iant]*mypl.sin(2*parang) |
166 199 | gains[0,0,:]/=mypl.sqrt(1.0+Qpsi) |
167 200 | gains[1,0,:]/=mypl.sqrt(1.0-Qpsi) |
168 201 | st.putcol('CPARAM',gains) |
169 202 | else: |
170 203 | st.close() |
171 204 | raise RuntimeError('Spurious fractional polarization!') |
172 - | |
173 205 | st.close() |
174 206 | |
175 207 | nantok=mypl.sum(antok) |
176 - | if nantok>0: |
208 + | if nantok==0: |
209 + | casalog.post('Found no good polarization solutions for Fld='+fldnames[ifld]+' Spw='+str(ispw),'WARN') |
210 + | mask[ispw,ifld]=False |
211 + | else: |
177 212 | Q[ispw,ifld]=mypl.sum(q)/nantok |
178 213 | U[ispw,ifld]=mypl.sum(u)/nantok |
179 214 | R[ispw,ifld]=mypl.sum(r)/nantok |
180 215 | mask[ispw,ifld]=True |
181 216 | |
182 217 | P=sqrt(Q[ispw,ifld]**2+U[ispw,ifld]**2) |
183 218 | X=0.5*atan2(U[ispw,ifld],Q[ispw,ifld])*180/pi |
184 219 | |
185 - | #casalog.post... 'Fld='+fldnames[ifld],'Spw='+str(ispw),'Ant=*', '(PA offset='+str(rang[iant,ispw]*180/pi+paoffset)+'deg)','Gx/Gy='+str(R[ispw,ifld]),'Q='+str(Q[ispw,ifld]),'U='+str(U[ispw,ifld]),'P='+str(P),'X='+str(X) |
186 - | |
187 - | casalog.post('Fld='+fldnames[ifld]+' Spw='+str(ispw)+' Ant=*'+' (PA offset='+str(rang[iant,ispw]*180/pi+paoffset)+'deg)'+' Q='+str(Q[ispw,ifld])+' U='+str(U[ispw,ifld])+' P='+str(P)+' X='+str(X)) |
220 + | casalog.post(' Ant=<*> '+ |
221 + | ' Q='+str(round(Q[ispw,ifld],4))+ |
222 + | ' U='+str(round(U[ispw,ifld],4))+ |
223 + | ' P='+str(round(P,4))+' X='+str(round(X,3))) |
188 224 | |
189 225 | IQUV[fldnames[ifld]]['Spw'+str(ispw)]=[1.0,Q[ispw,ifld],U[ispw,ifld],0.0] |
190 226 | |
191 - | else: |
192 - | mask[ispw,ifld]=False |
227 + | # if required, remove ant-averaged polarization from non-ok antennas (if any) |
228 + | if rempol and nantok<nant: |
229 + | badantlist=[i for i in range(nant) if not antok[i]] |
230 + | casalog.post(' (Correcting undersampled antennas ('+str(badantlist)+ |
231 + | ') with <*> solution.)') |
232 + | for iant in badantlist: |
233 + | if iant in parangbyant.keys(): |
234 + | qstring='FIELD_ID=='+str(ifld)+' && SPECTRAL_WINDOW_ID=='+str(ispw)+' && ANTENNA1=='+str(iant) |
235 + | st=mytb.query(query=qstring) |
236 + | if st.nrows()>0 and P<1.0: |
237 + | gains=st.getcol('CPARAM') |
238 + | parang=parangbyant[iant] |
239 + | Qpsi=Q[ispw,ifld]*mypl.cos(2*parang) + U[ispw,ifld]*mypl.sin(2*parang) |
240 + | gains[0,0,:]/=mypl.sqrt(1.0+Qpsi) |
241 + | gains[1,0,:]/=mypl.sqrt(1.0-Qpsi) |
242 + | st.putcol('CPARAM',gains) |
243 + | st.close() |
244 + | |
193 245 | |
194 246 | if sum(mask[:,ifld])>0: |
195 247 | casalog.post('For field='+fldnames[ifld]+' there are '+str(sum(mask[:,ifld]))+' good spws.') |
196 248 | Qm=mypl.mean(Q[mask[:,ifld],ifld]) |
197 249 | Um=mypl.mean(U[mask[:,ifld],ifld]) |
198 250 | IQUV[fldnames[ifld]]['SpwAve']=[1.0,Qm,Um,0.0] |
199 251 | Qe=mypl.std(Q[mask[:,ifld],ifld]) |
200 252 | Ue=mypl.std(U[mask[:,ifld],ifld]) |
201 253 | Pm=sqrt(Qm**2+Um**2) |
202 254 | Xm=0.5*atan2(Um,Qm)*180/pi |
203 - | casalog.post('Spw mean: Fld='+fldnames[ifld]+' Q='+str(Qm)+' U='+str(Um)+' P='+str(Pm)+' X='+str(Xm)) |
204 - | |
255 + | 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))) |
256 + | else: |
257 + | casalog.post('Found no good polarization solutions for Fld='+fldnames[ifld]+' in any spw.','WARN') |
258 + | |
205 259 | mytb.close() |
206 - | mymd.close() |
207 260 | |
208 261 | casalog.post("NB: Returning dictionary containing fractional Stokes results.") |
209 262 | return IQUV |
210 263 | |
211 264 | finally: |
212 - | mytb.close() |
213 265 | mymd.close() |
266 + | myme.done() |
267 + |