dustem_compute_sed.pro
5.43 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
FUNCTION dustem_compute_sed,p_dim,$
st=st,$
input_spec=input_spec, $
sed_spec=sed_spec, $
use_previous_fortran=use_previous_fortran, $
filters=filters, $
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
; filters = if set, use these filters to compute the SED (if not use filters in !dustem_data or all known filters if !dustem_data does not exist)
; 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 ptr_valid((*!dustem_fit).param_init_values) THEN BEGIN
use_param_init_values=(*(*!dustem_fit).param_init_values)
ENDIF ELSE BEGIN
use_param_init_values=p_dim
ENDELSE
;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/use_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=dustem_get_emission_conversion_factor()
SED_spec = st.sed.em_tot * fact
use_wavs = st.sed.wav
;===== Adding plugin output to spectrum ----------------
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 ELSE BEGIN
SED_spec=input_spec
use_wavs=dustem_get_wavelengths()
ENDELSE
;------------------------------------------
;This should not happen anyway
;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
IF ptr_valid((*!dustem_data).sed) THEN BEGIN
dustem_sed = (*(*!dustem_data).sed).values * 0.
use_filters=(*(*!dustem_data).sed).filt_names
use_filter_wavs=(*(*!dustem_data).sed).wav
ENDIF ELSE BEGIN
IF keyword_set(filters) THEN BEGIN
use_filters=filters
ENDIF ELSE BEGIN
use_filters=dustem_get_all_filter_names()
ENDELSE
Nfilters=n_elements(use_filters)
dustem_sed=(dustem_initialize_sed(Nfilters)).stokesI
use_filter_wavs=dustem_filter2wav(use_filters)
ENDELSE
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)
ind_sed=where(use_filters NE 'SPECTRUM' and use_wavs gt 0.,count_sed)
;ind_sed=where(use_filters NE 'SPECTRUM' and use_wavs gt 0.,count_sed,complement=ind_spec,Ncomplement=count_spec)
IF count_sed NE 0 THEN BEGIN
;filter_names=((*(*!dustem_data).sed).filt_names)[ind_sed]
filter_names=(use_filters)[ind_sed]
ssed=dustem_cc(use_wavs,SED_spec,filter_names,cc=cc)
dustem_sed[ind_sed]=ssed
ENDIF ELSE BEGIN
message,"**** WARNING **** NO FILTER defined in DATA for COLOUR CORRECTION",/continue
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)
ind_spec=where(use_filters 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]=interpol(SED_spec,use_wavs,use_filter_wavs[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:
RETURN,dustem_sed
END