dustem_compute_gb_sed.pro
3.14 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
FUNCTION dustem_compute_gb_sed,p_dim $
,_extra=extra $
,waves=waves $
,spec=spec $
,help=help
;+
; NAME:
; dustem_compute_gb_sed
; PURPOSE:
; Computes an SED for a modified blackbody
; CATEGORY:
; Dustem
; CALLING SEQUENCE:
; sed=dustem_compute_sed(p_dim[,st=][,_extra=][,/help])
; INPUTS:
; p_dim = [a0,a1,a2] modified blackbody parameter values
; a0 -- normalisation
; a1 -- dust temperature
; a2 -- emissivity
; OPTIONAL INPUT PARAMETERS:
; None
; OUTPUTS:
; sed = computed SED for filters in !dustem_data.sed
; OPTIONAL OUTPUT PARAMETERS:
; None
; ACCEPTED KEY-WORDS:
; help = if set, print this help
; COMMON BLOCKS:
; None
; SIDE EFFECTS:
; None
; RESTRICTIONS:
; The DustEMWrap IDL code must be installed
; PROCEDURE:
; None
; EXAMPLES
;
; MODIFICATION HISTORY:
; Written JP Bernard 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_compute_gb_sed'
dustem_sed=0.
goto,the_end
ENDIF
;==== Compute Black body spectrum
pp=double(p_dim)
NNw=20L
wwmin=50. & wwmax=10000.
ww=dindgen(NNw)/(1.*NNw-1)*(alog10(wwmax)-alog10(wwmin))+alog10(wwmin)
ww=10.^ww
waves=ww
;MAKE A COPY OF THE (EMPTY) OBSERVED VALUES
obs_sed = (*!dustem_data.sed).values
;COMPUTE THE MODEL SED
dustem_sed = obs_sed*0.D0
ind_sed=where((*!dustem_data.sed).filt_names NE 'SPECTRUM',count_sed)
IF !dustem_do_cc NE 0 AND !dustem_never_do_cc EQ 0 THEN BEGIN
message,'DOING color correction calculations',/info
ENDIF ELSE BEGIN
message,'SKIPPING color correction calculations',/info
ENDELSE
;stop
IF !dustem_do_cc NE 0 AND !dustem_never_do_cc EQ 0 THEN BEGIN
; message,'DOING color correction calculations',/info
IF count_sed NE 0 THEN BEGIN
spec=pp(0)*(ww)^(-1.*pp(2))*dustem_planck_function(pp(1),ww)
; filter_names=((*!dustem_data).filt_names)(ind_sed)
; FOR ii=0L,n_elements(filter_names)-1 DO BEGIN
ssed=dustem_cc(ww,spec,((*!dustem_data.sed).filt_names)(ind_sed),cc=cc)
;print,cc
; ENDFOR
dustem_sed[ind_sed]=ssed
ENDIF
ENDIF ELSE BEGIN
www=dustem_filter2wav(((*!dustem_data.sed).filt_names)[ind_sed])
dustem_sed[ind_sed]=pp(0)*(www)^(-1.*pp(2))*dustem_planck_function(pp(1),www)*(*!dustem_previous_cc)
message,'SKIPPING color correction calculations',/info
; dustem_sed(ind_sed)=interpol(spec,ww,1.D0*(((*!dustem_data).wav)(ind_sed)))*(*!dustem_previous_cc)
spec= dustem_sed[ind_sed]
ENDELSE
;stop
;For spectrum data points, interpolate in log-log
;Linear interpolation leads to wrong values, in particular where few
;wavelengths points exist in the model (long wavelengths).
ind_spec=where((*!dustem_data.sed).filt_names EQ 'SPECTRUM',count_spec)
IF count_spec NE 0 THEN BEGIN
dustem_sed[ind_spec]=interpol(spec,ww,(((*!dustem_data.sed).wav)[ind_spec]))
; dustem_sed(ind_spec)=10^interpol(alog10(spec),alog10(st.sed.wav),alog10((((*!dustem_data.sed).wav)(ind_spec))))
ENDIF
;clean pointers
heap_gc
the_end:
RETURN,dustem_sed
END