Blame view

src/idl/dustem_myisrf_example.pro 11.2 KB
e4285615   Ilyes Choubani   Uploading myisrf_...
1
PRO dustem_myisrf_example,model=model $
0dda35a7   Jean-Philippe Bernard   modified to inclu...
2
                         ,sed_file=sed_file $
b125205a   Annie Hughes   updated examples ...
3
                         ,isrf_file=isrf_file $
0dda35a7   Jean-Philippe Bernard   modified to inclu...
4
5
                         ,Nitermax=Nitermax $
                         ,postscript=postscript $
6735604d   Annie Hughes   making FITS read ...
6
                         ,fits_save=fits_save $
0dda35a7   Jean-Philippe Bernard   modified to inclu...
7
                         ,wait=wait $
a9349311   Annie Hughes   finished for tuto...
8
                         ,noobj=noobj $
0dda35a7   Jean-Philippe Bernard   modified to inclu...
9
10
                         ,verbose=verbose $
                         ,help=help
e4285615   Ilyes Choubani   Uploading myisrf_...
11
12
13

;+
; NAME:
0cca888a   Annie Hughes   added DPB90 as ex...
14
;    dustem_myisrf_example  
e4285615   Ilyes Choubani   Uploading myisrf_...
15
;
0cca888a   Annie Hughes   added DPB90 as ex...
16
17
18
19
; PURPOSE:
; This routine is an example of how to generate an alternative ISRF
; and then fit an observational SED (StokesI only) with DustEM and
; DustEMWrap using this non-default ISRF. The objective is to
b125205a   Annie Hughes   updated examples ...
20
21
22
; illustrate how to use dustem_plugin_modify_isrf and not to do
; science -- the fit obtained by running this example is likely to be
; poor.
e4285615   Ilyes Choubani   Uploading myisrf_...
23
;  
0cca888a   Annie Hughes   added DPB90 as ex...
24
; For this example, the code uses the SED in the file example_SED_1.xcat,
e4285615   Ilyes Choubani   Uploading myisrf_...
25
26
; which is distributed in the Data/EXAMPLE_OBSDATA/ directory 
;
e4285615   Ilyes Choubani   Uploading myisrf_...
27
28
29
30
; CATEGORY:
;    DustEMWrap, Distributed, High-Level, User Example
;
; CALLING SEQUENCE:
6735604d   Annie Hughes   making FITS read ...
31
;    dustem_myisrf_example[,model=][,sed_file=][,isrf_file=][,postscript=][,Nitermax=][,fits_save=][,/help,/wait,/verbose]
e4285615   Ilyes Choubani   Uploading myisrf_...
32
33
;
; INPUTS:
b125205a   Annie Hughes   updated examples ...
34
;  
e4285615   Ilyes Choubani   Uploading myisrf_...
35
36
37
38
39
40
41
42
43
44
; OPTIONAL INPUT PARAMETERS:
;    None
;
; OUTPUTS:
;    None
;
; OPTIONAL OUTPUT PARAMETERS:
;    Plots, results structure in binary FITS table format
;
; ACCEPTED KEY-WORDS:
6735604d   Annie Hughes   making FITS read ...
45
46
47
;    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.
e4285615   Ilyes Choubani   Uploading myisrf_...
48
49
;    sed_file = string naming the path to text file in .xcat format that
;          describes the observational SED. If not set, the file
0cca888a   Annie Hughes   added DPB90 as ex...
50
;          'Data/EXAMPLE_OBSDATA/example_SED_1.xcat' is used.  
5d4b92e3   Annie Hughes   finalised ISRF ex...
51
;    isrf_file = text file describing ISRF
e431104b   Annie Hughes   general updates f...
52
;    postscript = if set, final plot is saved as postscript file
e4285615   Ilyes Choubani   Uploading myisrf_...
53
;    Nitermax = maximum number of fit iterations. Default is 5.
6735604d   Annie Hughes   making FITS read ...
54
;    fits_save = if set, save the fit results in a binary
e431104b   Annie Hughes   general updates f...
55
;               FITS file. 
e4285615   Ilyes Choubani   Uploading myisrf_...
56
57
58
59
;    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
a9349311   Annie Hughes   finished for tuto...
60
;    noobj     = if set, runs with no object graphics
e4285615   Ilyes Choubani   Uploading myisrf_...
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
;
; 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
6735604d   Annie Hughes   making FITS read ...
76
;    dustem_myisrf_example,isrf_file='my_isrf.xcat',nitermax=10,fits_save='/tmp/mysavefile.fits'
e4285615   Ilyes Choubani   Uploading myisrf_...
77
78
79
80
81
82
83
84
;
; MODIFICATION HISTORY:
;    Evolution details on the DustEMWrap gitlab.
;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  
;-


