PRO dustem_fit_sed_readme,postcript=postcript,model=model,help=help,png=png,itermax=itermax ;This Readme describes how to fit SEDs with dustem ;It runs on the SED stored into the file sample_SED.xcat ;in this directory. Remember that the goal here is not necessarily to ;obtain a good fit in the end, but to illustrate the method. ;The provided SED has only photometric data points from SPITZER ;IRAC and MIPS and IRAS. No spectrum data points. ;SPECTRUM data points can be included and the corresponding filter ;filed must read SPECTRUM. Note that its is note necessary ;to use the .xcat file format, and data SED can be provided ;manually, but the observation structure must have the structure shown below. ;Obviously, the dustem package must have been installed succesfully (see ;dustem_cvs_readme.txt for install instructions). ;+ ; NAME: ; dustem_fit_sed_readme ; PURPOSE: ; This is an example of how to fit SEDs with the dustem wrapper. ; It is meant to be an example to follow when writing your own ; programs using the dustem IDL wrapper. ; It is not meant to reproduce the result in Compiegne et al 2010 ; The SED used here is in sample_SED.xcat ; CATEGORY: ; Dustem ; CALLING SEQUENCE: ; dustem_fit_sed_readme,postcript=postcript,model=model,help=help ; INPUTS: ; None ; OPTIONAL INPUT PARAMETERS: ; None ; OUTPUTS: ; None ; OPTIONAL OUTPUT PARAMETERS: ; None ; ACCEPTED KEY-WORDS: ; model = Selects one of the dust mixture used by dustem ; 'COMPIEGNE_ETAL2010' from Compiegne et al 2010 (default) ; 'DBP90' from Desert et al 1990 ; 'DL01' from Draine & Li 2001 ; 'DL07' from Draine & Li 2007 ; postcript = if set plot is done in DUSTEM/Docs/Figures/Last_dustem_fit.ps ; help = If set print this help ; COMMON BLOCKS: ; None ; SIDE EFFECTS: ; None ; RESTRICTIONS: ; The dustem fortran code must be installed ; The dustem idl wrapper must be installed ; PROCEDURE: ; None ; EXAMPLES ; dustem_fit_sed_readme,model='COMPIEGNE_ETAL2010' ; MODIFICATION HISTORY: ; Written by J.P. Bernard April 1st 2011 ; see evolution details on the dustem cvs maintained at CESR ; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems. ;- IF keyword_set(help) THEN BEGIN doc_library,'dustem_fit_sed_readme' goto,the_end ENDIF IF keyword_set(model) THEN BEGIN use_model=strupcase(model) ENDIF ELSE BEGIN use_model='THEMIS';G17_MODELA';'COMPIEGNE_ETAL2010' ;Default is last dustem model ENDELSE IF keyword_set(png) THEN BEGIN dir_png=!dustem_dat+'/Figures/' force_mkdir,dir_png ENDIF polarization=0. ;default is no polarization in models ;dustem_init,model=use_model,pol=polarization ;=== Set which parameters you want to fit CASE use_model OF 'DBP90':BEGIN pd = [ $ '(*!dustem_params).G0', $ ;G0 '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum iv = [1.0, 4.3e-4, 4.7e-4,6.4e-3,0.001] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) END 'DL01':BEGIN pd = [ $ '(*!dustem_params).G0', $ ;G0 '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction '(*!dustem_params).grains(2).mdust_o_mh', $ ;Gra '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum iv = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) END 'DL07':BEGIN pd = [ $ '(*!dustem_params).G0', $ ;G0 '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction '(*!dustem_params).grains(2).mdust_o_mh', $ ;Gra '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum ;] ; pd = [ $ ; '(*!dustem_params).gas.G0', $ ;G0 ; '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction ; '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction ; '(*!dustem_params).grains(2).mdust_o_mh', $ ;Gra ; '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra ; '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil ; 'dustem_create_continuum_2', $ ; 'dustem_create_stellar_population_O63',$ ; distance to the O3 stellar population ; 'dustem_create_stellar_population_O65',$ ; number of stars of the O3 stellar population ; 'dustem_create_stellar_population_B43',$ ; distance to the B4 stellar population ; 'dustem_create_stellar_population_B45'$ ; number of stars of the B4 stellar population ; ] iv = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001];,10,1.,10.,1.] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) ;llims =[1e-9 ,0. ,0. ,0. ,0. ,0. ,0. ,1. ,0. , 1. ,0. ] END 'COMPIEGNE_ETAL2010':BEGIN ;ORIGINAL LINES ;parameter description of parameters to be fitted pd = [ $ '(*!dustem_params).gas.G0', $ ;G0 'dustem_plugin_continuum_2', $ ;intensity of NIR continuum 'dustem_plugin_synchrotron_2',$ '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx '(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx '(*!dustem_params).grains(4).mdust_o_mh' $ ] ;initial parameter values for parameters to be fitted iv = [ $ 1.0, $ 0.002, $ ;intensity of NIR continuum 0.01,$ 7.8e-4, $ ;mass fraction of PAH0 7.8e-4, $ ;mass fraction of PAH1 1.65e-4, $ ;mass fraction of amCBEx 1.45e-3, $ ;mass fraction of amCBEx 7.8e-3 $ ;mass fraction of ] ;NEW LINES ;========COMMENTS======= ;amplitudes haven't been used so they're by default set to 1. ;using only the below stellar populations produces amplitudes that exceed 1. I do not understand the physical meaning of this. ; pd = [ $ ; '(*!dustem_params).gas.G0', $ ;G0 ; 'dustem_create_continuum_2', $ ;intensity of NIR continuum ; 'dustem_create_stellar_population_O63',$ ; distance to the O3 stellar population ; ;'dustem_create_stellar_population_O64',$ ; amplitude of the O3 stellar population ; 'dustem_create_stellar_population_O65',$ ; number of stars of the O3 stellar population ; 'dustem_create_stellar_population_B43',$ ; distance to the B4 stellar population ; ;'dustem_create_stellar_population_B44',$ ; amplitude of the B4 stellar population ; 'dustem_create_stellar_population_B45',$ ; number of stars of the B4 stellar population ; '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction ; '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction ; '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx ; '(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx ; '(*!dustem_params).grains(4).mdust_o_mh' $ ; ] ; pd = [ $ ; ;'(*!dustem_params).G0', $ ;G0 ; '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction ; '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction ; '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx ; '(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx ; '(*!dustem_params).grains(4).mdust_o_mh', $ ; ;'dustem_create_continuum_2', $ ;intensity of NIR continuum ; 'dustem_create_stellar_population_O63',$ ; distance to the O6 stellar population ; ;'dustem_create_stellar_population_O64',$ ; amplitude of the O6 stellar population ; ;'dustem_create_stellar_population_O65',$ ; number of stars of the O6 stellar population ; 'dustem_create_stellar_population_B43'$ ; distance to the B4 stellar population ; ;;'dustem_create_stellar_population_B44',$ ; amplitude of the B4 stellar population ; ;'dustem_create_stellar_population_B45'$ ; number of stars of the B4 stellar population ; ] ;initial parameter values for parameters to be fitted ; iv = [ $ ; 1., $ ;G0 ; 7.8e-4, $ ;mass fraction of PAH0 ; 7.8e-4, $ ;mass fraction of PAH1 ; 1.65e-4, $ ;mass fraction of amCBEx ; 1.45e-3, $ ;mass fraction of amCBEx ; 7.8e-3, $ ;mass fraction of big grains? ; 0.002, $ ;intensity of NIR continuum ; 10.,$ ; ;;1.,$ ; 1.,$ ; 10.,$ ; ;;1.,$ ; 1.$ ; ] ; ; iv = [ $ ; 7.8e-4, $ ;mass fraction of PAH0 ; 7.8e-4, $ ;mass fraction of PAH1 ; 1.65e-4, $ ;mass fraction of amCBEx ; 1.45e-3, $ ;mass fraction of amCBEx ; 7.8e-3, $ ;mass fraction of big grains? ; ;0.002, $ ;intensity of NIR continuum ; 10.,$ ; ;;1.,$ ; ;1.,$ ; 10.$ ; ;;1.,$ ; ;1.$ ; ] Npar=n_elements(pd) ;ulimed=[0 ,0 ,0 ,0 ,0 ,0 , 0 ,0 ,0 ,0 ,0 ,0 ,0] ;ulimed=[0 ,0 ,0 ,0 ,0 ,0 , 0 ,0 ,0 ,0] ;llimed=[1 , 1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1] ;llimed=[1 , 1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1] ;llims =[0. ,0. ,0. ,0. , 0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0.] ;llims =[1. ,0. ,10. ,0. , 0. ,0. ,0. ,0. ,0. ,0. ,0.] ;lower limiting the stellar distances to 1pc avoid singularities ;[1e-9 , ;llims =[1e-5 ,0. ,0. ,0. ,0. ,0. ,0. ,1. ,0. , 1. ,0. ] ;llims =[0. ,0. ,0. ,0. ,0. ,0. ,1. ,0. , 1. ,0. ] ;llims =[0. ,0. ,0. ,0. ,0. ,1. , 1. ] ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) ;=== Fixed parameters ;parameter description of parameters to be set to a non-default value ;fpd=[ $ ;'!dustem_isrf_star_add[0].amplitude', $ ;amplitiude of first stellar contribution to ISRF ;'!dustem_isrf_star_add[0].distance', $ ;amplitiude of first stellar contribution to ISRF ;'(*!dustem_params).gas.G0' $ ;multiplicative factor to total ISRF ;'dustem_create_continuum_2' $ ;intensity of NIR continuum ; ] ;initial parameter values for fixed parameters ;fiv=[ $ ;1.0, $ ;amplitiude of first stellar contribution to ISRF ;1.0, $ ;distance of first stellar contribution to ISRF ; 1. $ ;multiplicative factor to total ISRF ;1.e-6 $ ;intensity of NIR continuum ; ] END 'AJ13':BEGIN pd = [ $ '(*!dustem_params).G0', $ ;G0 '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx '(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum iv = [1.0, 7.8e-4, 7.8e-4,1.65e-4,1.45e-3,7.8e-3,0.001] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) END 'G17_MODELA':BEGIN pd = [ $ '(*!dustem_params).G0', $ ;G0 '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction '(*!dustem_params).grains(2).mdust_o_mh', $ ;Gra '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum iv = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) polarization=1 END 'G17_MODELB':BEGIN pd = [ $ '(*!dustem_params).G0', $ ;G0 '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction '(*!dustem_params).grains(2).mdust_o_mh', $ ;Gra '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum iv = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) polarization=1 END 'THEMIS':BEGIN pd = [ $ '(*!dustem_params).G0', $ ;G0 '(*!dustem_params).grains(0).mdust_o_mh',$ '(*!dustem_params).grains(1).mdust_o_mh',$ '(*!dustem_params).grains(2).mdust_o_mh', $ 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum iv = [1.0,5.4e-4,2.33e-3,8.27e-3,0.001] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) polarization=1 END ENDCASE dustem_init,model=use_model,pol=polarization !dustem_nocatch=1 !dustem_verbose=1 !dustem_show_plot=1 ;=== Read sample SED ;=== Composite SED from Compiegne et al 2010, gathered by C. Bot dir=!dustem_wrap_soft_dir+'/Data/SEDs/' file=dir+'Gal_composite_spectrum.xcat' ;file=dir+'Gal_composite_spectrum.xcat.old' spec=read_xcat(file,/silent) ;ind=where(spec.error EQ 0.,count) ;IF count NE 0 THEN spec(ind).error=0.2*spec(ind).spec ind=where(spec.sigmaII EQ 0.,count) IF count NE 0 THEN spec(ind).sigmaII=(0.2*spec(ind).StokesI)^2 ind=where(spec.instru EQ 'FIRAS',count) ;IF count NE 0 THEN spec(ind).error=0.2*spec(ind).spec IF count NE 0 THEN spec(ind).sigmaII=(0.2*spec(ind).StokesI)^2 ;=== SET THE OBSERVATION STRUCTURE st=dustem_set_data(sed=spec) ;== SET THE FITTED PARAMETERS dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims ;dustem_init_fixed_params,fpd,fiv dustem_init_plugins, pd ;initializing the scopes of the plugins to be read later ;=== RUN fit tol=1.e-16 ;use_Nitermax=50 ;maximum number of iteration. This is the criterium which will stop the fit procedure , NOTA BENE: 3 iterations aren't enough to fit the entirety of the newly added stellar population parameters. use_Nitermax=5 ;maximum number of iteration. This is the criterium which will stop the fit procedure , NOTA BENE: 3 iterations aren't enough to fit the entirety of the newly added stellar population parameters. IF keyword_set(itermax) THEN use_Nitermax=itermax t1=systime(0,/sec) yr=[1e-4,100] xr=[1,6e4] ;legend_xpos=0.6 ;legend_ypos=0.8 tit='Spectral Energy Distribution (Running)' ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') xtit=textoidl('\lambda (\mum)') 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) t2=systime(0,/sec) ;=== SAVE FIT RESULTS file_out='/tmp/DUSTEM_fit_example.sav' dustem_save_system_variables,file_out message,'Saved '+file_out,/continue ;stop ;====================================== ;====You can exit IDL here and re-enter ;====================================== file='/tmp/DUSTEM_fit_example.sav' dustem_restore_system_variables,file ;=== Plot best fit tit='Spectral Energy Distribution (Saved)' ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') xtit=textoidl('\lambda (\mum)') errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values) chi2=(*!dustem_fit).chi2 rchi2=(*!dustem_fit).rchi2 window,2 ;=== RESTORE FIT RESULTS res=*(*!dustem_fit).current_param_values chi2=(*!dustem_fit).chi2 rchi2=(*!dustem_fit).rchi2 errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values) ;=== Plot best fit loadct,13 IF keyword_set(postcript) THEN BEGIN set_plot,'PS' ; ps_file=getenv('DUSTEM_SOFT_DIR')+'/Docs/Figures/'+'Last_dustem_fit.ps' ps_file=!dustem_wrap_soft_dir+'Docs/Figures/'+'Last_dustem_fit.ps' device,filename=ps_file,/color ENDIF ;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 dustem_sed_plot,*(*!dustem_fit).current_param_values,ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog,legend_xpos=legend_xpos,legend_ypos=legend_ypos IF keyword_set(postcript) THEN BEGIN device,/close set_plot,'X' message,'wrote '+ps_file,/info ENDIF IF keyword_set(png) THEN BEGIN ;file_png=!dustem_wrap_soft_dir+'/Docs/Figures/'+'Last_dustem_fit.png' file_png=dir_png+'Last_dustem_fit_sed_readme.png' write_png,file_png,tvrd(/true) message,'Wrote '+file_png,/info ENDIF message,'dustem_mpfit_sed executed in '+strtrim(t2-t1,2)+' sec',/info the_end: END