PRO dustem_fit_sed_spinning_readme,postcript=postcript,mode=mode,help=help ;This is the version with spinning dust activated. Not very useful for the proposed SED ;This Readme describes how to fit SEDs with dustem with spinning dust turned on. ;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_spinning_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_spinning_readme,postcript=postcript,mode=mode,help=help ; INPUTS: ; None ; OPTIONAL INPUT PARAMETERS: ; None ; OUTPUTS: ; None ; OPTIONAL OUTPUT PARAMETERS: ; None ; ACCEPTED KEY-WORDS: ; mode = 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,mode='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_spinning_readme' goto,the_end ENDIF message,'Obsolete. use dustem_fit_spin_test_readme.pro' IF keyword_set(mode) THEN BEGIN use_mode=strupcase(mode) ENDIF ELSE BEGIN use_mode='COMPIEGNE_ETAL2010' ;Default is last dustem model ENDELSE dustem_init,mode=use_mode !dustem_verbose=1 !dustem_show_plot=1 ;===== Read the SED dir=!dustem_wrap_soft_dir+'/Data/SEDs/' file=dir+'Gal_composite_spectrum.xcat' 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.instru EQ 'FIRAS',count) IF count NE 0 THEN spec(ind).error=0.2*spec(ind).spec ;=== Set which parameters you want to fit CASE !dustem_which OF 'RELEASE':BEGIN ;=== This is to use Desert model ;=== must have done dustem_init,/DBP90 CASE use_mode 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_create_continuum_2'] ;Intensity of NIR continuum iv = [1.0, 4.3e-4, 4.7e-4,6.4e-3,0.001] ; iv = [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) ;=== Fixed parameters ; fpd=['(*!dustem_params).G0'] ;Fixed parameter description ; fiv=[1.5] ; dustem_init_fixed_params,fpd,fiv 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_create_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] ; iv = [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) ;=== Fixed parameters ; fpd=['(*!dustem_params).G0'] ;Fixed parameter description ; fiv=[1.0] ; dustem_init_fixed_params,fpd,fiv 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_create_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] ; iv = [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) ;=== Fixed parameters ; fpd=['(*!dustem_params).G0'] ;Fixed parameter description ; fiv=[1.0] ; dustem_init_fixed_params,fpd,fiv END 'COMPIEGNE_ETAL2010':BEGIN ; PUTTING G0 last DOES NOT WORK BECAUSE OF ORDER ISSUE 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_create_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] ; iv = [7.8e-4, 7.7e-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) ;=== Fixed parameters ; fpd=['(*!dustem_params).G0'] ;Fixed parameter description ; fiv=[1.0] ; dustem_init_fixed_params,fpd,fiv END ENDCASE END ELSE: BEGIN pd = ['(*!dustem_params).grains(0).propmass',$ ;PAH mass fraction '(*!dustem_params).grains(1).propmass',$ ;VSG mass fraction '(*!dustem_params).grains(2).propmass',$ ;BG mass fraction 'dustem_create_continuum_2', $ ;Intensity of NIR continuum 'dustem_create_isrf_1'] ;ISRF scaling factor ;=== set initial value of parameters you want to fit iv = [7e-5, 8e-4,6.6e-3,0.001,4.] ;in this example, we enforce parameters to be >0 and Xisrf >3 ;Caution, initial values must be within limits ... Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) END ENDCASE ;=== 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 ;stop ;help,(*!dustem_params).grains[0],/str (*!dustem_params).grains[0].TYPE_KEYWORDS='logn-spin' (*!dustem_params).grains[1].TYPE_KEYWORDS='logn-spin' ;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-chrg' ;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-chrg' ;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-mix-spin' ;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-mix-spin' ;dustem_sort_params ;hprint,*(*!dustem_fit).param_descs ;hprint,*(*!dustem_fit).param_init_values ;hprint,*(*!dustem_fit).param_func ;=== RUN fit tol=1.e-14 Nitermax=2 ;maximum number of iteration. This is the criterium which will stop the fit procedure loadct,13 !y.range=[1e-4,10] ;This is to ajust plot range from outside the routine t1=systime(0,/sec) res=dustem_mpfit_data(tol=tol,Nitermax=Nitermax,gtol=gtol,/xlog,/ylog,xrange=[1,20000]) t2=systime(0,/sec) ;=== SAVE FIT RESULTS ;file_out=getenv('DUSTEM_RES')+'DUSTEM_fit_example.sav' file_out=!dustem_res+'DUSTEM_fit_example.sav' dustem_save_sed_fit,file_out ;=== Plot best fit yr=[1e-4,10] xr=[1,20000] tit='DUSTEM Example' 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 loadct,13 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 ;stop ;====================================== ;====You can exit IDL here and re-enter ;====================================== window,1 ;=== RESTORE FIT RESULTS ;=== initialise dustem dustem_init,mode=use_mode ;CAUTION: MUST use the same mode as earlier !! file=!dustem_res+'DUSTEM_fit_example.sav' dustem_restore_sed_fit,file ;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-chrg' ;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-chrg' ;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-mix-spin' ;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-mix-spin' (*!dustem_params).grains[0].TYPE_KEYWORDS='logn-spin' (*!dustem_params).grains[1].TYPE_KEYWORDS='logn-spin' ;=== recover best fit values ;ind=where(*!dustem_fit_chi2 EQ min(*!dustem_fit_chi2)) ;res=(*!dustem_fit_params); res=*(*!dustem_fit).current_param_values chi2=(*!dustem_fit).chi2 rchi2=(*!dustem_fit).rchi2 errors=*(*!dustem_fit).current_param_errors ;=== Plot best fit yr=[1e-4,10] xr=[1,20000] tit='DUSTEM Example (Saved)' ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') xtit=textoidl('\lambda (\mum)') loadct,13 IF keyword_set(postcript) THEN BEGIN set_plot,'PS' ps_file=!dustem_wrap_soft_dir+'/Docs/Figures/'+'Last_dustem_fit.ps' device,filename=ps_file,/color ENDIF 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 IF keyword_set(postcript) THEN BEGIN device,/close set_plot,'X' message,'wrote '+ps_file,/info ENDIF message,'dustem_mpfit_sed executed in '+strtrim(t2-t1,2)+' sec',/info the_end: END