dustem_fit_spin_test_readme.pro
10.6 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
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
PRO dustem_fit_spin_test_readme,postcript=postcript,mode=mode,help=help
;This Readme describes how to fit SEDs with a modified Black-Body
;It runs on the SED stored into the file sample_SED.xcat
;Remember that the goal here is just to illustrate the method.
;Note that its is not necessary to use the .xcat file format, and data SED can be provided
;manually, but the observation structure must have the structure shown below.
;Obviously, the dustem package must have been installed succesfully
;see DustemWrap documentation for install instructions.
;+
; NAME:
; dustem_fit_modified_bb_readme
; PURPOSE:
; This is an example of how to fit SEDs with a modified Black-Body
; using the dustem wrapper.
; It is meant to be an example to follow when writing your own
; programs using the dustem IDL wrapper.
; The SED used here is in sample_SED.xcat
; CATEGORY:
; Dustem
; CALLING SEQUENCE:
; dustem_fit_spin_test_readme,postcript=postcript,help=help
; INPUTS:
; None
; OPTIONAL INPUT PARAMETERS:
; None
; OUTPUTS:
; None
; OPTIONAL OUTPUT PARAMETERS:
; None
; ACCEPTED KEY-WORDS:
; postcript = if set plot is done in DUSTEM/Docs/Figures/Last_dustem_fit.ps
; help = If set print this help
; COMMON BLOCKS:
; None
; SIDE EFFECTS:
; None
; RESTRICTIONS:
; The dustem idl wrapper must be installed
; PROCEDURE:
; None
; EXAMPLES
; dustem_fit_spin_test_readme
; MODIFICATION HISTORY:
; Written by J.P. Bernard June 9th 2011
; see evolution details on the dustem cvs maintained at CESR
; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-
IF keyword_set(help) THEN BEGIN
doc_library,'dustem_fit_spin_test_readme'
goto,the_end
ENDIF
;=== initialise dustem
dustem_init,mode='COMPIEGNE_ETAL2010'
!dustem_verbose=1
!dustem_show_plot=1
;=== make a fake one out for given filters
filters=['DIRBE1','DIRBE2','DIRBE3','DIRBE4','IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3']
Nfilt=n_elements(filters)
sed=dustem_initialize_sed(Nfilt)
sed.filter=filters
sed.wave=dustem_filter2wav(filters)
sed.instru=dustem_filter2instru(filters)
st=dustem_set_data(sed=sed)
;1.000000
;PAH0_MC10 25 logn-chrg-spin-zm 7.8000E-04 2.2400E+00 3.5000E-08 1.2000E-07 6.4000E-08 1.0000E-01
;PAH1_MC10 25 logn-chrg-spin-zm 7.8000E-04 2.2400E+00 3.5000E-08 1.2000E-07 6.4000E-08 1.0000E-01
;amCBEx 15 logn 1.6500E-04 1.8100E+00 6.0000E-08 2.0000E-06 2.0000E-07 3.5000E-01
;amCBEx 25 plaw-ed 1.4500E-03 1.8100E+00 4.0000E-07 2.0000E-04 -2.8000E+00 1.5000E-05 1.5000E-05 2.0000E+00
;aSilx 25 plaw-ed 6.7000E-03 3.0000E+00 4.0000E-07 2.0000E-04 -3.4000E+00 2.0000E-05 2.0000E-05 2.0000E+00
goto,freefree_tag
;=== Set which parameters you want to fit
pd = [ $
; '(*!dustem_params).G0', $ ;G0
'(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
'(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
'(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx
'(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx
'(*!dustem_params).grains(4).mdust_o_mh', $ ;aSilx
'(*!dustem_params).gas.G0', $ ;Gas G0
'(*!dustem_params).gas.Tgas', $ ;Gas temperature
'(*!dustem_params).gas.nH', $ ;Gas temperature
; 'dustem_create_freefree_1', $ ;Gas temperature for free-free emission
; 'dustem_create_freefree_2', $ ;H-alpha for free-free emission (R)
'dustem_create_continuum_2'] ;Intensity of NIR continuum
;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,100.,1000.,1000.,0.001]
;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,1.e4,0.01,0.1]
p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.1]
;=== set initial value of parameters you want to fit
;iv=p_truth ;exact solution
;iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,1.e4,0.01,0.001] ;shifted
iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.001] ;shifted
ulimed=[0 , 0 ,0,0,0,0,0,0,0]
llimed=[1 , 1 ,1,1,1,1,1,1,1]
llims =[0. ,0. ,0.,0.,0.,0.,0.,0.,0.]
goto,out_of_there
freefree_tag:
;=== Set which parameters you want to fit
pd = [ $
; '(*!dustem_params).G0', $ ;G0
'(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
'(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
'(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx
'(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx
'(*!dustem_params).grains(4).mdust_o_mh', $ ;aSilx
'(*!dustem_params).gas.G0', $ ;Gas G0
'(*!dustem_params).gas.Tgas', $ ;Gas temperature
'(*!dustem_params).gas.nH', $ ;Gas temperature
'dustem_create_freefree_1', $ ;Gas temperature for free-free emission
'dustem_create_freefree_2', $ ;H-alpha for free-free emission (R)
'dustem_create_continuum_2'] ;Intensity of NIR continuum
;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,100.,1000.,1000.,0.001]
p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,1.e4,1.e-11,0.1]
;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.1]
;=== set initial value of parameters you want to fit
;iv=p_truth ;exact solution
iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,1.e4,1.e-11,0.001] ;shifted
;iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.001] ;shifted
ulimed=[0 , 0 ,0,0,0,0,0,0,0,0,0]
llimed=[1 , 1 ,1,1,1,1,1,1,1,1,1]
llims =[0. ,0. ,0.,0.,0.,0.,0.,0.,0.,0,0]
synchrotron_tag:
;=== Set which parameters you want to fit
pd = [ $
; '(*!dustem_params).G0', $ ;G0
'(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
'(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
'(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx
'(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx
'(*!dustem_params).grains(4).mdust_o_mh', $ ;aSilx
'(*!dustem_params).gas.G0', $ ;Gas G0
'(*!dustem_params).gas.Tgas', $ ;Gas temperature
'(*!dustem_params).gas.nH', $ ;Gas temperature
'dustem_create_synchrotron_1', $ ;CR spectral index for synchrotron emission
'dustem_create_synchrotron_2', $ ;Amplitude for synchrotron emission (R)
'dustem_create_continuum_2'] ;Intensity of NIR continuum
;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,100.,1000.,1000.,0.001]
p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,3.,0.01,0.1]
;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.1]
;=== set initial value of parameters you want to fit
;iv=p_truth ;exact solution
iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,3.,0.01,0.1] ;shifted
;iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.001] ;shifted
ulimed=[0 , 0 ,0,0,0,0,0,0,0,0,0]
llimed=[1 , 1 ,1,1,1,1,1,1,1,1,1]
llims =[0. ,0. ,0.,0.,0.,0.,0.,0.,0.,0,0]
out_of_there:
;== SET THE FITTED PARAMETERS (now done in dustem_init_parinfo)
dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-chrg-spin'
(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-chrg-spin'
help,*!dustem_fit
;stop
sed.spec=dustem_compute_sed(p_truth,sst,_extra=extra)
sed.error=(sed.spec*0.1);>0.1
;sed.error[*]=100.
;=== SET THE OBSERVATION STRUCTURE
;ind=where(sed.wave GE 5.)
ind=where(sed.wave GE 0.1)
st=dustem_set_data(sed=sed[ind])
;=== RUN fit
;tol=1.e-20
tol=1.e-10
xtol=1.e-10
Nitermax=5 ;maximum number of iteration. This is the criterium which will stop the fit procedure
xrange=[1.,2.e4]
yrange=[1e-4,1.e3]
loadct,13
;!y.range=[1e-8,100] ;This is to ajust plot range from outside the routine
t1=systime(0,/sec)
res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=Nitermax,xrange=xrange,/xstyle,yrange=yrange,/ysty,/ylog,/xlog,xtit='wavelength [mic]',ytit='Brightness []')
t2=systime(0,/sec)
print,(res-p_truth)/p_truth*100.
;stop
;=== SAVE FIT RESULTS
file_out='/tmp/DUSTEM_spinning_fit_example.sav'
dustem_save_system_variables,file_out
message,'Saved '+file_out,/continue
stop
errors=(*(*!dustem_fit).current_param_errors)
chi2=(*!dustem_fit).chi2
rchi2=(*!dustem_fit).rchi2
loadct,13
dustem_sed_plot,(*(*!dustem_fit).current_param_values), $
ytit=ytit,xtit=xtit,title=tit,yr=yrange,xr=xrange,/ysty,/xsty, $
res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog
;stop
;======================================
;====You can exit IDL here and re-enter
;======================================
file='/tmp/DUSTEM_spinning_fit_example.sav'
dustem_restore_system_variables,file
window,1
xrange=[1.,2.e4]
yrange=[1e-4,1.e3]
;=== recover best fit values
res=*(*!dustem_fit).current_param_values
chi2=(*!dustem_fit).chi2
rchi2=(*!dustem_fit).rchi2
;errors=*(*!dustem_fit).current_param_errors
errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values)
;=== Plot best fit
tit='DUSTEM spinning dust Example (Saved)'
ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
xtit=textoidl('\lambda (\mum)')
loadct,13
IF keyword_set(postcript) THEN BEGIN
set_plot,'PS'
ps_file=!dustem_wrap_soft_dir+'/Docs/Figures/'+'Last_dustem_spinning_fit.ps'
device,filename=ps_file,/color
ENDIF
dustem_sed_plot,*(*!dustem_fit).current_param_values,ytit=ytit,xtit=xtit,title=tit,yr=yrange,xr=xrange,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog
;dustem_sed_plot,*(*!dustem_fit).current_param_values, $
; ytit=ytit,xtit=xtit,title=tit,yr=yrange,xr=xrange,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog
IF keyword_set(postcript) THEN BEGIN
device,/close
set_plot,'X'
message,'wrote '+ps_file,/info
ENDIF
print,'dustem_mpfit_sed executed in ',t2-t1,' sec'
the_end:
END