dustem_mpfit_run.pro
27.8 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
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
FUNCTION dustem_mpfit_run,x,p_min,_EXTRA=extra
;+
; NAME:
; dustem_mpfit_run_vg
; 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
; see evolution details on the dustem cvs maintained at CESR
; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-
;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
st = dustem_run(p_dim)
st_model=dustem_read_all(!dustem_soft_dir)
dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),res=res,st
chi2_sed = 0
chi2_ext = 0
chi2_polext = 0
chi2_polsed = 0
chi2_polfrac = 0
chi2_qsed =0
chi2_used =0
chi2_qext =0
chi2_uext =0
rchi2_sed = 0
rchi2_ext = 0
rchi2_polext = 0
rchi2_polsed = 0
rchi2_polfrac = 0
rchi2_qsed =0
rchi2_used =0
rchi2_qext =0
rchi2_uext =0
wrchi2_sed = 0
wrchi2_ext = 0
wrchi2_polext = 0
wrchi2_polsed = 0
wrchi2_polfrac = 0
wrchi2_qsed =0
wrich2_used =0
wrchi2_qext =0
wrich2_uext =0
n_sed = 0
n_ext = 0
n_polext = 0
n_polsed = 0
n_polfrac = 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)
;if exist('./noplot') then !dustem_show_plot = 0 else if !dustem_show_plot eq 0 then !dustem_show_plot = 1
FOR i_tag=0,count_data_tag-1 DO BEGIN
;stop
;print,i_tag,ind_data_tag(i_tag)
;print,tagnames(ind_data_tag(i_tag))
CASE tagnames(ind_data_tag(i_tag)) OF
'SED' : BEGIN
;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_SED
dustem_sed = dustem_compute_sed(p_dim,st);dustem_compute_sed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
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
(*!dustem_fit).chi2=chi2_sed
(*!dustem_fit).rchi2=rchi2_sed
; Plot if needed
IF !dustem_show_plot NE 0 THEN BEGIN
;title='DUSTEM WRAP('+strtrim(string(win),2)+')'
window,win,title='DUSTEM WRAP (SED)'
;cgwindow,WMULTI=[2,2,1],/CURRENT
;dustem_plot_fit_sed,st,dustem_sed,cont,freefree,synchrotron,_extra=extra,res=p_min*(*(*!dustem_fit).param_init_values),chi2=(*!dustem_fit).chi2,rchi2=(*!dustem_fit).rchi2
dustem_plot_fit_sed,st,dustem_sed,_extra=extra,res=p_min*(*(*!dustem_fit).param_init_values),chi2=(*!dustem_fit).chi2,rchi2=(*!dustem_fit).rchi2
;stop
ENDIF
dustem_out = [ dustem_out, dustem_sed ]
END
'EXT' : BEGIN
;ADDING PLUGIN TO SPECTRUM----------------
;if n_tags(!dustem_data.ext) gt 1 then begin
for i=0L,n_tags(*!dustem_scope)-1 do begin
if total(strsplit((*!dustem_scope).(i),'+',/extract) eq 'ADD_EXT') then st.ext.ext_tot+=(*(*!dustem_plugin).(i))
endfor
;endif
;------------------------------------------
dustem_ext = interpol(st.ext.ext_tot,st.ext.wav,(*!dustem_data.ext).wav)
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
; Plot if needed
IF !dustem_show_plot NE 0 THEN BEGIN
win+=1
;dustem_plot_ext,st,yr=[1.e-6,1.e0],/ysty,xr=[0.5,40],/xsty,win=win
xr=[0.3,40]
yr=[1.e-26,1.e-21]
dustem_plot_ext,st,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1
dustem_plot_ext,st,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2
win+=1
;dustem_plot_ext,st,/UV,yr=[0,2.5],/ysty,xr=[0,10],/xsty,win=win
xr=[0.2,10]
yr=[0,3e-21]
dustem_plot_ext,st,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1
dustem_plot_ext,st,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2
ENDIF
dustem_plot_ext,st,/print_Rv
dustem_out = [ dustem_out, dustem_ext ]
END
'POLEXT': BEGIN
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,dustem_qext,dustem_uext,Q_ext,U_ext)
chi2_polext = total(((*!dustem_data.polext).values-dustem_polext)^2/(*!dustem_data.polext).sigma^2)
n_polext = n_elements((*!dustem_data.polext).values)
rchi2_polext = chi2_polext / n_polext
wrchi2_polext = rchi2_polext * !fit_rchi2_weight.polext
print,"chi2 POLEXT = ",chi2_polext," rchi2 POLEXT = ",rchi2_polext;," weighted rchi2 POLEXT=",wrchi2_polext
;print,"chi2 per wl: ",((*!dustem_data.polext).values-dustem_polext)^2/(*!dustem_data.polext).sigma^2
; Fit in (pmax,lmax,K)
x=st.polext.wav
logx=alog(x)
logy=alog(st.polext.ext_tot)
;j=where(st.polext.ext_tot gt 0 and x gt 0.2 and x lt 1)
; Whittet+1992 (Table 1)
Uband = 0.36
Bband = 0.43
Vband = 0.55
Rband = 0.63
Iband = 0.78
Jband = 1.21
Hband = 1.64
Kband = 2.04
bands = [Uband,Bband,Vband,Rband,Iband,Jband,Hband]
nbands = n_elements(bands)
jband = indgen(nbands)*0
for k=0,nbands-1 do begin
xx=min(abs(st.polext.wav-bands(k)),j)
jband(k) = j
endfor
;j=where(st.polext.ext_tot gt 0 and x gt 0.2 and x lt 1)
res = poly_fit(logx(jband),logy(jband),2)
print
print,"Serkowski Fit with free K over ",string(st.polext(jband).wav,format='(100F5.2)')
print,"-------------------------"
K = -res(2)
lmax = exp(res(1)/(2*K))
pmax = exp(res(0)+K*alog(lmax)^2)
print,"* lmax = ",lmax," micron"
print,"* K = ",K
print,"* pmax = ",pmax," for NH=1e21"
print
; Plot if needed
IF !dustem_show_plot NE 0 THEN BEGIN
; ; Polarization fraction in the UV-IR
; win+=1
;
; dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,yr=[0,0.15],/ysty,xr=[0.1,50],/xsty,win=win;multi=0,ps=filename
; ;
; Polarisation curve (Serkowski)
win+=1
;xr=[0.3,1e4]
xr=[1e-4,3.5]
yr=[1e-30,5e-24]
dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,/donotclose,multi=1
dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,multi=2
; win+=1
; dustem_plot_fpol,win=win,xr=[10,1000],yr=[0,1.1]
;
ENDIF
;dustem_plot_polar,st,/print_ratio,cont=cont
dustem_out = [ dustem_out, dustem_polext ]
ENDIF
END
;USE OF POLFRAC DEPRECATED---------------------------------------
; 'POLFRAC': BEGIN
;
; ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_POLSED AND DUSTEM_COMPUTE_SED
;
; IF !run_lin then begin
;
; dustem_polsed = dustem_compute_polsed(p_dim,st);dustem_compute_polsed(p_dim,st.polsed,cont=cont,freefree=freefree,synchrotron=synchrotron)
; dustem_sed_tmp = dustem_compute_sed(p_dim,st);dustem_compute_sed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
; dustem_sed = dustem_polsed
; ind_filt=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_filt)
; for i = 0, count_filt-1 do begin
; j=where((*!dustem_data.polsed).filt_names(ind_filt(i)) eq (*!dustem_data.sed).filt_names,count)
; if count eq 1 then dustem_sed(i) = dustem_sed_tmp(j(0))
; endfor
; dustem_polfrac = dustem_polsed/dustem_sed
;
; ; We normalize at 353GHz
; ;x=min(((*!dustem_data.polfrac).wav-850)^2,j353)
; ;dustem_polfrac = dustem_polfrac/dustem_polfrac(j353)*(*!dustem_data.polfrac).values(j353)
;
; chi2_polfrac = total(((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2)
; n_polfrac = n_elements((*!dustem_data.polfrac).values)
; rchi2_polfrac = chi2_polfrac / n_polfrac
; wrchi2_polfrac = rchi2_polfrac * !fit_rchi2_weight.polfrac
; print,"chi2 POLFRAC = ",chi2_polfrac," rchi2 POLFRAC = ",rchi2_polfrac;," weighted rchi2 POLFRAC=",wrchi2_polfrac
; ;print,"chi2 per wl: ",((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2
;
; dustem_out = [ dustem_out, dustem_polfrac ]
;
; ; Plot if needed
; IF !dustem_show_plot NE 0 THEN BEGIN
;
; ; Also show prediction for polarized P/I
; win+=1
; dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/Pfrac,win=win,yr=[0,0.25],/ysty,xr=[10,5d3];dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/Pfrac,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,yr=[0,0.25],/ysty,xr=[10,5d3]
;
; ENDIF
;
; ENDIF
;
; END
;-----------------------------------------------------------------
'POLSED': BEGIN
;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_SED
IF !run_lin THEN BEGIN
dustem_polsed = dustem_compute_polsed(p_dim,st,dustem_qsed,dustem_used,Q_sed,U_sed);dustem_compute_polsed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
chi2_polsed = total(((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2)
n_polsed = n_elements((*!dustem_data.polsed).values)
rchi2_polsed = chi2_polsed / n_polsed
wrchi2_polsed = rchi2_polsed * !fit_rchi2_weight.polsed
message,"chi2 POLSED = "+strtrim(chi2_polsed,2)+" rchi2 POLSED = "+strtrim(rchi2_polsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
print,"chi2 per wl: ",((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2
;print,"SNR per wl: ",(*!dustem_data.polsed).values/(*!dustem_data.polsed).sigma
dustem_out = [ dustem_out, dustem_polsed ] ; classic command
;-----------------------------------------------
; Plot if needed
IF !dustem_show_plot NE 0 THEN BEGIN
win+=1
;
xr=[10,8e3]
yr=[1e-34,1e-28]
dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,xr=xr,yr=yr,/xsty,/donotclose,multi=1
dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,xr=xr,yr=yr,/xsty,multi=2
ENDIF
ENDIF
END
'QSED': BEGIN
;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_SED
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
if !dustem_show_plot NE 0 then begin
fact = (3.e10/(1.e-4*(st.polsed.wav)))/1.d40
fact1 = (3.e10/(1.e-4*((*!dustem_data.qsed).wav)))/1.d40 ; for data plotting
win+=1
window ,win ,title='DUSTEM WRAP (QSED)'
titstq='Stokes Q Polarization Intensity'
ytitstq=textoidl('\nuQ_\nu/N_H (W/m^2/sr/H)')
xtit=textoidl('\lambda (\mum)')
xr=[10,8e3]
normpol=Q_sed*fact
normpol1=dustem_qsed*fact1
;normpol1=(*!dustem_data.qsed).values*fact1
yr=[1e-34,1e-28]
indq=where(Q_sed lt 0, countq)
ylog=1
if countq ne 0 then begin
ylog=0
yr=[-5e-31,5e-31]
endif
cgplot,st.polsed.wav,Q_sed*fact,xtit='',ytit=ytitstq,tit=titstq,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)'
cgoplot,st.polsed.wav,Q_sed*fact,color='black'
cgoplot,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1,color='Tomato',psym=16,symsize=1,thick=2
err_bar,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1,yrms=3.*((*!dustem_data.qsed).sigma)/2.*fact1,color=cgColor('Tomato')
cgoplot,(*!dustem_data.qsed).wav,dustem_qsed*fact1,psym=6,color='red',symsize=2
cgplot,st.polsed.wav,Q_sed*fact/normpol,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1
cgoplot,st.polsed.wav,Q_sed*fact/normpol,color='black'
cgoplot,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1/normpol1,color='Tomato',psym=16,symsize=1,thick=2
err_bar,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1/normpol1,yrms=3.*((*!dustem_data.qsed).sigma)/2.*fact1/normpol1,color=cgColor('Tomato')
endif
END
'USED': BEGIN
;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_SED
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
if !dustem_show_plot ne 0 then begin
fact = (3.e10/(1.e-4*(st.polsed.wav)))/1.d40
fact1 = (3.e10/(1.e-4*((*!dustem_data.used).wav)))/1.d40 ; for data plotting
win+=1
window ,win ,title='DUSTEM WRAP (USED)'
titstu='Stokes U Polarization Intensity'
ytitstu=textoidl('\nuU_\nu/N_H (W/m^2/sr/H)')
xtit=textoidl('\lambda (\mum)')
xr=[10,8e3]
normpol2=U_sed*fact
normpol3=dustem_used*fact1
;yr=[min(normpol2)-min(normpol2)/10,max(normpol2)+max(normpol2)*10]
yr=[1e-34,1e-28]
indu=where(U_sed lt 0, countu)
ylog=1
if countu ne 0 then begin
ylog=0
yr=[-5e-31,5e-31]
endif
cgplot,st.polsed.wav,U_sed*fact,xtit='',ytit=ytitstu,tit=titstu,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',color='black'
cgoplot,st.polsed.wav,U_sed*fact,color='black'
cgoplot,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1,color='Teal',psym=16,symsize=1,thick=2
err_bar,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1,yrms=3.*((*!dustem_data.used).sigma)/2.*fact1,color=cgColor('Teal')
cgoplot,(*!dustem_data.used).wav,dustem_used*fact1,psym=6,color='red',symsize=2
cgplot,st.polsed.wav,U_sed*fact/normpol2,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1
cgoplot,st.polsed.wav,U_sed*fact/normpol2,color='black'
cgoplot,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1/normpol3,color='Teal',psym=16,symsize=1,thick=2
err_bar,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1/normpol3,yrms=3.*((*!dustem_data.used).sigma)/2.*fact1/normpol3,color=cgColor('Teal')
;cgwindow, 'plot',
endif
END
'QEXT': BEGIN ; Q - EXTINCTION CROSS SECTION
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
if !dustem_show_plot ne 0 then begin
win+=1
window ,win ,title='DUSTEM WRAP (QEXT)'
titqex='Q-polarized Cross-section'
ytitqex=textoidl('\sigma_{Q}/N_H (cm^2/H)')
xtit=textoidl('\lambda^{-1} (\mum^{-1})')
;xr=[0.3,1e4]
xr=[1e-4,3.5]
yr=[1e-34,5e-24]
indqex=where(Q_ext lt 0, countqex)
ylog=1
if countqex ne 0 then begin
ylog=0
yr=[-9e-27,5e-27]
endif
normpol = st.polext.ext_tot * 0. + 1.d21
cgplot,1/st.polext.wav,Q_ext,xtit='',ytit=ytitqex,tit=titqex,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',color='black'
cgoplot,1/st.polext.wav,Q_ext/normpol,color='black'
cgoplot,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/normpol,color='Indian Red',psym=16,symsize=1,thick=2
err_bar,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/normpol,yrms=3.*((*!dustem_data.qext).sigma)/2/normpol,color=cgColor('Indian Red')
cgoplot,1/(*!dustem_data.qext).wav,dustem_qext/normpol,psym=6,color='red',symsize=2
cgplot,1/st.polsed.wav,Q_ext/Q_ext,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1
cgoplot,1/st.polsed.wav,Q_ext/Q_ext,color='black'
cgoplot,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/dustem_qext,color='Indian Red',psym=16,symsize=1,thick=2
err_bar,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/dustem_qext,yrms=3.*((*!dustem_data.qext).sigma)/2.,color=cgColor('Indian Red')
endif
END
;
'UEXT': BEGIN
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
if !dustem_show_plot ne 0 then begin
win+=1
window ,win ,title='DUSTEM WRAP (UEXT)'
tituex='U-polarized Cross-section'
ytituex=textoidl('\sigma_{U}/N_H (cm^2/H)')
xtit=textoidl('\lambda^{-1} (\mum^{-1})')
;xr=[0.3,1e4]
xr=[1e-4,3.5]
yr=[1e-34,5e-24]
induex=where(U_ext lt 0, countuex)
ylog=1
if countuex ne 0 then begin
ylog=0
yr=[-9e-27,5e-27]
endif
normpol = st.polext.ext_tot * 0. + 1.d21
cgplot,1/st.polext.wav,U_ext,xtit='',ytit=ytituex,tit=tituex,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',color='black'
cgoplot,1/st.polext.wav,U_ext/normpol,color='black'
cgoplot,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/normpol,color='Steel Blue',psym=16,symsize=1,thick=2
err_bar,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/normpol,yrms=3.*((*!dustem_data.uext).sigma)/2/normpol,color=cgColor('Steel Blue')
cgoplot,1/(*!dustem_data.uext).wav,dustem_uext/normpol,psym=6,color='red',symsize=2
cgplot,1/st.polsed.wav,U_ext/U_ext,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1
cgoplot,1/st.polsed.wav,U_ext/U_ext,color='black'
cgoplot,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/dustem_uext,color='Steel Blue',psym=16,symsize=1,thick=2
err_bar,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/dustem_uext,yrms=3.*((*!dustem_data.uext).sigma)/2/normpol,color=cgColor('Steel Blue')
endif
END
ELSE : stop
ENDCASE
ENDFOR
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 !fit_rchi2_weight.polext ne 0 then begin
total_n += n_polext
total_chi2 += chi2_polext
total_rchi2 += rchi2_polext
endif
if !fit_rchi2_weight.polsed ne 0 then begin
total_n += n_polsed
total_chi2 += chi2_polsed
total_rchi2 += rchi2_polsed
endif
if !fit_rchi2_weight.polfrac ne 0 then begin
total_n += n_polfrac
total_chi2 += chi2_polfrac
total_rchi2 += rchi2_polfrac
endif
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
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
total_wchi2 = chi2_ext * !fit_rchi2_weight.ext $
+ chi2_sed * !fit_rchi2_weight.sed
if !run_pol then begin
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
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
message,'======================================================================',/continue
;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