dustem_mpfit_run.pro
16.7 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
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
FUNCTION dustem_mpfit_run,x,p_min,niter=niter_completed,_EXTRA=_extra
;+
; NAME:
; dustem_mpfit_run
;
; PURPOSE:
; function used to fit dustem model with dustem_mpfit_sed
;
; CATEGORY:
; Dustem
;
; CALLING SEQUENCE:
; res=dustem_mpfit_run_vg(x,params[,_EXTRA=extra][,cc_sed=cc_sed][,/help])
;
; INPUTS:
; x = necessary for interface but unused
; params = Dustem Model parameters normalized to initial values
;
; OPTIONAL INPUT PARAMETERS:
; cc_sed = SED color correction factors
; _extra = extra parameters for the plot routine
;
; OUTPUTS:
; None
;
; OPTIONAL OUTPUT PARAMETERS:
; None
;
; ACCEPTED KEY-WORDS:
; help = If set, print this help
;
; COMMON BLOCKS:
; None
;
; SIDE EFFECTS:
; SED and model are plotted if !dustem_show_plot
;
; RESTRICTIONS:
; The dustem idl wrapper must be installed
;
; PROCEDURE:
; None
;
; EXAMPLES
; res=dustem_mpfit_run_vg(x,params[,_EXTRA=extra][,cc_sed=cc_sed][,/help])
;
; MODIFICATION HISTORY:
; Adapted to extinction and polarization by Vincent Guillet (IAS)
; Written by J.-Ph. Bernard
; Evolution details on the DustEMWrap gitlab.
; See http://dustemwrap.irap.omp.eu/ for FAQ and help.
;-
;note that x is not used.
IF keyword_set(help) THEN BEGIN
doc_library,'dustem_mpfit_run'
sed=0.
goto,the_end
ENDIF
p_dim=p_min*(*(*!dustem_fit).param_init_values)
IF !run_ionfrac gt 0. THEN toto = dustem_create_ionfrac()
;RUN THE MODEL AND GET THE SPECTRUM
dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st=st ;st is an output...
st_model=dustem_read_all(!dustem_soft_dir)
chi2_sed = 0
chi2_ext = 0
chi2_qsed =0
chi2_used =0
chi2_qext =0
chi2_uext =0
rchi2_sed = 0
rchi2_ext = 0
rchi2_qsed =0
rchi2_used =0
rchi2_qext =0
rchi2_uext =0
wrchi2_sed = 0
wrchi2_ext = 0
wrchi2_qsed =0
wrich2_used =0
wrchi2_qext =0
wrich2_uext =0
n_sed = 0
n_ext = 0
n_qsed =0
n_used =0
n_qext =0
n_uext =0
win = 0
tagnames=tag_names((*!dustem_data))
dustem_out = [0]
ind_data_tag=dustem_data_tag(count=count_data_tag)
FOR i_tag=0,count_data_tag-1 DO BEGIN
CASE tagnames(ind_data_tag(i_tag)) OF
'SED' : BEGIN
dustem_sed = dustem_compute_sed(p_dim,st=st,SED_spec=SED_spec)
;help,*(*!dustem_fit).CURRENT_PARAM_VALUES
;stop
chi2_sed = total(((*(*!dustem_data).sed).values-dustem_sed)^2/(*(*!dustem_data).sed).sigma^2)
n_sed = n_elements((*(*!dustem_data).sed).values)
rchi2_sed = chi2_sed / n_sed
wrchi2_sed = rchi2_sed * (*!fit_rchi2_weight).sed
message,"chi2 SED = "+strtrim(chi2_sed,2)+" rchi2 SED = "+strtrim(rchi2_sed,2),/continue;," weighted rchi2 SED=",wrchi2_sed
;print,'wl ',(*(*!dustem_data).sed).wav
;print,"Delta SED per wl: ",(*(*!dustem_data).sed).values-dustem_sed
;print,"error per wl: ",(*(*!dustem_data).sed).sigma
;print,"chi2 per wl: ",((*(*!dustem_data).sed).values-dustem_sed)^2/(*(*!dustem_data).sed).sigma^2
alpha = total(((*(*!dustem_data).sed).values*dustem_sed)/(*(*!dustem_data).sed).sigma^2) $
/ total(((*(*!dustem_data).sed).values)^2/(*(*!dustem_data).sed).sigma^2)
beta = total(((*(*!dustem_data).sed).values*dustem_sed)/(*(*!dustem_data).sed).sigma^2) $
/ total(dustem_sed^2/(*(*!dustem_data).sed).sigma^2)
print,"Best SED chi2 obtained if f_HI = ",alpha, " or if we multiply model by ",beta
;These two lines were put here so that the chi2 is updated everytime we plot
;Moving it to the bottom so that all chi2 are taken into account
; (*!dustem_fit).chi2=chi2_sed
; (*!dustem_fit).rchi2=rchi2_sed
dustem_out = [ dustem_out, dustem_sed ]
END
'EXT' : BEGIN
dustem_ext = dustem_compute_ext(p_dim,st=st,EXT_spec=EXT_spec)
chi2_ext = total(((*(*!dustem_data).ext).values-dustem_ext)^2/(*(*!dustem_data).ext).sigma^2)
n_ext = n_elements((*(*!dustem_data).ext).values)
rchi2_ext = chi2_ext / n_ext
wrchi2_ext = rchi2_ext * (*!fit_rchi2_weight).ext
print,"chi2 EXT = ",chi2_ext," rchi2 EXT = ",rchi2_ext;," weighted rchi2 EXT=",wrchi2_ext
;print,"chi2 per wl: ",((*(*!dustem_data).ext).values-dustem_ext)^2/(*(*!dustem_data).ext).sigma^2
dustem_out = [ dustem_out, dustem_ext ]
END
'POLEXT': BEGIN
;Haven't investigated !run_circ=1 yet...
;---------------the procedure within this block has not been updated for a vert long time
IF !dustem_show_plot NE 0 and !run_circ then begin
win+=1
dustem_plot_circ,st,win=win,xr=[0.1,5],aligned=st_model.align.grains.aligned,yr=[-0.03,0.03]
endif
;--------------------------------------------------------------------------------------
IF !run_lin then begin
dustem_polext = dustem_compute_polext(p_dim,st=st,POLEXT_spec=POLEXT_spec,SPEXT_spec=SPEXT_spec,dustem_fpolext=dustem_fpolext,dustem_ext=dustem_ext)
ENDIF
END
'QSED': BEGIN
stokes = dustem_compute_stokes(p_dim,st=st,Q_spec=Q_spec,U_spec=U_spec,PSI_spec=PSI_spec,dustem_psi_em=dustem_psi_em)
dustem_qsed = stokes[0]
;taking care of the color correction for spass ; SHOULD NOT LEAVE THIS ON
;dustem_qsed(-1)=dustem_qsed(-2)
; testspass1 = ((*(*!dustem_data).qsed).filt_names)(-1) eq 'SPASS1'
; if testspass1 then dustem_qsed(-1)=dustem_qsed(-2)
;
chi2_qsed = total(((*(*!dustem_data).qsed).values-dustem_qsed)^2/(*(*!dustem_data).qsed).sigma^2)
n_qsed = n_elements((*(*!dustem_data).qsed).values)
rchi2_qsed = chi2_qsed / n_qsed
wrchi2_qsed = rchi2_qsed * (*!fit_rchi2_weight).qsed
dustem_out = [ dustem_out, dustem_qsed ] ; classic command
message,"chi2 QSED = "+strtrim(chi2_qsed,2)+" rchi2 QSED = "+strtrim(rchi2_qsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
print,"chi2 per wl: ",((*(*!dustem_data).qsed).values-dustem_qsed)^2/(*(*!dustem_data).qsed).sigma^2
END
'USED': BEGIN
dustem_used = stokes[1]
chi2_used = total(((*(*!dustem_data).used).values-dustem_used)^2/(*(*!dustem_data).used).sigma^2)
n_used = n_elements((*(*!dustem_data).used).values)
rchi2_used = chi2_used / n_used
wrchi2_used = rchi2_used * (*!fit_rchi2_weight).used
dustem_out = [ dustem_out, dustem_used] ; classic command
message,"chi2 USED = "+strtrim(chi2_used,2)+" rchi2 USED = "+strtrim(rchi2_used,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
print,"chi2 per wl: ",((*(*!dustem_data).used).values-dustem_used)^2/(*(*!dustem_data).used).sigma^2
END
'POLSED': BEGIN
IF !run_pol and !run_lin THEN BEGIN
;the keyword for dustem_sed here is used so that the function uses it to compute the polarization fraction instead of
;running dustem_compute_sed again in the procedure. Which it used to do.
dustem_polsed = dustem_compute_polsed(p_dim,st=st,P_spec=P_spec,SP_spec=SP_spec,dustem_polfrac=dustem_polfrac,dustem_sed=dustem_sed)
ENDIF
END
'QEXT': BEGIN ;UPDATED
;stop
Stokext = dustem_compute_stokext(p_dim,st=st,QEXT_spec=QEXT_spec,UEXT_spec=UEXT_spec,PSIEXT_spec=PSIEXT_spec,dustem_psi_ext=dustem_psi_ext)
dustem_qext = stokext[0]
chi2_qext =total(((*(*!dustem_data).qext).values-dustem_qext)^2/(*(*!dustem_data).qext).sigma^2)
n_qext = n_elements((*(*!dustem_data).qext).values)
rchi2_qext = chi2_qext / n_qext
wrchi2_qext = rchi2_qext * (*!fit_rchi2_weight).qext
dustem_out = [ dustem_out, dustem_qext] ; classic command
message,"chi2 QEXT = "+strtrim(chi2_qext,2)+" rchi2 QEXT = "+strtrim(rchi2_qext,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
print,"chi2 per wl: ",((*(*!dustem_data).qext).values-dustem_qext)^2/(*(*!dustem_data).qext).sigma^2
END
;
'UEXT': BEGIN ;UPDATED
dustem_uext = stokext[1]
chi2_uext =total(((*(*!dustem_data).uext).values-dustem_uext)^2/(*(*!dustem_data).uext).sigma^2)
n_uext = n_elements((*(*!dustem_data).uext).values)
rchi2_uext = chi2_uext / n_uext
wrchi2_uext = rchi2_uext * (*!fit_rchi2_weight).uext
dustem_out = [ dustem_out, dustem_uext] ; classic command
message,"chi2 UEXT = "+strtrim(chi2_uext,2)+" rchi2 UEXT = "+strtrim(rchi2_uext,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
print,"chi2 per wl: ",((*(*!dustem_data).uext).values-dustem_uext)^2/(*(*!dustem_data).uext).sigma^2
END
ELSE : ;Do nothing
ENDCASE
ENDFOR
;=== computing chi^2
DOF = n_elements(p_min)
total_chi2 = 0
total_rchi2 = 0
total_n = 0
if (*!fit_rchi2_weight).ext ne 0 then begin
total_n += n_ext
total_chi2 += chi2_ext
total_rchi2 += rchi2_ext
endif
if (*!fit_rchi2_weight).sed ne 0 then begin
total_n += n_sed
total_chi2 += chi2_sed
total_rchi2 += rchi2_sed
endif
if !run_pol then begin
if tag_exist(*!fit_rchi2_weight,'qsed') then begin ;takes into account the activation$
;of the !run_pol keyword without the user knowing it.
;ie: the use of a model with polarization and without polarization data or the user wanting to fit polarization.
if (*!fit_rchi2_weight).qsed ne 0 then begin
total_n += n_qsed
total_chi2 += chi2_qsed
total_rchi2 += rchi2_qsed
endif
if (*!fit_rchi2_weight).used ne 0 then begin
total_n += n_used
total_chi2 += chi2_used
total_rchi2 += rchi2_used
endif
endif
if tag_exist(*!fit_rchi2_weight,'qext') then begin
if (*!fit_rchi2_weight).qext ne 0 then begin
total_n += n_qext
total_chi2 += chi2_qext
total_rchi2 += rchi2_qext
endif
if (*!fit_rchi2_weight).uext ne 0 then begin
total_n += n_uext
total_chi2 += chi2_uext
total_rchi2 += rchi2_uext
endif
endif
endif
total_wchi2 = chi2_ext * (*!fit_rchi2_weight).ext $
+ chi2_sed * (*!fit_rchi2_weight).sed
if !run_pol then begin
if tag_exist(*!fit_rchi2_weight,'qsed') THEN total_wchi2+=chi2_qsed * (*!fit_rchi2_weight).qsed $
+chi2_used * (*!fit_rchi2_weight).used
if tag_exist(*!fit_rchi2_weight,'qext') THEN total_wchi2+=chi2_qext * (*!fit_rchi2_weight).qext $
+chi2_used * (*!fit_rchi2_weight).used
; total_wchi2 += $;chi2_polext * !fit_rchi2_weight.polext $
; ; + chi2_polsed * !fit_rchi2_weight.polsed $
; ; + chi2_polfrac * !fit_rchi2_weight.polfrac $
; + chi2_qsed * (*!fit_rchi2_weight).qsed $
; + chi2_used * (*!fit_rchi2_weight).used $
; + chi2_qext * (*!fit_rchi2_weight).qext $
; + chi2_uext * (*!fit_rchi2_weight).uext
endif
;Using this expression for now for the weighting.
(*!dustem_fit).chi2 = total_wchi2
(*!dustem_fit).rchi2 = total_wchi2/(total_n-DOF)
message,'======================================================================',/continue
print,"Total weighted chi2 driving mpfitfun = ", total_wchi2
print,"Nb of data points = ", total_n
print,"DOF = ", DOF
print,"Reduced chi2 EXT/SED/POLEXT/POLSED/POLFRAC",rchi2_ext,rchi2_sed;,rchi2_polext;,rchi2_polsed,rchi2_polfrac
print,"Weighted reduced chi2 EXT/SED/POLEXT/POLSED/POLFRAC",wrchi2_ext,wrchi2_sed;,wrchi2_polext;,wrchi2_polsed,wrchi2_polfrac
if !run_pol then print,"Weights",(*!fit_rchi2_weight).ext,(*!fit_rchi2_weight).sed $;,!fit_rchi2_weight.polext $;,!fit_rchi2_weight.polsed,!fit_rchi2_weight.polfrac $
else print,"Weights",(*!fit_rchi2_weight).ext,(*!fit_rchi2_weight).sed
print,"Mean weighted reduced chi2 = ", total_wchi2/(total_n-DOF)
print,"Unweighted reduced Total chi2 = ", total_chi2/(total_n-DOF)
print,"Unweighted mean reduced chi2 = ", total_rchi2/4*total_n/(total_n-DOF);total_rchi2/4*total_n/(total_n-DOF) original ; I do not understand this line. Maybe redo some computation.
message,'======================================================================',/continue
;the following is needed under fawlty (not IDL) which message,/reset does not work like in IDL
!Error_State.msg=''
If !dustem_show_plot THEN BEGIN
;=== plotting the fit
if !dustem_iter.act NE !dustem_iter.prv or !dustem_iter.prv EQ 1 then begin
IF !dustem_noobj EQ 0 THEN BEGIN
dustemwrap_plot, p_dim, $
st = st, $
dustem_sed = dustem_sed, $
SED_spec = SED_spec, $
dustem_qsed = dustem_qsed, $
Q_spec = Q_spec, $
dustem_used = dustem_used, $
U_spec = U_spec, $
dustem_polsed = dustem_polsed, $
P_spec = P_spec, $
dustem_polfrac = dustem_polfrac, $
SP_spec = SP_spec, $
dustem_psi_em = dustem_psi_em, $
PSI_spec = PSI_spec, $
dustem_ext = dustem_ext, $
EXT_spec = EXT_spec, $
dustem_qext = dustem_qext, $
QEXT_spec = QEXT_spec, $
dustem_uext = dustem_uext, $
UEXT_spec = UEXT_spec, $
dustem_polext = dustem_polext, $
POLEXT_spec = POLEXT_spec, $
dustem_fpolext = dustem_fpolext, $
SPEXT_spec = SPEXT_spec, $
dustem_psi_ext = dustem_psi_ext, $
PSIEXT_spec = PSIEXT_spec, $
_extra=_extra
ENDIF ELSE BEGIN
dustemwrap_plot_noobj, p_dim,$
st = st,$
dustem_sed = dustem_sed, $
SED_spec = SED_spec, $
dustem_qsed = dustem_qsed, $
Q_spec = Q_spec, $
dustem_used = dustem_used, $
U_spec = U_spec, $
dustem_polsed = dustem_polsed, $
P_spec = P_spec, $
dustem_polfrac = dustem_polfrac, $
SP_spec = SP_spec, $
dustem_psi_em = dustem_psi_em, $
PSI_spec = PSI_spec, $
dustem_ext = dustem_ext, $
EXT_spec = EXT_spec, $
dustem_qext = dustem_qext, $
QEXT_spec = QEXT_spec, $
dustem_uext = dustem_uext, $
UEXT_spec = UEXT_spec, $
dustem_polext = dustem_polext, $
POLEXT_spec = POLEXT_spec, $
dustem_fpolext = dustem_fpolext, $
SPEXT_spec = SPEXT_spec, $
dustem_psi_ext = dustem_psi_ext, $
PSIEXT_spec = PSIEXT_spec, $
_extra = _extra
ENDELSE
if !dustem_iter.act EQ 1 then begin
!dustem_iter.prv=0
!dustem_iter.act=0
endif else !dustem_iter.prv=!dustem_iter.act
endif
ENDIF
;print,"Total reduced chi2 = ", (wrchi2_sed+wrchi2_ext+wrchi2_polext+wrchi2_polsed+wrchi2_polfrac)/(n_sed+n_ext+n_polext+n_polsed+n_polfrac-DOF)
; Reduced chi2 is like if there were the same nb of data points for each curve
;print,"Mean reduced chi2 = ", (wrchi2_sed+wrchi2_ext+wrchi2_polext+wrchi2_polsed+wrchi2_polfrac) * total_n/(total_n-DOF)
; Plot size distribution
win+=1
xr=[0.3,5e3]
yr=[0.1,100]
;IF !dustem_show_plot NE 0 THEN dustem_plot_sdist,xr=xr,yr=yr,ps=filename,win=win;,/donotclose
; Suppress the artificial [0] at the beginning
dustem_out = dustem_out[1:*]
;IF defined(extra) THEN BEGIN
; IF keyword_set(extra.plot_it) THEN BEGIN
;; message,'params='+strtrim(p_min,2),/info,/continue
;; stop
; print,'dustem_mpfit_run_vg: Current parameters values:',p_dim
; ENDIF
;ENDIF
the_end:
RETURN,dustem_out
END