IF keyword_set(help) THEN BEGIN
b125205a   Annie Hughes   updated examples ...
85
  doc_library,'dustem_myisrf_example'
e4285615   Ilyes Choubani   Uploading myisrf_...
86
87
88
89
90
91
  goto,the_end
END

IF keyword_set(model) THEN BEGIN
  use_model=strupcase(model)
ENDIF ELSE BEGIN
fabcd2dd   Annie Hughes   code now conforms...
92
  use_model='DBP90'    ;default dust model for this example
e4285615   Ilyes Choubani   Uploading myisrf_...
93
94
ENDELSE

6735604d   Annie Hughes   making FITS read ...
95
96
97
exists=dustem_test_model_exists(use_model)
if exists ne 1 then $
   message,'Unknown dust model'
e4285615   Ilyes Choubani   Uploading myisrf_...
98
99

use_verbose=0
b125205a   Annie Hughes   updated examples ...
100
101
use_polarization=0   ; initialize Dustemwrap in no polarization mode 
use_window=2          ; default graphics window number to use for plotting the results
e4285615   Ilyes Choubani   Uploading myisrf_...
102
103
104
105
106
107
108
109
110
111
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
;=== Refer to the DustEM and DustEMWrap userguides for an explanation
;    of the different grain types

fabcd2dd   Annie Hughes   code now conforms...
112
113
114
115
116
117
118
119
120
121
122
123
124
;=== AN EXAMPLE FOR DBP90
;=== Here we fit the PAH abundance of the model, and the amplitude
;=== of a user-supplied ISRF. A synchrotron component is included
;=== The other grain parameters are fixed to their default values.
;=== The free parameters are all lower-bounded at zero.
;=== Note that if '(*!dustem_params).G0' is not set explicitly as a fixed parameter, then
;=== the ISRF used by DustEM will include a standard Mathis field with
;=== amplitude, G0=1.
use_model='DBP90' ; you should specify this above, or in the command line
pd = [ $
     '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
     'dustem_plugin_modify_isrf_1', $           ; Amplitude of user-defined ISRF
     'dustem_plugin_synchrotron_2']             ;Synchrotron amplitude at 10 mm
a9349311   Annie Hughes   finished for tuto...
125
iv =   [5.4e-4, 0.1, 0.01]
5eebdedb   Jean-Philippe Bernard   changed plugin st...
126
iv =   [5.4e-4, 0.222222, 0.333333]
fabcd2dd   Annie Hughes   code now conforms...
127
128
129
130
Npar=n_elements(pd)
ulimed=replicate(0,Npar)
llimed=replicate(1,Npar)
llims=replicate(1.e-15,Npar)
a9349311   Annie Hughes   finished for tuto...
131
; To suppress the Mathis field altogether
fabcd2dd   Annie Hughes   code now conforms...
132
;fpd=['(*!dustem_params).G0']
a9349311   Annie Hughes   finished for tuto...
133
;fiv=[1.e-12]
fabcd2dd   Annie Hughes   code now conforms...
134
135


a9349311   Annie Hughes   finished for tuto...
136
;; ;=== AN EXAMPLE FOR DL07 
b125205a   Annie Hughes   updated examples ...
137
;; ;=== Here we fit the PAH abundances of the model, and the amplitude
0cca888a   Annie Hughes   added DPB90 as ex...
138
;; ;=== of a user-supplied ISRF. A synchrotron component is included
b125205a   Annie Hughes   updated examples ...
139
140
141
;; ;=== The other grain parameters are fixed to their default values.
;; ;=== The free parameters are all lower-bounded at zero.
;; use_model='DL07' ; you should specify this above, or in the command line
a9349311   Annie Hughes   finished for tuto...
142
143
144
145
146
147
148
149
150
151
152
;; pd = [ $
;;      '(*!dustem_params).grains(0).mdust_o_mh',$     ;PAH0 mass fraction
;;      '(*!dustem_params).grains(1).mdust_o_mh',$         ;PAH1 mass fraction
;;      'dustem_plugin_modify_isrf_1', $                   ; Ampitude of user-defined ISRF
;;      'dustem_plugin_synchrotron_2']                     ;Synchrotron amplitude at 10 mm
;; iv =   [5.4e-4, 5.4e-4, 0.1, 0.01]
;; Npar=n_elements(pd)
;; ulimed=replicate(0,Npar)
;; llimed=replicate(1,Npar)
;; llims=replicate(1.e-15,Npar)
;; ; To suppress the Mathis field altogether
6735604d   Annie Hughes   making FITS read ...
153
;; ;fpd=['(*!dustem_params).G0']
a9349311   Annie Hughes   finished for tuto...
154
;; ;fiv=[1.e-12]
e4285615   Ilyes Choubani   Uploading myisrf_...
155
156
157
158
159
160
161


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

