improved comments and removed unimplemented keywords

Annie Hughes
1 parent 79bbc7a3
Showing 1 changed file with 103 additions and 110 deletions Show diff stats
src/idl/dustem_fit_intensity_example.pro
 PRO dustem_fit_intensity_example,model=model $
                                 ,sed_file=sed_file $
                                 ,Nitermax=Nitermax $
-                                ,postscript=postscript $
                                 ,fits_save=fits_save $
                                 ,help=help $
                                 ,wait=wait $
@@ -31,7 +30,7 @@ PRO dustem_fit_intensity_example,model=model $
 ;    DustEMWrap, Distributed, High-Level, User Example
 ;
 ; CALLING SEQUENCE:
-;    dustem_fit_intensity_example[,model=][sed_file=][,postscript=][,Nitermax=][,fits_save=][,/help,/wait,/verbose]
+;    dustem_fit_intensity_example[,model=][sed_file=][,Nitermax=][,fits_save=][,/help,/wait,/verbose]
 ;
 ; INPUTS:
 ;    None
@@ -45,16 +44,15 @@ PRO dustem_fit_intensity_example,model=model $
 ; OPTIONAL OUTPUT PARAMETERS:
 ;    Plots, results structure in binary FITS table format
 ;
-; ACCEPTED KEY-WORDS:
+; ACCEPTED KEY-WORDS
 ;    model = specifies the interstellar dust mixture used by
 ;            DustEM. See userguide or dustem_test_model_exists.pro
 ;            for more details about available models in current release.
 ;    sed_file = string naming the path to text file in .xcat format that
 ;          describes the observational SED. If not set, the file
 ;          'Data/EXAMPLE_OBSDATA/example_SED_1.xcat' is used.  
-;    postscript = if set, final plot is saved as postscript file
 ;    Nitermax = maximum number of fit iterations. Default is 5.
-;    fits_save = if set, save the fit results in a binary
+;    fits_save = if set, save the DustEMWrap fitting results in a binary
 ;               FITS file. 
 ;    help      = if set, print this help
 ;    wait      = if set, wait this many seconds between each step of
@@ -90,23 +88,27 @@ IF keyword_set(help) THEN BEGIN
   goto,the_end
 END
+;; ;===================================
+;; ;=== SET THE DUST MODEL HERE
+;; ;===================================
 IF keyword_set(model) THEN BEGIN
   use_model=strupcase(model)
 ENDIF ELSE BEGIN
   use_model='DBP90'    ;Example with default keywords uses the DBP90 model
 ENDELSE
-
-exists=dustem_test_model_exists(use_model)
+use_polarization=0   ; here we initialize DustEMWrap in no polarization mode since we are only fitting Stokes I
+exists=dustem_test_model_exists(use_model,/silent)
 if exists ne 1 then $
    message,'Unknown dust model'
-use_polarization=0   ; initialize Dustemwrap in no polarization mode 
-use_window=2          ; default graphics window number to use for plotting the results
+;; ;===================================
+;; ;=== PARSE OTHER INPUTS AND SET SOME BASIC RUN PARAMETERS
+;; ;===================================
 use_verbose=0
 if keyword_set(verbose) then use_verbose=1
 use_Nitermax=5        ; maximum number of iterations for the fit
 IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax
-
+use_window=0          ; default graphics window number to use for plotting the results
 dustem_define_la_common
 ;=== Set the (model-dependent) parameters that you want to fit (pd),
@@ -120,22 +122,27 @@ dustem_define_la_common
 ;=== Examples are provided for some of the dust models. 
 ;=== To try them, uncomment the model that you want to try and re-run
-;=== AN EXAMPLE FOR DBP90
+;; ;===================================
+;; ;=== AN EXAMPLE FOR DBP90
+;; ;===================================
 ;=== Here we fit the dust abundances of the DBP90 model, the
-;=== intensity of the dust-heating radiation field, as well as a plug-in
-;=== for synchrotron emission
+;=== intensity of the dust-heating radiation field. We also invoke a plug-in
+;=== for that fits parameters of synchrotron emission
 ;=== The free parameters are all lower-bounded at zero.
