first commit

Annie Hughes
1 parent 52a3c8dc
Showing 2 changed files with 476 additions and 0 deletions Show diff stats
src/idl/dustem_fit_sed_ext_stokesi_example.pro
src/idl/dustem_fitsio_example.pro
@@ -0,0 +1,300 @@
+PRO dustem_fit_sed_ext_stokesi_example,model=model $
+                                       ,sed_file=sed_file $
+                                       ,ext_file=ext_file $
+                                       ,Nitermax=Nitermax $
+                                       ,fits_save=fits_save $
+                                       ,wait=wait $
+                                       ,noobj=noobj $
+                                       ,help=help $
+                                       ,verbose=verbose
+
+;+
+; NAME:
+;    dustem_fit_sed_ext_stokesi_example  
+;
+; PURPOSE: This routine is an example of how to fit an observational
+; data in emission + extinction SED (StokesI only) with DustEM and
+; DustEMWrap.
+;  
+; See the DustEMWrap User Guide for more information.
+;
+; CATEGORY:
+;    DustEMWrap, Distributed, High-Level, User Example
+;
+; CALLING SEQUENCE:
+;    dustem_fit_sed_ext_stokesi_example[,model=][sed_file=][ext_file=][,Nitermax=][,fits_save=][,/help,/wait,/verbose,/noobj]
+;
+; 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. See userguide or dustem_test_model_exists.pro
+;            for more details about available models in current
+;            release. Default for this example is J13.
+;    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_3.xcat' is used to
+;          identify the filters, and the model is used to generate a
+;          predicted SED.  
+;    ext_file = string naming the path to text file in .xcat format that
+;          describes the observational extinction SED. If not set, the file
+;          'Data/EXAMPLE_OBSDATA/Mathis90Fitz99_DISM_NH20.xcat' is
+;          used to identify the wavelength range, and the dust model
+;          is used to generate the corresponding extinction. 
+;    Nitermax = maximum number of fit iterations. Default is 5.
+;    fits_save = 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
+;    noobj     = if set, runs with no object graphics
+;
+; 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_sed_ext_stokesi_example
+;    dustem_fit_sed_ext_stokesi_example,Nitermax=1,fits_save='/tmp/mysavefile.fits'
+;    dustem_fit_sed_ext_stokesi_example,model='MC10'
+;
+; MODIFICATION HISTORY:
+;    Written by JPB Apr-2011
+;    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_fit_sed_ext_stokesi_example'
+  goto,the_end
+END
+
+IF keyword_set(model) THEN BEGIN
+  use_model=strupcase(model)
+ENDIF ELSE BEGIN
+  use_model='J13'    
+ENDELSE
+
+exists=dustem_test_model_exists(use_model)
+if exists ne 1 then $
+   message,'Unknown dust model: '+use_model
+
+use_polarization=0   ;default is no polarization in models
+use_window=2       ; default graphics window number to use for plotting the results
+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
+
+;=== Set the (model-dependent) parameters that you want to fit (pd),
+;=== their initial values (iv)
+;=== and whether they are bounded (ulimed,llimed,llims,ulims).
+;=== Fixed parameters (fpd) and their values (fiv) are also set here.
+;=== Refer to the DustEM and DustEMWrap User guides for an explanation
+;=== of the physical meaning of dust model and plug-in parameters, and
+;=== how to specify them.
+
+;; ;=== 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:
+;; ;=== (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.
+
+;=== If you choose a different dust model, you will need to adjust the
+;=== parameter information accordingly.
+
+pd = [ $
+     '(*!dustem_params).G0' $                                   ;G0
+     ,'(*!dustem_params).grains(0).mdust_o_mh'$                  ;CM20 -- power law size distribution
+     ,'(*!dustem_params).grains(1).mdust_o_mh'$                 ;CM20 -- logN size distribution
+     ,'(*!dustem_params).grains(2).mdust_o_mh' $                ;aPyM5
+     ,'(*!dustem_params).grains(3).mdust_o_mh' $                ;aOlM5
+     ,'dustem_plugin_continuum_2']                              ;Intensity at peak of the BB continuum
+true_vals =   [1.5,1.8e-3, 6.e-4, 3.5e-3, 3.55e-3, 0.003]         
+iv =   [1.2,1.7e-3, 6.3e-4, 2.55e-3, 2.55e-3, 0.001]         
+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.]                       
+    
+
+
+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
+!dustem_nocatch=1
+!dustem_verbose=use_verbose
+IF keyword_set(noobj) THEN !dustem_noobj=1
+!EXCEPT=2 ; for debugging
+
+
+;=== READ TEMPLATE EXTINCTION DATA 
+;=== IF THE USER DOES NOT PROVIDE AN EXTINCTION CURVE FILE , WE READ MATHIS TO INITIALISE THE OBSERVATIONAL STRUCTURE.
+dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'
+use_ext_file=dir+'Mathis90Fitz99_DISM_NH20.xcat' 
+IF keyword_set(ext_file) THEN use_ext_file=ext_file
+ext=read_xcat(use_ext_file,/silent)
+IF not keyword_set(ext_file) THEN BEGIN
+   message, 'You have not provided an extinction curve. I will generate one from the '+use_model+' dust model',/info
+   ext.EXT_I=1.
+ENDIF   
+
+
+
+;=== READ TEMPLATE EMISSION DATA 
+;=== IF THE USER DOES NOT PROVIDE A SED FILE , WE READ AN EXAMPLE_SED FILE TO INITIALISE THE OBSERVATIONAL STRUCTURE.
+dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'
+file=dir+'example_SED_2.xcat'
+if keyword_set(sed_file) then file=sed_file 
+sed=read_xcat(file,/silent)
+IF not keyword_set(sed_file) THEN BEGIN
+   message, 'You have not provided an SED. I will generate one from the '+use_model+' dust model',/info
+   sed.StokesI=1.
+ENDIF   
+
+
+;=== HERE WE SET THE DATA SO THAT THE COMPUTE_ FUNCTIONS BELOW HAVE
+;=== THE NECESSARY FIELDS AND INFORMATION
+dustem_set_data,m_fit=sed,m_show=sed,x_fit=ext,x_show=ext 
+
+;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS FOR THE FIT 
+;== AND ACTIVATE ANY PLUGINS (HENCE WE NEED TO DO THIS BEFORE THE COMPUTE_ CALLS)
+dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims,polarization=use_polarization
+
+if keyword_set(wait) then begin
+   message,'Finished initializing DustEMWrap, including plugins and fixed parameters',/info
+   wait,wait
+end
+
+;== GENERATE THE EXTINCTION DATA USING THE MODEL AND TRUE VALUES
+if not keyword_set(ext_file) then $
+   ext.EXT_I = dustem_compute_ext(true_vals,st=st)
+
+;== GENERATE THE EMISSION DATA USING THE MODEL AND TRUE VALUES
+IF not keyword_set(sed_file) THEN $
+   sed.StokesI = dustem_compute_sed(true_vals,st=st)
+
+;=== RESET THE OBSERVATIONAL STRUCTURE IF WE HAVE GENERATED DATA USING
+;=== A MODEL
+if not keyword_set(ext_file) or  not keyword_set(sed_file) THEN BEGIN
+   dustem_set_data, m_fit=sed,m_show=sed,x_fit=ext,x_show=ext
+   if keyword_set(wait) then begin
+      message,'Finished generating data from the model '+use_model,/info
+      wait,wait
+   end
+end
+
+
+;== INFORMATION TO RUN THE FIT
+tol=1.e-10  ;fit tolerence
+
+;=== INFORMATION TO MAKE THE PLOT
+;=== _x/_extinction means extinction plots, _m/_emissions means emission plots
+xr_x = [0.01,30]
+yr_x = [5.00E-8,10]
+xr_m = [1.,5e5]
+yr_m = [5e-8,1.00e6]
+tit_m='Spectral Energy Distribution' 
+tit_x='Dust Optical Depth' 
+
+;===  RUN THE FIT
+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,xr_x=xr_x,yr_x=yr_x,xtit=xtit,ytit=ytit,tit_m=tit_m,tit_x=tit_x $
+                      ,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),st=dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,tit_m=tit_m+' (Final fit)',tit_x=tit_x+' (Final fit)' 
+ENDIF ELSE BEGIN
+  dustemwrap_plot,(*(*!dustem_fit).CURRENT_PARAM_VALUES),st=dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,tit_m=tit_m+' (Final fit)',tit_x=tit_x+' (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) 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
+
+;;    dustem_read_fits_table,filename=fits_save,dustem_st=dustem_spectra_st
+;; ;==== plot result taken from the saved fits table
+;;    res=*(*!dustem_fit).CURRENT_PARAM_VALUES
+;;   ;If the user prefers non-OOP plotting then reset !dustem_noobj=1 here since dustem_init is called in reads_fits_table.
+;;    IF !dustem_noobj THEN BEGIN
+;;       dustemwrap_plot_noobj,res,st=dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,tit_m=tit_m+' (From Saved FITS file)',tit_x=tit_x+' (From Saved FITS file)'
+;;    ENDIF ELSE BEGIN
+;;       dustemwrap_plot,res,st=dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,tit_m=tit_m+' (From Saved FITS file)',tit_x=tit_x+' (From Saved FITS file)'
+;;    ENDELSE
+
+   IF keyword_set(wait) THEN BEGIN
+      message,'Saved the results as FITS in the file: '+fits_save,/info
+      wait,wait
+   ENDIF
+ENDIF
+
+the_end:
+message,'Finished dustem_fit_sed_ext_stokesi_example',/info
+
+END
@@ -0,0 +1,176 @@
+PRO dustem_fitsio_example,help=help
+
+;+
+; NAME:
+;    dustem_fitsio_example
+;
+; PURPOSE:
+;    Testing the save/restore FITS files
+;
+; CATEGORY:
+;    DustEMWrap, Distributed, Low-Level, Development
+;
+; CALLING SEQUENCE:
+;    dustem_fitsio_example
+;
+; INPUTS:
+;    None
+;
+; OPTIONAL INPUT PARAMETERS:
+;
+; OUTPUTS:
+;
+; OPTIONAL OUTPUT PARAMETERS:
+;
+; ACCEPTED KEY-WORDS:
+;    help = print this help
+;  
+; 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
+; 
+; MODIFICATION HISTORY:
+;    Written by AH Apr-2023
+;    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_fitsio_example'
+  goto,the_end
+END
+
+
+dustem_fit_intensity_example,Nitermax=2,fits_save='/tmp/example_fit_intensity_only.fits'
+;exit/.reset
+dustem_read_fits_table,filename='/tmp/example_fit_intensity_only.fits',dustem_st=dustem_spectra_st
+
+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
+res=*(*!dustem_fit).CURRENT_PARAM_VALUES
+dustemwrap_plot,res,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
+
+
+stop
+dustem_fit_spectro_example,Nitermax=2,fits_save='/tmp/example_fit_spectro.fits'
+;exit/.reset
+dustem_read_fits_table,filename='/tmp/example_fit_spectro_only.fits',dustem_st=dustem_spectra_st
+xr = [1.,5e5]
+yr = [5e-8,1.00e6]
+tit='FIT SPECTRO EXAMPLE' ; plot title
+res=*(*!dustem_fit).CURRENT_PARAM_VALUES
+dustemwrap_plot,res,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
+
+
+stop
+
+dustem_fit_intensity_mbb_example,Nitermax=2,fits_save='/tmp/example_fit_mbb.fits'
+;exit/.reset
+dustem_read_fits_table,filename='/tmp/example_fit_mbb.fits',dustem_st=dustem_spectra_st
+yr=[1.00e-4,1.00E2] ; y-axis limits
+xr=[1.00E0,6.00e4] ; x-axis limits
+tit='MBB FIT 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
+res=*(*!dustem_fit).CURRENT_PARAM_VALUES
+dustemwrap_plot,res,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
+
+
+stop
+dustem_fit_polarization_example,Nitermax=2,fits_save='/tmp/example_fit_polar.fits'
+;exit/.reset
+dustem_read_fits_table,filename='/tmp/example_fit_polar.fits',dustem_st=dustem_spectra_st
+yr=[1.00e-4,1.00E2] ; y-axis limits
+xr=[1.00E0,6.00e4] ; x-axis limits
+tit='FIT POLARISATION SED EXAMPLE' ; plot title
+res=*(*!dustem_fit).CURRENT_PARAM_VALUES
+dustemwrap_plot,res,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
+
+
+stop
+dustem_fit_ext_example,Nitermax=2,fits_save='/tmp/example_fit_ext.fits'
+;exit/.reset
+dustem_read_fits_table,filename='/tmp/example_fit_ext.fits',dustem_st=dustem_spectra_st
+xr=[0.01,30]
+yr=[1.00E-10,10]
+tit='FIT EXTINCTION EXAMPLE' ; plot title
+res=*(*!dustem_fit).CURRENT_PARAM_VALUES
+dustemwrap_plot,res,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
+
+
+stop
+dustem_fit_ext_pol_example,Nitermax=2,fits_save='/tmp/example_fit_ext_polar.fits'
+;exit/.reset
+dustem_read_fits_table,filename='/tmp/example_fit_ext_polar.fits',dustem_st=dustem_spectra_st
+xr=[0.01,30]
+yr=[1.00E-10,50]
+tit='FIT POLARISED EXTINCTION EXAMPLE' ; plot title
+res=*(*!dustem_fit).CURRENT_PARAM_VALUES
+dustemwrap_plot,res,st=dustem_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
+
+stop
+dustem_fit_sed_ext_stokesi_example,Nitermax=2,fits_save='/tmp/example_fit_sed_ext_stokesi.fits'
+;exit/.reset
+dustem_read_fits_table,filename='/tmp/example_fit_sed_ext_stokesi.fits',dustem_st=dustem_spectra_st
+xr_x = [0.01,30]
+yr_x = [5.00E-8,10]
+xr_m = [1.,5e5]
+yr_m = [5e-8,1.00e6]
+tit_m='Spectral Energy Distribution' 
+tit_x='Dust Optical Depth' 
+tit='FIT STOKESI EXTINCTION AND EMISSION EXAMPLE' ; plot title
+res=*(*!dustem_fit).CURRENT_PARAM_VALUES
+dustemwrap_plot,res,st=dustem_st,xr_x=xr_x,xr_m=xr_m,/xstyle,yr_x=yr_x,yr_m=yr_m,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
+
+
+stop
+dustem_fit_sed_ext_pol_example,Nitermax=2,fits_save='/tmp/example_fit_sed_ext_polar.fits'
+;exit/.reset
+dustem_read_fits_table,filename='/tmp/example_fit_sed_ext_polar.fits',dustem_st=dustem_spectra_st
+xr_x = [0.01,30]
+yr_x = [5.00E-8,10]
+xr_m = [1.,5e5]
+yr_m = [5e-8,1.00e6]
+tit_m='Spectral Energy Distribution' 
+tit_x='Dust Optical Depth' 
+tit='FIT POLARISED EXTINCTION AND EMISSION EXAMPLE' ; plot title
+res=*(*!dustem_fit).CURRENT_PARAM_VALUES
+dustemwrap_plot,res,st=dustem_st,xr_x=xr_x,xr_m=xr_m,/xstyle,yr_x=yr_x,yr_m=yr_m,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
+
+stop
+dustem_myisrf_example,Nitermax=2,fits_save='/tmp/example_fit_isrf.fits'
+;exit/.reset
+dustem_read_fits_table,filename='/tmp/example_fit_isrf.fits',dustem_st=dustem_spectra_st
+xr_m = [1.,5e5]
+yr_m = [5e-8,1.00e6]
+tit_m='Spectral Energy Distribution' 
+tit='FIT WITH USER ISRF EXAMPLE' ; plot title
+res=*(*!dustem_fit).CURRENT_PARAM_VALUES
+dustemwrap_plot,res,st=dustem_st,xr_m=xr_m,/xstyle,yr_m=yr_m,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
+
+
+
+
+
+;=== fails when run with only the minimal dustemwrap installed (stellar . txt file is missing)
+dustem_stellarpopisrf_example,Nitermax=3
+
+
+
+
+; END LIST OF TOP-LEVEL ROUTINES
+
+end