5eebdedb   Jean-Philippe Bernard   changed plugin st...
162
dustem_init,model=use_model,polarization=use_polarization,/show_plot
e4285615   Ilyes Choubani   Uploading myisrf_...
163
!dustem_nocatch=1
b125205a   Annie Hughes   updated examples ...
164
!dustem_verbose=use_verbose
a9349311   Annie Hughes   finished for tuto...
165
IF keyword_set(noobj) THEN !dustem_noobj=1
b125205a   Annie Hughes   updated examples ...
166
167
168
169
170
171
172

;== Fill a structure with default inputs to the model. This
;== includes the default ISRF, which we want here just for plotting a
;comparison to the ISRF that we construct.
dir_in=!dustem_soft_dir  
st_model=dustem_read_all(dir_in)

0cca888a   Annie Hughes   added DPB90 as ex...
173
;=== CREATE A DIFFERENT ISRF USING DUSTEM_CREATE_RFIELD AND OUTPUT TO FILE
b125205a   Annie Hughes   updated examples ...
174
175
176
if not keyword_set(isrf_file) then begin
   mywaves=[]
   my_isrf_file='./myisrf_3bb_8000K.dat'
a9349311   Annie Hughes   finished for tuto...
177
   myisrf=dustem_create_rfield([8000,8000,8000],wdil=[1.d-11,1.d-11,1.d-11],isrf=0,fname=my_isrf_file,x=mywaves)
6735604d   Annie Hughes   making FITS read ...
178
;   myisrf=dustem_create_rfield([8000],wdil=[1.d-12],isrf=0,fname=my_isrf_file,x=mywaves)
b125205a   Annie Hughes   updated examples ...
179
180
181
182
183
184
185
186
187
188
;myisrf=dustem_create_rfield([3500,4000,5500,6500,8000],wdil=[1.d-13,1.d-13,1.d-13,1.d-13,1.d-13],isrf=1,fname='./myisrf_mathis_5bb.dat',x=mywaves)
;myisrf=dustem_create_rfield([3500,4000,5500,6500,8000],wdil=[1.d-13,1.d-13,1.d-13,1.d-13,1.d-13],isrf=2,fname='./myisrf_habing_5bb.dat',x=mywaves)
;myisrf=dustem_create_rfield([3500,4000,5500,6500,8000],wdil=[1.d-13,1.d-13,1.d-13,1.d-13,1.d-13],isrf=0,fname='./myisrf_5bb.dat',x=mywaves)
   isrf_file=my_isrf_file
endif


;=== PLOT AS A SANITY CHECK
xtit=textoidl('\lambda (\mum)')
ytit=textoidl('log ISRF [4 \pi I_\nu (erg/cm^2/s/Hz)]')
e431104b   Annie Hughes   general updates f...
189
tit='MYISRF EXAMPLE'
b125205a   Annie Hughes   updated examples ...
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
yr=[min(st_model.isrf.isrf) < min(myisrf),5.*(max(st_model.isrf.isrf)>max(myisrf))]
xr=[0.001,1e6]
window,use_window,xs=600,ys=400,tit=tit & use_window=use_window+2
cgplot,st_model.isrf.lambisrf,st_model.isrf.isrf $
       ,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog $
       ,title=tit,xtit=xtit,ytit=ytit,color=cgcolor('black'),/nodata
oplot,st_model.isrf.lambisrf,st_model.isrf.isrf,col=cgcolor('black'),thick=2
oplot,mywaves,myisrf,col=cgcolor('red'),thick=2
al_legend,/top,/right,clear=0,box=0 $
          ,['default ISRF','my ISRF'] $
          ,linsize=0.5,lines=0 $
          ,color=[cgcolor('black'),cgcolor('red')] $
          ,thick=2,charsize=1.3

;=== READ THE ISRF
a9349311   Annie Hughes   finished for tuto...
205
206
207
208
if keyword_set(isrf_file) then begin
   message,'Setting ISRF component from file: '+isrf_file,/info
   !dustem_isrf_file=ptr_new(isrf_file)
end
b125205a   Annie Hughes   updated examples ...
209
210

;=== READ EXAMPLE SED DATA: 
e4285615   Ilyes Choubani   Uploading myisrf_...
211
dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'
b125205a   Annie Hughes   updated examples ...
212
file=dir+'example_SED_1.xcat'
0dda35a7   Jean-Philippe Bernard   modified to inclu...
213
IF keyword_set(sed_file) THEN file=sed_file
e4285615   Ilyes Choubani   Uploading myisrf_...
214
215
216
217
218
219
220
221
sed=read_xcat(file,/silent)

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