-;=== use_model='DBP90' ; you should specify this above!
+;=== use_model='DBP90' ; you should specify this above (line 93)
+
+;=== pd is the structure of free parameters
 pd = [ $
      '(*!dustem_params).G0', $                           ;G0
      '(*!dustem_params).grains(0).mdust_o_mh',$          ;PAH0 mass fraction
      '(*!dustem_params).grains(1).mdust_o_mh',$          ;VSG mass fraction
-     '(*!dustem_params).grains(2).mdust_o_mh'];,$          ;BG mass fraction
-;     'dustem_plugin_synchrotron_1', $                    ;Spectral index of CREs
-;     'dustem_plugin_synchrotron_2']                      ;Synchrotron amplitude at 10 mm
+     '(*!dustem_params).grains(2).mdust_o_mh'];,$        ;BG mass fraction
+;     'dustem_plugin_synchrotron_1', $                   ;Spectral index of CREs
+;     'dustem_plugin_synchrotron_2']                     ;Synchrotron amplitude at 10 mm
-; initial values vector for run without synchrotron plugin
-; [G0,PAH0,VSG,BG]
+;=== iv is the vector of initial values for the free parameters
+;=== following lines are for a run without fitting the synchrotron
+;=== i.e. [ initial values for G0,PAH0,VSG,BG]
 iv =   [1.6, 2.2e-4, 5.7e-4, 3.4e-3] 
 ; initial values vector for run including the synchrotron plugin
@@ -143,13 +150,15 @@ iv =   [1.6, 2.2e-4, 5.7e-4, 3.4e-3]
 ;iv =   [1.6, 2.2e-4, 5.7e-4, 3.4e-3, 2.7,0.01]
 ; initial values vector for run including the synchrotron plugin and
-; fixing the G0 and BB continuum parameters
+; fixing the G0 and BB continuum parameters (see below)
+; [PAH0,VSG,BG,alpha_CR,Amp_syn]
 ;iv =   [2.2e-4, 5.7e-4, 3.4e-3, 2.7,0.01]
 Npar=n_elements(pd)
-ulimed=replicate(0,Npar)
-llimed=replicate(1,Npar) 
-llims=replicate(1.e-15,Npar)
+ulimed=replicate(0,Npar) ; flag ON=1, OFF=0
+llimed=replicate(1,Npar) ; flag ON=1, OFF=0
+llims=replicate(1.e-15,Npar) ; lower limit value for each free parameter
+; ulims=replicate(1.e15,Npar) ; upper limit value for each free parameter
 fpd=[] & fiv=[]
 ; example using fixed parameters (ISRF and a NIR continuum)
@@ -164,40 +173,51 @@ fpd=[] &amp; fiv=[]
 ;fiv=[0.5,750,1.e-2]
 ;; ;==================================
-;; ;=== EXAMPLES FOR OTHER DUST MODELS
+;; ;=== EXAMPLES FOR OTHER PHYSICAL DUST MODELS
 ;; ;=== START BELOW HERE
 ;; ;==================================
-
+;; ;===================================
 ;; ;=== AN EXAMPLE FOR DL07
+;; ;===================================
 ;; ;=== Here we fit the dust abundances of the model and the
 ;; ;=== intensity of the ISRF (via gas.G0 since the model
 ;; ;=== includes spinning dust)
 ;; ;=== The free parameters are all lower-bounded at zero.
