Source
xxxxxxxxxx
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])))
# enable local tools
from taskinit import *
import os
from math import pi,floor,atan2,sin,cos,sqrt
import pylab as mypl
mytb=tbtool()
myme=metool()
mymd=msmdtool()
def polfromgain(vis,tablein,caltable,paoffset):
casalog.origin('polfromgain')
casalog.post("Deriving calibrator linear polarization from gain ratios.")
try:
if ((type(vis)==str) & (os.path.exists(vis))):
mymd.open(vis)
else:
raise Exception, 'Visibility data set not found - please verify the name'
nant=mymd.nantennas()
nfld=mymd.nfields()
fldnames=mymd.fieldnames()
nspw=mymd.nspw()
rempol=False
if ((type(tablein)==str) & (os.path.exists(tablein))):
if type(caltable)==str and len(caltable)>0:
if os.path.exists(caltable):
raise Exception, 'Output caltable='+caltable+' exists. Choose another name or delete it.'
casalog.post("New caltable, "+caltable+", corrected for linear polarization, will be generated.")
mytb.open(tablein)
myout=mytb.copy(newtablename=caltable,deep=True)
mytb.close()
myout.close()
rempol=True
else:
casalog.post("No new caltable will be generated")
caltable=tablein
else:
raise Exception, 'input calibration table not found - please verify the name'
if paoffset!=0.0:
casalog.post("NB: default band position angle will be offset by "+str(paoffset)+"deg.")
# Field coords
mytb.open(caltable+'/FIELD')
dirs=mytb.getcol('DELAY_DIR')[:,0,:]
mytb.close()
# Must retrieve nominal feed angles from MS.FEED!
mytb.open(vis+'/FEED')
nfeed=mytb.nrows()
fang=mytb.getcol('RECEPTOR_ANGLE')
fspw=mytb.getcol('SPECTRAL_WINDOW_ID')
fant=mytb.getcol('ANTENNA_ID')
rang=mypl.zeros((nant,nspw));
for ifeed in range(nfeed):
rang[fant[ifeed],fspw[ifeed]]=fang[0,ifeed]
mytb.close()
R=mypl.zeros((nspw,nfld))
Q=mypl.zeros((nspw,nfld))
U=mypl.zeros((nspw,nfld))
mask=mypl.zeros((nspw,nfld),dtype=bool)
IQUV={}
nomod=not rempol
mytb.open(caltable,nomodify=nomod)
uflds=mypl.unique(mytb.getcol('FIELD_ID'))
uspws=mypl.unique(mytb.getcol('SPECTRAL_WINDOW_ID'))
for ifld in uflds:
rah=dirs[0,ifld]*12.0/pi
decr=dirs[1,ifld]
IQUV[fldnames[ifld]]={}
for ispw in uspws:
r=mypl.zeros(nant)
q=mypl.zeros(nant)
u=mypl.zeros(nant)
antok=mypl.zeros(nant,dtype=bool)
for iant in range(nant):
qstring='FIELD_ID=='+str(ifld)+' && SPECTRAL_WINDOW_ID=='+str(ispw)+' && ANTENNA1=='+str(iant)
st=mytb.query(query=qstring)