updated examples to conform with webpages

Annie Hughes
1 parent 553fad92
Showing 4 changed files with 92 additions and 216 deletions Show diff stats
src/idl/dustem_fit_ext_pol_example.pro
src/idl/dustem_fit_intensity_example.pro
src/idl/dustem_myisrf_example.pro
src/idl/dustem_run_example.pro
@@ -216,7 +216,7 @@ if not keyword_set(ext_file) then begin
    dustem_set_data, m_fit,m_show,ext,ext
  
    IF keyword_set(wait) THEN BEGIN
-      message,'Finished generating a synthetic EXT from the model',/info
+      message,'Finished generating a synthetic EXT from the model '+use_model,/info
      wait,wait
   ENDIF
  
@@ -138,10 +138,8 @@ dustem_define_la_common
  
 ;=== 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 several plug-ins:
-;=== (i) continuum due to a blackbody  
-;=== (ii) free-free emission
-;=== (iii) synchrotron emission
+;=== intensity of the dust-heating radiation field, as well as a plug-in
+;=== for synchrotron emission
 ;=== The free parameters are all lower-bounded at zero.
 ;=== use_model='DBP90' ; you should specify this above!
 pd = [ $
 PRO dustem_myisrf_example,model=model $
                          ,sed_file=sed_file $
+                         ,isrf_file=isrf_file $
                          ,Nitermax=Nitermax $
                          ,postscript=postscript $
                          ,fits_save_and_restore=fits_save_and_restore $
@@ -11,11 +12,11 @@ PRO dustem_myisrf_example,model=model $
 ; NAME:
 ;    dustem_fit_myisrf_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 modify_isrf plugin and not to do science -- the fit
-; obtained by running this example is likely to be poor.
+; 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 dustem_plugin_modify_isrf 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 
@@ -30,11 +31,12 @@ PRO dustem_myisrf_example,model=model $
 ;    DustEMWrap, Distributed, High-Level, User Example
 ;
 ; CALLING SEQUENCE:
-;    dustem_fit_myisrf_example[,model=][sed_file=][,postscript=][,Nitermax=][,fits_save_and_restore=][,/help,/wait,/verbose]
+;    dustem_fit_myisrf_example[,model=][,sed_file=][,isrf_file=][,postscript=][,Nitermax=][,fits_save_and_restore=][,/help,/wait,/verbose]
 ;
 ; INPUTS:
-;    None
-;
+;    sed_file = file in .xcat format containing observational SED
+;    isrf_file = text file describing ISRF
+;  
 ; OPTIONAL INPUT PARAMETERS:
 ;    None
 ;
@@ -92,7 +94,7 @@ PRO dustem_myisrf_example,model=model $
  
  
 IF keyword_set(help) THEN BEGIN
-  doc_library,'dustem_stellarpopisrf_example'
+  doc_library,'dustem_myisrf_example'
   goto,the_end
 END
  
@@ -104,6 +106,8 @@ ENDELSE
  
  
 use_verbose=0
+use_polarization=0   ; initialize Dustemwrap in no polarization mode 
+use_window=2          ; default graphics window number to use for plotting the results
 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
@@ -114,155 +118,25 @@ dustem_define_la_common
 ;=== 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_modify_isrf_1']  ;amplitude of user-defined ISRF (here G0 because mathis is chosen)
-        rv =   [4.3e-4, 4.7e-4,6.4e-3,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,-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_modify_isrf_1'] ;Amplitude of user-defined ISRF
-             
-        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 ;takes at least 16 iterations to find minima.
-          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_modify_isrf_1',   $ ;amplitude of user-defined ISRF
-              'dustem_plugin_stellar_population_O5V3'   $ ;distance to stellar pop O5
-              ] 
-                     
-        rv =   [5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,5.,1.];,10,1.,10.,1.]
-        iv =   rv+ [5.e-4, 5.e-4,5.e-4,5.e-3,5.e-3,5.,2]
-        
-        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
-            '(*!dustem_params).G0'                  $
-            ]
-        ;initial parameter values for fixed parameters
-        fiv=[ $
-          30., $      ;fixed stellar population number
-          10.  $
-            ]
-        
-      END
-      'MC10':BEGIN
-        
-        pd = [ $
-              'dustem_plugin_modify_isrf_1', $  ;amplitude of user-defined ISRF
-              '(*!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_modify_isrf_1']                 ;amplitude of user-defined ISRF
-                           
-        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
+
+;; ;=== AN EXAMPLE FOR DL07 (or DL01 or WD01_RV5P5B)
+;; ;=== Here we fit the PAH abundances of the model, and the amplitude
+;; ;=== of a user-supplied ISRF.
+;; ;=== The other grain parameters are fixed to their default values.
+;; ;=== 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).grains(0).mdust_o_mh',$     ;PAH0 mass fraction
+     '(*!dustem_params).grains(1).mdust_o_mh',$         ;PAH1 mass fraction
+     'dustem_plugin_modify_isrf_1', $                   ; Ampitude of user-defined ISRF
+     'dustem_plugin_synchrotron_2']                     ;Synchrotron amplitude at 10 mm
+    iv =   [5.4e-4, 5.4e-4,0.1, 0.01]
+   Npar=n_elements(pd)
+   ulimed=replicate(0,Npar)
+   llimed=replicate(1,Npar)
+   llims=replicate(0.,Npar)
+;   fpd=['(*!dustem_params).G0']
+;   fiv=[1.e-10]
  
  
 if keyword_set(wait) then begin
