OV-GSO-DC / dustem-wrapper_idl

Blame view

src/idl/dustem_fit_sed_extsed_readme.pro 5.13 KB
Edit Raw Normal View History
5b76b765   Jean-Philippe Bernard   First commit 
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
 
PRO dustem_fit_sed_extsed_readme,postcript=postcript,model=model,help=help

;This Readme describes how to fit total intensity and polarization SEDs jointly
;Obviously, the dustem package must have been installed succesfully
;see DustemWrap documentation for install instructions.

;+
; NAME:
;    dustem_fit_sed_polsed_readme
; PURPOSE:
;    This is an example of how to fit total intensity and polarization SEDs jointly
;    It is meant to be an example to follow when writing your own
;    programs using the dustem IDL wrapper.
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    dustem_fit_sed_polsed_readme,postcript=postcript,mode=mode,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_sed_polsed_readme
; MODIFICATION HISTORY:
;    Written by J.P. Bernard June 29th 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_sed_extsed_readme'
  goto,the_end
ENDIF

;=== initialise dustem
;dustem_init,mode='COMPIEGNE_ETAL2010',/pol
use_model='COMPIEGNE_ETAL2010'
IF keyword_set(model) THEN use_model=model

dustem_init,mode=use_model

!dustem_verbose=1
!dustem_show_plot=1

;stop

;==== read emission SED

dir=!dustem_wrap_soft_dir+'/Data/SEDs/'
file=dir+'Gal_composite_spectrum.xcat'
;file=dir+'Mathis1990_DISM.xcat'

;file=dir+'Gal_composite_spectrum_wAKARI.xcat'
;stop
spec=read_xcat(file,/silent)
ind=where(spec.error EQ 0.,count)
IF count NE 0 THEN spec(ind).error=0.2*spec(ind).spec
ind=where(spec.instru EQ 'FIRAS',count)
IF count NE 0 THEN spec(ind).error=0.2*spec(ind).spec

;==== read extinction SED

dir=!dustem_wrap_soft_dir+'/Data/EXTs/'
file=dir+'Mathis1990_DISM.xcat'
extsed=read_xcat(file,/silent)
ind=where(extsed.error EQ 0.,count)
IF count NE 0 THEN extsed[ind].error=0.2*spec(ind).spec


st=dustem_set_data(ext=extsed,sed=spec,rchi2_weight=rchi2_weight,f_HI=f_HI)

;=== Set which parameters you want to fit
pd = [ $
             '(*!dustem_params).gas.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
              ]                   

p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04]    

;=== set initial value of parameters you want to fit
;iv=p_truth   ;exact solution
iv=p_truth+[0.,3.e-5,1.e-5,-2.e-5]    ;shifted from solution

ulimed=[0   , 0   ,0   ,0 ]
llimed=[1   , 1   ,1   ,1 ]
llims =[0.  ,0.   ,0.  ,0.]

;== 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

;=== RUN fit
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)

yr=[1e-4,10]
xr=[1,1.e5]
legend_xpos=0.6
legend_ypos=0.8
tit='DUSTEMWrap Intensity SED (Running)'
ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
xtit=textoidl('\lambda (\mum)')
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

;help,(*!dustem_fit),/STRUCTURE

;=== SAVE FIT RESULTS
file_out='/tmp/DUSTEM_sedextsed_fit_example.sav'
dustem_save_system_variables,file_out
message,'Saved '+file_out,/continue

stop

;======================================
;====You can exit IDL here and re-enter
;======================================
file='/tmp/DUSTEM_sedextsed_fit_example.sav'
dustem_restore_system_variables,file

;=== Plot best fit
win=1
window,win & win=win+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 polarization 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,/pol

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