Blame view

src/idl/dustem_fit_intensity_phangs_example.pro 10.2 KB
6ff7d98f   Jean-Philippe Bernard   First commit
1
2
3
4
5
6
7
8
9
10
11
12
PRO dustem_fit_intensity_phangs_example,model=model $
                                ,sed_file=sed_file $
                                ,Nitermax=Nitermax $
                                ,postscript=postscript $
                                ,fits_save=fits_save $
                                ,help=help $
                                ,wait=wait $
                                ,noobj=noobj $
                                ,verbose=verbose

;+
; NAME:
88f84272   Jean-Philippe Bernard   changed
13
;    dustem_fit_intensity_phangs_example  
6ff7d98f   Jean-Philippe Bernard   First commit
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
;
; PURPOSE:
; This routine is an example of how to fit an observational SED
; (StokesI only) with DustEM and DustEMWrap. The objective is to
; illustrate how to use DustEMWrap and not to do science -- the fit
; obtained by running this example is likely to be poor.
;  
; For this example, the code uses the SED in the file example_SED_1.xcat,
; which is distributed in the Data/EXAMPLE_OBSDATA/ directory 
;
; The example SED has Stokes I photometric data points from
; IRAC, MIPS and IRAS. Examples illustrating running DustEMWrap to
; fit spectral data, polarisation data and extinction data
; are provided in other _example routines in the src/idl/
; directory. See the DustEMWrap User Guide for more information.
;
; CATEGORY:
;    DustEMWrap, Distributed, High-Level, User Example
;
; CALLING SEQUENCE:
;    dustem_fit_intensity_example[,model=][sed_file=][,postscript=][,Nitermax=][,fits_save=][,/help,/wait,/verbose]
;
; 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. See userguide or dustem_test_model_exists.pro
;            for more details about available models in current release.
;    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_1.xcat' is used.  
;    postscript = if set, final plot is saved as postscript file
;    Nitermax = maximum number of fit iterations. Default is 5.
;    fits_save = if set, save the fit results in a binary
;               FITS file. 
;    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
;    noobj     = if set, runs with no object graphics
;
; 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:
;    
; EXAMPLES
;    dustem_fit_intensity_example
;    dustem_fit_intensity_example,Nitermax=1,fits_save='/tmp/mysavefile.fits'
;    dustem_fit_intensity_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.  
;-

IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_fit_intensity_phangs_example'
  goto,the_end
END

IF keyword_set(model) THEN BEGIN
  use_model=strupcase(model)
ENDIF ELSE BEGIN
  use_model='DBP90'    ;Example with default keywords uses the DBP90 model
ENDELSE

exists=dustem_test_model_exists(use_model)
if exists ne 1 then $
   message,'Unknown dust model'

use_polarization=0   ; initialize Dustemwrap in no polarization mode 
use_window=2          ; default graphics window number to use for plotting the results
use_verbose=0
if keyword_set(verbose) then use_verbose=1
use_Nitermax=5        ; maximum number of iterations for the fit
IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax

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.

;=== Examples are provided for some of the dust models. 
;=== To try them, uncomment the model that you want to try and re-run

;=== AN EXAMPLE FOR DBP90
;=== Here we fit the dust abundances of the DBP90 model, the
;=== intensity of the dust-heating radiation field, as well as a plug-in
;=== for synchrotron emission
;=== The free parameters are all lower-bounded at zero.
;=== use_model='DBP90' ; you should specify this above!

age_values=[0.03, 0.05, 0.08, 0.15, 0.25, 0.40, 0.60, 1.0, 1.75, 3.0, 5.0, 8.5, 13.5]
metalicity_values=[-1.49, -0.96, -0.35, +0.06, +0.26, +0.4]
Nage=n_elements(age_values)
Nmetalicity=n_elements(metalicity_values)
ii=where(age_values EQ 0.2500,counta)
jj=where(metalicity_values EQ -0.35,countb)
ij=lonarr(1,2)
ij[0,0]=ii
ij[0,1]=jj
ind=ij2index(ij,[Nage,Nmetalicity])+1
pd = [ $
     '(*!dustem_params).G0', $                           ;G0
     '(*!dustem_params).grains(0).mdust_o_mh',$          ;PAH0 mass fraction
     '(*!dustem_params).grains(1).mdust_o_mh',$          ;VSG mass fraction
     '(*!dustem_params).grains(2).mdust_o_mh',$          ;BG mass fraction
     'dustem_plugin_phangs_stellar_continuum_0']         ;phangs tau
iv =   [1.6, 2.2e-4, 5.7e-4, 3.4e-3,1.] 

pd = [ $
     '(*!dustem_params).G0', $                           ;G0
     '(*!dustem_params).grains(0).mdust_o_mh',$          ;PAH0 mass fraction
     '(*!dustem_params).grains(1).mdust_o_mh',$          ;VSG mass fraction
     '(*!dustem_params).grains(2).mdust_o_mh',$          ;BG mass fraction
     'dustem_plugin_phangs_stellar_continuum_'+strtrim(ind[0],2)]


; initial values vector for run without synchrotron plugin
iv =   [1.6, 2.2e-4, 5.7e-4, 3.4e-3,100.] 