@@ -271,26 +145,53 @@ if keyword_set(wait) then begin
 end
  
  
-dustem_init,model=use_model
-
+dustem_init,model=use_model,polarization=use_polarization
 !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.
-
+!dustem_verbose=use_verbose
+
+;== Fill a structure with default inputs to the model. This
+;== includes the default ISRF, which we want here just for plotting a
+;comparison to the ISRF that we construct.
+dir_in=!dustem_soft_dir  
+st_model=dustem_read_all(dir_in)
+
+;=== CREATE A DIFFERENT ISRF AND DUMP TO FILE
+if not keyword_set(isrf_file) then begin
+   mywaves=[]
+   my_isrf_file='./myisrf_3bb_8000K.dat'
+   myisrf=dustem_create_rfield([8000,8000,8000],wdil=[1.d-11,1.d-11,1.d-11],isrf=0,fname=my_isrf_file,x=mywaves)
+;myisrf=dustem_create_rfield([3500,4000,5500,6500,8000],wdil=[1.d-13,1.d-13,1.d-13,1.d-13,1.d-13],isrf=1,fname='./myisrf_mathis_5bb.dat',x=mywaves)
+;myisrf=dustem_create_rfield([3500,4000,5500,6500,8000],wdil=[1.d-13,1.d-13,1.d-13,1.d-13,1.d-13],isrf=2,fname='./myisrf_habing_5bb.dat',x=mywaves)
+;myisrf=dustem_create_rfield([3500,4000,5500,6500,8000],wdil=[1.d-13,1.d-13,1.d-13,1.d-13,1.d-13],isrf=0,fname='./myisrf_5bb.dat',x=mywaves)
+   isrf_file=my_isrf_file
+endif
+
+
+;=== PLOT AS A SANITY CHECK
+xtit=textoidl('\lambda (\mum)')
+ytit=textoidl('log ISRF [4 \pi I_\nu (erg/cm^2/s/Hz)]')
+tit='DUSTEM MYISRF EXAMPLE'
+yr=[min(st_model.isrf.isrf) < min(myisrf),5.*(max(st_model.isrf.isrf)>max(myisrf))]
+xr=[0.001,1e6]
+window,use_window,xs=600,ys=400,tit=tit & use_window=use_window+2
+cgplot,st_model.isrf.lambisrf,st_model.isrf.isrf $
+       ,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog $
+       ,title=tit,xtit=xtit,ytit=ytit,color=cgcolor('black'),/nodata
+oplot,st_model.isrf.lambisrf,st_model.isrf.isrf,col=cgcolor('black'),thick=2
+oplot,mywaves,myisrf,col=cgcolor('red'),thick=2
+al_legend,/top,/right,clear=0,box=0 $
+          ,['default ISRF','my ISRF'] $
+          ,linsize=0.5,lines=0 $
+          ,color=[cgcolor('black'),cgcolor('red')] $
+          ,thick=2,charsize=1.3
+
+;=== READ THE ISRF
+if keyword_set(isrf_file) then !dustem_isrf_file=ptr_new(isrf_file)
+
+
+;=== READ EXAMPLE SED DATA: 
 dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'
-file=dir+'example_SED_2.xcat'
+file=dir+'example_SED_1.xcat'
 IF keyword_set(sed_file) THEN file=sed_file
 sed=read_xcat(file,/silent)
  
@@ -300,36 +201,15 @@ if keyword_set(wait) then begin
 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
-
+;== SET THE OBSERVATIONAL STRUCTURE
+;== sed is passed twice -- the first occurrence is the SED that you
+;== wish to fit, the second occurrence is the SED that you wish to visualise. 
 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
@@ -337,23 +217,19 @@ if keyword_set(wait) then begin
 end
  
 ;== RUN THE FIT
-tol=1.e-18
+tol=1.e-16
 use_Nitermax=5  ;maximum number of iterations. 
 IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax
  
-xr_m = [1.,5e5]
-yr_m = [5e-8,1.00e6]
-
-
-tit='Spectral Energy Distribution'
+xr = [1.,5e5]
+yr = [5e-8,1.00e6]
+tit='MYISRF EXAMPLE'
 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 $
+                      ,/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,show_plot=show_plot)
 t2=systime(0,/sec)
@@ -411,7 +287,7 @@ IF keyword_set(fits_save_and_restore) THEN BEGIN
   ENDIF
 ENDIF
  
-message,'Finished dustem_fit_intensity_example',/info
+message,'Finished dustem_myisrf_example',/info
  
  
 the_end:
@@ -25,6 +25,7 @@ PRO dustem_run_example, model $
 ;                 polarisation. See Tables 2 and 3 of that paper for details.
 ;           'G17_ModelB' model B from Guillet et al (2018)
 ;           'G17_ModelC' model C from Guillet et al (2018)
+;           'G17_ModelD' model D from Guillet et al (2018)
 ; OPTIONAL INPUT PARAMETERS:
 ;       show = vector specifying the type of plot(s) to show. Possible
 ;       values (for all dust models) are
@@ -97,6 +98,7 @@ dir_in=!dustem_soft_dir
 ;== Fill the following structure contains with default inputs to the model
 st_model=dustem_read_all(dir_in)
  
+
 ;=== You could modify the model inputs here by directly editing the values in the
 ;=== st_model structure. To inspect contents, type IDL> help,st_model,/str
 ;=== e.g.