Blame view

src/idl/dustem_make_sed_table.pro 7.82 KB
2b6fd350   Jean-Philippe Bernard   added pre-computi...
1
2
3
4
5
6
7
PRO dustem_make_sed_table,model, $
											    parameters_description, $
											    iv_min, $
											    iv_max, $
											    iv_Nvalues, $
						  						fpd=fpd, $
						  						fiv=fiv, $
4ddd208e   Jean-Philippe Bernard   chnaged to store ...
8
						  						grain_keywords=grain_keywords, $
2b6fd350   Jean-Philippe Bernard   added pre-computi...
9
10
						  						filename=filename, $
						  						filters=filters, $
e2ce1dee   Jean-Philippe Bernard   defined grid stru...
11
						  						plog=plog, $
2b6fd350   Jean-Philippe Bernard   added pre-computi...
12
						  						show_seds=show_seds, $
54583f40   Jean-Philippe Bernard   improved
13
						  						print_params=print_params, $
2b6fd350   Jean-Philippe Bernard   added pre-computi...
14
						  						help=help
eb263d61   Jean-Philippe Bernard   commented
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

;+
; NAME:
;    dustem_make_sed_table
; PURPOSE:
;    makes an SED lookup table using dustemwrap for a given dust model and free parameters
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    dustem_make_sed_table,model,parameters_description,parvalues_min,parvalues_max,parvalues_Nvalues, $
;                     	   [,fpd=][,fiv=][,filename=][,filters=][,log=][,/show_seds]
; INPUTS:
;    model                  : dustemwrap model name to be used
;    parameters_description : dustemwrap parameter description array
;    parvalues_min          : minimum values for each parameter
;    parvalues_max          : maximum values for each parameter
2b6fd350   Jean-Philippe Bernard   added pre-computi...
31
;    parvalues_Nvalues      : number of parameter values for each parameter
eb263d61   Jean-Philippe Bernard   commented
32
; OPTIONAL INPUT PARAMETERS:
54583f40   Jean-Philippe Bernard   improved
33
;    filters                : name of dustemwrap filters to be included in the SED calculations (default = ALL known filters)
eb263d61   Jean-Philippe Bernard   commented
34
35
;    fpd                    : fixed parameter description
;    fiv                    : fixed parameter values
eb263d61   Jean-Philippe Bernard   commented
36
;    filename               : output fits file name for the table (default ='/tmp/dustem_seds.fits')
e2ce1dee   Jean-Philippe Bernard   defined grid stru...
37
;    plog                   : array indicating if a given parameter is to be sampled in linear (0) or log scale (1)
eb263d61   Jean-Philippe Bernard   commented
38
39
40
41
42
43
;    show_seds              : if set, plots SEDs as they are computed.
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
eb263d61   Jean-Philippe Bernard   commented
44
45
46
47
48
49
50
51
52
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    A file is written
; RESTRICTIONS:
;    The DustEM fortran code must be installed
;    The DustEMWrap IDL code must be installed
; PROCEDURE:
54583f40   Jean-Philippe Bernard   improved
53
;    Filters in the table are ordered by increasing wavelength
eb263d61   Jean-Philippe Bernard   commented
54
; EXAMPLES
2b6fd350   Jean-Philippe Bernard   added pre-computi...
55
;    dustem_make_sed_table,'DBP90',['(*!dustem_params).G0','(*!dustem_params).grains(0).mdust_o_mh'], $
e2ce1dee   Jean-Philippe Bernard   defined grid stru...
56
;																[0.1,1.e-3],[100.,1.e-1],[10,3],filename='/tmp/IRAS_GO.fits',plog=[1,0],/show_seds, $
2b6fd350   Jean-Philippe Bernard   added pre-computi...
57
;                              	filters=[dustem_instru2filters('IRAC'),dustem_instru2filters('SPIRE')]
eb263d61   Jean-Philippe Bernard   commented
58
59
60
61
62
63
64
65
66
67
; MODIFICATION HISTORY:
;    Written by J.-Ph. Bernard (2023)
;    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_make_sed_table'
  goto,the_end
