dustem_cc.pro
8.29 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
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
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:
; DustEMWrap, Mid-Level, Distributed
; 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',
; 'FLAMBDA=cste', 'HFI', 'LABOCA'
; help = If set, print this help
; COMMON BLOCKS:
; None
; SIDE EFFECTS:
; None
; RESTRICTIONS:
; The DustEM fortran code must be installed
; The DustEMWrap IDL code must be installed
; PROCEDURE:
; None
; EXAMPLES
;
; MODIFICATION HISTORY:
; Written JPB Apr-2011
; Evolution details on the DustEMWrap gitlab.
; See http://dustemwrap.irap.omp.eu/ for FAQ and help.
;-
IF keyword_set(help) THEN BEGIN
doc_library,'dustem_cc'
sed=0.
goto,the_end
ENDIF
;eps=1.e-30
eps=0.
;=== constants needed for CMB flux convention
Tcmb=2.725D0
;=== constants 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 BEGIN
message,'Instrument '+instrus(i)+' not found',/continue
stop
ENDIF
st=(*!dustem_filters).(iinstru)
; stop
ii=(where(st.filter_names EQ filter_names[i]))[0]
spec0=interpol(spec2,wavein2,st.central_wavelengths[ii])
IF spec0 EQ 0 THEN BEGIN
sed[i]=spec0
cc[i]=1
goto,next_filter
ENDIF
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
'FLAMBDA=CSTE': BEGIN
; stop
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_wmin[ii],st.use_wmax[ii],/double)
; print,st.use_wmin[ii],st.use_wmax[ii]
; print,st.central_wavelengths[ii]
; t=*(st.use_transmissions[ii])
; w=*(st.use_wavelengths[ii])
; xx=where(t gt 0)
; print,min(w[xx]),max(w[xx])
; stop
; print,spec0*den
; stop
cc[i]=(num*(st.central_wavelengths[ii])^2)/(spec0*den)
; print,filter_names[i],st.central_wavelengths[ii],cc[i],st.use_wmin[ii],st.use_wmax[ii],minmax(*(st.use_wavelengths[ii]))
; stop
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
ELSE:BEGIN
message,strupcase(fluxconvs[i])+' not recognized',/info
cc[i]=1.
;stop
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
IF finite(cc[i]) THEN sed[i]=spec0*cc[i] ELSE sed[i]=spec0
next_filter:
ENDFOR
the_end:
RETURN,sed
END