Blame view

src/idl/dustem_cc.pro 8.34 KB
427f1205   Jean-Michel Glorian   version 4.2 merged
1
2
3
4
5
FUNCTION dustem_cc,wavein,spec,filter_names,cc=cc,fluxconv=fluxconv,help=help

;+
; NAME:
;    dustem_cc
1e005c88   Annie Hughes   inlie help
6
;
427f1205   Jean-Michel Glorian   version 4.2 merged
7
; PURPOSE:
1e005c88   Annie Hughes   inlie help
8
9
;    Computes observed SED and color correction for a given spectrum in a given filter
;
427f1205   Jean-Michel Glorian   version 4.2 merged
10
; CATEGORY:
1e005c88   Annie Hughes   inlie help
11
12
;    DustEMWrap, Mid-Level, Distributed, Color corrections
;
427f1205   Jean-Michel Glorian   version 4.2 merged
13
14
; CALLING SEQUENCE:
;    sed=dustem_cc(wave,spec,filter_names,cc=cc,fluxconv=fluxconv,help=help)
1e005c88   Annie Hughes   inlie help
15
;
427f1205   Jean-Michel Glorian   version 4.2 merged
16
17
18
; INPUTS:
;    wave: array of wavelengths for spec
;    spec: spectrum (must be in brightness units)
1e005c88   Annie Hughes   inlie help
19
20
21
;    filter_names: names of filters for which color correction is
;    needed
;  
427f1205   Jean-Michel Glorian   version 4.2 merged
22
23
; OPTIONAL INPUT PARAMETERS:
;    None
1e005c88   Annie Hughes   inlie help
24
;
427f1205   Jean-Michel Glorian   version 4.2 merged
25
26
; OUTPUTS:
;    sed: SED as observed in filters provided in filter_names
1e005c88   Annie Hughes   inlie help
27
;
427f1205   Jean-Michel Glorian   version 4.2 merged
28
29
; OPTIONAL OUTPUT PARAMETERS:
;    cc        = color correction coefficients (spec*cc=sed)
1e005c88   Annie Hughes   inlie help
30
;
427f1205   Jean-Michel Glorian   version 4.2 merged
31
32
; ACCEPTED KEY-WORDS:
;    fluxconv  = if set, these are taken for flux conventions
997bd805   Annie Hughes   comments and gene...
33
34
;              possible values: 'nuInu=cste', 'IRAC', 'MIPS', 'CMB',
;              'FLAMBDA=cste', 'HFI', 'LABOCA'
427f1205   Jean-Michel Glorian   version 4.2 merged
35
;    help      = If set, print this help
1e005c88   Annie Hughes   inlie help
36
;
427f1205   Jean-Michel Glorian   version 4.2 merged
37
38
; COMMON BLOCKS:
;    None
1e005c88   Annie Hughes   inlie help
39
;
427f1205   Jean-Michel Glorian   version 4.2 merged
40
41
; SIDE EFFECTS:
;    None
1e005c88   Annie Hughes   inlie help
42
;
427f1205   Jean-Michel Glorian   version 4.2 merged
43
; RESTRICTIONS:
997bd805   Annie Hughes   comments and gene...
44
45
;    The DustEM fortran code must be installed
;    The DustEMWrap IDL code must be installed
1e005c88   Annie Hughes   inlie help
46
;
427f1205   Jean-Michel Glorian   version 4.2 merged
47
48
; PROCEDURE:
;    None
1e005c88   Annie Hughes   inlie help
49
;
427f1205   Jean-Michel Glorian   version 4.2 merged
50
51
52
; EXAMPLES
;    
; MODIFICATION HISTORY:
997bd805   Annie Hughes   comments and gene...
53
54
55
;    Written JPB Apr-2011
;    Evolution details on the DustEMWrap gitlab.
;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  
427f1205   Jean-Michel Glorian   version 4.2 merged
56
57
;-

427f1205   Jean-Michel Glorian   version 4.2 merged
58
59
60
61
62
63
64
65
IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_cc'
  sed=0.
  goto,the_end
ENDIF

;eps=1.e-30
eps=0.
997bd805   Annie Hughes   comments and gene...
66
67

;=== constants needed for CMB flux convention
427f1205   Jean-Michel Glorian   version 4.2 merged
68
69
Tcmb=2.725D0

997bd805   Annie Hughes   comments and gene...
70
;=== constants needed for MIPS flux convention
427f1205   Jean-Michel Glorian   version 4.2 merged
71
72
73
74
75
76
77
78
79
80
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)
997bd805   Annie Hughes   comments and gene...
81
82
spec2 = [eps,spec[order]]
wavein2 = [eps,wavein[order]]
427f1205   Jean-Michel Glorian   version 4.2 merged
83
84
85
86
87
88
89
90
91
92
93

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
36e8b879   Jean-Michel Glorian   update from deborah
94
95
;j=where(instrus eq 'HFI',count)
;if count gt 0 then instrus(j) = 'HFI_V3P02'
427f1205   Jean-Michel Glorian   version 4.2 merged
96
97
98
99
100
101
102
103
104

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)
b5ccb706   Jean-Philippe Bernard   improved to fit p...
105
106
107
108
  IF countinstr EQ 0 THEN BEGIN
    message,'Instrument '+instrus(i)+' not found',/continue
    stop
  ENDIF