ENDIF
6ff7d98f   Jean-Philippe Bernard   First commit
68

6ff7d98f   Jean-Philippe Bernard   First commit
69
70
71
72
Nparams=n_elements(parameters_description)

;== INITIALISE DUSTEM
use_polarization=0
4ddd208e   Jean-Philippe Bernard   chnaged to store ...
73
dustem_init,model=model,pol=use_polarisation,grain_keywords=grain_keywords
0297394d   Jean-Philippe Bernard   test
74
!quiet=1
f28cb4f4   Jean-Philippe Bernard   improved
75
76
;!dustem_verbose=1
;!dustem_show_plot=1
6ff7d98f   Jean-Philippe Bernard   First commit
77
78
!dustem_which='RELEASE'

429a6075   Jean-Philippe Bernard   improved
79
;=== define default filters
e2ce1dee   Jean-Philippe Bernard   defined grid stru...
80
use_filters=dustem_get_all_filter_names(all_instrus=all_instrus)
6ff7d98f   Jean-Philippe Bernard   First commit
81
IF keyword_set(filters) THEN use_filters=filters
e2ce1dee   Jean-Philippe Bernard   defined grid stru...
82
83
84
85
86
87
88
89
ind=where(use_filters NE 'SPECTRUM',count)
IF count NE 0 THEN BEGIN
	use_filters=use_filters[ind]
ENDIF ELSE BEGIN
	message,'No usable filter found',/info
	stop
ENDELSE

eb263d61   Jean-Philippe Bernard   commented
90
91
92
wavs=dustem_filter2wav(use_filters)
order=sort(wavs)
use_filters=use_filters[order]
6ff7d98f   Jean-Philippe Bernard   First commit
93
94
95
Nfilters=n_elements(use_filters)

;==== compute all combinations of model parameter values
e2ce1dee   Jean-Philippe Bernard   defined grid stru...
96
parameter_values=dustem_param_range2param_values(iv_min,iv_max,iv_Nvalues,Nc=Nc,log=plog)
6ff7d98f   Jean-Philippe Bernard   First commit
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

;=== initialize filters for the SED
Nfilt=n_elements(use_filters)
;stop
sed=dustem_initialize_sed(Nfilt)

sed.filter=use_filters
sed.wave=dustem_filter2wav(use_filters)
sed.instru=dustem_filter2instru(use_filters)

;=== initialize IQU and associated errors to avoid problems when checking SED in dustem_set_data.pro
sed[*].StokesI=1.e-10
sed.StokesQ=sed.StokesI/100.
sed.StokesU=sed.StokesI/100.
sed.SigmaII=sed.StokesI/100.
sed.SigmaQQ=sed.StokesI/100.
sed.SigmaUU=sed.StokesI/100.
sed.SigmaIQ=sed.StokesI/100.
sed.SigmaIU=sed.StokesI/100.
sed.SigmaQU=sed.StokesI/100.

;=== initialize wavelengths for the EXTINCTION
use_next=100 ; number of extinction measurements
use_wave=10^(alog10(0.01)+alog10(50)*findgen(use_next)/float(use_next)) ; wavelengths at which extinction is defined
if keyword_set(next) then use_next=next
if keyword_set(wave) then $
   use_wave=range_gen(use_next,[wave[0],wave[1]],/log)

ext=dustem_initialize_ext(use_next)
ext.instru='EXTINCTION'
ext.filter='SPECTRUM'
ext.wave=use_wave

;=== initialize IQU and associated errors to avoid problems when checking EXT in dustem_set_data.pro
ext[*].EXT_I=1.e-10
ext.EXT_Q=ext.EXT_I/100.
ext.EXT_U=ext.EXT_I/100.
ext.SIGEXTII=ext.EXT_I/1000.
ext.SIGEXTQQ=ext.EXT_I/1000.
ext.SIGEXTUU=ext.EXT_I/1000.
ext.SIGEXTIQ=ext.EXT_I/1000.
ext.SIGEXTIU=ext.EXT_I/1000.
ext.SIGEXTQU=ext.EXT_I/1000.

