forked from rvweeren/lofar_facet_selfcal
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathlin2circ.py
146 lines (118 loc) · 5.21 KB
/
lin2circ.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
#!/usr/bin/env python
"""
lin2circ.py
Converts linear XX,XY,YX,YY correlations to circular RR,RL,LR,LL
assuming that the definitions follow the IAU ones, which are
described by Hamaker & Bregman (1996, A&AS, 117, 161).
In particular, I use the coordinate transformation given in Section 3,
1 / 1 +i \
C_A = ------- | |
sqrt(2) \ 1 -i /
The Hermitian conjugate of this is,
1 / 1 1 \
C+_A = ------- | |
sqrt(2) \ -i +i /
So V_RL = C_A * V_XY * C+_A
where V_XY is the visibilities in linear coordinates,
V_RL is the visibilities in circular coordinates.
This reduces to:
RR = XX - iXY + iYX + YY
RL = XX + iXY + iYX - YY
LR = XX - iXY - iYX - YY
LL = XX + iXY - iYX + YY
Version 1.0 written 15 April 2010 by George Heald
Version 1.1 (18 March 2012) adds update of POLARIZATION table
Version 1.2 (July 6, 2012) - added inverse capability, convert from circular
to linear polarization.
Version 1.3 (2014, Reinout van Weeren) - corrected the flux error, factor of 0.5 (/2) was missing. Added option to convert back from circular to linear
Since the transformation matrices are hermitian, C_A ^-1 = C+_A and C+_A^-1 = C_A.
So, we have:
V_XY = C+_A*V_RL*C_A
which is:
XX = RR + RL + LR + LL
XY = iRR - iRL + iLR - iLL
YX = -iRR - iRL + iLR + iLL
YY = RR - RL - LR + LL
"""
import optparse
import casacore.tables as pt
import numpy
def main(options):
stepsize = 200000
cI = numpy.complex(0.,1.)
inms = options.inms
if inms == '':
print('Error: you have to specify an input MS, use -h for help')
return
column = options.column
outcol = options.outcol
t = pt.table(inms, readonly=False, ack=True)
if options.back:
lincol = options.lincol
if lincol not in t.colnames():
print('Adding the output linear polarization column',lincol,'to',inms)
desc = t.getcoldesc(column)
newdesc = pt.makecoldesc(lincol, desc)
newdmi = t.getdminfo(column)
newdmi['NAME'] = 'Dysco' + lincol
t.addcols(newdesc, newdmi)
### RVW EDIT 2012
print('Reading the input column (circular)', column)
if column not in t.colnames():
print('Error: Input column does not exist')
return
for row in range(0,t.nrows(),stepsize):
print('Doing row', row, 'out of', t.nrows())
### RVW EDIT 2012 Input column with the -c switch
cirdata = t.getcol(column, startrow=row, nrow=stepsize, rowincr=1)
#cirdata = t.getcol(outcol)
#print 'SHAPE ARRAY', numpy.shape(cirdata)
print('Computing the linear polarization terms...')
lindata = numpy.transpose(numpy.array([
0.5*(cirdata[:,:,0]+cirdata[:,:,1]+cirdata[:,:,2]+cirdata[:,:,3]),
0.5*(cI*cirdata[:,:,0]-cI*cirdata[:,:,1]+cI*cirdata[:,:,2]-cI*cirdata[:,:,3]),
0.5*(-cI*cirdata[:,:,0]-cI*cirdata[:,:,1]+cI*cirdata[:,:,2]+cI*cirdata[:,:,3]),
0.5*(cirdata[:,:,0]-cirdata[:,:,1]-cirdata[:,:,2]+cirdata[:,:,3])]),
(1,2,0))
print('Finishing up...')
t.putcol(lincol, lindata, startrow=row, nrow=stepsize, rowincr=1)
t.close()
else:
if outcol not in t.colnames():
print('Adding the output column',outcol,'to',inms)
desc = t.getcoldesc(column)
newdesc = pt.makecoldesc(outcol, desc)
newdmi = t.getdminfo(column)
newdmi['NAME'] = 'Dysco' + outcol
t.addcols(newdesc, newdmi)
print('Reading the input column (linear)', column)
for row in range(0,t.nrows(),stepsize):
print('Doing row', row, 'out of', t.nrows(), row+stepsize)
data = t.getcol(column, startrow=row, nrow=stepsize, rowincr=1)
print('Computing the output circular column')
outdata = numpy.transpose(numpy.array([
0.5*(data[:,:,0]-cI*data[:,:,1]+cI*data[:,:,2]+data[:,:,3]),
0.5*(data[:,:,0]+cI*data[:,:,1]+cI*data[:,:,2]-data[:,:,3]),
0.5*(data[:,:,0]-cI*data[:,:,1]-cI*data[:,:,2]-data[:,:,3]),
0.5*(data[:,:,0]+cI*data[:,:,1]-cI*data[:,:,2]+data[:,:,3])]),
(1,2,0))
print('Finishing up...')
t.putcol(outcol, outdata, startrow=row, nrow=stepsize, rowincr=1)
t.close()
if options.poltable:
print('Updating the POLARIZATION table...')
tp = pt.table(inms+'/POLARIZATION',readonly=False,ack=True)
### RVW EDIT 2012
if options.back:
tp.putcol('CORR_TYPE',numpy.array([[9,10,11,12]],dtype=numpy.int32)) # FROM CIRC-->LIN
else:
tp.putcol('CORR_TYPE',numpy.array([[5,6,7,8]],dtype=numpy.int32)) # FROM LIN-->CIRC
opt = optparse.OptionParser()
opt.add_option('-i','--inms',help='Input MS [no default]',default='')
opt.add_option('-c','--column',help='Input column [default DATA]',default='DATA')
opt.add_option('-o','--outcol',help='Output column [default DATA_CIRC]',default='DATA_CIRC')
opt.add_option('-p','--poltable',help='Update POLARIZATION table? [default False]',default=False,action='store_true')
opt.add_option('-b','--back',help='Go back to linear polarization [default False]',default=False,action='store_true')
opt.add_option('-l','--lincol',help='Output linear polarization column, if the -b switch is used [default DATA_LIN]; we want to keep the original DATA column',default='DATA_LIN')
options, arguments = opt.parse_args()
main(options)