Blame view

src/idl/dustem_fit_sed_ext_pol_example.pro 10.1 KB
1d10fa6d   Ilyes Choubani   extinction exampl...
1
PRO dustem_fit_sed_ext_pol_example,model=model $
9e027842   Jean-Philippe Bernard   updated to new fu...
2
3
4
5
6
                                  ,sed_file=sed_file $
                                  ,ext_file=ext_file $
                                  ,Nitermax=Nitermax $
                                  ,fits_save_and_restore=fits_save_and_restore $
                                  ,wait=wait $
8f9df3fb   Annie Hughes   tidied up code, a...
7
8
9
                                   ,noobj=noobj $
                                   ,help=help $
                                  ,verbose=verbose
9e027842   Jean-Philippe Bernard   updated to new fu...
10
11
12
13
14

;+
; NAME:
;    dustem_fit_sed_ext_pol_example  
;
8f9df3fb   Annie Hughes   tidied up code, a...
15
16
17
; PURPOSE: This routine is an example of how to fit an observational
; data in emission + extinction SED (StokesI,Q,U) with DustEM and
; DustEMWrap.
9e027842   Jean-Philippe Bernard   updated to new fu...
18
19
20
21
22
23
24
;  
; See the DustEMWrap User Guide for more information.
;
; CATEGORY:
;    DustEMWrap, Distributed, High-Level, User Example
;
; CALLING SEQUENCE:
8f9df3fb   Annie Hughes   tidied up code, a...
25
;    dustem_fit_ext_pol_example[,model=][sed_file=][ext_file=][,Nitermax=][,fits_save_and_restore=][,/help,/wait,/verbose]
9e027842   Jean-Philippe Bernard   updated to new fu...
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
;
; INPUTS:
;    None
;
; OPTIONAL INPUT PARAMETERS:
;    None
;
; OUTPUTS:
;    None
;
; OPTIONAL OUTPUT PARAMETERS:
;    Plots, results structure in binary FITS table format
;
; ACCEPTED KEY-WORDS:
;    model = specifies the interstellar dust mixture used by DustEM
;           'MC10' model from Compiegne et al 2010 
;           'DBP90' model from Desert et al 1990
;           'DL01' model from Draine & Li 2001
8f9df3fb   Annie Hughes   tidied up code, a...
44
;           'WD01_RV5P5B' model from Weingartner & Draine 2002 with Rv=5.5
9e027842   Jean-Philippe Bernard   updated to new fu...
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
;           'DL07' model from Draine & Li 2007
;           'J13' model from Jones et al 2013, as updated in
;                 Koehler et al 2014
;           'G17_ModelA' model A from Guillet et al (2018). Includes
;                 polarisation. See Tables 2 and 3 of that paper for details.
;           'G17_ModelB' model B from Guillet et al (2018)
;           'G17_ModelC' model C from Guillet et al (2018)
;           'G17_ModelD' model A from Guillet et al (2018)
;    sed_file = string naming the path to text file in .xcat format that
;          describes the observational SED. If not set, the file
;          'Data/EXAMPLE_OBSDATA/example_SED_2.xcat' is used.  
;    ext_file = string naming the path to text file in .xcat format that
;          describes the observational extinction SED. If not set, the file
;          'Data/EXAMPLE_OBSDATA/Mathis1990_DISM.xcat' is used.  
;    Nitermax = maximum number of fit iterations. Default is 5.
;    fits_save_and_restore = if set, save the fit results in a binary
;               FITS file. The code then restore this file and plots
;               the results using the saved results information.
;    help      = if set, print this help
;    wait      = if set, wait this many seconds between each step of
;                the code (for illustration purposes)
;    verbose      = if set, subroutines will run in verbose mode
8f9df3fb   Annie Hughes   tidied up code, a...
67
;    noobj     = if set, runs with no object graphics
9e027842   Jean-Philippe Bernard   updated to new fu...
68
69
70
71
72
73
74
75
76
77
78
79
;
; COMMON BLOCKS:
;    None
;
; SIDE EFFECTS:
;    None
;
; RESTRICTIONS:
;    The DustEM fortran code must be installed
;    The DustEMWrap IDL code must be installed
;
; PROCEDURES AND SUBROUTINES USED:
9e027842   Jean-Philippe Bernard   updated to new fu...
80
81
82
83
84
85
86
87
88
89
90
;
; EXAMPLES
;    dustem_fit_sed_ext_pol_example
;    dustem_fit_sed_ext_pol_example,Nitermax=1,fits_save_and_restore='/tmp/mysavefile.fits'
;    dustem_fit_sed_ext_pol_example,model='DBP90'
;
; MODIFICATION HISTORY:
;    Written by JPB Apr-2011
;    Evolution details on the DustEMWrap gitlab.
;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  
;-
1d10fa6d   Ilyes Choubani   extinction exampl...
91
92