b125205a   Annie Hughes   updated examples ...
222
223
224
;== SET THE OBSERVATIONAL STRUCTURE
;== sed is passed twice -- the first occurrence is the SED that you
;== wish to fit, the second occurrence is the SED that you wish to visualise. 
88872c38   Ilyes Choubani   small update. Cha...
225
dustem_set_data,m_fit=sed,m_show=sed
e4285615   Ilyes Choubani   Uploading myisrf_...
226
227
228

; ;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS THAT WILL BE
; ;== ADJUSTED DURING THE FIT
e4285615   Ilyes Choubani   Uploading myisrf_...
229
230
dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims

e4285615   Ilyes Choubani   Uploading myisrf_...
231
232
233
234
235
236
237

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

;== RUN THE FIT
b125205a   Annie Hughes   updated examples ...
238
tol=1.e-16
0dda35a7   Jean-Philippe Bernard   modified to inclu...
239
240
use_Nitermax=5  ;maximum number of iterations. 
IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax
e4285615   Ilyes Choubani   Uploading myisrf_...
241

b125205a   Annie Hughes   updated examples ...
242
243
xr = [1.,5e5]
yr = [5e-8,1.00e6]
e431104b   Annie Hughes   general updates f...
244
tit='FIT WITH MYISRF EXAMPLE'
e4285615   Ilyes Choubani   Uploading myisrf_...
245
246
ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
xtit=textoidl('\lambda (\mum)')
b125205a   Annie Hughes   updated examples ...
247

e4285615   Ilyes Choubani   Uploading myisrf_...
248
t1=systime(0,/sec)
5eebdedb   Jean-Philippe Bernard   changed plugin st...
249
250
251
252
;help,(*!dustem_plugin),/str
;help,(*!dustem_plugin).(0),/str
;help,(*!dustem_plugin).(1).scope
;stop
e4285615   Ilyes Choubani   Uploading myisrf_...
253
res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $
b125205a   Annie Hughes   updated examples ...
254
                      ,/xlog,/ylog,xr=xr,yr=yr,xtit=xtit,ytit=ytit,title=tit $
e4285615   Ilyes Choubani   Uploading myisrf_...
255
256
257
258
259
260
261
262
263
264
265
                      ,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

e4285615   Ilyes Choubani   Uploading myisrf_...
266
267
;=== MAKE THE FINAL PLOT
IF keyword_set(postscript) THEN BEGIN
e431104b   Annie Hughes   general updates f...
268
269
270
   mydevice=!d.name
   set_plot,'PS'
    ps_file=postscript
e4285615   Ilyes Choubani   Uploading myisrf_...
271
272
273
274
    device,filename=ps_file,/color
ENDIF

IF !dustem_noobj THEN BEGIN
6735604d   Annie Hughes   making FITS read ...
275
  dustemwrap_plot_noobj,*(*!dustem_fit).CURRENT_PARAM_VALUES,st=dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (Final fit)'
e4285615   Ilyes Choubani   Uploading myisrf_...
276
 ENDIF ELSE BEGIN
6735604d   Annie Hughes   making FITS read ...
277
  dustemwrap_plot,*(*!dustem_fit).CURRENT_PARAM_VALUES,st=dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (Final fit)'
e4285615   Ilyes Choubani   Uploading myisrf_...
278
279
280
ENDELSE

IF keyword_set(postscript) THEN BEGIN
e4285615   Ilyes Choubani   Uploading myisrf_...
281
  device,/close
e431104b   Annie Hughes   general updates f...
282
  set_plot,mydevice
e4285615   Ilyes Choubani   Uploading myisrf_...
283
284
285
286
287
288
289
290
  message,'Wrote '+ps_file,/info
ENDIF

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

6735604d   Annie Hughes   making FITS read ...
291
292
293
IF keyword_set(fits_save) THEN BEGIN
   message,'Writing out results structure: '+fits_save,/info
   dustem_write_fits_table,filename=fits_save,help=help
e4285615   Ilyes Choubani   Uploading myisrf_...
294
295
;=== 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
6735604d   Annie Hughes   making FITS read ...
296
297
298
299
300
301
302
303
304
305
306
307
308

; Moved following to dustem_fitsio_example
  ;; dustem_read_fits_table,filename=fits_save,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,title=tit+' (From Saved FITS file)'
  ;; ENDIF ELSE BEGIN
  ;;   dustemwrap_plot,res,dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
  ;; ENDELSE

   IF keyword_set(wait) THEN BEGIN
     message,'Saved the results as FITS in the file: '+fits_save,/info
e4285615   Ilyes Choubani   Uploading myisrf_...
309
310
311
312
     wait,wait
  ENDIF
ENDIF

e4285615   Ilyes Choubani   Uploading myisrf_...
313
314
315


the_end:
e431104b   Annie Hughes   general updates f...
316
message,'Finished dustem_myisrf_example',/info
e4285615   Ilyes Choubani   Uploading myisrf_...
317
318

END