/
m_mkaimg.pro
167 lines (143 loc) · 4.78 KB
/
m_mkaimg.pro
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
165
166
;+
; NAME:
; x_mkaimg
; Version 1.1
;
; PURPOSE:
; Given the 2D solution for the slope of the lines as a function of
; position this code creates a wavelength image (i.e. assigns a
; unique wavelength to each pixel in each order). The xoffset is
; input in order to properly determine the edges of each order. A
; simple spline interpolation is used to determine the values.
;
; Note, you should have shifted the order structure as necessary
; prior to calling this routine.
;
; CALLING SEQUENCE:
; x_mkaimg, arc_fil, ordr_str, arc2d_fil, fil_fittrc, $
; out_fil, /CHK, /CLOBBER, BAD_ANLY=
;
; INPUTS:
; arc_fil -- Name of arc file
; ordr_str -- Order strucure describing the echelle footprint
; arc2d_fil -- Name of 2D wavelength solution
; fil_fittrc -- Name of FITS file containing the 2D fit to the
; tilted arc lines
;
; RETURNS:
;
; OUTPUTS:
; 2D wavelength image with name like 'Arcs/Arc_mb0439I.fits'
;
; OPTIONAL KEYWORDS:
; /CLOBBER - Overwrite previous image
; /CHK - Display the final image
;
; OPTIONAL OUTPUTS:
; BAD_ANLY= - Set to 1 if the code finds double valued solutions.
;
; COMMENTS:
;
; EXAMPLES:
;
;
; PROCEDURES/FUNCTIONS CALLED:
;
; REVISION HISTORY:
; 15-May-2003 Written by SB
;-
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function m_mkaimg, arc_fil, ordr_str, arc2d_fil, fil_fittrc, $
out_fil, CHK=chk, CLOBBER=clobber, $
BAD_ANLY=bad_anly
;
if N_params() LT 5 then begin
print,'Syntax - ' + $
'rslt = x_mkaimg( arc_fil, ordr_str, arc2d_fil, fil_fittrc, ' + $
'out_fil, /CHK, ' + $
'/CLOBBER, SHFTPRM= ) [v1.1]'
return, -1
endif
;; Optional keywords
bad_anly = 0
;; Order structure
nordr = n_elements(ordr_str)
;; Grab Arc 2D Fit
arc_2dfit = xmrdfits(arc2d_fil,1,/silent)
;; Check arc_fil for size
head = xheadfits(arc_fil)
sz = lonarr(2)
sz[0] = sxpar(head, 'NAXIS1')
sz[1] = sxpar(head, 'NAXIS2')
ximage = lindgen(sz[0]) # replicate(1,sz[1])
yimage = lindgen(sz[1]) ## replicate(1,sz[0])
;; Grab Slope fit
mkaimg_str = xmrdfits(fil_fittrc,1,/silent)
;; Create Final image
aimg = dblarr(sz[0],sz[1])
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; LOOPING
;; Create the ordermask image
ordermask = x_ordermask(sz[0], sz[1], ordr_str, trim=0.1)
print, 'x_mkaimg: Looping on Orders'
bad_anly = 0
for q=0L,nordr-1 do begin
print, 'x_mkaimg: Order', ordr_str[q].order
;; Set order
mkaimg_ordr = ordr_str[q].order
;; xcen
xcen = (ordr_str[q].lhedg + ordr_str[q].rhedg)/2.
inorder = where(ordermask EQ ordr_str[q].order, nin)
if nin GT 300000. then stop
ystart = 1.0d*yimage[inorder]
xstart = (ximage - xcen ## replicate(1,sz[0]))[inorder]
ycol = dindgen(sz[1])
slope_spline = spl_init(ycol, double(ordr_str[q].arc_m))
;
; Four iterations is plenty, converges quickly after the second call
;
ywave = ystart
for islope=1,4 do begin
slope = spl_interp(ycol, ordr_str[q].arc_m, slope_spline, ywave)
ywave = ystart - slope*xstart
oldslope = slope
endfor
p = 2.0d * (ywave - arc_2dfit.nrm[0])/arc_2dfit.nrm[1]
t = (2.0d *(ordr_str[q].order - arc_2dfit.nrmt[0]) $
/arc_2dfit.nrmt[1])
worky = flegendre(p, arc_2dfit.ny)
workt = flegendre(t[0], arc_2dfit.no)
; aimg[inorder] = $
; worky # transpose(reform(arc_2dfit.res, arc_2dfit.no, $
; arc_2dfit.ny) ## workt)
aimg[inorder] = alog10($
worky # transpose(reform(arc_2dfit.res, arc_2dfit.no, $
arc_2dfit.ny) ## workt) / ordr_str[q].order)
ys = sort(ywave)
nss = n_elements(ys)
multi_forw = total( aimg[inorder[ys[1:*]]] - aimg[inorder[ys[0:nss-2]]] GT 0)
multi_back = total( aimg[inorder[ys[1:*]]] - aimg[inorder[ys[0:nss-2]]] LT 0)
if keyword_set(chk) then print, 'x_mkaimg: mono? ', multi_forw, multi_back
if (multi_forw GT 0) AND (multi_back GT 0) then begin
print, 'x_mkaimg: This order is not single-valued in wavelength.'
print, 'The 2D wavelength fit is likely bad!'
print, 'Setting flg_anly to 0'
stop
bad_anly = 1
endif
endfor
if keyword_set( CHK ) then begin
xatv, aimg, /block, min=3.5, max=4.0
endif
if bad_anly GT 0 then begin
print, 'x_mkaimg: Bad wavelength solution, not writing to disk'
return, out_fil
endif
;; Output
print, 'x_mkaimg_work: Writing: ', out_fil
mwrfits, aimg, out_fil, /create
print, 'x_mkaimg_work: Compressing..'
spawn, 'gzip -f '+out_fil
print, 'x_mkaimg: All done..'
return, out_fil
end