OV-GSO-DC / dustem-wrapper_idl

Commit 5135b1945f75d60ee8a5b94854739b2343938e38

Authored by Jean-Philippe Bernard
1 parent 26361fb3
Exists in master

cleaned

Showing 1 changed file with 32 additions and 58 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_fit_sed_readme.pro
  Diff comments   View file @5135b19
1   -PRO dustem_fit_sed_readme,postcript=postcript,mode=mode,help=help,png=png
  1 +PRO dustem_fit_sed_readme,postcript=postcript,model=model,help=help,png=png,itermax=itermax
2 2  
3 3 ;This Readme describes how to fit SEDs with dustem
4 4 ;It runs on the SED stored into the file sample_SED.xcat
... ... @@ -25,7 +25,7 @@ PRO dustem_fit_sed_readme,postcript=postcript,mode=mode,help=help,png=png
25 25 ; CATEGORY:
26 26 ; Dustem
27 27 ; CALLING SEQUENCE:
28   -; dustem_fit_sed_readme,postcript=postcript,mode=mode,help=help
  28 +; dustem_fit_sed_readme,postcript=postcript,model=model,help=help
29 29 ; INPUTS:
30 30 ; None
31 31 ; OPTIONAL INPUT PARAMETERS:
... ... @@ -35,7 +35,7 @@ PRO dustem_fit_sed_readme,postcript=postcript,mode=mode,help=help,png=png
35 35 ; OPTIONAL OUTPUT PARAMETERS:
36 36 ; None
37 37 ; ACCEPTED KEY-WORDS:
38   -; mode = Selects one of the dust mixture used by dustem
  38 +; model = Selects one of the dust mixture used by dustem
39 39 ; 'COMPIEGNE_ETAL2010' from Compiegne et al 2010 (default)
40 40 ; 'DBP90' from Desert et al 1990
41 41 ; 'DL01' from Draine & Li 2001
... ... @@ -52,7 +52,7 @@ PRO dustem_fit_sed_readme,postcript=postcript,mode=mode,help=help,png=png
52 52 ; PROCEDURE:
53 53 ; None
54 54 ; EXAMPLES
55   -; dustem_fit_sed_readme,mode='COMPIEGNE_ETAL2010'
  55 +; dustem_fit_sed_readme,model='COMPIEGNE_ETAL2010'
56 56 ; MODIFICATION HISTORY:
57 57 ; Written by J.P. Bernard April 1st 2011
58 58 ; see evolution details on the dustem cvs maintained at CESR
... ... @@ -64,10 +64,10 @@ IF keyword_set(help) THEN BEGIN
64 64 goto,the_end
65 65 ENDIF
66 66  
67   -IF keyword_set(mode) THEN BEGIN
68   - use_mode=strupcase(mode)
  67 +IF keyword_set(model) THEN BEGIN
  68 + use_model=strupcase(model)
69 69 ENDIF ELSE BEGIN
70   - use_mode='COMPIEGNE_ETAL2010' ;Default is last dustem model
  70 + use_model='COMPIEGNE_ETAL2010' ;Default is last dustem model
71 71 ENDELSE
72 72  
73 73 IF keyword_set(png) THEN BEGIN
... ... @@ -79,7 +79,7 @@ ENDIF
79 79 ;=== initialise dustem
80 80 ;dustem_init,/wrap
81 81 ;dustem_init
82   -dustem_init,mode=use_mode
  82 +dustem_init,model=use_model
83 83 !dustem_verbose=1
84 84 !dustem_show_plot=1
85 85  
... ... @@ -99,21 +99,13 @@ ind=where(spec.error EQ 0.,count)
99 99 IF count NE 0 THEN spec(ind).error=0.2*spec(ind).spec
100 100 ind=where(spec.instru EQ 'FIRAS',count)
101 101 IF count NE 0 THEN spec(ind).error=0.2*spec(ind).spec
102   -;=== Draine Model SED for U=1. by C. Bot
103   -;;dir=getenv('DUSTEM_SOFT_DIR')+'/src/dustem3.8_web/SEDs/'
104   -;dir=!dustem_wrap_soft_dir+'/Data/SEDs/'
105   -;file=dir+'SED_DraineModel_U1.0.xcat'
106   -;spec=read_xcat(file,/silent)
107   -;=== Fake SED made from Dustem outputs
108   -;spec=dustem_compute_sed(p_dim,st=st,cont=cont,_extra=extra)
109 102  
110 103 ;=== Set which parameters you want to fit
111   -;CASE getenv('DUSTEM_WHICH') OF
112 104 CASE !dustem_which OF
113 105 'WEB3p8':BEGIN
114 106 ;=== This is to use Desert model
115 107 ;=== must have done dustem_init,/DBP90
116   - CASE use_mode OF
  108 + CASE use_model OF