dustem_set_data,m_fit=sed,m_show=sed,x_fit=ext,x_show=ext

;=== initialize at least one parameter as 'free'
;=== this is only to keep dustemwrap happy in its management of variables
;=== no fitting is going to be done
25370043   Jean-Philippe Bernard   src/idl/dustem_ac...
146
147
148
;pd = ['(*!dustem_params).gas.G0'] ; G0
;pval = [-1.e-20]
pd=parameters_description
6ff7d98f   Jean-Philippe Bernard   First commit
149
150
151
152
153

;=== declare some other parameters as fixed
;=== this shows how to add any non-dust plugins to the observed spectra/SED
;=== for the dust-model parameters, it is not necessary, but allows the
;=== user to (i) change default values and (2) see the adopted values in the graphical output windows
25370043   Jean-Philippe Bernard   src/idl/dustem_ac...
154
;fpd=parameters_description
6ff7d98f   Jean-Philippe Bernard   First commit
155
seds=ptrarr(Nc)
e2ce1dee   Jean-Philippe Bernard   defined grid stru...
156
157
em_spectra=ptrarr(Nc)
ext_spectra=ptrarr(Nc)
6ff7d98f   Jean-Philippe Bernard   First commit
158
159
yrange=[1.e-5,1.e5]
colors=long(range_gen(Nc,[0,255]))
14c90b64   Jean-Philippe Bernard   removed loadct co...
160
IF keyword_set(show_seds) THEN loadct,13
6ff7d98f   Jean-Philippe Bernard   First commit
161
162
163
FOR i=0L,Nc-1 DO BEGIN
	;dustem_init,model=model,pol=use_polarisation
	;dustem_set_data,m_fit=sed,m_show=sed,x_fit=ext,x_show=ext
e2ce1dee   Jean-Philippe Bernard   defined grid stru...
164
  st=0 ;otherwise model is not recomputed !!
6ff7d98f   Jean-Philippe Bernard   First commit
165

25370043   Jean-Philippe Bernard   src/idl/dustem_ac...
166
167
	;fpval = *parameter_values[i]
  pval = *parameter_values[i]
f28cb4f4   Jean-Philippe Bernard   improved
168
	message,'= computing model '+strtrim(i,2)+'/'+strtrim(Nc,2)+' ('+strtrim(1.*i/Nc*100.)+'%)',/continue
c2457c07   Jean-Philippe Bernard   improved
169
170
	message,'free parameter values=',/continue
	print,'free_parameter_values=',pval
ea1acf9c   Jean-Philippe Bernard   improved
171
172
173
174
	IF keyword_set(fiv) THEN BEGIN
		message,'fixed parameter values=',/continue
		print,'fixed_parameter_values=',fiv
	ENDIF
dc10161d   Jean-Philippe Bernard   unstopped
175
  ;stop
6ff7d98f   Jean-Philippe Bernard   First commit
176
	;=== initialise all this information into dustemwrap's brain
d339db1c   Jean-Philippe Bernard   improved
177
	dustem_init_params,model,pd,pval,fpd=fpd,fiv=fiv,pol=use_polarisation,/force_reset
54583f40   Jean-Philippe Bernard   improved
178

e2ce1dee   Jean-Philippe Bernard   defined grid stru...
179
180
181
182
  ;IF keyword_set(print_params) THEN BEGIN
 ; 	dustem_print_params,0
 ; 	;stop
 ; ENDIF
6ff7d98f   Jean-Philippe Bernard   First commit
183
184

	;=== compute the predicted SED and extinction curve for the dust-model
e2ce1dee   Jean-Philippe Bernard   defined grid stru...
185
186
187
	;=== and any plugins
	dustem_Ised=dustem_compute_sed(pval,st=st)   ;is in MJy/sr/1e20 H/cm2
	
