OV-GSO-DC / dustem-wrapper_idl

Commit 906f3c27c86fd44cb347e92baeef51ef0eac525c

Authored by Jean-Philippe Bernard
1 parent 0d9e4f2f
Exists in master

First commit

Showing 1 changed file with 202 additions and 0 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_fit_spin_test_readme.pro 0 → 100644
  Diff comments   View file @906f3c2
... ... @@ -0,0 +1,202 @@
  1 +PRO dustem_fit_spin_test_readme,postcript=postcript,mode=mode,help=help
  2 +
  3 +;This Readme describes how to fit SEDs with a modified Black-Body
  4 +;It runs on the SED stored into the file sample_SED.xcat
  5 +;Remember that the goal here is just to illustrate the method.
  6 +;Note that its is not necessary to use the .xcat file format, and data SED can be provided
  7 +;manually, but the observation structure must have the structure shown below.
  8 +;Obviously, the dustem package must have been installed succesfully
  9 +;see DustemWrap documentation for install instructions.
  10 +
  11 +;+
  12 +; NAME:
  13 +; dustem_fit_modified_bb_readme
  14 +; PURPOSE:
  15 +; This is an example of how to fit SEDs with a modified Black-Body
  16 +; using the dustem wrapper.
  17 +; It is meant to be an example to follow when writing your own
  18 +; programs using the dustem IDL wrapper.
  19 +; The SED used here is in sample_SED.xcat
  20 +; CATEGORY:
  21 +; Dustem
  22 +; CALLING SEQUENCE:
  23 +; dustem_fit_modified_bb_readme,postcript=postcript,help=help
  24 +; INPUTS:
  25 +; None
  26 +; OPTIONAL INPUT PARAMETERS:
  27 +; None
  28 +; OUTPUTS:
  29 +; None
  30 +; OPTIONAL OUTPUT PARAMETERS:
  31 +; None
  32 +; ACCEPTED KEY-WORDS:
  33 +; postcript = if set plot is done in DUSTEM/Docs/Figures/Last_dustem_fit.ps
  34 +; help = If set print this help
  35 +; COMMON BLOCKS:
  36 +; None
  37 +; SIDE EFFECTS:
  38 +; None
  39 +; RESTRICTIONS:
  40 +; The dustem idl wrapper must be installed
  41 +; PROCEDURE:
  42 +; None
  43 +; EXAMPLES
  44 +; dustem_fit_modified_bb_readme
  45 +; MODIFICATION HISTORY:
  46 +; Written by J.P. Bernard June 9th 2011
  47 +; see evolution details on the dustem cvs maintained at CESR
  48 +; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
  49 +;-
  50 +
  51 +
  52 +IF keyword_set(help) THEN BEGIN
  53 + doc_library,'dustem_fit_spin_test_readme'
  54 + goto,the_end
  55 +ENDIF
  56 +
  57 +;=== initialise dustem
  58 +dustem_init,mode='COMPIEGNE_ETAL2010'
  59 +!dustem_verbose=1
  60 +!dustem_show_plot=1
  61 +
  62 +;=== make a fake one out for given filters
  63 +filters=['IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3']
  64 +Nfilt=n_elements(filters)
  65 +sed=dustem_initialize_sed(Nfilt)
  66 +sed.filter=filters
  67 +sed.wave=dustem_filter2wav(filters)
  68 +sed.instru=dustem_filter2instru(filters)
  69 +st=dustem_set_data(sed=sed)
  70 +
  71 +;1.000000
  72 +;PAH0_MC10 25 logn-chrg-spin-zm 7.8000E-04 2.2400E+00 3.5000E-08 1.2000E-07 6.4000E-08 1.0000E-01
  73 +;PAH1_MC10 25 logn-chrg-spin-zm 7.8000E-04 2.2400E+00 3.5000E-08 1.2000E-07 6.4000E-08 1.0000E-01
  74 +;amCBEx 15 logn 1.6500E-04 1.8100E+00 6.0000E-08 2.0000E-06 2.0000E-07 3.5000E-01
  75 +;amCBEx 25 plaw-ed 1.4500E-03 1.8100E+00 4.0000E-07 2.0000E-04 -2.8000E+00 1.5000E-05 1.5000E-05 2.0000E+00
  76 +;aSilx 25 plaw-ed 6.7000E-03 3.0000E+00 4.0000E-07 2.0000E-04 -3.4000E+00 2.0000E-05 2.0000E-05 2.0000E+00
  77 +
  78 +;=== Set which parameters you want to fit
  79 +pd = [ $
  80 +; '(*!dustem_params).G0', $ ;G0
  81 + '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
  82 + '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
  83 + '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx
  84 + '(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx
  85 + '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSilx
  86 + '(*!dustem_params).gas.G0', $ ;Gas G0
  87 + '(*!dustem_params).gas.Tgas', $ ;Gas temperature
  88 + '(*!dustem_params).gas.nH', $ ;Gas temperature
  89 + 'dustem_create_continuum_2'] ;Intensity of NIR continuum
  90 +
  91 +p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,100.,1000.,1000.,0.001]
  92 +
  93 +;=== set initial value of parameters you want to fit
  94 +iv=p_truth ;exact solution
  95 +iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,90.,500.,600.,0.001] ;shifted
  96 +
  97 +ulimed=[0 , 0 ,0,0,0,0,0,0,0]
  98 +llimed=[1 , 1 ,1,1,1,1,1,1,1]
  99 +llims =[0. ,0. ,0.,0.,0.,0.,0.,0.,0.]
  100 +
  101 +;== SET THE FITTED PARAMETERS (now done in dustem_init_parinfo)
  102 +dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
  103 +
  104 +(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-spin'
  105 +(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-spin'
  106 +
  107 +help,*!dustem_fit
  108 +
  109 +;stop
  110 +
  111 +sed.spec=dustem_compute_sed(p_truth,sst,_extra=extra)
  112 +sed.error=(sed.spec*0.1);>0.1
  113 +;sed.error[*]=100.
  114 +
  115 +;=== SET THE OBSERVATION STRUCTURE
  116 +
  117 +ind=where(sed.wave GE 5.)
  118 +st=dustem_set_data(sed=sed[ind])
  119 +
  120 +
  121 +;=== RUN fit
  122 +;tol=1.e-20
  123 +tol=1.e-10
  124 +xtol=1.e-10
  125 +Nitermax=5 ;maximum number of iteration. This is the criterium which will stop the fit procedure
  126 +
  127 +loadct,13
  128 +;!y.range=[1e-8,100] ;This is to ajust plot range from outside the routine
  129 +t1=systime(0,/sec)
  130 +res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=Nitermax,xrange=[1.,2.e4],/xstyle,yrange=[1e-4,1.e3],/ysty,/ylog,/xlog,xtit='wavelength [mic]',ytit='Brightness []')
  131 +t2=systime(0,/sec)
  132 +
  133 +print,(res-p_truth)/p_truth*100.
  134 +
  135 +stop
  136 +
  137 +;=== SAVE FIT RESULTS
  138 +;file_out=getenv('DUSTEM_RES')+'DUSTEM_fit_example.sav'
  139 +file_out=!dustem_res+'DUSTEM_bb_fit_example.sav'
  140 +dustem_save_sed_fit,file_out
  141 +
  142 +;=== Plot best fit
  143 +yr=[1e-3,100]
  144 +xr=[50,4000]
  145 +tit='DUSTEM Modified BB Example'
  146 +ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
  147 +xtit=textoidl('\lambda (\mum)')
  148 +;errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values)
  149 +errors=(*(*!dustem_fit).current_param_errors)
  150 +chi2=(*!dustem_fit).chi2
  151 +rchi2=(*!dustem_fit).rchi2
  152 +
  153 +loadct,13
  154 +dustem_sed_plot,(*(*!dustem_fit).current_param_values), $
  155 + ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty, $
  156 + res=res,errors=errors,chi2=chi2,rchi2=rchi2,function_name='dustem_greybody_mpfit'
  157 +
  158 +;stop
  159 +
  160 +;======================================
  161 +;====You can exit IDL here and re-enter
  162 +;======================================
  163 +window,1
  164 +
  165 +;=== RESTORE FIT RESULTS
  166 +;=== initialise dustem
  167 +dustem_init
  168 +file=!dustem_res+'DUSTEM_bb_fit_example.sav'
  169 +dustem_restore_sed_fit,file
  170 +
  171 +;=== recover best fit values
  172 +res=*(*!dustem_fit).current_param_values
  173 +chi2=(*!dustem_fit).chi2
  174 +rchi2=(*!dustem_fit).rchi2
  175 +errors=*(*!dustem_fit).current_param_errors
  176 +
  177 +;=== Plot best fit
  178 +yr=[1e-3,100]
  179 +xr=[50,4000]
  180 +tit='DUSTEM Modified BB Example (Saved)'
  181 +ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
  182 +xtit=textoidl('\lambda (\mum)')
  183 +
  184 +loadct,13
  185 +IF keyword_set(postcript) THEN BEGIN
  186 + set_plot,'PS'
  187 + ps_file=!dustem_wrap_soft_dir+'/Docs/Figures/'+'Last_dustem_bb_fit.ps'
  188 + device,filename=ps_file,/color
  189 +ENDIF
  190 +dustem_sed_plot,*(*!dustem_fit).current_param_values, $
  191 + ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,function_name='dustem_greybody_mpfit'
  192 +IF keyword_set(postcript) THEN BEGIN
  193 + device,/close
  194 + set_plot,'X'
  195 + message,'wrote '+ps_file,/info
  196 +ENDIF
  197 +
  198 +print,'dustem_mpfit_sed executed in ',t2-t1,' sec'
  199 +
  200 +the_end:
  201 +
  202 +END
... ...