dustem_cc.pro
7.13 KB
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
FUNCTION dustem_cc,wavein,spec,filter_names,cc=cc,fluxconv=fluxconv,help=help
;+
; NAME:
; dustem_cc
; PURPOSE:
; Computes observed sed and color correction for a given spectrum in a given filter
; CATEGORY:
; Dustem
; CALLING SEQUENCE:
; sed=dustem_cc(wave,spec,filter_names,cc=cc,fluxconv=fluxconv,help=help)
; INPUTS:
; wave: array of wavelengths for spec
; spec: spectrum (must be in brightness units)
; filter_names: names of filters for which color correction is needed
; OPTIONAL INPUT PARAMETERS:
; None
; OUTPUTS:
; sed: SED as observed in filters provided in filter_names
; OPTIONAL OUTPUT PARAMETERS:
; cc = color correction coefficients (spec*cc=sed)
; ACCEPTED KEY-WORDS:
; fluxconv = if set, these are taken for flux conventions
; possible values: 'nuInu=cste', 'IRAC', 'MIPS', 'CMB'
; help = If set, print this help
; COMMON BLOCKS:
; None
; SIDE EFFECTS:
; None
; RESTRICTIONS:
; The dustem idl wrapper must be installed
; PROCEDURE:
; None
; EXAMPLES
;
; MODIFICATION HISTORY:
; Written by J.P. Bernard April 18th 2011
; replaces individual cc routines for various instruments
; see evolution details on the dustem cvs maintained at CESR
; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-
;For MIPS color corrections, see
;http://irsa.ipac.caltech.edu/data/SPITZER/docs/mips/mipsinstrumenthandbook/51/
IF keyword_set(help) THEN BEGIN
doc_library,'dustem_cc'
sed=0.
goto,the_end
ENDIF
;eps=1.e-30
eps=0.
;=== needed for CMB flux convention
Tcmb=2.725D0
;=== needed for MIPS flux convention
h=6.62606876*1d-34 ; J s
c=2.99792458*1d8 ; m*s-1
K=1.380658*1d-23 ; J/K
T_0=1d4 ; calibration star (K units)
fact=h*c
;fact1=h*c
fact3=K*T_0
; sort input
order =sort(wavein)
spec2 = [eps,spec(order)]
wavein2 = [eps,wavein(order)]
IF not keyword_set(fluxconv) THEN BEGIN
fluxconvs=dustem_filter2fluxconv(filter_names)
ENDIF ELSE BEGIN
fluxconvs=fluxconv
ENDELSE
instrus=dustem_filter2instru(filter_names)
inst_tags=tag_names((*!dustem_filters))
;DEBUG VG
j=where(instrus eq 'HFI',count)
if count gt 0 then instrus(j) = 'HFI_V3P02'
Nbands=n_elements(filter_names)
sed=dblarr(Nbands)
cc=dblarr(Nbands)
cc(*)=1.d0
FOR i=0L,Nbands-1 DO BEGIN
iinstru=where(inst_tags EQ instrus(i),countinstr)
IF countinstr EQ 0 THEN stop
st=(*!dustem_filters).(iinstru)
ii=(where(st.filter_names EQ filter_names(i)))(0)
spec0=interpol(spec2,wavein2,st.central_wavelengths(ii))
IF !dustem_do_cc NE 0 and not !dustem_never_do_cc THEN BEGIN ;Color corrections are computed
; print,'DOING color correction calculations'
CASE strupcase(fluxconvs(i)) OF
'HFI0': BEGIN ;This is purely for test (stupid since reads filters each time)
filter_dir=!dustem_wrap_soft_dir+'/Data/FILTERS/'
dir_hfi=filter_dir+'HFI/OFFICIAL'+'/'
; CAUTION: This below should be consistent with what is in dustem_read_filters.pro
; hfi_filter_version='v101'
hfi_filter_version='v201'
bp = hfi_read_bandpass(hfi_filter_version, /rimo,path_rimo=dir_hfi)
spec_int=interpol(spec2,wavein2,*(st.use_wavelengths(ii)))
nu=*(st.use_frequencies(ii))/1.d9
order=sort(nu)
spec_int=spec_int(order)
nu=nu(order)
channel=strtrim(round(dustem_filter2freq(filter_names(i))),2)
;stop
ccc=hfi_color_correction_dustem(bp.name,bp.freq/1.d9,bp.trans,channel, bolo_cc,nu, spec_int)
; ccc=hfi_color_correction(bp,channel, bolo_cc,nu, spec_int)
cc(i)=ccc.cc
END
'HFI': BEGIN
spec_int=interpol(spec2,wavein2,*(st.use_wavelengths(ii)))
nu=*(st.use_frequencies(ii))/1.d9
trans=*(st.use_transmissions(ii))
order=sort(nu)
spec_int=spec_int(order)
nu=nu(order)
trans=trans(order)
channel=strtrim(round(dustem_filter2freq(filter_names(i))),2)
ccc=hfi_color_correction_dustem(channel,nu,trans,channel, bolo_cc,nu, spec_int)
cc(i)=ccc.cc
END
'NUINU=CSTE': BEGIN
spec_int=interpol(spec2,wavein2,*(st.use_wavelengths(ii)))
num=integral(*(st.use_wavelengths(ii)),spec_int*(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii)))^2.,st.use_wmin(ii),st.use_wmax(ii),/double)
den=integral(*(st.use_wavelengths(ii)),(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii))),st.use_wmin(ii),st.use_wmax(ii),/double)
cc(i)=num/den/spec0*(st.central_wavelengths(ii))
END
'MIPS': BEGIN
expon=exp(fact/(*(st.use_wavelengths(ii))*1e-6*fact3))-1.
spec_int=interpol(spec2,wavein2,*(st.use_wavelengths(ii)))
num=integral(*(st.use_wavelengths(ii)),spec_int*(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii)))^2.,st.use_wmin(ii),st.use_wmax(ii),/double)
den=integral(*(st.use_wavelengths(ii)),*(st.use_transmissions(ii))/(*(st.use_wavelengths(ii)))^5/expon,st.use_wmin(ii),st.use_wmax(ii),/double)
ffact=st.central_wavelengths(ii)*1e-6*K*T_0
cc(i)=num/den/spec0/(st.central_wavelengths(ii))^3/(exp(fact/ffact)-1.)
END
'IRAC': BEGIN
spec_int=interpol(spec2,wavein2,*(st.use_wavelengths(ii)))
num=integral(*(st.use_wavelengths(ii)),spec_int*(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii))),st.use_wmin(ii),st.use_wmax(ii),/double)
den=integral(*(st.use_wavelengths(ii)),*(st.use_transmissions(ii)),st.use_wmin(ii),st.use_wmax(ii),/double)
cc(i)=num/den/spec0*(st.central_wavelengths(ii))
END
'CMB': BEGIN
spec_int=interpol(spec2,wavein2,*(st.use_wavelengths(ii)))
num=integral(*(st.use_wavelengths(ii)),spec_int*(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii)))^2.,st.use_wmin(ii),st.use_wmax(ii),/double)
den=integral(*(st.use_wavelengths(ii)),(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii))^2)*dustem_planck_function(Tcmb,*(st.use_wavelengths(ii))),st.use_wmin(ii),st.use_wmax(ii),/double)
cc(i)=num/den*dustem_planck_function(Tcmb,st.central_wavelengths(ii))/spec0
END
'LABOCA': BEGIN
; spec_int=interpol(spec2,wavein2,*(st.use_wavelengths(ii)))
; num=integral(*(st.use_wavelengths(ii)),spec_int*(*(st.use_transmissions(ii))),st.use_wmin(ii),st.use_wmax(ii),/double)
; den=integral(*(st.use_wavelengths(ii)),*(st.use_transmissions(ii))/(*(st.use_wavelengths(ii))^2.5),st.use_wmin(ii),st.use_wmax(ii),/double)
; cc(i)=num/den/spec0/(st.central_wavelengths(ii))^2.5
alpha_new=0.
spec_int=interpol(spec2,wavein2,*(st.use_wavelengths(ii)))
num=integral(*(st.use_wavelengths(ii)),spec_int*(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii))),st.use_wmin(ii),st.use_wmax(ii),/double)
den=integral(*(st.use_wavelengths(ii)),*(st.use_transmissions(ii))/(*(st.use_wavelengths(ii))^(alpha_new+2.)),st.use_wmin(ii),st.use_wmax(ii),/double)
cc(i)=num/den/spec0*(st.central_wavelengths(ii))^(-(alpha_new+1.))
END
ENDCASE
!dustem_previous_cc=ptr_new(cc)
ENDIF ELSE BEGIN ;When no new cc calculation is required
IF !dustem_never_do_cc EQ 0 THEN BEGIN
cc=(*!dustem_previous_cc)
ENDIF
ENDELSE
sed(i)=spec0*cc(i)
ENDFOR
the_end:
RETURN,sed
END