6ff7d98f   Jean-Philippe Bernard   First commit
188
189
	dustem_Qsed=dustem_Ised*0.
	dustem_Used=dustem_Ised*0.
6ff7d98f   Jean-Philippe Bernard   First commit
190
191

	;== Now we have a predicted SED (I), we fill up a structure that we
6ff7d98f   Jean-Philippe Bernard   First commit
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
	sed.stokesI=dustem_Ised
	sed.stokesQ=dustem_Qsed
	sed.stokesU=dustem_Used
	;=== set uncertainties for I,Q,U to 0.
	sed.sigmaII=0.
	sed.sigmaQQ=0.
	sed.sigmaUU=0.
	;=== set covariances to 0.
	sed.sigmaIQ=0.
	sed.sigmaIU=0.
	sed.sigmaQU=0.

	;==== fill in dependent columns of the SED.
	sed=dustem_fill_sed_dependent_columns(sed)
	seds[i]=ptr_new(sed)
e2ce1dee   Jean-Philippe Bernard   defined grid stru...
207
208
	em_spectra[i]=ptr_new(st.sed)
	ext_spectra[i]=ptr_new(st.ext)
25370043   Jean-Philippe Bernard   src/idl/dustem_ac...
209
  IF keyword_set(show_seds) THEN BEGIN
e2ce1dee   Jean-Philippe Bernard   defined grid stru...
210
  	  xrange=[0.1,1.e3]
6ff7d98f   Jean-Philippe Bernard   First commit
211
	    IF i EQ 0 THEN BEGIN
25370043   Jean-Philippe Bernard   src/idl/dustem_ac...
212
213
214
		    tit='Grid SEDs model '+model
		    xtit='Wavelength [mic]'
		    ytit='SED [MJy/sr for NH=1.e20 H/cm^2]'
e2ce1dee   Jean-Philippe Bernard   defined grid stru...
215
216
217
		    wavs=(*seds[i]).wave
		    Is=(*seds[i]).stokesI
	    	cgplot,(*seds[i]).wave,(*seds[i]).stokesI,psym=-4,/ylog,yr=yrange,xtit=xtit,ytit=ytit,tit=tit,/xlog,xrange=xrange,/xsty
6ff7d98f   Jean-Philippe Bernard   First commit
218
	    ENDIF ELSE BEGIN
eb263d61   Jean-Philippe Bernard   commented
219
	    	cgoplot,(*seds[i]).wave,(*seds[i]).stokesI,psym=-4,color=colors[i]
6ff7d98f   Jean-Philippe Bernard   First commit
220
	    ENDELSE
e2ce1dee   Jean-Philippe Bernard   defined grid stru...
221
	    ;stop
6ff7d98f   Jean-Philippe Bernard   First commit
222
223
224
	ENDIF
	;stop
ENDFOR
54583f40   Jean-Philippe Bernard   improved
225
;stop
6ff7d98f   Jean-Philippe Bernard   First commit
226

e2ce1dee   Jean-Philippe Bernard   defined grid stru...
227
dustem_write_sed_grid,model,parameters_description,iv_min,iv_max,iv_Nvalues,plog $
907a3ade   Jean-Philippe Bernard   improved
228
229
										 ,parameter_values,seds,em_spectra=em_spectra,ext_spectra=ext_spectra,fpd=fpd,fiv=fiv $
										 ,filename=filename,use_polarisation=use_polarisation,grain_keywords=grain_keywords
6ff7d98f   Jean-Philippe Bernard   First commit
230

6ff7d98f   Jean-Philippe Bernard   First commit
231

e2ce1dee   Jean-Philippe Bernard   defined grid stru...
232
;toto=mrdfits(fits_file,1,header)
6ff7d98f   Jean-Philippe Bernard   First commit
233
234
235
;hprint,header

;stop
eb263d61   Jean-Philippe Bernard   commented
236
the_end:
6ff7d98f   Jean-Philippe Bernard   First commit
237
238

END