-
Notifications
You must be signed in to change notification settings - Fork 1
/
flux_calib2.pro
141 lines (123 loc) · 5.66 KB
/
flux_calib2.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
;;
;; version 2 of flux_calib: split into "get_CF" and "flux_calib"
@$SINFOTOOLS/get_atmo
;;
;;
;; FUNCTION get_CF
;;
;; PURPOSE
;; given H,K magnitude of a star and the reduced cube (units: counts),
;; extract 1D spectrum, derive conversion factor to be applied to
;; science cube
;;
function get_CF, Hmag,Kmag,cube_std_atmo,spectrum_std_atmo,std_name, grat1=grat1
;;
;; extract spectrum from cube_std_atmo in large aperture to include all flux
spec = cube2spec(cube_std_atmo, spectrum_std_atmo, 10, wave=wave)
stdspec = read_fits_wave(spectrum_std_atmo,naxis=1,hdr=hdr)
dateobs=sxpar(hdr,"DATE-OBS")
cmd='sqlite3 $OBSDB "select grat1 from sinfo where dateobs=\"'+dateobs+'\";"'
spawn, cmd, grat1
spawn, 'dfits ' + cube_std_atmo + " | grep EXPTIME | awk -F ' ' '{ print $3 }'", exptime_cal
exptime_cal = double(exptime_cal[0])
params=get_sinfo_drs_parameters()
H0=params[where(params.name eq 'H0')].value ; H band zero point in W/(m^2 um)
K0=params[where(params.name eq 'K0')].value ; K band zero point in W/(m^2 um)
H=[params[where(params.name eq 'H_lower')].value,params[where(params.name eq 'H_upper')].value]
K=[params[where(params.name eq 'K_lower')].value,params[where(params.name eq 'K_upper')].value]
twomass_k_lam=float(params[where(params.name eq 'twomass_k_lam')].value)
twomass_k_dlam=float(params[where(params.name eq 'twomass_k_dlam')].value)
twomass_k0=double(params[where(params.name eq 'twomass_k0')].value)
twomass_h_lam=float(params[where(params.name eq 'twomass_h_lam')].value)
twomass_h_dlam=float(params[where(params.name eq 'twomass_h_dlam')].value)
twomass_h0=double(params[where(params.name eq 'twomass_h0')].value)
; should be using 2MASS wavelength, bandwidth and zeropoint, but it makes almost no difference, so I will stick to the standard K band for the time being.
; H0=twomass_h0
; K0=twomass_k0
; H=[twomass_h_lam-twomass_h_dlam/2,twomass_h_lam+twomass_h_dlam/2]
; K=[twomass_k_lam-twomass_k_dlam/2,twomass_k_lam+twomass_k_dlam/2]
; H/K flux in 1e-17 W/(m^2 um)
Hflux = 1.e17 * H0 * 10^(-Hmag/2.5)
Kflux = 1.e17 * K0 * 10^(-Kmag/2.5)
if grat1 eq 'H+K' then begin
inband_H=where(stdspec.w ge H[0] and stdspec.w le H[1])
inband_K=where(stdspec.w ge K[0] and stdspec.w le K[1])
cts_H = median(stdspec.f[inband_H])
cts_K = median(stdspec.f[inband_K])
; print, "standard cts_K:"
; print, cts_K
; conversion factor from counts/s to W/(m^2 um)
CF_H = Hflux * exptime_cal / cts_H
CF_K = Kflux * exptime_cal / cts_K
CF_avg = (CF_H+CF_K)/2.
spawn, 'sinfolog "Conversion factors for ' + std_name + ' at ' + dateobs + ' in H/K/avg: ' + string(CF_H, CF_K, CF_avg) + '"'
endif else if grat1 eq 'K' then begin
inband_K=where(stdspec.w ge K[0] and stdspec.w le K[1])
cts_K = median(stdspec.f[inband_K])
; conversion factor from counts/s to W/(m^2 um)
CF_K = Kflux * exptime_cal / cts_K
CF_avg = CF_K
spawn, 'sinfolog "Conversion factor for ' + std_name + ' at ' + dateobs + ' in K: ' + string(CF_K) + '"'
endif else stop
return, CF_avg
end
;;
;; given a certain H and K band magnitude, return conversion factors
;;
;; INPUT values
;;
;; Hmag, Kmag -- magnitudes of star
;; cube_std_atmo -- cube of standard star with atmospheric transmission spectrum removed
;; spectrum_std_atmo -- spectrum extracted from cube of standard star with atmospheric transmission removed
;; scifile -- optional: apply the average conversion factor to science file, taking exposure time into account
;;
;; USAGE (from bash command line)
;; idl -e flux_calib -args "
;;
;pro flux_calib, Hmag, Kmag, cube_std_atmo, spectrum_std_atmo, scifile=scifile
pro flux_calib2
args=command_line_args()
Hmag=float(args[0])
Kmag=float(args[1])
cube_std_atmo=args[2]
spectrum_std_atmo=args[3]
scifile=args[4]
outfile=args[5]
std_name=args[6]
CF_avg = get_CF(Hmag,Kmag,cube_std_atmo,spectrum_std_atmo,std_name, grat1=grat1)
params=get_sinfo_drs_parameters()
H0=params[where(params.name eq 'H0')].value ; H band zero point in W/(m^2 um)
K0=params[where(params.name eq 'K0')].value ; K band zero point in W/(m^2 um)
H=[params[where(params.name eq 'H_lower')].value,params[where(params.name eq 'H_upper')].value]
K=[params[where(params.name eq 'K_lower')].value,params[where(params.name eq 'K_upper')].value]
if keyword_set(scifile) then begin
spawn, 'dfits ' + scifile + " | grep EXPTIME | awk -F ' ' '{ print $3 }'", exptime_sci
exptime_sci = double(exptime_sci[0])
; print, exptime_sci, format='("Exposure time SCIENCE is: ", f8.3, " s")'
cube = readfits(scifile,hdr)
CF_sci = CF_avg / exptime_sci
spawn, 'sinfolog "' + string(CF_sci, format='("Multiplying cube by ", e10.4)') + '"'
cube_cal = cube * CF_sci
writefits, outfile, cube_cal, hdr
spawn, 'sinfolog "Wrote ' + outfile + '"'
;;
;; extract flux from science cube and determine magnitude of nucleus
spectrum_obj_cal='/tmp/spectrum_obj_cal.fits'
spec = cube2spec(outfile, spectrum_obj_cal, 10, wave=wave)
spec_obj_cal = read_fits_wave(spectrum_obj_cal,naxis=1,hdr=hdr)
if grat1 eq 'H+K' then begin
inband_H=where(spec_obj_cal.w ge H[0] and spec_obj_cal.w le H[1])
inband_K=where(spec_obj_cal.w ge K[0] and spec_obj_cal.w le K[1])
cts_H = median(spec_obj_cal.f[inband_H])
cts_K = median(spec_obj_cal.f[inband_K])
Hmag_obj=2.5*alog10(H0/(1.d-17 * cts_H))
Kmag_obj=2.5*alog10(K0/(1.d-17 * cts_K))
spawn, 'sinfolog "Nuclear magnitude of target in H/K: ' + string(Hmag_obj, Kmag_obj) + '"'
endif else if grat1 eq 'K' then begin
inband_K=where(spec_obj_cal.w ge K[0] and spec_obj_cal.w le K[1])
cts_K = median(spec_obj_cal.f[inband_K])
Kmag_obj=2.5*alog10(K0/(1.d-17 * cts_K))
spawn, 'sinfolog "Nuclear magnitude of target in K: ' + string(Kmag_obj) + '"'
endif else stop
endif
end