-;; use_model='DL07' ; you should specify this above, or in the command line
-;; pd = [  '(*!dustem_params).gas.G0'], $             
+;; ;=== use_model='DL07' ; you should specify this above (line 93), or in the command line when you run the example
+;; 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
-;;    iv =   [1.5 ,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3]
-;;    Npar=n_elements(pd)
-;;    ulimed=replicate(0,Npar)
-;;    llimed=replicate(1,Npar)
-;;    llims=replicate(1.e-15,Npar)
-;; ;   fpd=['(*!dustem_params).gas.G0']
-;; ;   fiv=[2.]
-   
-
-;; ;; ;=== AN EXAMPLE FOR MC10
+;; iv =   [1.5 ,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3]
+;; 
+;; ;=== Uncomment the following lines if you want to fix the ISRF instead of leaving it as a free parameter
+;; ;pd = ['(*!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
+;; ;iv =   [5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3]
+;; ;fpd=['(*!dustem_params).gas.G0']
+;; ;fiv=[2.]
+;; 
+;; Npar=n_elements(pd)
+;; ulimed=replicate(0,Npar)
+;; llimed=replicate(1,Npar)
+;; llims=replicate(1.e-15,Npar)
+
+;; ;===================================
+;; ;=== AN EXAMPLE FOR MC10
+;; ;===================================
 ;; ;; ;=== Here we fit the dust abundances of the MC10 model, the
 ;; ;; ;=== intensity of the dust-heating radiation field as well as a plug-in:
 ;; ;; ;===  (i) continuum due to a blackbody 
 ;; ;; ;=== The intensity of the dust-heating radiation field is fixed to
 ;; ;; ;=== 1.5*G0 and the tmperature of the blackbody is fixed to 1200K
 ;; ;; ;=== The free parameters in the fit are lower-bounded at zero.
-;; use_model='MC10' ; you should specify this above, or in the command line
+;; ;; ;=== use_model='MC10' ; you should specify this above, or in the command line
 ;; pd = [ $
 ;;      '(*!dustem_params).grains(0).mdust_o_mh'$      ;PAH0 mass fraction
 ;;      ,'(*!dustem_params).grains(1).mdust_o_mh' $    ;PAH1 mass fraction
@@ -217,16 +237,17 @@ fpd=[] &amp; fiv=[]
 ;;     ]
 ;; fiv=[1.5, 1200.]
-
-
+;; ;===================================
 ;; ;=== AN EXAMPLE FOR J13
+;; ;===================================
 ;; ;=== Here we fit the dust abundances of the J13 model, the
-;; ;=== intensity of the dust-heating radiation field as well as two plug-ins:
+;; ;=== intensity of the dust-heating radiation field as well as the free 
+;; ;=== parameters of two plug-ins:
 ;; ;=== (i) free-free emission
 ;; ;=== (ii)continuum due to a blackbody 
 ;; ;=== The temperature of the blackbody is fixed to 1000K
 ;; ;=== The free parameters in the fit are all lower-bounded at zero.
-;; use_model='J13' ; you should specify this above, or in the command line
+;; ;=== use_model='J13' ; you should specify this above, or in the command line
 ;; pd = [ $
 ;;      '(*!dustem_params).G0' $                                   ;G0
 ;;      ,'(*!dustem_params).grains(0).mdust_o_mh'$                  ;CM20 -- power law size distribution
@@ -235,36 +256,27 @@ fpd=[] &amp; fiv=[]
 ;;      ,'(*!dustem_params).grains(3).mdust_o_mh' $                ;aOlM5
 ;;      ,'dustem_plugin_freefree_1' $                              ;ionized gas temperature
 ;;      ,'dustem_plugin_freefree_2' $                              ;free-free amplitude
-;;      ,'dustem_plugin_continuum_2']                              ;Intensity at peak of the BB continuum
+;;      ,'dustem_plugin_continuum_2']                              ;intensity at peak of the BB continuum
 ;; iv =   [1.2,1.7e-3, 6.3e-4, 2.55e-3, 2.55e-3, 7500., 0.4, 0.001]         
+;; fpd=[ 'dustem_plugin_continuum_1'] ; Temperature of the BB
+;; fiv=[1000.]                       
 ;; Npar=n_elements(pd)
 ;; ulimed=replicate(0,Npar)
 ;; llimed=replicate(1,Npar)
 ;; llims=replicate(1.e-15,Npar)
-;; fpd=[ 'dustem_plugin_continuum_1'] ; Temperature of the BB
-;; fiv=[1000.]                       
-
-
-;;Nfix=n_elements(fpd)
-;;if n_elements(fiv) ne Nfix then begin
-;;   message,'Number of fixed parameters (fpd) does not equal number of initial values of fixed parameters (fiv)',/info
-;;   stop
-;;end
-
-if keyword_set(wait) then begin
-   message,'Finished setting dust model and plug-in parameters: '+use_model,/info
-   wait,wait
-end
-;== INITIALISE DUSTEM
-dustem_init,model=use_model,polarization=use_polarization,/show_plot
-;stop
+;;===================================
+;;=== INITIALISE DUSTEM
+;;===================================
+dustem_init,model=use_model,polarization=use_polarization
 !dustem_nocatch=1
 !dustem_verbose=use_verbose
 IF keyword_set(noobj) THEN !dustem_noobj=1
 !EXCEPT=2 ; for debugging
-;=== READ EXAMPLE SED DATA
+;;===================================
+;;=== READ EXAMPLE SED DATA
+;;===================================
 dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'
 file=dir+'example_SED_1.xcat'
 IF keyword_set(sed_file) THEN file=sed_file
@@ -275,17 +287,22 @@ if keyword_set(wait) then begin
    wait,wait
 end
+;;===================================
 ;;=== ADJUST THE UNCERTAINTIES FROM WITHIN THE CODE (FOR ILLUSTRATION)
+;;===================================
 ind=where(sed.sigmaII LT (0.2*sed.StokesI)^2,count)
 IF count NE 0 THEN sed[ind].sigmaII=(0.2*sed[ind].StokesI)^2
