Uploading stellarpopisrf_example and some minor debugging

Ilyes Choubani
1 parent 3d0f6b67
Showing 5 changed files with 442 additions and 27 deletions Show diff stats
src/idl/dustem_fit_ext_pol_example.pro
src/idl/dustem_fit_intensity_example.pro
src/idl/dustem_plot_dataset.pro
src/idl/dustem_plugin_modify_isrf.pro
src/idl/dustem_stellarpopisrf_example.pro
@@ -388,7 +388,6 @@ dustem_init_plugins, pd,fpd=fpd
  
 dustem_init_fixed_params,fpd,fiv
  
-
 ;generation of the fake data (since dustem_set_data takes ONE structure containing all the info)
 !dustem_psi_ext=25
  
@@ -396,8 +395,6 @@ ext.EXT_I = dustem_compute_ext(rv,st=st)
 toto = dustem_compute_stokext(rv,st=st,qext_spec=dustem_qext,uext_spec=dustem_uext)
 ext.EXT_Q = toto[0]
 ext.EXT_U = toto[1]
-;stop
-
  
 ;==== JPB
 ;below m_fit does not exist, so there is no emission data to plot, and therefore, no emission data to plot.
@@ -420,6 +417,7 @@ yr_x=[1.00E-10,10]
 tit='Spectral Energy Distribution'
 ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
 xtit=textoidl('\lambda (\mum)')
+
 ;Set show_plot to 0 to hide plot
 ;Commented or set to 1 is the same since !dustem_show_plot (existing sysvar) is initialized to 1 in dustem_init
 ;show_plot = 0 
@@ -336,16 +336,6 @@ if keyword_set(wait) then begin
 end
  
  
-
-;IC: -
-;The final plot is automatically plotted
-;and a (Final run) string is automatically added
-;to modify this go to dustemwrap_plot and search for the string '(FINAL_RUN)' 
-;This is true for OOP plotting.
-;However, even thought the last iteration for noobj plotting is plotted, the '(FINAL_RUN)' string isn't applied
-;to modify this, got to dustemwrap_plot_noobj and search for the '(RUNNING)' string. 
-;an extra condition using !dustem_end (signaling end of fit) will have to be used. 
-
 ;=== MAKE THE FINAL PLOT
 IF keyword_set(postscript) THEN BEGIN
     dir_ps='./'
@@ -353,6 +343,7 @@ IF keyword_set(postscript) THEN BEGIN
     ps_file=dir_ps+postscript
     device,filename=ps_file,/color
 ENDIF
+
 IF !dustem_noobj THEN BEGIN
   dustemwrap_plot_noobj,*(*!dustem_fit).CURRENT_PARAM_VALUES,dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (Final fit)'
  ENDIF ELSE BEGIN
@@ -370,13 +361,6 @@ if keyword_set(wait) then begin
    wait,wait
 end
  
-;IC: Last iteration is automatically plotted.
-;I think modifications will have to be applied to dustemwrap_plot...
-;idea: create a system variable that signals if the fits file has been successfully saved and use it in the dustemwrap_plot(obj/noobj) routine(s) 
-;to add the '(From saved FITS file)' string to the title
-
-;stop
-
 IF keyword_set(fits_save_and_restore) THEN BEGIN
    message,'Writing out results structure: '+fits_save_and_restore,/info
    dustem_write_fits_table,filename=fits_save_and_restore,help=help
@@ -3218,13 +3218,25 @@ if keyword_set(dataset) then begin
  
                     Nfpar = n_elements((*(*!dustem_fit).fixed_param_descs))
                     res = (*(*!dustem_fit).fixed_param_init_values)
+                    
+                    iii=0
+                    FOR i=0L,Nfpar-1 DO BEGIN
+                        
+                        parameter_description = (*(*!dustem_fit).fixed_param_descs)[i]     
+                        parameter_type = dustem_parameter_description2type(parameter_description,string_name=string_name)
+                        if parameter_type NE 'PLUGIN' then iii+=1 
+                        
+                    ENDFOR
+                    
                     k = 0