IF keyword_set(help) THEN BEGIN
9e027842   Jean-Philippe Bernard   updated to new fu...
93
  doc_library,'dustem_fit_sed_ext_pol_example'
1d10fa6d   Ilyes Choubani   extinction exampl...
94
95
96
97
98
99
  goto,the_end
END

IF keyword_set(model) THEN BEGIN
  use_model=strupcase(model)
ENDIF ELSE BEGIN
8f9df3fb   Annie Hughes   tidied up code, a...
100
  use_model='G17_MODELA'    
1d10fa6d   Ilyes Choubani   extinction exampl...
101
102
ENDELSE

8f9df3fb   Annie Hughes   tidied up code, a...
103
104
105
106
107
108
109
known_mdls=['MC10','DBP90','DL01','WD01_RV5P5B','DL07','J13','G17_MODELA','G17_MODELB','G17_MODELC','G17_MODELD'] 
test_model = where(known_mdls eq use_model,ct)
if ct eq 0 then begin
   message,'ISM dust model '+use_model+' unknown',/continue
   message,'Known models are MC10,DBP90,DL01,WD01_RV5P5B,DL07,J13,G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD',/continue
   stop
end
dc84fd94   Ilyes Choubani   Small corrections...
110

8f9df3fb   Annie Hughes   tidied up code, a...
111
use_polarization=0.   ;default is no polarization in models
1d10fa6d   Ilyes Choubani   extinction exampl...
112
use_window=2       ; default graphics window number to use for plotting the results
9e027842   Jean-Philippe Bernard   updated to new fu...
113
114
use_verbose=0
if keyword_set(verbose) then use_verbose=1
8f9df3fb   Annie Hughes   tidied up code, a...
115
116
use_Nitermax=5        ; maximum number of iterations for the fit
IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax
1d10fa6d   Ilyes Choubani   extinction exampl...
117

8f9df3fb   Annie Hughes   tidied up code, a...
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
dustem_define_la_common

;=== Set the (model-dependent) parameters that you want to fit (pd),
;=== their initial values (iv)
;=== and whether they are bounded (ulimed,llimed,llims,ulims).
;=== Fixed parameters (fpd) and their values (fiv) are also set here.
;=== Refer to the DustEM and DustEMWrap User guides for an explanation
;=== of the physical meaning of dust model and plug-in parameters, and
;=== how to specify them.

;=== Example is provided for G17_MODELA
;=== If you choose a different dust model, you will need to adjust the
;=== parameter information accordingly.
pd = [ $
     '(*!dustem_params).grains(0).mdust_o_mh',$         ;PAH0 mass fraction
     '(*!dustem_params).grains(1).mdust_o_mh',$         ;amCBE mass fraction
     '(*!dustem_params).grains(2).mdust_o_mh',$         ;aSil
     'dustem_plugin_modify_dust_polx_2' ]

true_vals = [5.4e-4, 5.4e-4,1.8e-4,50] ; these are the true values that will be used to generate the data
iv = true_vals+[10.00E-4,10.00E-4,10.00E-4,-30] ; this is the vector of starting guesses for the fit

Npar=n_elements(pd)
ulimed=replicate(0,Npar)
llimed=replicate(1,Npar) & llimed[3]=0
llims=replicate(0.,Npar)
use_polarization=1

;=== Fixed parameters
fpd=[ $
    '(*!dustem_params).G0' ,  $     ;multiplicative factor to total ISRF
    'dustem_plugin_continuum_2' $       ;intensity of a BB continuum
    ]

;=== Initial parameter values for fixed parameters
fiv=[ $
    1. , $                      ;multiplicative factor to total ISRF
    3.e-3$                      ;intensity of BB continuum
]
1d10fa6d   Ilyes Choubani   extinction exampl...
157
        
1d10fa6d   Ilyes Choubani   extinction exampl...
158
159


8f9df3fb   Annie Hughes   tidied up code, a...
160
161
162
163
if keyword_set(wait) then begin
   message,'Finished setting dust model and plug-in parameters: '+use_model,/info
   wait,wait
end
1d10fa6d   Ilyes Choubani   extinction exampl...
164
165

