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 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 ] Npar=n_elements(pd) 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,kwords=['plaw-ed','logn','logn'];,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(st_sed_fit=spec);sed=spec) ;stop ;== 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=[1.00e-4,1.00E2] xr=[1.00E0,6.00e4] ;legend_xpos=0.6 ;legend_ypos=0.8 tit='Spectral Energy Distribution' 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 tit='Spectral Energy Distribution (Final run)' dustemwrap_plot,res,stp,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit message,'dustem_mpfit_sed executed in '+strtrim(t2-t1,2)+' sec',/info the_end: END