OV-GSO-DC / dustem-wrapper_idl

Commit 4caf33df772d49ca8f9b84dcf973bcfd620db688

Authored by Jean-Philippe Bernard
1 parent b0873df1
Exists in master

First commit

Showing 1 changed file with 235 additions and 0 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_fit_redshift_readme.pro 0 → 100644
  Diff comments   View file @4caf33d
... ... @@ -0,0 +1,235 @@
  1 +PRO dustem_fit_redshift_readme,postcript=postcript,model=model,help=help,png=png,itermax=itermax
  2 +
  3 +;This Readme describes how to fit SEDs with dustem
  4 +;It runs on the SED stored into the file sample_SED.xcat
  5 +;in this directory. Remember that the goal here is not necessarily to
  6 +;obtain a good fit in the end, but to illustrate the method.
  7 +;The provided SED has only photometric data points from SPITZER
  8 +;IRAC and MIPS and IRAS. No spectrum data points.
  9 +;SPECTRUM data points can be included and the corresponding filter
  10 +;filed must read SPECTRUM. Note that its is note necessary
  11 +;to use the .xcat file format, and data SED can be provided
  12 +;manually, but the observation structure must have the structure shown below.
  13 +;Obviously, the dustem package must have been installed succesfully (see
  14 +;dustem_cvs_readme.txt for install instructions).
  15 +
  16 +;+
  17 +; NAME:
  18 +; dustem_fit_sed_readme
  19 +; PURPOSE:
  20 +; This is an example of how to fit SEDs with the dustem wrapper.
  21 +; It is meant to be an example to follow when writing your own
  22 +; programs using the dustem IDL wrapper.
  23 +; It is not meant to reproduce the result in Compiegne et al 2010
  24 +; The SED used here is in sample_SED.xcat
  25 +; CATEGORY:
  26 +; Dustem
  27 +; CALLING SEQUENCE:
  28 +; dustem_fit_sed_readme,postcript=postcript,model=model,help=help
  29 +; INPUTS:
  30 +; None
  31 +; OPTIONAL INPUT PARAMETERS:
  32 +; None
  33 +; OUTPUTS:
  34 +; None
  35 +; OPTIONAL OUTPUT PARAMETERS:
  36 +; None
  37 +; ACCEPTED KEY-WORDS:
  38 +; model = Selects one of the dust mixture used by dustem
  39 +; 'COMPIEGNE_ETAL2010' from Compiegne et al 2010 (default)
  40 +; 'DBP90' from Desert et al 1990
  41 +; 'DL01' from Draine & Li 2001
  42 +; 'DL07' from Draine & Li 2007
  43 +; postcript = if set plot is done in DUSTEM/Docs/Figures/Last_dustem_fit.ps
  44 +; help = If set print this help
  45 +; COMMON BLOCKS:
  46 +; None
  47 +; SIDE EFFECTS:
  48 +; None
  49 +; RESTRICTIONS:
  50 +; The dustem fortran code must be installed
  51 +; The dustem idl wrapper must be installed
  52 +; PROCEDURE:
  53 +; None
  54 +; EXAMPLES
  55 +; dustem_fit_sed_readme,model='COMPIEGNE_ETAL2010'
  56 +; MODIFICATION HISTORY:
  57 +; Written by J.P. Bernard April 1st 2011
  58 +; see evolution details on the dustem cvs maintained at CESR
  59 +; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
  60 +;-
  61 +
  62 +IF keyword_set(help) THEN BEGIN
  63 + doc_library,'dustem_fit_redshift_readme'
  64 + goto,the_end
  65 +ENDIF
  66 +
  67 +IF keyword_set(model) THEN BEGIN
  68 + use_model=strupcase(model)
  69 +ENDIF ELSE BEGIN
  70 + use_model='COMPIEGNE_ETAL2010' ;Default is last dustem model
  71 +ENDELSE
  72 +
  73 +IF keyword_set(png) THEN BEGIN
  74 + dir_png=!dustem_dat+'/Figures/'
  75 + force_mkdir,dir_png
  76 +ENDIF
  77 +
  78 +;dbp=1
  79 +;=== initialise dustem
  80 +;dustem_init,/wrap
  81 +;dustem_init
  82 +dustem_init,model=use_model
  83 +!dustem_verbose=1
  84 +!dustem_show_plot=1
  85 +
  86 +;=== Read sample SED
  87 +;;dir=getenv('DUSTEM_SOFT_DIR')+'/src/'
  88 +;dir=!dustem_wrap_soft_dir+'/Data/SEDs/'
  89 +;file=dir+'sample_SED.xcat'
  90 +;spec=read_xcat(file,/silent)
  91 +;=== Composite SED from Compiegne et al 2010, gathered by C. Bot
  92 +;dir=getenv('DUSTEM_SOFT_DIR')+'/src/dustem3.8_web/SEDs/'
  93 +dir=!dustem_wrap_soft_dir+'/Data/SEDs/workshop/'
  94 +file=dir+'GN26_SED.xcat'
  95 +;file=dir+'Gal_composite_spectrum_wAKARI.xcat'
  96 +;stop
  97 +spec=read_xcat(file,/silent)
  98 +ind=where(spec.error EQ 0.,count)
  99 +IF count NE 0 THEN spec(ind).error=0.2*spec(ind).spec
  100 +ind=where(spec.instru EQ 'FIRAS',count)
  101 +IF count NE 0 THEN spec(ind).error=0.2*spec(ind).spec
  102 +
  103 +;=== Set which parameters you want to fit
  104 +CASE !dustem_which OF
  105 + 'WEB3p8':BEGIN
  106 +;=== This is to use Desert model
  107 +;=== must have done dustem_init,/DBP90
  108 + CASE use_model OF
  109 + 'COMPIEGNE_ETAL2010':BEGIN
  110 +; PUTTING G0 last DOES NOT WORK BECAUSE OF ORDER ISSUE
  111 + pd = ['!dustem_redshift', $
  112 + '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
  113 + '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
  114 + '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx
  115 + '(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx
  116 + '(*!dustem_params).grains(4).mdust_o_mh'] ;aSil
  117 + iv = [1.0, 7.8e-4, 7.8e-4,1.65e-4,1.45e-3,7.8e-3]
  118 + Npar=n_elements(pd)
  119 + ulimed=replicate(0,Npar)
  120 + llimed=replicate(1,Npar)
  121 + llims=replicate(0.,Npar)
  122 + ;=== Fixed parameters
  123 + fpd=['(*!dustem_params).gas.G0'] ;Fixed parameter description
  124 + fiv=[53.1]
  125 + dustem_init_fixed_params,fpd,fiv
  126 + END
  127 + ENDCASE
  128 + END
  129 + ELSE: BEGIN
  130 + pd = ['(*!dustem_params).grains(0).propmass',$ ;PAH mass fraction
  131 + '(*!dustem_params).grains(1).propmass',$ ;VSG mass fraction
  132 + '(*!dustem_params).grains(2).propmass',$ ;BG mass fraction
  133 + 'dustem_create_continuum_2', $ ;Intensity of NIR continuum
  134 + 'dustem_create_isrf_1'] ;ISRF scaling factor
  135 + ;=== set initial value of parameters you want to fit
  136 + iv = [7e-5, 8e-4,6.6e-3,0.001,4.]
  137 + ;in this example, we enforce parameters to be >0 and Xisrf >3
  138 + ;Caution, initial values must be within limits ...
  139 + Npar=n_elements(pd)
  140 + ulimed=replicate(0,Npar)
  141 + llimed=replicate(1,Npar)
  142 + llims=replicate(0.,Npar)
  143 + END
  144 +ENDCASE
  145 +
  146 +!dustem_redshift=1.22
  147 +
  148 +;=== SET THE OBSERVATION STRUCTURE
  149 +
  150 +st=dustem_set_data(sed=spec)
  151 +
  152 +;== SET THE FITTED PARAMETERS
  153 +dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
  154 +;==== below is to reset the defaults keywords for PAH which in newer version of the fortran include spin and mix, ...
  155 +;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn'
  156 +;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn'
  157 +
  158 +;=== RUN fit
  159 +tol=1.e-14
  160 +use_Nitermax=2 ;maximum number of iteration. This is the criterium which will stop the fit procedure
  161 +IF keyword_set(itermax) THEN use_Nitermax=itermax
  162 +
  163 +loadct,13
  164 +yrange=[0.1,1.e4]
  165 +xrange=[10,10000]
  166 +tit='DUSTEM Intensity SED fit Example'
  167 +ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
  168 +xtit=textoidl('\lambda (\mum)')
  169 +legend_xpos=0.6
  170 +legend_ypos=0.8
  171 +t1=systime(0,/sec)
  172 +res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol,/xlog,/ylog,yrange=yrange,xrange=xrange,/xsty,/ysty,title=tit,xtitle=xtit,ytitle=ytit,legend_xpos=legend_xpos,legend_ypos=legend_ypos)
  173 +t2=systime(0,/sec)
  174 +
  175 +;=== SAVE FIT RESULTS
  176 +file_out='/tmp/DUSTEM_fit_example.sav'
  177 +dustem_save_system_variables,file_out
  178 +message,'Saved '+file_out,/continue
  179 +
  180 +;stop
  181 +
  182 +;======================================
  183 +;====You can exit IDL here and re-enter
  184 +;======================================
  185 +
  186 +file='/tmp/DUSTEM_fit_example.sav'
  187 +dustem_restore_system_variables,file
  188 +
  189 +;=== Plot best fit
  190 +tit='DUSTEM Intensity SED fit Example (Saved)'
  191 +ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
  192 +xtit=textoidl('\lambda (\mum)')
  193 +errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values)
  194 +chi2=(*!dustem_fit).chi2
  195 +rchi2=(*!dustem_fit).rchi2
  196 +
  197 +window,2
  198 +
  199 +;=== RESTORE FIT RESULTS
  200 +res=*(*!dustem_fit).current_param_values
  201 +chi2=(*!dustem_fit).chi2
  202 +rchi2=(*!dustem_fit).rchi2
  203 +errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values)
  204 +
  205 +;=== Plot best fit
  206 +loadct,13
  207 +IF keyword_set(postcript) THEN BEGIN
  208 + set_plot,'PS'
  209 +; ps_file=getenv('DUSTEM_SOFT_DIR')+'/Docs/Figures/'+'Last_dustem_fit.ps'
  210 + ps_file=!dustem_wrap_soft_dir+'/Docs/Figures/'+'Last_dustem_fit.ps'
  211 + device,filename=ps_file,/color
  212 +ENDIF
  213 +legend_xpos=1000
  214 +legend_ypos=3
  215 +
  216 +;dustem_sed_plot,(*!dustem_fit_params),ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2
  217 +dustem_sed_plot,*(*!dustem_fit).current_param_values,ytit=ytit,xtit=xtit,title=tit,yrange=yrange,xrange=xrange,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog, $
  218 + legend_xpos=legend_xpos,legend_ypos=legend_ypos
  219 +IF keyword_set(postcript) THEN BEGIN
  220 + device,/close
  221 + set_plot,'X'
  222 + message,'wrote '+ps_file,/info
  223 +ENDIF
  224 +IF keyword_set(png) THEN BEGIN
  225 + ;file_png=!dustem_wrap_soft_dir+'/Docs/Figures/'+'Last_dustem_fit.png'
  226 + file_png=dir_png+'Last_dustem_fit_sed_readme.png'
  227 + write_png,file_png,tvrd(/true)
  228 + message,'Wrote '+file_png,/info
  229 +ENDIF
  230 +
  231 +message,'dustem_mpfit_sed executed in '+strtrim(t2-t1,2)+' sec',/info
  232 +
  233 +the_end:
  234 +
  235 +END
... ...