;== INITIALISE DUSTEM
8f9df3fb   Annie Hughes   tidied up code, a...
166
dustem_init,model=use_model,polarization=use_polarization
1d10fa6d   Ilyes Choubani   extinction exampl...
167
!dustem_nocatch=1
9e027842   Jean-Philippe Bernard   updated to new fu...
168
!dustem_verbose=use_verbose
8f9df3fb   Annie Hughes   tidied up code, a...
169
170
IF keyword_set(noobj) THEN !dustem_noobj=1
!EXCEPT=2 ; for debugging
1d10fa6d   Ilyes Choubani   extinction exampl...
171

1d10fa6d   Ilyes Choubani   extinction exampl...
172

8f9df3fb   Annie Hughes   tidied up code, a...
173
174
175
;=== READ EXAMPLE DATA IN EXTINCTION AND EMISSION TO INITIALISE
;=== THE OBSERVATIONAL STRUCTURES. NOTE THAT WE GENERATE THE
;=== ACTUAL DATA USED IN THE FIT USING THE MODEL AND TRUE_VALS
1d10fa6d   Ilyes Choubani   extinction exampl...
176
177

dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'
8f9df3fb   Annie Hughes   tidied up code, a...
178
179
file=dir+'Mathis1990_DISM.xcat' 
if keyword_set(ext_file) then file=ext_file 
1d10fa6d   Ilyes Choubani   extinction exampl...
180
181
182
183
ext=read_xcat(file,/silent)

dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'
file=dir+'example_SED_2.xcat'
8f9df3fb   Annie Hughes   tidied up code, a...
184
if keyword_set(sed_file) then file=sed_file 
1d10fa6d   Ilyes Choubani   extinction exampl...
185
186
sed=read_xcat(file,/silent)

8f9df3fb   Annie Hughes   tidied up code, a...
187
;;=== CREATE THE UNCERTAINTIES (NEEDED OTHERWISE DUSTEM_SET_DATA WILL COMPLAIN)
1d10fa6d   Ilyes Choubani   extinction exampl...
188
189
190
191
192
ext.EXT_I=1.
for i=4l,n_tags(ext)-1 do begin
    ext.(i) = ext.EXT_I/1E13
endfor

1d10fa6d   Ilyes Choubani   extinction exampl...
193
194
195
196
197
sed.StokesI=1.
for i=4l,n_tags(sed)-1 do begin
    sed.(i) = sed.StokesI/1E10
endfor

8f9df3fb   Annie Hughes   tidied up code, a...
198
199
;=== HERE WE SET THE DATA SO THAT THE COMPUTE_ FUNCTIONS BELOW HAVE
;=== THE NECESSARY FIELDS AND INFORMATION
1d10fa6d   Ilyes Choubani   extinction exampl...
200
201
dustem_set_data,sed,sed,ext,ext 

8f9df3fb   Annie Hughes   tidied up code, a...
202
203
204
;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS FOR THE FIT 
;== AND ACTIVATE ANY PLUGINS (HENCE WE NEED TO DO THIS BEFORE THE COMPUTE_ CALLS)
dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims,polarization=use_polarization
1d10fa6d   Ilyes Choubani   extinction exampl...
205

8f9df3fb   Annie Hughes   tidied up code, a...
206
207
208
209
if keyword_set(wait) then begin
   message,'Finished initializing DustEMWrap, including plugins and fixed parameters',/info
   wait,wait
end
1d10fa6d   Ilyes Choubani   extinction exampl...
210

8f9df3fb   Annie Hughes   tidied up code, a...
211
212
213
;== GENERATE THE DATA USING THE MODEL AND TRUE VALUES
ext.EXT_I = dustem_compute_ext(true_vals,st=st)
toto = dustem_compute_stokext(true_vals,st=st,Qext_spec=dustem_qext,Uext_spec=dustem_uext)
9e027842   Jean-Philippe Bernard   updated to new fu...
214
215
ext.EXT_Q = toto[0]
ext.EXT_U = toto[1]
1d10fa6d   Ilyes Choubani   extinction exampl...
216

8f9df3fb   Annie Hughes   tidied up code, a...
217
218
sed.StokesI = dustem_compute_sed(true_vals,st=st)
toto = dustem_compute_stokes(true_vals,st=st,Q_spec=dustem_qsed,U_spec=dustem_used)
9e027842   Jean-Philippe Bernard   updated to new fu...
219
220
sed.StokesQ = toto[0]
sed.StokesU = toto[1]
1d10fa6d   Ilyes Choubani   extinction exampl...
221

8f9df3fb   Annie Hughes   tidied up code, a...
222
223
;== SET THE OBSERVATIONAL STRUCTURE, NOW WITH THE DATA GENERATED USING
;== THE MODEL
dc84fd94   Ilyes Choubani   Small corrections...
224
dustem_set_data, sed,sed,ext,ext
1d10fa6d   Ilyes Choubani   extinction exampl...
225