-;== SET THE OBSERVATIONAL STRUCTURE
-;== sed is passed twice in the call -- the first occurrence is the SED that you
-;== wish to fit, the second occurrence is the SED that you wish to visualise. 
+;;===================================
+;=== SET THE OBSERVATIONAL STRUCTURE
+;;===================================
+;== sed is passed twice in the call -- the first occurrence (m_sed) is the SED that you
+;== wish to fit, the second occurrence (m_show) is the SED that you wish to visualise. 
 dustem_set_data,m_fit=sed,m_show=sed
-;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS THAT WILL BE
-;== ADJUSTED DURING THE FIT
+;;===================================
+;=== INITIALISE DUSTEM WITH INITIAL VALUES AND LIMITS OF ALL PARAMETERS 
+;;===================================
 dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims
@@ -294,18 +311,24 @@ if keyword_set(wait) then begin
    wait,wait
 end
-;=== INFORMATION TO RUN THE FIT
+;;===================================
+;;=== INFORMATION TO RUN THE FIT
+;;===================================
 tol=1.e-16     ;fit tolerence
-;=== INFORMATION TO MAKE THE PLOT
+;;===================================
+;;=== INFORMATION TO MAKE THE PLOT
+;;===================================
 yr=[1.00e-4,1.00E2] ; y-axis limits
 xr=[1.00E0,6.00e4] ; x-axis limits
 tit='FIT INTENSITY EXAMPLE' ; plot title
 ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') ; y-axis title
 xtit=textoidl('\lambda (\mum)') ; x-axis title
-;stop
-;===  RUN THE FIT
+
+;;===================================
+;;===  RUN THE FIT
+;;===================================
 t1=systime(0,/sec)
 res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $
                       ,/xlog,/ylog,xr=xr,yr=yr,xtit=xtit,ytit=ytit,title=tit $
@@ -319,58 +342,28 @@ if keyword_set(wait) then begin
    wait,wait
 end
-;=== MAKE THE FINAL PLOT
-IF keyword_set(postscript) THEN BEGIN
-;    dir_ps='./'
-    mydevice=!d.name
-    set_plot,'PS'
-;    ps_file=dir_ps+postscript
-    ps_file=postscript
-    device,filename=ps_file,/color
-ENDIF
-
-;stop
-
+;;===================================
+;;=== MAKE THE FINAL PLOT
+;;===================================
 IF !dustem_noobj THEN BEGIN
-  dustemwrap_plot_noobj,*(*!dustem_fit).CURRENT_PARAM_VALUES,st=dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (Final fit)'
+  dustem_plot_noobj,*(*!dustem_fit).CURRENT_PARAM_VALUES,st=dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (Final fit)'
  ENDIF ELSE BEGIN
   dustemwrap_plot,*(*!dustem_fit).CURRENT_PARAM_VALUES,st=dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (Final fit)'
 ENDELSE
-
-IF keyword_set(postscript) THEN BEGIN
-  device,/close
-  set_plot,mydevice
-  message,'Wrote '+ps_file,/info
-ENDIF
-
 if keyword_set(wait) then begin
    message,'Made the plot of the final results',/info
    wait,wait
 end
+;;===================================
+;;=== WRITE OUT THE DUSTEMWRAP FITTING RESULTS IN A BINARY FITS FILE
+;;===================================
 IF keyword_set(fits_save) THEN BEGIN
    message,'Writing out results structure: '+fits_save,/info
    dustem_write_fits_table,filename=fits_save,help=help
-;=== At this point, you could erase all dustem system variables, or exit idl... all the
-;=== information needed to recover the results and remake the plots has been saved in the FITS table
-
-;=== Moved the following to dustem_fits_io_example
-  ;; dustem_read_fits_table,filename=fits_save,dustem_st=dustem_spectra_st
-  ;; ;stop
-  ;; ;!dustem_show=ptr_new(*!dustem_data)   ;test
-  ;; !dustem_noobj=1
-  ;; ;==== plot result taken from the saved fits table
-  ;; res=*(*!dustem_fit).CURRENT_PARAM_VALUES
-  ;; IF !dustem_noobj THEN BEGIN
-  ;;   ;dustemwrap_plot_noobj,res,st=dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
-  ;;   dustemwrap_plot_noobj,res,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
-  ;; ENDIF ELSE BEGIN
-  ;;   ;dustemwrap_plot,res,st=dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
-  ;;   dustemwrap_plot,res,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
-  ;; ENDELSE
-  ;; ;stop
-
-  IF keyword_set(wait) THEN BEGIN
+;=== At this point, you could erase all dustem system variables, or exit IDL... all the
+;=== information needed to recover the fitting results and remake any plots has been saved in the FITS table
+IF keyword_set(wait) THEN BEGIN
      message,'Saved the results as FITS in the file: '+fits_save,/info
      wait,wait
   ENDIF