Commits

Neal Schweighart authored 3240a9cc903 
Changed the pylab import back to what it was originally
No tags
c7c21837a95
master

gcwrap/python/scripts/task_polfromgain.py

Modified
1 1 # enable local tools
2 2 import os
3 3 from math import pi,floor,atan2,sin,cos,sqrt
4 -import pylab
4 +import pylab as mypl
5 5 ​
6 6 try:
7 7 import casatools
8 8 from casatasks import casalog
9 9
10 10 mytb=casatools.table()
11 11 myme=casatools.measures()
12 12 mymd=casatools.msmetadata()
13 13 except ImportError:
14 14 from taskinit import *
62 62 mytb.open(caltable+'/FIELD')
63 63 dirs=mytb.getcol('DELAY_DIR')[:,0,:]
64 64 mytb.close()
65 65
66 66 # Must retrieve nominal feed angles from MS.FEED!
67 67 mytb.open(vis+'/FEED')
68 68 nfeed=mytb.nrows()
69 69 fang=mytb.getcol('RECEPTOR_ANGLE')
70 70 fspw=mytb.getcol('SPECTRAL_WINDOW_ID')
71 71 fant=mytb.getcol('ANTENNA_ID')
72 - rang=pylab.zeros((nant,nspw));
72 + rang=mypl.zeros((nant,nspw));
73 73 for ifeed in range(nfeed):
74 74 rang[fant[ifeed],fspw[ifeed]]=fang[0,ifeed]
75 75 mytb.close()
76 76 ​
77 - R=pylab.zeros((nspw,nfld))
78 - Q=pylab.zeros((nspw,nfld))
79 - U=pylab.zeros((nspw,nfld))
80 - mask=pylab.zeros((nspw,nfld),dtype=bool)
77 + R=mypl.zeros((nspw,nfld))
78 + Q=mypl.zeros((nspw,nfld))
79 + U=mypl.zeros((nspw,nfld))
80 + mask=mypl.zeros((nspw,nfld),dtype=bool)
81 81 ​
82 82 IQUV={}
83 83 nomod=not rempol
84 84 mytb.open(caltable,nomodify=nomod)
85 - uflds=pylab.unique(mytb.getcol('FIELD_ID'))
86 - uspws=pylab.unique(mytb.getcol('SPECTRAL_WINDOW_ID'))
85 + uflds=mypl.unique(mytb.getcol('FIELD_ID'))
86 + uspws=mypl.unique(mytb.getcol('SPECTRAL_WINDOW_ID'))
87 87 for ifld in uflds:
88 88 rah=dirs[0,ifld]*12.0/pi
89 89 decr=dirs[1,ifld]
90 90 IQUV[fldnames[ifld]]={}
91 91 for ispw in uspws:
92 92 ​
93 - r=pylab.zeros(nant)
94 - q=pylab.zeros(nant)
95 - u=pylab.zeros(nant)
96 - antok=pylab.zeros(nant,dtype=bool)
93 + r=mypl.zeros(nant)
94 + q=mypl.zeros(nant)
95 + u=mypl.zeros(nant)
96 + antok=mypl.zeros(nant,dtype=bool)
97 97 ​
98 98 for iant in range(nant):
99 99 qstring='FIELD_ID=='+str(ifld)+' && SPECTRAL_WINDOW_ID=='+str(ispw)+' && ANTENNA1=='+str(iant)
100 100 st=mytb.query(query=qstring)
101 101 nrows=st.nrows()
102 102 if nrows > 0:
103 103 ​
104 104 times=st.getcol('TIME')
105 105 gains=st.getcol('CPARAM')
106 106 flags=st.getcol('FLAG')
107 - flags=pylab.logical_or(flags[0,0,:],flags[1,0,:]) # 1D
107 + flags=mypl.logical_or(flags[0,0,:],flags[1,0,:]) # 1D
108 108 ​
109 109 # Escape if insufficient data
110 - if (nrows-pylab.sum(flags))<3:
110 + if (nrows-mypl.sum(flags))<3:
111 111 antok[iant]=False
112 112 st.close()
113 113 continue
114 114 ​
115 115 ​
116 116 # parang
117 - parang=pylab.zeros(len(times))
117 + parang=mypl.zeros(len(times))
118 118
119 119 apos=mymd.antennaposition(iant)
120 120 latr=myme.measure(apos,'WGS84')['m1']['value']
121 121 myme.doframe(apos)
122 - har=pylab.zeros(nrows)
122 + har=mypl.zeros(nrows)
123 123 for itim in range(len(times)):
124 124 tm=myme.epoch('UTC',str(times[itim])+'s')
125 125 last=myme.measure(tm,'LAST')['m0']['value']
126 126 last-=floor(last) # days
127 127 last*=24.0 # hours
128 128 ha=last-rah # hours
129 129 har[itim]=ha*2.0*pi/24.0 # radians
130 130
131 - parang=pylab.arctan2( (cos(latr)*pylab.sin(har)),
132 - (sin(latr)*cos(decr)-cos(latr)*sin(decr)*pylab.cos(har)) )
131 + parang=mypl.arctan2( (cos(latr)*mypl.sin(har)),
132 + (sin(latr)*cos(decr)-cos(latr)*sin(decr)*mypl.cos(har)) )
133 133 ​
134 134 parang+=rang[iant,ispw]
135 135 parang+=(paoffset*pi/180.) # manual feed pa offset
136 136
137 137 # indep var matrix
138 - A=pylab.ones((nrows,3))
139 - A[:,1]=pylab.cos(2*parang)
140 - A[:,2]=pylab.sin(2*parang)
138 + A=mypl.ones((nrows,3))
139 + A[:,1]=mypl.cos(2*parang)
140 + A[:,2]=mypl.sin(2*parang)
141 141 A[flags,:]=0.0 # zero flagged rows
142 142
143 143 # squared gain amplitude ratio
144 - amps=pylab.absolute(gains)
144 + amps=mypl.absolute(gains)
145 145 amps[amps==0.0]=1.0
146 - gratio2=pylab.square(amps[0,0,:]/amps[1,0,:])
146 + gratio2=mypl.square(amps[0,0,:]/amps[1,0,:])
147 147 gratio2[flags]=0.0 # zero flagged samples
148 148
149 - fit=pylab.lstsq(A,gratio2)
149 + fit=mypl.lstsq(A,gratio2)
150 150 ​
151 151 r[iant]=fit[0][0]
152 152 q[iant]=fit[0][1]/r[iant]/2.0
153 153 u[iant]=fit[0][2]/r[iant]/2.0
154 154 r[iant]=sqrt(r[iant])
155 155 p=sqrt(q[iant]**2+u[iant]**2)
156 156 x=0.5*atan2(u[iant],q[iant])*180/pi
157 157
158 158 antok[iant]=True;
159 159 ​
160 160 #print '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 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])))
162 162 ​
163 163 if rempol:
164 164 if p<1.0:
165 - Qpsi=q[iant]*pylab.cos(2*parang) + u[iant]*pylab.sin(2*parang)
166 - gains[0,0,:]/=pylab.sqrt(1.0+Qpsi)
167 - gains[1,0,:]/=pylab.sqrt(1.0-Qpsi)
165 + Qpsi=q[iant]*mypl.cos(2*parang) + u[iant]*mypl.sin(2*parang)
166 + gains[0,0,:]/=mypl.sqrt(1.0+Qpsi)
167 + gains[1,0,:]/=mypl.sqrt(1.0-Qpsi)
168 168 st.putcol('CPARAM',gains)
169 169 else:
170 170 st.close()
171 171 raise Exception('Spurious fractional polarization!')
172 172 ​
173 173 st.close()
174 174 ​
175 - nantok=pylab.sum(antok)
175 + nantok=mypl.sum(antok)
176 176 if nantok>0:
177 - Q[ispw,ifld]=pylab.sum(q)/nantok
178 - U[ispw,ifld]=pylab.sum(u)/nantok
179 - R[ispw,ifld]=pylab.sum(r)/nantok
177 + Q[ispw,ifld]=mypl.sum(q)/nantok
178 + U[ispw,ifld]=mypl.sum(u)/nantok
179 + R[ispw,ifld]=mypl.sum(r)/nantok
180 180 mask[ispw,ifld]=True
181 181
182 182 P=sqrt(Q[ispw,ifld]**2+U[ispw,ifld]**2)
183 183 X=0.5*atan2(U[ispw,ifld],Q[ispw,ifld])*180/pi
184 184 ​
185 185 #print '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 186 ​
187 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))
188 188 ​
189 189 IQUV[fldnames[ifld]]['Spw'+str(ispw)]=[1.0,Q[ispw,ifld],U[ispw,ifld],0.0]
190 190 ​
191 191 else:
192 192 mask[ispw,ifld]=False
193 193 ​
194 194 if sum(mask[:,ifld])>0:
195 195 casalog.post('For field='+fldnames[ifld]+' there are '+str(sum(mask[:,ifld]))+' good spws.')
196 - Qm=pylab.mean(Q[mask[:,ifld],ifld])
197 - Um=pylab.mean(U[mask[:,ifld],ifld])
196 + Qm=mypl.mean(Q[mask[:,ifld],ifld])
197 + Um=mypl.mean(U[mask[:,ifld],ifld])
198 198 IQUV[fldnames[ifld]]['SpwAve']=[1.0,Qm,Um,0.0]
199 - Qe=pylab.std(Q[mask[:,ifld],ifld])
200 - Ue=pylab.std(U[mask[:,ifld],ifld])
199 + Qe=mypl.std(Q[mask[:,ifld],ifld])
200 + Ue=mypl.std(U[mask[:,ifld],ifld])
201 201 Pm=sqrt(Qm**2+Um**2)
202 202 Xm=0.5*atan2(Um,Qm)*180/pi
203 203 casalog.post('Spw mean: Fld='+fldnames[ifld]+' Q='+str(Qm)+' U='+str(Um)+' P='+str(Pm)+' X='+str(Xm))
204 204
205 205 mytb.close()
206 206 mymd.close()
207 207 ​
208 208 casalog.post("NB: Returning dictionary containing fractional Stokes results.")
209 209 return IQUV
210 210 ​

Everything looks good. We'll let you know here if there's anything you should know about.

Add shortcut