plot_lastrun_vg.pro
3.21 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
PRO plot_lastrun_vg,run=run
if keyword_set(run) then dustem_init,/wraptest
loadct,13
defsysv,'!dustem_inputs',ptr_new()
dustem_inputs={GRAIN:'GRAIN.DAT',ISRF:'ISRF.DAT'}
!dustem_inputs=ptr_new(dustem_inputs)
;=== Read model ===
st_model=dustem_read_all(!dustem_soft_dir)
defsysv,'!dustem_params' , ptr_new(st_model)
st=dustem_read_all_res(!dustem_res,/silent)
;=== Read observational data ===
data_files = { $
sed: !dustem_wrap_soft_dir + 'Data/SEDs/Gal_composite_spectrum.xcat', $
ext: !dustem_wrap_soft_dir + 'Data/EXTs/Mathis1990_DISM.xcat', $
polext: !dustem_wrap_soft_dir + 'Data/POLEXTs/Serkowski_DISM.xcat', $
polsed: '' $
}
defsysv, '!dustem_data', { $ ;Data to fit
sed: ptr_new(), $
ext: ptr_new(), $
polext: ptr_new() $
; polsed: ptr_new() $
}
tagnames = tag_names(!dustem_data)
for i = 0, n_elements(tagnames)-1 do begin
instr='!dustem_data.('+strtrim(i,2)+') = dustem_set_data(read_xcat(data_files.'+tagnames(i)+',/silent))'
toto=execute(instr)
endfor
if tag_exist(st,"sed") then begin
window,1
;data_sed = read_xcat(data_files.sed,/silent)
fact = 1./4./!pi*1.e4/1.e7*1.e20
dustem_plot_nuinu_em,st,yr=[1e-11,1.e-7],/ysty,xr=[1,5e3],/xsty
dustem_sed = interpol(st.sed.em_tot*fact,st.sed.wav,(*!dustem_data.sed).wav)
chi2_sed = total(((*!dustem_data.sed).values*(3.e8/1.e-6/(*!dustem_data.sed).wav)/1.d20-dustem_sed)^2/((*!dustem_data.sed).sigma*(3.e8/1.e-6/(*!dustem_data.sed).wav)/1.d20)^2)
n_sed = n_elements((*!dustem_data.sed).values)
rchi2_sed = chi2_sed / n_sed
;wrchi2_sed = rchi2_sed * !fit_rchi2_weight.sed
print,"chi2 SED = ",chi2_sed," rchi2 SED = ",rchi2_sed
endif
if tag_exist(st,"ext") then begin
window,2
dustem_plot_ext,st,st_model,yr=[1.e-6,1.e0],/ysty,xr=[0.3,400],/xsty
window,3
dustem_plot_ext,st,st_model,/UV,yr=[0,2.5],/ysty,xr=[0,10],/xsty
dustem_ext = interpol(st.ext.ext_tot,st.sed.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
endif
if tag_exist(st,"polext") and tag_exist(st,"polsed") then begin
; Polarization fraction in the UV-IR
window,4
dustem_plot_polar,st,st_model,yr=[0,0.08],/ysty,xr=[0.1,4],/xsty
; Polarisation curve (Serkowski)
window,5
dustem_plot_polar,st,st_model,/UV,yr=[1e-3,3e-2],/ysty,xr=[0.1,10],/xsty
; Polarization fraction from the UV to the submm, extinction and emission
window,6
dustem_plot_polar,st,st_model,yr=[0,0.30],/ysty,xr=[0.1,5e3],/xsty
dustem_polext = interpol(st.polext.ext_tot,st.polext.wav,(*!dustem_data.polext).wav)
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_sed = rchi2_sed * !fit_rchi2_weight.polext
print,"chi2 POLEXT = ",chi2_polext," rchi2 POLEXT = ",rchi2_polext
endif
dustem_plot_sdist
stop
END