forked from rkotulla/SALT_RSS_longslit_pipeline
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathwlmodel.py
executable file
·155 lines (111 loc) · 4.81 KB
/
wlmodel.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
147
148
149
150
151
152
153
154
155
#!/usr/bin/env python
import os, sys
import pyfits
import numpy
import scipy
import pysalt
datadir="/work/salt/sandbox_official/polSALT/polsalt/data/"
# based on http://www.sal.wisc.edu/PFIS/docs/rss-vis/archive/protected/pfis/3170/3170AM0010_Spectrograph_Model_Draft_2.pdf
# and https://github.com/saltastro/SALTsandbox/blob/master/polSALT/polsalt/specpolmap.py
def rssmodelwave(#grating,grang,artic,cbin,refimg,
header, img,
xbin=1, ybin=1):
# compute wavelengths from model (this can probably be done using pyraf spectrograph model)
ncols = img.shape[0]
nrows = img.shape[1]
# compute X/Y position for each pixel
y,x = numpy.indices(img.shape)
# also account for binning
y *= ybin
x *= xbin
#
#
# Load spectrograph parameters
#
spec=numpy.loadtxt(datadir+"spec.txt",usecols=(1,))
grating_rotation_home_error = spec[0]
Grat0,Home0,ArtErr,T2Con,T3Con=spec[0:5]
FCampoly=spec[5:11]
grating_names=numpy.loadtxt(datadir+"gratings.txt",dtype=str,usecols=(0,))
#grname=numpy.loadtxt(datadir+"gratings.txt",dtype=str,usecols=(0,))
grlmm,grgam0=numpy.loadtxt(datadir+"gratings.txt",usecols=(1,2),unpack=True)
#
# Load all necessary information from FITS header
#
grating_angle = header['GR-ANGLE'] # alpha_C
articulation_angle = header['CAMANG'] #GRTILT'] # A_C
grating_name = header['GRATING']
print "grating-angle:", grating_angle
print "articulation angle:", articulation_angle
print "grating name:", grating_name
# get grating data: lines per mm
#grnum = numpy.where(grname==grating)[0][0]
#lmm = grlmm[grnum]
grnum = numpy.where(grating_names==grating_name)[0][0]
grating_lines_per_mm = grlmm[grating_name == grating_names][0]
print "grating lines/mm:", grating_lines_per_mm
#alpha_r = numpy.radians(grang+Grat0)
alpha_r = numpy.radians(grating_angle+Grat0)
#beta0_r = numpy.radians(artic*(1+ArtErr)+Home0)-alpha_r
beta0_r = numpy.radians(articulation_angle*(1+ArtErr)+Home0)-alpha_r
gam0_r = numpy.radians(grgam0[grnum])
print "alpha-r:", alpha_r
print "beta_r :", beta0_r
print "gamma_r:", gam0_r
# compute reference wavelength at center of focal plane
#lam0 = 1e7*numpy.cos(gam0_r)*(numpy.sin(alpha_r) + numpy.sin(beta0_r))/lmm
lam0 = 1e7*numpy.cos(gam0_r)*(numpy.sin(alpha_r) + numpy.sin(beta0_r))/grating_lines_per_mm
print "reference wavelength:", lam0
# compute camera focal length
ww = (lam0-4000.)/1000.
fcam = numpy.polyval(FCampoly,ww)
print "camera focal length @ 4000A:", fcam,"mm"
# compute dispersion per pixel
disp = (1e7*numpy.cos(gam0_r)*numpy.cos(beta0_r)/grating_lines_per_mm) / (fcam/.015)
#disp = (1e7*numpy.cos(gam0_r)*numpy.cos(beta0_r)/lmm)/(fcam/.015)
print "dispersion:", disp," angstroems/pixel"
#
# Iteratively compute a lambda for each pixel, refine the focal length as
# fct of lambda, and recompute lambda
#
_x = (x - 3162) * 0.015 #/ 3162.
_y = (y - 2048) * 0.015 # / 2048.
print numpy.min(x), numpy.max(x)
print numpy.min(y), numpy.max(y)
print numpy.min(_x), numpy.max(_x)
print numpy.min(_y), numpy.max(_y)
alpha = numpy.ones(img.shape) * alpha_r
for iteration in range(4):
beta = _x/fcam + beta0_r
gamma = _y/fcam + gam0_r
print beta.shape, gamma.shape
# compute lambda (1e7 = angstroem/mm)
_lambda = 1e7 * numpy.cos(gamma) * (numpy.sin(beta) + numpy.sin(alpha)) / grating_lines_per_mm
L = (_lambda - 4000.) / 1000.
fcam = numpy.polyval(FCampoly,L)
print "ITER", iteration, fcam.shape
pyfits.PrimaryHDU(data=_lambda).writeto("lambda_%d.fits" % (iteration+1), clobber=True)
return _lambda
# now compute F_cam as function of lambda_0
# use polynomial fit from ZEMAX camera model (that's step from Ken's docu)
dfcam = 3.162*disp*numpy.polyval([FCampoly[x]*(5-x) for x in range(5)],ww)
#T2 = -0.25*(1e7*numpy.cos(gam0_r)*numpy.sin(beta0_r)/lmm)/(fcam/47.43)**2 + T2Con*disp*dfcam
T2 = -0.25*(1e7*numpy.cos(gam0_r)*numpy.sin(beta0_r)/grating_lines_per_mm)/(fcam/47.43)**2 + T2Con*disp*dfcam
T3 = (-1./24.)*3162.*disp/(fcam/47.43)**2 + T3Con*disp
T0 = lam0 + T2
T1 = 3162.*disp + 3*T3
return
# compute normalized X-position (range [-1,1], X=0 is center of middle chip)
X = (numpy.array(range(cols))+1-cols/2)*cbin/3162.
lam_X = T0+T1*X+T2*(2*X**2-1)+T3*(4*X**3-3*X)
return lam_X
if __name__ == "__main__":
fn = sys.argv[1]
hdulist = pyfits.open(fn)
binning = pysalt.get_binning(hdulist)
xbin, ybin = 4,4
print "binning", binning, xbin, ybin
rssmodelwave(#grating,grang,artic,cbin,refimg,
header=hdulist[0].header,
img=hdulist['SCI'].data,
xbin=xbin, ybin=ybin)