-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathchain.pro
executable file
·361 lines (317 loc) · 11.8 KB
/
chain.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
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
pro chain,chaindir=chaindir,bin=bin,logfile=logfile
;basic dir/files
home=getenv('HOME')
if not keyword_set(chaindir) then chaindir=home+'/idl/chain'
if n_elements(bin) eq 0 then bin=1
if not keyword_set(logfile) then logfile='logfile'
;params
scat='yes' ; activate scattered light subtraction
la_minexptime=9.e10; min. exp. time (seconds) to activate la_cosmic removal
xmformat='vo' ; format for order-merged output: '', 'vo', 'three-col' or 'rana'
if bin eq 1 then begin
flatname='flat.fits'
ap1name='ap1.fits'
apname='ap.fits'
endif else begin
flatname='uflat.fits'
ap1name='uap1.fits'
apname='uap.fits'
endelse
;collect input files
files=file_search('00*fits')
;rebin unbinned data and rebin
if bin eq 0 then begin
for i=0,n_elements(files)-1 do begin
newfile=files[i]
strput,newfile,'01',0
sbin,files[i],newfile
endfor
files=file_search('01*fits')
endif
;make data inventory
inventory,files,st,/bin
if n_elements(st) eq 0 then begin
if keyword_set(bin) then print,'% CHAIN: no binned (2x8) data' else $
print,'% CHAIN: no unbinned (1x1) data'
return
endif
openw,10,logfile
print,'make inventory of fits files ...'
print,' filename object type RA[deg] DEC[deg] t[s]'
print,'-------------------------------------------------------------------------------'
printf,10,'make inventory of fits files ...'
printf,10,' filename object type RA[deg] DEC[deg] t[s]'
printf,10,'-------------------------------------------------------------------------------'
for i=0,n_elements(st.obstype)-1 do begin
sfilename=strmid(st[i].filename,strpos(st[i].filename,'/',/reverse_search)+1,1000)
printf,10,sfilename,st[i].object,strmid(st[i].obstype,0,3),$
st[i].ra,st[i].dec,st[i].exptime,$
format='(a42,x,a10,x,a3,x,f8.4,x,f7.3,x,i4)'
print,sfilename,st[i].object,strmid(st[i].obstype,0,3),$
st[i].ra,st[i].dec,st[i].exptime,$
format='(a42,x,a10,x,a3,x,f8.4,x,f7.3,x,i4)'
endfor
print,'-------------------------------------------------------------------------------'
printf,10,'-------------------------------------------------------------------------------'
wcal=where(strmid(st.obstype,0,3) eq 'Cal')
wspe=where(strmid(st.obstype,0,3) eq 'Spe')
wfla=where(strmid(st.obstype,0,3) eq 'Fla')
wdar=where(strmid(st.obstype,0,3) eq 'Dar')
wob=[wcal,wspe]
;have got enough data?
if max(wcal) lt 0 then begin
print,'% CHAIN: no Cal images available, I quit'
return
endif
if max(wfla) lt 0 then begin
print,'% CHAIN: no Fla images available, I quit'
return
endif
if max(wspe) lt 0 then begin
print,'% CHAIN: no Spe images available, only Cal frames will collapsed'
endif
if (max(wdar) lt 0) then begin
print,'% CHAIN: no Dar image available -- the dark current will not be subtracted'
endif else begin
if (n_elements(wdar) eq 1) then begin
dark = readfits(st[wdar].filename)
endif else begin
print,'% CHAIN: Warning -- multiple Dar images available -- the first one will be used'
dark = readfits(st[wdar[0]].filename)
endelse
endelse
;average binned flats and extract average
printf,10,'creating average binned flat ...'
print,'creating average flat ...'
imadd,st[wfla].filename,flatname,/av
f = readfits(flatname,header)
;find order information from flat...
printf,10,'find order information from flat spectrum...'
dispdir,f,idisp
hbias,f,/bin
;we either find the apertures on the fly from the flat
if (1 eq 0) then begin
hfind,f,idisp,ap1,delta1,fwhm=4,width=0.65,smoothinglength=2,xorder=2
if n_elements(ap1) ne 27 then begin
print,'% CHAIN: ERROR -- did not find exactly 27 orders!'
print,' -- this is a requirement for xcal.'
printf,10,'% CHAIN: ERROR -- did not find exactly 27 orders!'
printf,10,' -- this is a requirement for xcal.'
close,10
endif
endif else begin
; or we take it from a reference image and use a ref flat to find the appropropriate offset
ap1=readfits(chaindir+'/rap1.fits')
ap=readfits(chaindir+'/rap.fits')
;derive spatial-direction offset of orders using flat
rflat=readfits(chaindir+'/rflat.fits')
trflat=total(rflat,2)
tflat=total(f,2)
;if not keyword_set(bin) then tflat=rebin(tflat,514)
xc,trflat,tflat,trflat/100.,tflat/100.,fshift,efshift
ap1=ap1-fshift
;ap1=ap1/3.-fshift ;old ref. flat
ap=ap-fshift
width=0.72
endelse
;compute delta1 (order width)
;it must be done consistently in the inspect/collapse (1/2) routines
delta1=median(ap1-shift(ap1,1))/2.*width
np=n_elements(ap[0,*])
delta=median(ap[*,np/2]-shift(ap[*,np/2],1))/2.*width
;save aperture info
writefits,strcompress(ap1name,/rem),ap1
writefits,apname,ap
printf,10,'delta1=',delta1
print,'delta1=',delta1
printf,10,'delta=',delta
print,'delta=',delta
printf,10,'fshift=',fshift
print,'fshift=',fshift
;find actual left/right limits for each aperture
xwindows,ap,delta,left,right
gain = sxpar(header,'GAIN')
rdnoise = sxpar(header,'RDNOISE')
f = f * gain
vf = f + rdnoise^2
collapse, f, idisp, left,right, xf, vf= vf, vs=xfv ,/clean
ws,'x'+flatname,xf,xfv, hd=header
;remove cosmics, scattered light and extract spes
printf,10,'removing cosmics and scattered light and extracting spes ...'
print,'removing cosmics and scattered light and extracting spes ...'
for i=0,n_elements(wspe)-1 do begin
print,i
j=wspe[i]
filename=st[j].filename
frame = float(readfits(filename,header))
if st[j].exptime gt la_minexptime then begin
gain=sxpar(header,'gain')
rdnoise=sxpar(header,'rdnoise')
writefits,'tmp.fits',frame,header
la_cosmic,'tmp.fits',gain=gain,readn=rdnoise,sigfrac=20.,sigclip=4.5
frame = readfits('tmp-out.fits')
file_delete,'tmp.fits','tmp-out.fits','tmp-mask.fits'
endif
if n_elements(dark) gt 0 then begin
print,'subtracting dark to image ',filename
frame = frame - dark/mean(dark[0:100,*])*mean(frame[0:100,*])
rdn=stddev(frame[*,1:16])
endif else begin
print,'subtracting bias to image ',filename
hbias,frame,rdn=rdn,/bin
endelse
if scat eq 'yes' then begin
scatter,frame,idisp,delta1,back
frame = frame - back
endif
gain=sxpar(header,'GAIN')
rdnoise=sxpar(header,'RDNOISE')
printf,10,filename,' rdnoise=',rdnoise/gain,rdn,' (counts)'
print,filename,' rdnoise=',rdnoise/gain,rdn,' (counts)'
frame = frame * gain
vframe = frame + rdnoise^2
collapse, frame, idisp, left, right, xframe, vf = vframe, vs = xvframe ,/clean
xfilename='x'+strmid(filename,strpos(filename,'/',/reverse_search)+1)
if strmid(xfilename,strlen(xfilename)-3,3) eq '.gz' then $
xfilename=strmid(xfilename,0,strlen(xfilename)-3)
;upgrade header with HORuS coords (HORUSRA/HORUSDEC) and reduction time stamp
hheader,header
;write extracted file
ws, xfilename, xframe, xvframe, hd=header
;check memory usage
;mem,used,tot,avail
;print,'memory used-total-avail (MB) = ',used,tot,avail
;printf,10,'memory used-total-avail (MB) = ',used,tot,avail
endfor
;extract cals
printf,10,'extracing cals ...'
print,'extracing cals ...'
for i=0,n_elements(wcal)-1 do begin
j=wcal[i]
filename=st[j].filename
frame = readfits(filename,header)
hbias,frame,rdn=rdn,/bin
gain=sxpar(header,'GAIN')
rdnoise=sxpar(header,'RDNOISE')
printf,10,filename,' rdnoise=',rdnoise/gain,rdn,' (counts)'
print,filename,' rdnoise=',rdnoise/gain,rdn,' (counts)'
frame = frame * gain
vframe = frame + rdnoise^2
collapse, frame, idisp, left, right, xframe, vf = vframe, vs = xvframe ,/clean
xfilename='x'+strmid(filename,strpos(filename,'/',/reverse_search)+1)
if strmid(xfilename,strlen(xfilename)-3,3) eq '.gz' then $
xfilename=strmid(xfilename,0,strlen(xfilename)-3)
;upgrade header with HORuS coords (HORUSRA/HORUSDEC) and reduction time stamp
hheader,header
;write extratcted file
ws, xfilename, xframe, xvframe, hd=header
;check memory usage
;mem,used,tot,avail
;print,'memory used-total-avail (MB) = ',used,tot,avail
;printf,10,'memory used-total-avail (MB) = ',used,tot,avail
endfor
;wave cal. cals
printf,10,'wavelength calibration of cals ...'
print,'wavelength calibration of cals ...'
for i=0,n_elements(wcal)-1 do begin
j=wcal[i]
filename=st[j].filename
xfilename='x'+strmid(filename,strpos(filename,'/',/reverse_search)+1)
if strmid(xfilename,strlen(xfilename)-3,3) eq '.gz' then $
xfilename=strmid(xfilename,0,strlen(xfilename)-3)
printf,10,'calibrating ... '+xfilename
print,'calibrating ... '+xfilename
rs, xfilename, xframe, xvframe, hd=header
xcal2, xframe, wframe,calstats=cs, chaindir=chaindir
ratio=cs[6,*]/cs[2,*]
minratio=min(ratio)
maxratio=max(ratio)
medratio=median(ratio)
print,'min-max-median of rms/dispersion =',minratio,maxratio,medratio,format='(a30,3(1x,f10.5))'
printf,10,'min-max-median of rms/dispersion =',minratio,maxratio,medratio,format='(a30,3(1x,f10.5))'
ws, xfilename, xframe, xvframe, w = wframe, hd=header
if i eq 0 then begin
calfiles = xfilename
calmjd = st[j].mjd0
endif else begin
calfiles = [ calfiles, xfilename ]
calmjd = [ calmjd, st[j].mjd0 ]
endelse
endfor
;average flatfield (approx. cal. using last cal)
rs,'x'+flatname, xframe, xvframe, hd=header
ws, 'x'+flatname, xframe, xvframe, w = wframe, hd=header
print,calfiles
print,calmjd
;spec
for i=0,n_elements(wspe)-1 do begin
j=wspe[i]
filename=st[j].filename
xfilename='x'+strmid(filename,strpos(filename,'/',/reverse_search)+1)
if strmid(xfilename,strlen(xfilename)-3,3) eq '.gz' then $
xfilename=strmid(xfilename,0,strlen(xfilename)-3)
nfilename='n'+strmid(filename,strpos(filename,'/',/reverse_search)+1)
if strmid(nfilename,strlen(nfilename)-3,3) eq '.gz' then $
nfilename=strmid(nfilename,0,strlen(nfilename)-3)
printf,10,'calibrating ... '
print,'calibrating ... '
printf,10,xfilename+' mjd=',st[j].mjd0
print,xfilename+' mjd=',st[j].mjd0
rs, xfilename, xframe, xvframe, norder=norder, hd=header
creject, xframe, xframe2
xframe = xframe2
wpre = max(where(st[j].mjd0-calmjd gt 0.))
wpos = min(where(st[j].mjd0-calmjd lt 0.))
if max(wpre) gt -1 then begin
printf,10,' 1-'+calfiles[wpre]+' mjd=',calmjd[wpre]
print,' 1-'+calfiles[wpre]+' mjd=',calmjd[wpre]
rs, calfiles[wpre], x1, xv1, w = w1
sxaddpar,header,'FILECAL1',calfiles[wpre],'Cal. spectrum preceding observation'
sxaddpar,header,'MJDCAL1',calmjd[wpre],'MJD Ca. spectrum preceding'
endif
if max(wpos) gt -1 then begin
printf,10,' 2-'+calfiles[wpos]+' mjd=',calmjd[wpos]
print,' 2-'+calfiles[wpos]+' mjd=',calmjd[wpos]
rs, calfiles[wpos], x2, xv2, w = w2
sxaddpar,header,'FILECAL2',calfiles[wpos],'Cal. spectrum following observation'
sxaddpar,header,'MJDCAL2',calmjd[wpos],'MJD Ca. spectrum following'
endif
if max(wpre) gt -1 and max(wpos) gt -1 then begin
d = (st[j].mjd0 - calmjd[wpre] ) / (calmjd[wpos] - calmjd[wpre] )
endif else begin
if max(wpre) gt -1 and max(wpos) lt 0 then d=0.
if max(wpre) lt 0 and max(wpos) gt -1 then d=1.
if max(wpre) lt 0 and max(wpos) lt 0 then d=-1.
endelse
if d lt 0. then begin
printf,10,'ERROR -- cannot wave calibrate!'
print,'ERROR -- cannot wave calibrate!'
endif else begin
if abs(d) lt 1.e-7 then w2=w1
if abs(d-1.) lt 1.e-7 then w1=w2
wframe = (1. - d) * w1 + d * w2
ws, xfilename, xframe, xvframe, w = wframe, hd=header
;normalize and order merge to create the m* files
if xmformat ne '' then xm, xfilename,xmformat=xmformat
;normalize by the flatfield
for k=0,norder-1 do begin
ratio = 1.0d0 ; median(xf[k,*])/median(xframe[k,*])
xframe[k,*] = xframe[k,*] / xf[k,*] * ratio
xvframe[k,*] = xvframe[k,*] / xf[k,*]^2 * ratio^2
endfor
;trim
margin=150
np=n_elements(xframe[0,*])
xframe = xframe[*,margin:np-1-margin]
xvframe = xvframe[*,margin:np-1-margin]
wframe = wframe[*,margin:np-1-margin]
ws, nfilename, xframe, xvframe, w = wframe, hd=header
endelse
endfor
close,10
;create plots
;mkpng
cmd = 'python3 '+chaindir+'/horus.py'
spawn,cmd,back
print,back
end