427f1205   Jean-Michel Glorian   version 4.2 merged
109
  st=(*!dustem_filters).(iinstru)
4caf5d4a   Annie Hughes   changes for JWST ...
110
111
;  stop
  ii=(where(st.filter_names EQ filter_names[i]))[0]
997bd805   Annie Hughes   comments and gene...
112
  spec0=interpol(spec2,wavein2,st.central_wavelengths[ii])
cc7d8a09   Jean-Philippe Bernard   changed
113
114
115
116
117
  IF spec0 EQ 0 THEN BEGIN
    sed[i]=spec0
    cc[i]=1
    goto,next_filter
  ENDIF
997bd805   Annie Hughes   comments and gene...
118
  IF !dustem_do_cc NE 0 and not !dustem_never_do_cc THEN BEGIN ;Color corrections are computed
427f1205   Jean-Michel Glorian   version 4.2 merged
119
;    print,'DOING color correction calculations'
997bd805   Annie Hughes   comments and gene...
120
    CASE strupcase(fluxconvs[i]) OF
427f1205   Jean-Michel Glorian   version 4.2 merged
121
122
123
124
125
126
127
      '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)
997bd805   Annie Hughes   comments and gene...
128
        spec_int=interpol(spec2,wavein2,*(st.use_wavelengths[ii]))
427f1205   Jean-Michel Glorian   version 4.2 merged
129
130
131
132
133
134
135
136
        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)
cc7d8a09   Jean-Philippe Bernard   changed
137
        cc[i]=ccc.cc
427f1205   Jean-Michel Glorian   version 4.2 merged
138
139
140
141
142
143
144
145
146
147
148
      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)
cc7d8a09   Jean-Philippe Bernard   changed
149
        cc[i]=ccc.cc
427f1205   Jean-Michel Glorian   version 4.2 merged
150
151
      END
      'NUINU=CSTE': BEGIN
109b7cd3   Annie Hughes   added flambda=cste
152
153
154
         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)               
cc7d8a09   Jean-Philippe Bernard   changed
155
         cc[i]=num/den/spec0*(st.central_wavelengths(ii))
109b7cd3   Annie Hughes   added flambda=cste
156
157
      END
      'FLAMBDA=CSTE': BEGIN
4caf5d4a   Annie Hughes   changes for JWST ...
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
;         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
427f1205   Jean-Michel Glorian   version 4.2 merged
174
175
176
177
178
179
180
      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
cc7d8a09   Jean-Philippe Bernard   changed
181
        cc[i]=num/den/spec0/(st.central_wavelengths(ii))^3/(exp(fact/ffact)-1.)
427f1205   Jean-Michel Glorian   version 4.2 merged
182
183
184
185
186
      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)
cc7d8a09   Jean-Philippe Bernard   changed
187
        cc[i]=num/den/spec0*(st.central_wavelengths(ii))
427f1205   Jean-Michel Glorian   version 4.2 merged
188
189
190
191
192
      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)
cc7d8a09   Jean-Philippe Bernard   changed
193
        cc[i]=num/den*dustem_planck_function(Tcmb,st.central_wavelengths(ii))/spec0
427f1205   Jean-Michel Glorian   version 4.2 merged
194
195
196
197
198
199
200
     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

b5ccb706   Jean-Philippe Bernard   improved to fit p...
201
        alpha_new=0.
427f1205   Jean-Michel Glorian   version 4.2 merged
202
203
204
        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)
cc7d8a09   Jean-Philippe Bernard   changed
205
        cc[i]=num/den/spec0*(st.central_wavelengths(ii))^(-(alpha_new+1.))
427f1205   Jean-Michel Glorian   version 4.2 merged
206
207

      END
b5ccb706   Jean-Philippe Bernard   improved to fit p...
208
      ELSE:BEGIN
cc7d8a09   Jean-Philippe Bernard   changed
209
        message,strupcase(fluxconvs[i])+' not recognized',/info
b5ccb706   Jean-Philippe Bernard   improved to fit p...
210
211
212
        cc[i]=1.
        ;stop
      END
427f1205   Jean-Michel Glorian   version 4.2 merged
213
214
215
216
217
218
219
    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
9074d6cd   Jean-Philippe Bernard   fixed a bug of Na...
220
  IF finite(cc[i]) THEN sed[i]=spec0*cc[i] ELSE sed[i]=spec0
cc7d8a09   Jean-Philippe Bernard   changed
221
  next_filter:
427f1205   Jean-Michel Glorian   version 4.2 merged
222
223
224
225
226
227
ENDFOR

the_end:
RETURN,sed

END