dustem_compute_sed.pro
4.67 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
FUNCTION dustem_compute_sed,p_dim,$
st=st,$
input_spec=input_spec, $
sed_spec=sed_spec, $
use_previous_fortran=use_previous_fortran, $
help=help
;+
; NAME:
; dustem_compute_sed
;
; PURPOSE:
; Computes a SED for a given Dustem model
;
; CATEGORY:
; Dustem
;
; CALLING SEQUENCE:
; sed=dustem_compute_sed(p_dim[,st][,input_spec][,sed_spec][,/help])
;
; INPUTS:
; p_dim = parameter values. This is actually optional, and
; will not be used if st or input_spec are provided, but
; leaving like this for backwards compatbility.
;
; OPTIONAL INPUT PARAMETERS:
; st = Dustem Fortran output structure. Takes precedence
; over p_dim. NB THIS IS ALSO THE OUTPUT VARIABLE.
; input_spec = input spectrum that will be used to calculate
; SED. Expected units are MJy/sr. Takes precedence
; over p_dim and st
;
; OUTPUTS:
; sed = computed SED for filters in !dustem_data
;
; OPTIONAL OUTPUT PARAMETERS:
; st = Dustem fortran output structure NB THIS IS ALSO AN INPUT VARIABLE.
; sed_spec = final emission spectrum
;
; ACCEPTED KEY-WORDS:
; use_previous_fortran = if set, uses the output of the previous fortran run, and does not run the fortran.
; 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 by J.-Ph. Bernard 2008
; 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_sed'
dustem_sed=0.
goto,the_end
ENDIF
;stop
IF not keyword_set(st) and not keyword_set(input_spec) THEN BEGIN
;==== JPB: This division leads to NaN for initial values at 0 !!
dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st=st,use_previous_fortran=use_previous_fortran
ENDIF
;stop
IF not keyword_set(input_spec) THEN BEGIN
; Convert into MJy/sr/1d20
fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/st.sed.wav)*1.e20/1.e7 ;this is correct despite the way it is presented (1.e4/1e.7*1e.20=1.e17 which is the factor to convert from ergs to Mega Janskys)
SED_spec = st.sed.em_tot * fact
use_wavs = st.sed.wav
;===== Adding plugin output to spectrum ----------------
;Nplugins=n_elements(*!dustem_plugin)
;scopes=(*!dustem_plugin).scope
plugin_tags=(*!dustem_plugin).name
ind=where(plugin_tags NE 'NONE',Nplugins)
IF Nplugins NE 0 THEN BEGIN
FOR i=0L,Nplugins-1 DO BEGIN
IF total(strsplit((*!dustem_plugin)[i].scope,'+',/extract) eq 'ADD_SED') THEN BEGIN
SED_spec+=(*(*!dustem_plugin)[i].spec)[*,0]
ENDIF
ENDFOR
ENDIF
ENDIF
if keyword_set(input_spec) then begin
SED_spec=input_spec
use_wavs=dustem_get_wavelengths()
endif
;------------------------------------------
;stop
IF not isarray(SED_spec) THEN stop
;COMPUTE THE MODEL SED
;JPB: problem there: we should be able to compute an SED even when there is no data in !dustem_data
dustem_sed = (*(*!dustem_data).sed).values * 0.
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: we take previous values of cc',/info
ENDELSE
ind_sed=where((*(*!dustem_data).sed).filt_names NE 'SPECTRUM' and (*(*!dustem_data).sed).wav gt 0.,count_sed)
IF count_sed NE 0 THEN BEGIN
filter_names=((*(*!dustem_data).sed).filt_names)[ind_sed]
; ssed=dustem_cc(st.sed.wav,SED_spec,filter_names,cc=cc)
ssed=dustem_cc(use_wavs,SED_spec,filter_names,cc=cc)
dustem_sed[ind_sed]=ssed
ENDIF ELSE BEGIN
print,"**** WARNING **** NO FILTER defined in DATA for COLOUR CORRECTION"
ENDELSE
;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(SED_spec,st.sed.wav,(((*(*!dustem_data).sed).wav)[ind_spec]))
dustem_sed[ind_spec]=interpol(SED_spec,use_wavs,(((*(*!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
;out_st=st
;clean pointers
heap_gc
the_end:
;VG : error in cc IRAS1 when no PAH
; j=where(dustem_sed ne dustem_sed,count)
; if count gt 0 then dustem_sed[j]=0.
;stop
RETURN,dustem_sed
END