added

Jean-Philippe Bernard
1 parent 6c9d0328
Showing 2 changed files with 373 additions and 0 deletions Show diff stats
src/idl/dustem_fit_sed_spinning_readme.pro
src/idl/dustem_plot_extinction.pro
@@ -0,0 +1,287 @@
+PRO dustem_fit_sed_spinning_readme,postcript=postcript,mode=mode,help=help
+
+;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
+
+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
+  'WEB3p8':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='mix-logn-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)
+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,2000]
+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
+
+;stop
+
+;======================================
+;====You can exit IDL here and re-enter
+;======================================
+window,1
+
+;=== RESTORE FIT RESULTS
+;=== initialise dustem
+dustem_init,mode=mode
+;file=getenv('DUSTEM_RES')+'DUSTEM_fit_example.sav'
+file=!dustem_res+'DUSTEM_fit_example.sav'
+dustem_restore_sed_fit,file
+
+;=== 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
+;model_sed=*!dustem_fit_sed
+;errors=*!dustem_fit_params_error;
+errors=*(*!dustem_fit).current_param_errors
+
+;=== Plot best fit
+yr=[1e-4,10]
+xr=[1,2000]
+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=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
+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
@@ -0,0 +1,86 @@
+PRO dustem_plot_extinction,st,st_model,_Extra=extra,help=help
+
+;+
+; NAME:
+;     dustem_plot_extinction
+; PURPOSE:
+;    Plots a Dustem extinction
+; CATEGORY:
+;    Dustem
+; CALLING SEQUENCE:
+;    dustem_plot_extinction,st,st_model[,/help][,_extra=]
+; INPUTS:
+;    st        = Dustem model output structure
+;    st_model  = Dustem model parameters structure
+; OPTIONAL INPUT PARAMETERS:
+;    _extra    = extra parameters for the plot routine
+; OUTPUTS:
+;    None
+; OPTIONAL OUTPUT PARAMETERS:
+;    None
+; ACCEPTED KEY-WORDS:
+;    help      = If set, print this help
+; COMMON BLOCKS:
+;    None
+; SIDE EFFECTS:
+;    model is plotted
+; RESTRICTIONS:
+;    The dustem idl wrapper must be installed
+; PROCEDURE:
+;    None
+; EXAMPLES
+;    dustem_init,/wrap
+;    dir_in=getenv('DUSTEM_SOFT_DIR')+'/src/dustem3.8_web/'
+;    st_model=dustem_read_all(dir_in)
+;    dir=getenv('DUSTEM_RES')
+;    st=dustem_read_all_res(dir_in)
+;    loadct,13
+;    dustem_plot_extinction,st,st_model,yr=[1.e-6,1.e0],/ysty,xr=[1.,400],/xsty
+; MODIFICATION HISTORY:
+;    Written by J.-Ph. Bernard
+;    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_plot_extinction'
+  goto,the_end
+ENDIF
+
+;stop
+;Ngrains=n_tags(st.dustem)-2
+Ngrains=st_model.Ngrains
+
+colors=dustem_grains_colors(Ngrains)
+;colmax=210 & colmin=30
+;colors=findgen(Ngrains)/(1.*Ngrains-1)*(colmax-colmin)+colmin
+
+xtit=textoidl('\lambda (\mum)')
+;ytit=textoidl('extinction (cm^2/g)') & fact=1.
+ytit=textoidl('\sigma_{ext} (1e-21 cm^2/H)') & mH=1.66e-24 & fact=mH*1.e21
+tit='Dustem extinction'
+
+;stop
+
+i=0L
+;plot_oo,st.ext.wav,st.ext.abs_grain(i)*fact,/nodata,_Extra=extra,xtit=xtit,ytit=ytit,tit=tit
+ffact=fact*st_model.grains(i).MDUST_O_MH
+plot_oo,st.ext.wav,(st.ext.abs_grain(i)+st.ext.sca_grain(i))*ffact,/nodata,_Extra=extra,xtit=xtit,ytit=ytit,tit=tit
+abs_tot=st.ext.abs_grain(i)*0.
+sca_tot=st.ext.sca_grain(i)*0.
+FOR i=0L,Ngrains-1 DO BEGIN
+  ffact=fact*st_model.grains(i).MDUST_O_MH
+  oplot,st.ext.wav,(st.ext.abs_grain(i)+st.ext.sca_grain(i))*ffact,color=colors(i),psym=psym
+;  oplot,st.ext.wav,st.ext.abs_grain(i)*fact,color=colors(i-1),psym=psym
+;  oplot,st.ext.wav,st.ext.sca_grain(i)*fact,color=colors(i-1),psym=psym,linestyle=2
+  abs_tot=abs_tot+st.ext.abs_grain(i)*ffact
+  sca_tot=sca_tot+st.ext.sca_grain(i)*ffact
+ENDFOR
+oplot,st.ext.wav,(abs_tot+sca_tot),color=255,psym=psym
+;oplot,st.dustem.wav,(abs_tot+sca_tot),color=255,psym=psym
+;oplot,st.dustem.wav,abs_tot*fact,color=255,psym=psym
+;oplot,st.dustem.wav,sca_tot*fact,color=255,psym=psym,linestyle=2
+
+the_end:
+
+END