First commit

Jean-Philippe Bernard
1 parent 666cba76
Showing 2 changed files with 373 additions and 0 deletions Show diff stats
src/idl/dustem_define_isrf_stellar_modifyer_variable.pro
src/workshop/dustem_fit_sed_n26_twice.pro
@@ -0,0 +1,28 @@
+FUNCTION dustem_define_isrf_stellar_modifyer_variable,stellar_values=stellar_values
+
+one_st={radius:0.,temperature:0.,distance:0.,amplitude:0.}
+
+IF not keyword_set(stellar_values) THEN BEGIN
+  ;The defaults is one sun at 1au, but with zero amplitude
+  Nstar=1
+  st=replicate(one_st,Nstar)
+  one_pc=3.08e16
+  one_au=1.49e11
+  st[0].radius=1.          ;Rsun
+  st[0].temperature=5770.    ;K
+  st[0].distance=one_au/one_pc         ;pc
+  st[0].amplitude=0.         ;--
+  status=1
+ENDIF ELSE BEGIN
+  Nstars=n_elements(stellar_values)
+  ;st=replicate(one_st,Nstar)
+  st=stellar_values
+  status=2
+ENDELSE
+
+defsysv,'!dustem_isrf_star_add',st
+message,'defined !dustem_isrf_star_add',/continue
+
+RETURN,status
+
+END
 \ No newline at end of file
