-
Notifications
You must be signed in to change notification settings - Fork 1
/
mage_quicklook.pro
146 lines (119 loc) · 5.14 KB
/
mage_quicklook.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
pro mage_quicklook, fitsfile, psout, SNR=snr
ARCHIVE = getenv('MAGE_DIR')+'/ARCHIVE/'
slitfile = strtrim(ARCHIVE,2)+'mage_archive_orders.fits'
sensfile = strtrim(ARCHIVE,2)+'GD108_mar09_std.sav'
wavefile = strtrim(ARCHIVE,2)+'arcimg.fits.gz'
pixfile = strtrim(ARCHIVE,2)+'Piximg.fits'
ostruct = strtrim(ARCHIVE,2)+'OStr_mage.fits'
filstd = strtrim(ARCHIVE,2)+'ObjStr_GD108.fits'
; Generate a mask of the order numbers
tset_slits = mrdfits(slitfile, 1)
slitmask = long_slits2mask(tset_slits)
ordermask = 0*slitmask
ordermask[WHERE(slitmask GT 0)] = -slitmask[WHERE(slitmask GT 0)] + 21L
; Read in the files
exptime=0
filenames = [fitsfile]
nfiles = n_elements(filenames)
IF nfiles GT 0 THEN BEGIN
sciimg = 0
var_tot = 0
FOR ifile = 0L, nfiles-1L DO BEGIN
mage_proc, filenames[ifile], sciimg1, scivar1, hdr=header ;, pixflatfile=pixflatfile
exptime+=fxpar(header,'EXPTIME')
sciimg = sciimg + sciimg1
var_tot = var_tot + 1.0D/(scivar1 + (scivar1 EQ 0))
var_tot = var_tot*(scivar1 GT 0)
ENDFOR
scivar = (var_tot GT 0)/(var_tot + (var_tot EQ 0))
ENDIF ELSE begin
mage_proc, filenames, sciimg1, scivar, hdr=header ; pixflatfile=pixflatfile
exptime+=fxpar(header,'EXPTIME')
endelse
ximg = long_slits2x(tset_slits, edgmask = edgmask,TOL_EDG=3)
waveimg = xmrdfits(wavefile, 0, /fscale)
piximg = xmrdfits(pixfile, 0)
; Mask out the object for first pass sky subtraction
bsp = 0.6
FWHM = 3.0
print, "Masking out the object for first pass sky subtraction..."
objstruct1 = long_objfind(sciimg, tset_slits = tset_slits, invvar = sciivar $
, skymask = skymask1, objmask = objmask1 $
, nperslit = 1L, peakthresh = reduxthresh $
, fwhm = FWHM, ISLIT = ISLIT, /silent)
print, "Performing 2D sky subtraction"
skyimage = long_skysub(sciimg, sciivar, piximg, slitmask, skymask1, edgmask $
, bsp = bsp, ISLIT = ISLIT)
box_rad = 5L
velpix = 22.0d
;; The order struct is output by running mage_traceorders.pro
ordr_strct = xmrdfits(ostruct, 1)
obj_strct = m_mkobjstr(15)
for i=0,14 do obj_strct[i].exp=exptime
sciivar = 0.0*sciivar
sciivar[where(skyimage GT 0)] = 1.0d/(skyimage[where(skyimage GT 0)])
; Object Finding
print, "Fancy object finding..."
;m_fntobj, obj_strct, ordr_strct, sciimg-skyimage, sciivar, qafil
obj_strct=mage_findobj(sciimg-skyimage,sciivar,waveimg,tset_slits $
,filstd=filstd,chk=chk)
; Sanity check
;xatv, (sciimg-skyimage)*(slitmask GT 0), min = -20.0, max = 200.0 $
; , wv = waveimg, /block
; for qq=0L,14 do xatvplot, obj_strct[qq].trace[0:2048-1], findgen(2048)
; Quick boxcar extract
print, "Boxcar extracting object..."
;m_extechopt, sciimg, sciimg-skyimage, sciivar, ordr_strct, obj_strct, $
; velpix, img_arc=waveimg, skyfil=skyimage, helio=0.0, $
; obj_name="Quicklook", ordermask=ordermask, /boxonly
outfil='ql_obj.fits'
mage_echextobj,sciimg,sciivar,header,skyimage $
,piximg,waveimg,tset_slits $
,obj_strct,outfil=outfil $
,/boxcar, box_rad=box_rad
obj_strct.exp=float(sxpar(header,'EXPTIME'))
obj_strct.wave = obj_strct.box_wv
obj_strct.fx = obj_strct.box_fx
obj_strct.var = obj_strct.box_var
mage_flux, sensfile, obj_strct, rej=0.05
mage_combspec, obj_strct, fspec
set_plot, "ps"
device, filename=psout, /landscape, /color
colors=getcolor(/load)
res = 299792.458/4100.*0.7
print, "Combining to 1D spectrum..."
mage_1dspec, fspec, "ql.fits", "ql_sig.fits", "ql_comb.fits", resvel= res ;, /rebin
t1 = mrdfits('ql.fits',0,hdr)
t2 = mrdfits('ql_sig.fits')
wv=10^(sxpar(hdr,"CRVAL1")+dindgen(n_elements(t1))*sxpar(hdr, "CDELT1"))
!p.multi=[0,1,2]
sf = sort(t1)
ninety = t1[sf[0.97*n_elements(sf)]]
plot, wv, t1, xrange=[3200,6500], yrange=[-0.1,1.7*ninety], /xsty, /ysty, xtitle="Wavelength (A)", ytitle="Relative Flux"
oplot, wv, t2, color=colors.green
oplot, wv, t1-t1, color=colors.red
plot, wv, t1, xrange=[6500,10000], yrange=[-0.1,1.5*ninety], /xsty, /ysty, xtitle="Wavelength (A)", ytitle="Relative Flux"
oplot, wv, t2, color=colors.green
oplot, wv, t1-t1, color=colors.red
!p.multi=[0,2,2]
for i=0, 14 do begin
if not keyword_set(SNR) then begin
sf = sort(obj_strct[i].fx)
ninety = obj_strct[i].fx[sf[0.95*n_elements(sf)]]
medfx = median(obj_strct[i].fx)
plot, obj_strct[i].wave, obj_strct[i].fx, yrange=[-0.1*medfx, 1.4*ninety], xtitle="Wavelength (A)", ytitle="Counts"
oplot, obj_strct[i].wave, sqrt(obj_strct[i].var), color=colors.green
endif else begin
sf = sort(obj_strct[i].fx)
ninety = obj_strct[i].fx[sf[0.95*n_elements(sf)]]
medfx = median(obj_strct[i].fx)
mednoise = median(sqrt(obj_strct[i].var))
plot, obj_strct[i].wave, obj_strct[i].fx/sqrt(obj_strct[i].var), yrange=[-0.1*medfx/mednoise, 1.4*ninety/mednoise]
endelse
endfor
device, /close
set_plot, "x"
!p.multi=0
x_specplot, "ql.fits", "ql_sig.fits", /block
; stop
end