117 109 'DBP90':BEGIN
118 110 pd = [ $
119 111 '(*!dustem_params).G0', $ ;G0
... ... @@ -193,6 +185,25 @@ CASE !dustem_which OF
193 185 ; fiv=[1.0]
194 186 ; dustem_init_fixed_params,fpd,fiv
195 187 END
  188 + 'AJ13':BEGIN
  189 + pd = [ $
  190 + '(*!dustem_params).G0', $ ;G0
  191 + '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
  192 + '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
  193 + '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx
  194 + '(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx
  195 + '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil
  196 + 'dustem_create_continuum_2'] ;Intensity of NIR continuum
  197 + iv = [1.0, 7.8e-4, 7.8e-4,1.65e-4,1.45e-3,7.8e-3,0.001]
  198 + Npar=n_elements(pd)
  199 + ulimed=replicate(0,Npar)
  200 + llimed=replicate(1,Npar)
  201 + llims=replicate(0.,Npar)
  202 + ;=== Fixed parameters
  203 +; fpd=['(*!dustem_params).G0'] ;Fixed parameter description
  204 +; fiv=[1.0]
  205 +; dustem_init_fixed_params,fpd,fiv
  206 + END
196 207 ENDCASE
197 208 END
198 209 ELSE: BEGIN
... ... @@ -217,54 +228,29 @@ ENDCASE
217 228  
218 229 st=dustem_set_data(sed=spec)
219 230  
220   -;VG : dustem_set_data turned into a function to be used indifferently for sed, ext, polext
221   -;dustem_set_data,spec
222   -;defsysv, '!dustem_data', {sed:ptr_new()} ;Data to fit
223   -
224 231 ;== SET THE FITTED PARAMETERS
225 232 dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
226 233 ;==== below is to reset the defaults keywords for PAH which in newer version of the fortran include spin and mix, ...
227 234 ;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn'
228 235 ;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn'
229   -(*!dustem_params).grains[0].TYPE_KEYWORDS='logn'
230   -(*!dustem_params).grains[1].TYPE_KEYWORDS='logn'
231 236 ;stop
232   -;dustem_sort_params
233 237 ;hprint,*(*!dustem_fit).param_descs
234 238 ;hprint,*(*!dustem_fit).param_init_values
235 239 ;hprint,*(*!dustem_fit).param_func
  240 +
236 241 ;=== RUN fit
237 242 tol=1.e-14
238   -;tol=1.e-2
239   -;Nitermax=20 ;maximum number of iteration. This is the criterium which will stop the fit procedure
240   -Nitermax=2 ;maximum number of iteration. This is the criterium which will stop the fit procedure
  243 +use_Nitermax=2 ;maximum number of iteration. This is the criterium which will stop the fit procedure
  244 +IF keyword_set(itermax) THEN use_Nitermax=itermax
241 245  
242 246 loadct,13
243 247 !y.range=[1e-4,10] ;This is to ajust plot range from outside the routine
244 248 t1=systime(0,/sec)
245   -res=dustem_mpfit_data(tol=tol,Nitermax=Nitermax,gtol=gtol,/xlog,/ylog)
246   -;res=dustem_mpfit_sed(tol=tol,Nitermax=Nitermax) ;,func_name=func_name,cf_min=cf_min)
  249 +res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol,/xlog,/ylog)
247 250 t2=systime(0,/sec)
248 251  
249   -;=== With FXD Model
250   -;Final parameters : 0.00012425071 0.00022732277 0.0021897348 0.00047272697 4.4916225
251   -;Iter 6 CHI-SQUARE = 76.345245 DOF = 6
252   -;=== With MC Model
253   -;Final parameters : 0.00012292305 0.00020672300 0.0020245400 0.00076911024 4.9751991
254   -;Iter 6 CHI-SQUARE = 76.047516 DOF = 6
255   -;=== With LV Model
256   -;Final parameters : 0.00013302325 0.0010706767 0.0039392709 0.00074149811 4.7549752
257   -;Iter 10 CHI-SQUARE = 90.478020 DOF = 6
258   -;=== With LV Model using GRAIN_DBP.DAT
259   -;Iter 7 CHI-SQUARE = 117.27229 DOF = 6
260   -;Final parameters : 0.00011163016 0.00015163894 0.0017526653 0.00096639410 5.8391213
261   -
262   -;stop
263   -
264 252 ;=== SAVE FIT RESULTS
265   -;file_out=getenv('DUSTEM_RES')+'DUSTEM_fit_example.sav'
266 253 file_out='/tmp/DUSTEM_fit_example.sav'
267   -;dustem_save_sed_fit,file_out
268 254 dustem_save_system_variables,file_out
269 255 message,'Saved '+file_out,/continue
270 256  
... ... @@ -277,9 +263,6 @@ message,'Saved '+file_out,/continue
277 263 file='/tmp/DUSTEM_fit_example.sav'
278 264 dustem_restore_system_variables,file
279 265  
280   -;file=!dustem_res+'DUSTEM_fit_example.sav'
281   -;dustem_restore_sed_fit,file
282   -
283 266 ;=== Plot best fit
284 267 yr=[1e-4,10]
285 268 xr=[1,2000]
... ... @@ -293,15 +276,6 @@ rchi2=(*!dustem_fit).rchi2
293 276 window,2
294 277  
295 278 ;=== RESTORE FIT RESULTS
296   -;=== initialise dustem
297   -;dustem_init,mode=use_mode
298   -;file=getenv('DUSTEM_RES')+'DUSTEM_fit_example.sav'
299   -
300   -;==== below is to reset the defaults keywords for PAH which in newer version of the fortran include spin and mix, ...
301   -;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn'
302   -;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn'
303   -
304   -;=== recover best fit values
305 279 res=*(*!dustem_fit).current_param_values
306 280 chi2=(*!dustem_fit).chi2
307 281 rchi2=(*!dustem_fit).rchi2
... ...