@@ -0,0 +1,345 @@
+PRO dustem_fit_sed_n26_twice,mode=mode,help=help,png=png,postscript=postscript,nitermax=nitermax
+
+;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_n26
+; 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_n26,postscript=postscript,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
+;    postscript = 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_n26,mode='COMPIEGNE_ETAL2010',nitermax=10
+; 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_n26_twice'
+  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
+
+IF keyword_set(png) THEN BEGIN
+  dir_png='/tmp/Dustem/Figures/'
+  force_mkdir,dir_png
+ENDIF
+IF keyword_set(postscript) THEN BEGIN
+  dir_ps='/tmp/Dustem/Figures/'
+  force_mkdir,dir_ps
+ENDIF
+
+;dbp=1
+;=== initialise dustem
+;dustem_init,/wrap
+;dustem_init
+dustem_init,mode=use_mode
+!dustem_verbose=1
+!dustem_show_plot=1
+
+;=== Read sample SED
+;=== N26 SED from Longji
+dir=!dustem_wrap_soft_dir+'/Data/SEDs/workshop/'
+file=dir+'GN26_SED.xcat'
+!dustem_redshift=1.22
+
+spec=read_xcat(file,/silent)
+;stop
+
+ind=where(spec.error EQ 0.,count)
+IF count NE 0 THEN spec(ind).error=0.2*spec(ind).spec
+;=== Draine Model SED for U=1. by C. Bot
+;;dir=getenv('DUSTEM_SOFT_DIR')+'/src/dustem3.8_web/SEDs/'
+;dir=!dustem_wrap_soft_dir+'/Data/SEDs/'
+;file=dir+'SED_DraineModel_U1.0.xcat'
+;spec=read_xcat(file,/silent)
+;=== Fake SED made from Dustem outputs
+;spec=dustem_compute_sed(p_dim,st=st,cont=cont,_extra=extra)
+
+;=== Set which parameters you want to fit
+;CASE getenv('DUSTEM_WHICH') OF
+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-3, 5.4e-3,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
+        ; pd = [ $
+        ;      '(*!dustem_params).gas.G0', $      ;G0
+        ;      '(*!dustem_params).grains(0).mdust_o_mh',$     ;PAH1 mass fraction
+        ;      '(*!dustem_params).grains(1).mdust_o_mh',$  ;PAH0 mass fraction
+        ;      '(*!dustem_params).grains(2).mdust_o_mh', $    ;amCBEx
+        ;      '(*!dustem_params).grains(2).alpha_o_a0',$     ;amCBEx slope of size distribution (logn)
+        ;      '(*!dustem_params).grains(3).mdust_o_mh', $    ;amCBEx
+        ;      '(*!dustem_params).grains(4).mdust_o_mh']      ;aSil
+        ; iv =   [45., 7.8e-5, 7.8e-5,1.65e-4,2.e-7,1.45e-3,7.8e-3]
+        ; Npar=n_elements(pd)
+        ; ulimed=replicate(0,Npar)
+        ; llimed=replicate(1,Npar)
+        ; llims=[0.,0.,0.,0.,0.,0.,0.]
+
+        pd = [ $
+             '(*!dustem_params).gas.G0', $      ;G0
+             '(*!dustem_params).grains(0).mdust_o_mh',$     ;PAH1 mass fraction
+             '(*!dustem_params).grains(1).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(3).alpha_o_a0',$     ;amCBEx slope of size distribution (logn)
+             '(*!dustem_params).grains(4).mdust_o_mh']      ;aSil
+        iv =   [45., 7.8e-5, 7.8e-5,1.65e-4,1.45e-3,-2.,7.8e-3]
+        Npar=n_elements(pd)
+        ulimed=replicate(0,Npar)
+        llimed=replicate(0,Npar)
+        llims=[0.,0.,0.,0.,0.,0.,0.]
+
+        ;=== 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
+;=== RUN fit
+tol=1.e-14
+use_Nitermax=2               ;maximum number of iteration. This is the criterium which will stop the fit procedure
+IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax
+
+loadct,13
+;!y.range=[1e-4,10]              ;This is to ajust plot range from outside the routine
+t1=systime(0,/sec)
+yr=[1e-1,1.e3]
+xr=[10,10000]
+legend_xpos=0.6
+legend_ypos=0.8
+tit='DUSTEMWrap Intensity SED GN26 (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=getenv('DUSTEM_RES')+'DUSTEM_fit_example.sav'
+file_out='/tmp/DUSTEM_fit_example1.sav'
+;dustem_save_sed_fit,file_out
+dustem_save_system_variables,file_out
+message,'Saved '+file_out,/continue
+
+;stop
+
+;==== Doing a second fit
+dir=!dustem_wrap_soft_dir+'/Data/SEDs/workshop/'
+file=dir+'GN26_SPEC.xcat'
+
+        pd = [ $
+             '(*!dustem_params).grains(0).mdust_o_mh',$     ;PAH1 mass fraction
+             '(*!dustem_params).grains(1).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(3).alpha_o_a0',$     ;amCBEx slope of size distribution (logn)
+             '(*!dustem_params).grains(4).mdust_o_mh']      ;aSil
+        iv =   [7.8e-5, 7.8e-5,1.65e-4,1.45e-3,-2.,7.8e-3]
+        Npar=n_elements(pd)
+        ulimed=replicate(0,Npar)
+        llimed=replicate(0,Npar)
+        llims=[0.,0.,0.,0.,0.,0.,0.]
+        fpd=['(*!dustem_params).gas.G0'] ;Fixed parameter description
+        fiv=[30.]
+        dustem_init_fixed_params,fpd,fiv
+
+dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
+
+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
+
+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
+;=== RUN fit
+tol=1.e-14
+use_Nitermax=2               ;maximum number of iteration. This is the criterium which will stop the fit procedure
+IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax
+
+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)
+
+file_out='/tmp/DUSTEM_fit_example2.sav'
+;dustem_save_sed_fit,file_out
+dustem_save_system_variables,file_out
+message,'Saved '+file_out,/continue
+
+;======================================
+;====You can exit IDL here and re-enter
+;======================================
+
+file='/tmp/DUSTEM_fit_example1.sav'
+dustem_restore_system_variables,file
+
+p_min1=*(*!dustem_fit).current_param_values
+dustem_sed1=dustem_compute_sed(p_min1,st=st,cont=cont,freefree=freefree,synchrotron=synchrotron,out_st=out_st1)
+st1=out_st1
+z1=!dustem_redshift
+help,st1.sed,/str
+
+file='/tmp/DUSTEM_fit_example2.sav'
+dustem_restore_system_variables,file
+
+p_min2=*(*!dustem_fit).current_param_values
+dustem_sed2=dustem_compute_sed(p_min2,st=st,cont=cont,freefree=freefree,synchrotron=synchrotron,out_st=out_st2)
+st2=out_st2
+z2=!dustem_redshift
+help,st2.sed,/str
+
+;This is a magic factor to bring SED to the right units
+fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/st1.sed.wav)*1.e20/1.e7
+
+yr=[1e-1,1.e3]
+xr=[10,10000]
+tit='DUSTEMWrap Intensity SED GN26 (Comparison)'
+ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
+xtit=textoidl('\lambda (\mum)')
+
+cgplot,st1.sed.wav*(1.+z1),st1.sed.em_tot*fact,/xlog,/ylog,yr=yr,/ysty,xr=xr,/xsty,xtit=xtit,ytit=ytit,tit=tit
+cgoplot,st1.sed.wav*(1.+z1),st1.sed.em_grain_4*fact,color='red'
+
+cgoplot,st2.sed.wav*(1.+z2),st2.sed.em_tot*fact,linestyle=2
+cgoplot,st2.sed.wav*(1.+z2),st2.sed.em_grain_4*fact,color='pink',linestyle=2
+
+stop
+
+the_end:
+
+END