8f9df3fb   Annie Hughes   tidied up code, a...
226
227
228
229
if keyword_set(wait) then begin
   message,'Finished generating data from the model',/info
   wait,wait
end
1d10fa6d   Ilyes Choubani   extinction exampl...
230

1d10fa6d   Ilyes Choubani   extinction exampl...
231

8f9df3fb   Annie Hughes   tidied up code, a...
232
233
234
235
236
;== INFORMATION TO RUN THE FIT
tol=1.e-10  ;fit tolerence

;=== INFORMATION TO MAKE THE PLOT
;=== _x/_extinction means extinction plots, _m/_emissions means emission plots
1d10fa6d   Ilyes Choubani   extinction exampl...
237
xr_x = [0.01,30]
dc84fd94   Ilyes Choubani   Small corrections...
238
yr_x = [5.00E-8,10]
1d10fa6d   Ilyes Choubani   extinction exampl...
239
240
xr_m = [1.,5e5]
yr_m = [5e-8,1.00e6]
8f9df3fb   Annie Hughes   tidied up code, a...
241
242
tit_m='Spectral Energy Distribution' 
tit_x='Dust Optical Depth' 
3d0f6b67   Ilyes Choubani   rest of plotting ...
243
244

;These two are not needed anymore for any of the plottings.... 
8f9df3fb   Annie Hughes   tidied up code, a...
245
246
;ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
;xtit=textoidl('\lambda (\mum)')
3d0f6b67   Ilyes Choubani   rest of plotting ...
247

8f9df3fb   Annie Hughes   tidied up code, a...
248
;===  RUN THE FIT
1d10fa6d   Ilyes Choubani   extinction exampl...
249
250
t1=systime(0,/sec)
res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $
3d0f6b67   Ilyes Choubani   rest of plotting ...
251
                      ,/xlog,/ylog,xr_m=xr_m,yr_m=yr_m,xr_x=xr_x,yr_x=yr_x,xtit=xtit,ytit=ytit,tit_m=tit_m,tit_x=tit_x $
1d10fa6d   Ilyes Choubani   extinction exampl...
252
                      ,legend_xpos=legend_xpos,legend_ypos=legend_ypos $
dc84fd94   Ilyes Choubani   Small corrections...
253
                      ,errors=errors,chi2=chi2,rchi2=rchi2,show_plot=show_plot)
1d10fa6d   Ilyes Choubani   extinction exampl...
254
255
t2=systime(0,/sec)

8f9df3fb   Annie Hughes   tidied up code, a...
256
257
258
259
260
261
262
263
if keyword_set(wait) then begin
   message,'Finished running DustEMWrap, using Niters: '+strtrim(string(use_Nitermax),2),/info
   message,'Time taken [sec]: '+sigfig(t2-t1,2,/sci),/info
   wait,wait
end



9e027842   Jean-Philippe Bernard   updated to new fu...
264
265
266
IF keyword_set(fits_save_and_restore) THEN BEGIN
   message,'Writing out results structure: '+fits_save_and_restore,/info
   dustem_write_fits_table,filename=fits_save_and_restore,help=help
8f9df3fb   Annie Hughes   tidied up code, a...
267
268
269
270
271
272
273
274
275
276
277
278
   dustem_read_fits_table,filename=fits_save_and_restore,dustem_st=dustem_spectra_st
                                ;==== plot result taken from the saved fits table
   res=*(*!dustem_fit).CURRENT_PARAM_VALUES
   IF !dustem_noobj THEN BEGIN
      dustemwrap_plot_noobj,res,dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,tit_m=tit_m+' (From Saved FITS file)',tit_x=tit_x+' (From Saved FITS file)'
   ENDIF ELSE BEGIN
      dustemwrap_plot,res,dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,tit_m=tit_m+' (From Saved FITS file)',tit_x=tit_x+' (From Saved FITS file)'
   ENDELSE
   IF keyword_set(wait) THEN BEGIN
      message,'Saved the results as FITS in the file: '+fits_save_and_restore+', and made a plot using the data in this file',/info
      wait,wait
   ENDIF
1d10fa6d   Ilyes Choubani   extinction exampl...
279
280
ENDIF

9e027842   Jean-Philippe Bernard   updated to new fu...
281
message,'Finished dustem_fit_sed_ext_pol_example',/info
1d10fa6d   Ilyes Choubani   extinction exampl...
282
283
284
285

the_end:

END