OV-GSO-DC / dustem-wrapper_idl

Blame view

src/idl/dustem_cc.pro 7.14 KB
Edit Raw Normal View History
427f1205   Jean-Michel Glorian   version 4.2 merged 
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
 
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
36e8b879   Jean-Michel Glorian   update from deborah 
81
82
 
;j=where(instrus eq 'HFI',count)
;if count gt 0 then instrus(j) = 'HFI_V3P02'
427f1205   Jean-Michel Glorian   version 4.2 merged 
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
 

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