6ff7d98f   Jean-Philippe Bernard   First commit
159
160
161
162
163
164
Npar=n_elements(pd)
ulimed=replicate(0,Npar)
llimed=replicate(1,Npar) 
llims=replicate(1.e-15,Npar)
fpd=[] & fiv=[]

6ff7d98f   Jean-Philippe Bernard   First commit
165
166
167
168
169
170
171

if keyword_set(wait) then begin
   message,'Finished setting dust model and plug-in parameters: '+use_model,/info
   wait,wait
end

;== INITIALISE DUSTEM
88f84272   Jean-Philippe Bernard   changed
172
173
174
;!dustem_show_plot=1
dustem_init,model=use_model,polarization=use_polarization,/show_plot
;stop
6ff7d98f   Jean-Philippe Bernard   First commit
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
!dustem_nocatch=1
!dustem_verbose=use_verbose
IF keyword_set(noobj) THEN !dustem_noobj=1
!EXCEPT=2 ; for debugging

;=== READ EXAMPLE SED DATA
dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'
file=dir+'example_SED_5.xcat'
IF keyword_set(sed_file) THEN file=sed_file
sed=read_xcat(file,/silent)

if keyword_set(wait) then begin
   message,'Finished reading SED data: '+file,/info
   wait,wait
end

;;=== ADJUST THE UNCERTAINTIES FROM WITHIN THE CODE (FOR ILLUSTRATION)
ind=where(sed.sigmaII LT (0.2*sed.StokesI)^2,count)
IF count NE 0 THEN sed[ind].sigmaII=(0.2*sed[ind].StokesI)^2

;== SET THE OBSERVATIONAL STRUCTURE
;== sed is passed twice in the call -- the first occurrence is the SED that you
;== wish to fit, the second occurrence is the SED that you wish to visualise. 
dustem_set_data,m_fit=sed,m_show=sed

;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS THAT WILL BE
;== ADJUSTED DURING THE FIT
dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims


if keyword_set(wait) then begin
   message,'Finished initializing DustEMWrap, including plugins and fixed parameters',/info
   wait,wait
end

;=== INFORMATION TO RUN THE FIT
tol=1.e-16     ;fit tolerence

;=== INFORMATION TO MAKE THE PLOT
yr=[1.00e-4,1.00E2] ; y-axis limits
xr=[1.00E-1,6.00e4] ; x-axis limits
tit='FIT INTENSITY EXAMPLE' ; plot title
ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') ; y-axis title
xtit=textoidl('\lambda (\mum)') ; x-axis title


;===  RUN THE FIT
t1=systime(0,/sec)
res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $
                      ,/xlog,/ylog,xr=xr,yr=yr,xtit=xtit,ytit=ytit,title=tit $
                      ,legend_xpos=legend_xpos,legend_ypos=legend_ypos $
                      ,errors=errors,chi2=chi2,rchi2=rchi2,show_plot=show_plot)
t2=systime(0,/sec)

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

;=== MAKE THE FINAL PLOT
IF keyword_set(postscript) THEN BEGIN
;    dir_ps='./'
    mydevice=!d.name
    set_plot,'PS'
;    ps_file=dir_ps+postscript
    ps_file=postscript
    device,filename=ps_file,/color
ENDIF

;stop

IF !dustem_noobj THEN BEGIN
  dustemwrap_plot_noobj,*(*!dustem_fit).CURRENT_PARAM_VALUES,st=dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (Final fit)'
 ENDIF ELSE BEGIN
  dustemwrap_plot,*(*!dustem_fit).CURRENT_PARAM_VALUES,st=dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (Final fit)'
ENDELSE

IF keyword_set(postscript) THEN BEGIN
  device,/close
  set_plot,mydevice
  message,'Wrote '+ps_file,/info
ENDIF

if keyword_set(wait) then begin
   message,'Made the plot of the final results',/info
   wait,wait
end

IF keyword_set(fits_save) THEN BEGIN
   message,'Writing out results structure: '+fits_save,/info
   dustem_write_fits_table,filename=fits_save,help=help
;=== At this point, you could erase all dustem system variables, or exit idl... all the
;=== information needed to recover the results and remake the plots has been saved in the FITS table

;=== Moved the following to dustem_fits_io_example
  ;; dustem_read_fits_table,filename=fits_save,dustem_st=dustem_spectra_st
  ;; ;stop
  ;; ;!dustem_show=ptr_new(*!dustem_data)   ;test
  ;; !dustem_noobj=1
  ;; ;==== plot result taken from the saved fits table
  ;; res=*(*!dustem_fit).CURRENT_PARAM_VALUES
  ;; IF !dustem_noobj THEN BEGIN
  ;;   ;dustemwrap_plot_noobj,res,st=dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
  ;;   dustemwrap_plot_noobj,res,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
  ;; ENDIF ELSE BEGIN
  ;;   ;dustemwrap_plot,res,st=dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
  ;;   dustemwrap_plot,res,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
  ;; ENDELSE
  ;; ;stop

  IF keyword_set(wait) THEN BEGIN
     message,'Saved the results as FITS in the file: '+fits_save,/info
     wait,wait
  ENDIF
ENDIF


the_end:
message,'Finished dustem_fit_intensity_example',/info

END