-                    xxpos = 0.0
-                    yypos = 0.9
-                    ;xxpos = 0.0 ;might need to be modified
-                    cgtext, xxpos, yypos, '-FROZEN-', charsize=-1, width=thiswidth, /normal  
-                    xxpos=(1-thiswidth)/2
-                    xyouts, xxpos, yypos, '-FROZEN-', charsize=1.2, /normal, color=0
+                    if iii ne 0 then begin
+                        xxpos = 0.0
+                        yypos = 0.9
+                        ;xxpos = 0.0 ;might need to be modified
+                        cgtext, xxpos, yypos, '-FROZEN-', charsize=-1, width=thiswidth, /normal  
+                        xxpos=(1-thiswidth)/2
+                        xyouts, xxpos, yypos, '-FROZEN-', charsize=1.2, /normal, color=0
+                    endif
  
                     FOR i=0L,Nfpar-1 DO BEGIN
  
@@ -65,7 +65,7 @@ ENDIF
  
 if isa(!dustem_isrf_file) then user_ISRF = dustem_read_isrf(*!dustem_isrf_file)
  
-out=amp*User_ISRF     
+out=amp*User_ISRF.isrf     
  
 the_end:
 scope='USER_ISRF'
@@ -0,0 +1,421 @@
+PRO dustem_stellarpopisrf_example,model=model $
+                          ,sed=sed $
+                          ,Nitermax=Nitermax $
+                          ,postscript=postscript $
+                          ,fits_save_and_restore=fits_save_and_restore $
+                          ,wait=wait $
+                          ,verbose=verbose $
+                          ,help=help
+
+;+
+; NAME:
+;    dustem_fit_stellarpopisrf_example  
+;
+; PURPOSE:
+; This routine is an example of how to fit an observational SED
+; (StokesI only) with DustEM and DustEMWrap. The objective is to
+; illustrate how to use the stellar population plugin and not to do science -- the fit
+; obtained by running this example is likely to be poor.
+;  
+; For this example, the code uses the SED in the file example_SED_2.xcat,
+; which is distributed in the Data/EXAMPLE_OBSDATA/ directory 
+;
+; The example SED has Stokes I photometric data points from
+; IRAC, MIPS and IRAS. Examples illustrating running DustEMWrap to
+; fit spectral data, polarisation data and extinction data
+; are provided in other _example routines in the src/idl/
+; directory. See the DustEMWrap User Guide for more information.
+;
+; CATEGORY:
+;    DustEMWrap, Distributed, High-Level, User Example
+;
+; CALLING SEQUENCE:
+;    dustem_fit_stellarpopisrf_example[,model=][sed_file=][,postscript=][,Nitermax=][,fits_save_and_restore=][,/help,/wait,/verbose]
+;
+; INPUTS:
+;    None
+;
+; OPTIONAL INPUT PARAMETERS:
+;    None
+;
+; OUTPUTS:
+;    None
+;
+; OPTIONAL OUTPUT PARAMETERS:
+;    Plots, results structure in binary FITS table format
+;
+; ACCEPTED KEY-WORDS:
+;    model = specifies the interstellar dust mixture used by DustEM
+;           'MC10' model from Compiegne et al 2010 
+;           'DBP90' model from Desert et al 1990
+;           'DL01' model from Draine & Li 2001
+;           'WD01_RV5p5B' model from Weingartner & Draine 2002 with Rv=5.5
+;           'DL07' model from Draine & Li 2007
+;           'J13' model from Jones et al 2013, as updated in
+;                 Koehler et al 2014
+;    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_2.xcat' is used.  
+;    postscript = if set, final plot is saved as postscript in the
+;                 current working directory
+;    Nitermax = maximum number of fit iterations. Default is 5.
+;    fits_save_and_restore = if set, save the fit results in a binary
+;               FITS file. The code then restore this file and plots
+;               the results using the saved results information.
+;    help      = if set, print this help
+;    wait      = if set, wait this many seconds between each step of
+;                the code (for illustration purposes)
+;    verbose      = if set, subroutines will run in verbose mode
+;
+; COMMON BLOCKS:
+;    None
+;
+; SIDE EFFECTS:
+;    None
+;
+; RESTRICTIONS:
+;    The DustEM fortran code must be installed
+;    The DustEMWrap IDL code must be installed
+;
+; PROCEDURES AND SUBROUTINES USED:
+;    
+;
+; EXAMPLES
+;    dustem_fit_stellarpopisrf_example
+;    dustem_fit_stellarpopisrf_example,Nitermax=1,fits_save_and_restore='/tmp/mysavefile.fits'
+;    dustem_fit_stellarpopisrf_example,model='DBP90'
+;
+; MODIFICATION HISTORY:
+;    Evolution details on the DustEMWrap gitlab.
+;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  
+;-
+
+
+IF keyword_set(help) THEN BEGIN
+  doc_library,'dustem_stellarpopisrf_example'
+  goto,the_end
+END
+
+IF keyword_set(model) THEN BEGIN
+  use_model=strupcase(model)
+ENDIF ELSE BEGIN
+  use_model='J13'    ;Default is last dustem model
+ENDELSE
+
+
+
+
+
+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
+
+dustem_define_la_common
+
+;=== Set the (model-dependent) parameters that you want to fit
+;=== Refer to the DustEM and DustEMWrap userguides for an explanation
+;    of the different grain types
+
+CASE use_model OF
+      'DBP90':BEGIN
+         pd = [ $
+             '(*!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', $    ;amCBEx
+             'dustem_plugin_stellar_population_O7V3',$ ;distance to stellar population
+             'dustem_plugin_stellar_population_O5V4']  ;number of stars
+        rv =   [4.3e-4, 4.7e-4,6.4e-3,1.5,100.] ;rv because we're generating fake data and try to retrieve the original values
+        
+        iv = rv+[5.00E-4,5.00E-4,5.00E-4,-1.,-30] ;###setting the initial parameter vector
+        
+        Npar=n_elements(pd)
+        ulimed=replicate(0,Npar)
+        llimed=replicate(1,Npar)
+        llims=replicate(0.,Npar)         
+        
+        
+        fpd=[ $
+          'dustem_plugin_stellar_population_O7V4',$ ;number of stars of stellar population' ,  $     ;multiplicative factor to total ISRF
+          'dustem_plugin_stellar_population_O5V3'$ ; stellar distance
+            ]
+        ;initial parameter values for fixed parameters
+        fiv=[ $
+          2. , $        ;number of stars of pop O7
+          45$      ;dustance of pop O5
+            ]
+
+     END
+      'DL01':BEGIN
+         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
+             'dustem_plugin_stellar_population_O5V4'] ;number of stars of pop O5
+             
+        rv =   [5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,50]
+        iv =   rv+[5.e-4, 5.e-4,5.e-5,5.e-3,5.e-3,100]
+        
+        fpd=[ $
+            'dustem_plugin_stellar_population_O5V3'$ ; stellar distance
+            ]
+        ;initial parameter values for fixed parameters
+        fiv=[ $
+          45$      ;fixed stellar distance
+            ]
+
+        
+        Npar=n_elements(pd)
+        ulimed=replicate(0,Npar)
+        llimed=replicate(1,Npar)
+        llims=replicate(0.,Npar)
+      END
+      'WD01_RV5p5B':BEGIN
+;;  ***COMMENT AH***
+;; we need to implement this, or remove         
+         message, 'WD01 model not yet implemented in DustEMWrap',/info
+;;  ***END COMMENT AH***
+       END
+      'DL07':BEGIN
+          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
+              'dustem_plugin_stellar_population_O5V3'   $ ;distance to stellar pop O5V
+              ] 
+                     
+        rv =   [5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,5.];,10,1.,10.,1.]
+        iv =   rv+ [5.e-4, 5.e-4,5.e-4,5.e-3,5.e-3,25.]
+        
+        Npar=n_elements(pd)
+        ulimed=replicate(0,Npar)
+        llimed=replicate(1,Npar)
+        llims=replicate(0.,Npar)
+        
+        fpd=[ $
+            'dustem_plugin_stellar_population_O5V4'$ ; stellar number of pop O5V
+            ]
+        ;initial parameter values for fixed parameters
+        fiv=[ $
+          30$      ;fixed stellar population number
+            ]
+        
+      END
+      'MC10':BEGIN
+        
+        pd = [ $
+              'dustem_plugin_stellar_population_O5V4', $  ;stellar number of pop O5V
+              '(*!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' $     ;aSilx
+              ]
+        ;initial parameter values for parameters to be fitted   
+        rv =   [ $  
+          50, $    
+          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 aSilx
+          ]
+          
+          iv = rv + [-30,5.E-4,5.E-4,5.E-4,5.E-3,5.E-3]
+          
+          fpd = ['dustem_plugin_stellar_population_O5V3'] ;distance to stellar population O5V
+          fiv  = [25]
+              
+     
+        Npar=n_elements(pd)
+        
+        ulimed=replicate(0,Npar)
+        llimed=replicate(1,Npar)
+        llims=replicate(0.,Npar)
+        
+      END
+     'J13':BEGIN
+        pd = [ $
+             '(*!dustem_params).grains(1).mdust_o_mh',$     ;CM20 mass fraction
+             '(*!dustem_params).grains(0).mdust_o_mh',$     ;CM20 mass fraction
+             '(*!dustem_params).grains(2).mdust_o_mh', $    ;aPyM5
+             '(*!dustem_params).grains(3).mdust_o_mh', $    ;aPyM5
+             'dustem_plugin_stellar_population_O5V4']         ;stellar number of pop O5V
+                           
+        rv =   [7.8e-4, 7.8e-4,1.65e-4,1.45e-3,35]
+        iv = rv+ [5.e-4,5.e-4,5.e-4,5.e-3,35]
+        
+        Npar=n_elements(pd)
+        ulimed=replicate(0,Npar)
+        llimed=replicate(1,Npar)
+        llims=replicate(0.,Npar)
+        
+        fpd = ['dustem_plugin_stellar_population_O5V3'] ;distance to stellar population O5V
+        fiv  = [5]
+      
+      END
+      'ELSE':BEGIN
+         message,'model '+model+' unknown',/continue
+         message,'Known models are MC10,DBP90,DL01,DL07,J13,G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD',/continue
+         stop
+      END
+ENDCASE
+
+
+if keyword_set(wait) then begin
+   message,'Finished setting dust model and plug-in parameters: '+use_model,/info
+   wait,wait
+end
+
+
+dustem_init,model=use_model
+
+!dustem_nocatch=1
+!dustem_verbose=1
+;!dustem_show_plot=1
+!EXCEPT=2 ;so I can locate the plotting error...
+;!dustem_dim=1 ; this option is to dim the stellar population ISRF with the current Dustem extinction
+;We're fitting total optical depths so the dimmed-ISRF scenario is a reasonable assumption.
+;Because of this the final parameter values are slightly lower than the initial (real) values.
+;But this is because the first run was not extinct with the current dustem extinction 
+;since it comes prior to the dustem run. Meaning the first ISRF has to be extinct to see if the paramters are retrieved.
+
+
+!dustem_isrf_file=ptr_new(!dustem_soft_dir+'data/ISRF_MATHIS.DAT')
+
+
+;=== READ EXAMPLE DATA: EXTINCTION
+
+;NB:HERE WE ARE NOT FITTING THE DATA. WE'RE JUST USING THE EXISTING FILTERS.
+
+dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'
+file=dir+'example_SED_2.xcat'
+sed=read_xcat(file,/silent)
+
+if keyword_set(wait) then begin
+   message,'Finished reading SED data: '+file,/info
+   wait,wait
+end
+
+
+;=== GENERATE EXAMPLE DATA: EXTINCTION
+
+;=== initializing IQU and associated errors to avoid problems when checking SED in dustem_set_data.pro
+;sed.StokesI=1.
+;setting the errors (well equating all the tags to the sigmaii array but this is an initialization)
+for i=4l,n_tags(sed)-1 do begin
+    sed.(i) = sed.sigmaii
+endfor
+
+;###first call to dustem_set_data so that the !dustem_data tags that are used in the 'compute_' procedures are set
+
+dustem_set_data,sed,sed
+
+; ;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS THAT WILL BE
+; ;== ADJUSTED DURING THE FIT
+; dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
+; 
+; ;== INITIALIZE ANY PLUGINS
+; dustem_init_plugins, pd,fpd=fpd 
+; 
+; dustem_init_fixed_params,fpd,fiv
+
+dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims
+
+;Generating the fake emission data
+sed.StokesI = dustem_compute_sed(rv,st=st)
+
+;== SET THE OBSERVATIONAL STRUCTURE
+dustem_set_data, sed,sed
+
+
+if keyword_set(wait) then begin
+   message,'Finished initializing DustEMWrap, including plugins and fixed parameters',/info
+   wait,wait
+end
+
+;== RUN THE FIT
+tol=1.e-18
+use_Nitermax=20  ;maximum number of iterations. 
+IF keyword_set(itermax) THEN use_Nitermax=itermax
+
+xr_m = [1.,5e5]
+yr_m = [5e-8,1.00e6]
+
+
+tit='Spectral Energy Distribution'
+ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
+xtit=textoidl('\lambda (\mum)')
+;Set show_plot to 0 to hide plot
+;Commented or set to 1 is the same since !dustem_show_plot (existing sysvar) is initialized to 1 in dustem_init
+;show_plot = 0 
+t1=systime(0,/sec)
+res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $
+                      ,/xlog,/ylog,xr_m=xr_m,yr_m=yr_m,xtit=xtit,ytit=ytit,title=tit $
+                      ,legend_xpos=legend_xpos,legend_ypos=legend_ypos $
+                      ,errors=errors,chi2=chi2,rchi2=rchi2,show_plot=show_plot)
+t2=systime(0,/sec)
+
+
+if keyword_set(wait) then begin
+   message,'Finished running DustEMWrap, using Niters: '+strtrim(string(use_Nitermax),2),/info
+   message,'Time taken [sec]: '+sigfig(t2-t1,2,/sci),/info
+   wait,wait
+end
+
+
+;=== MAKE THE FINAL PLOT
+IF keyword_set(postscript) THEN BEGIN
+    dir_ps='./'
+    set_plot,'PS'
+    ps_file=dir_ps+postscript
+    device,filename=ps_file,/color
+ENDIF
+
+IF !dustem_noobj THEN BEGIN
+  dustemwrap_plot_noobj,*(*!dustem_fit).CURRENT_PARAM_VALUES,dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (Final fit)'
+ ENDIF ELSE BEGIN
+  dustemwrap_plot,*(*!dustem_fit).CURRENT_PARAM_VALUES,dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (Final fit)'
+ENDELSE
+
+IF keyword_set(postscript) THEN BEGIN
+  set_plot,'X'
+  device,/close
+  message,'Wrote '+ps_file,/info
+ENDIF
+
+if keyword_set(wait) then begin
+   message,'Made the plot of the final results',/info
+   wait,wait
+end
+
+IF keyword_set(fits_save_and_restore) THEN BEGIN
+   message,'Writing out results structure: '+fits_save_and_restore,/info
+   dustem_write_fits_table,filename=fits_save_and_restore,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
+  ;stop
+  dustem_read_fits_table,filename=fits_save_and_restore,dustem_spectra_st=dustem_spectra_st
+  ;==== plot result taken from the saved fits table
+  res=*(*!dustem_fit).CURRENT_PARAM_VALUES
+  IF !dustem_noobj THEN BEGIN
+    dustemwrap_plot_noobj,res,dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
+  ENDIF ELSE BEGIN
+    dustemwrap_plot,res,dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
+  ENDELSE
+  IF keyword_set(wait) THEN BEGIN
+     message,'Saved the results as FITS in the file: '+fits_save_and_restore+', and made a plot using the data in this file',/info
+     wait,wait
+  ENDIF
+ENDIF
+
+message,'Finished dustem_fit_intensity_example',/info
+
+
+the_end:
+
+END