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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
|
import datetime
import numpy
from numpy import linalg
#from pytools import wcsutil
from makewcs import MakeWCS
from dgeo import DGEO
class TDDCorr(object):
"""
Purpose
=======
Apply time dependent distortion correction to SIP coefficients and basic
WCS keywords. Atthi stime it is applicable only to ACS/WFC data.
Ref: Jay Anderson, ACS ISR 2007-08, Variation of the Distortion
Solution of the WFC
"""
def __init__(self, owcs=None, refwcs=None):
"""
:Parameters:
`owcs`: HSTWCS object
An extension HSTWCS object to be modified
`refwcs`: HSTWCS object
A reference HSTWCS object
"""
self.wcsobj = owcs
if self.wcsobj.DOTDDCorr == 'PERFORM':
self.updateWCS()
self.wcsobj.DOTDDCorr = 'COMPLETE'
else:
pass
def updateWCS(self):
"""
- Calculates alpha and beta for ACS/WFC data.
- Writes 2 new kw to the extension header: TDDALPHA and TDDBETA
"""
if self.wcsobj.instrument == 'ACS' and self.wcsobj.detector == 'WFC':
newkw = ['TDDALPHA', 'TDDBETA']
self.set_alpha_beta()
self.updatehdr(newkeywords=newkw)
else:
pass
def set_alpha_beta(self):
# Compute the time dependent distortion skew terms
# default date of 2004.5 = 2004-7-1
datedefault = datetime.datetime(2004,7,1)
year,month,day = self.wcsobj.date_obs.split('-')
rdate = datetime.datetime(int(year),int(month),int(day))
self.TDDALPHA = 0.095+0.090*((rdate-datedefault).days/365.25)/2.5
self.TDDBETA = -0.029-0.030*((rdate-datedefault).days/365.25)/2.5
self.TDDALPHA = 0.0
self.TDDBETA = 0.0
def updatehdr(self, newkeywords=None):
for kw in newkeywords:
self.wcsobj.hdr.update(kw, self.__getattribute__(kw))
class VACorr(object):
"""
Purpose
=======
Apply velocity aberation correction to WCS keywords.
Modifies the CD matrix and CRVAL1/2
"""
def __init__(self, owcs=None, refwcs=None):
self.wcsobj = owcs
self.refwcs = refwcs
if self.wcsobj.DOVACorr == 'PERFORM':
self.updateWCS()
self.wcsobj.DOVACorr == 'COMPLETE'
else:
pass
def updateWCS(self):
if self.wcsobj.vafactor != 1:
self.wcsobj.cd = self.wcsobj.cd * self.wcsobj.vafactor
#self.wcsobj.crval[0] = self.wcsobj.ra_targ + self.wcsobj.vafactor*(self.wcsobj.crval[0] - self.wcsobj.ra_targ)
#self.wcsobj.crval[1] = self.wcsobj.dec_targ + self.wcsobj.vafactor*(self.wcsobj.crval[1] - self.wcsobj.dec_targ)
ref_backup = self.refwcs.restore()
crval1 = ref_backup['CRVAL1']
crval2 = ref_backup['CRVAL2']
self.wcsobj.crval[0] = crval1 + self.wcsobj.vafactor*diff_angles(self.wcsobj.crval[0], crval1)
self.wcsobj.crval[1] = crval2 + self.wcsobj.vafactor*diff_angles(self.wcsobj.crval[1], crval2)
kw2update={'CD1_1': self.wcsobj.cd[0,0], 'CD1_2':self.wcsobj.cd[0,1],
'CD2_1':self.wcsobj.cd[1,0], 'CD2_2':self.wcsobj.cd[1,1],
'CRVAL1':self.wcsobj.crval[0], 'CRVAL2':self.wcsobj.crval[1]}
self.updatehdr(newkeywords=kw2update)
else:
pass
def updatehdr(self, newkeywords=None):
for kw in newkeywords:
self.wcsobj.hdr.update(kw, newkeywords[kw])
class CompSIP(object):
"""
Purpose
=======
Compute SIP coefficients from idc table coefficients.
"""
def __init__(self, owcs, refwcs):
self.wcsobj = owcs
self.updateWCS()
self.DOCOMPSIP = 'COMPLETE'
def updateWCS(self):
kw2update = {}
order = self.wcsobj.idcmodel.norder
kw2update['A_ORDER'] = order
kw2update['B_ORDER'] = order
pscale = self.wcsobj.idcmodel.refpix['PSCALE']
cx = self.wcsobj.idcmodel.cx
cy = self.wcsobj.idcmodel.cy
matr = numpy.array([[cx[1,1],cx[1,0]], [cy[1,1],cy[1,0]]], dtype=numpy.float)
imatr = linalg.inv(matr)
akeys1 = numpy.zeros((order+1,order+1), dtype=numpy.float)
bkeys1 = numpy.zeros((order+1,order+1), dtype=numpy.float)
for n in range(order+1):
for m in range(order+1):
if n >= m and n>=2:
idcval = numpy.array([[cx[n,m]],[cy[n,m]]])
sipval = numpy.dot(imatr, idcval)
akeys1[m,n-m] = sipval[0]
bkeys1[m,n-m] = sipval[1]
Akey="A_%d_%d" % (m,n-m)
Bkey="B_%d_%d" % (m,n-m)
kw2update[Akey] = sipval[0,0]
kw2update[Bkey] = sipval[1,0]
self.updatehdr(newkw=kw2update)
def updatehdr(self, newkw=None):
if not newkw: return
for kw in newkw.keys():
self.wcsobj.hdr.update(kw, newkw[kw])
def diff_angles(a,b):
"""
Perform angle subtraction a-b taking into account
small-angle differences across 360degree line.
"""
diff = a - b
if diff > 180.0:
diff -= 360.0
if diff < -180.0:
diff += 360.0
return diff
|
