OV-GSO-DC / dustem-wrapper_idl

Commit 9d93919edb9cf2d741bebe9bc6cf1640fdd85935

Authored by Annie Hughes
1 parent bd6149e4
Exists in master

first commit

Showing 2 changed files with 475 additions and 0 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_fit_ext_example.pro 0 → 100644
  Diff comments   View file @9d93919
... ... @@ -0,0 +1,234 @@
  1 +PRO dustem_fit_ext_example,model=model $
  2 + ,ext_file=ext_file $
  3 + ,Nitermax=Nitermax $
  4 + ,fits_save_and_restore=fits_save_and_restore $
  5 + ,wait=wait $
  6 + ,noobj=noobj $
  7 + ,help=help $
  8 + ,verbose=verbose
  9 +
  10 +;+
  11 +; NAME:
  12 +; dustem_fit_ext_example
  13 +;
  14 +; PURPOSE:
  15 +; This routine is an example of how to fit an observational extinction SED
  16 +; (StokesI only) with DustEM and DustEMWrap.
  17 +;
  18 +; See the DustEMWrap User Guide for more information.
  19 +;
  20 +; CATEGORY:
  21 +; DustEMWrap, Distributed, High-Level, User Example
  22 +;
  23 +; CALLING SEQUENCE:
  24 +; dustem_fit_ext_example[,model=][ext_file=][,Nitermax=][,fits_save_and_restore=][,/help,/wait,/verbose,/noobj]
  25 +;
  26 +; INPUTS:
  27 +; None
  28 +;
  29 +; OPTIONAL INPUT PARAMETERS:
  30 +; None
  31 +;
  32 +; OUTPUTS:
  33 +; None
  34 +;
  35 +; OPTIONAL OUTPUT PARAMETERS:
  36 +; Plots, results structure in binary FITS table format
  37 +;
  38 +; ACCEPTED KEY-WORDS:
  39 +; model = specifies the interstellar dust mixture used by DustEM
  40 +; 'MC10' model from Compiegne et al 2010
  41 +; 'DBP90' model from Desert et al 1990
  42 +; 'DL01' model from Draine & Li 2001
  43 +; 'WD01_RV5P5B' model from Weingartner & Draine 2002 with Rv=5.5
  44 +; 'DL07' model from Draine & Li 2007
  45 +; 'J13' model from Jones et al 2013, as updated in
  46 +; Koehler et al 2014
  47 +; 'G17_ModelA' model A from Guillet et al (2018). Includes
  48 +; polarisation. See Tables 2 and 3 of that paper for details.
  49 +; 'G17_ModelB' model B from Guillet et al (2018)
  50 +; 'G17_ModelC' model C from Guillet et al (2018)
  51 +; 'G17_ModelD' model A from Guillet et al (2018)
  52 +; ext_file = string naming the path to text file in .xcat format that
  53 +; describes the observational extinction SED. If not set, the file
  54 +; 'Data/EXAMPLE_OBSDATA/Mathis90Fitz99_DISM_NH20.xcat' is used.
  55 +; Nitermax = maximum number of fit iterations. Default is 5.
  56 +; fits_save_and_restore = if set, save the fit results in a binary
  57 +; FITS file. The code then restore this file and plots
  58 +; the results using the saved results information.
  59 +; help = if set, print this help
  60 +; wait = if set, wait this many seconds between each step of
  61 +; the code (for illustration purposes)
  62 +; verbose = if set, subroutines will run in verbose mode
  63 +; noobj = if set, runs with no object graphics
  64 +;
  65 +; COMMON BLOCKS:
  66 +; None
  67 +;
  68 +; SIDE EFFECTS:
  69 +; None
  70 +;
  71 +; RESTRICTIONS:
  72 +; The DustEM fortran code must be installed
  73 +; The DustEMWrap IDL code must be installed
  74 +;
  75 +; PROCEDURES AND SUBROUTINES USED:
  76 +;
  77 +;
  78 +; EXAMPLES
  79 +; dustem_fit_ext_example
  80 +; dustem_fit_ext_example,Nitermax=1,fits_save_and_restore='/tmp/mysavefile.fits'
  81 +; dustem_fit_ext_example,model='DBP90'
  82 +;
  83 +; MODIFICATION HISTORY:
  84 +; Written by AH Oct 2022
  85 +; Evolution details on the DustEMWrap gitlab.
  86 +; See http://dustemwrap.irap.omp.eu/ for FAQ and help.
  87 +;-
  88 +
  89 +
  90 +IF keyword_set(help) THEN BEGIN
  91 + doc_library,'dustem_fit_ext_example'
  92 + goto,the_end
  93 +END
  94 +
  95 +IF keyword_set(model) THEN BEGIN
  96 + use_model=strupcase(model)
  97 +ENDIF ELSE BEGIN
  98 + use_model='J13'
  99 +ENDELSE
  100 +
  101 +known_mdls=['MC10','DBP90','DL01','WD01_RV5P5B','DL07','J13','G17_MODELA','G17_MODELB','G17_MODELC','G17_MODELD']
  102 +test_model = where(known_mdls eq use_model,ct)
  103 +if ct eq 0 then begin
  104 + message,'ISM dust model '+use_model+' unknown',/continue
  105 + message,'Known models are MC10,DBP90,DL01,WD01_RV5P5B,DL07,J13,G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD',/continue
  106 + stop
  107 +end
  108 +
  109 +use_polarization=0. ;default is no polarization in models
  110 +use_window=2 ; default graphics window number to use for plotting the results
  111 +use_verbose=0
  112 +if keyword_set(verbose) then use_verbose=1
  113 +use_Nitermax=5 ; maximum number of iterations for the fit
  114 +IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax
  115 +
  116 +dustem_define_la_common
  117 +
  118 +;=== Set the (model-dependent) parameters that you want to fit
  119 +;=== Refer to the DustEM and DustEMWrap userguides for an explanation
  120 +; of the different grain types
  121 +
  122 +;=== Set the (model-dependent) parameters that you want to fit (pd),
  123 +;=== their initial values (iv)
  124 +;=== and whether they are bounded (ulimed,llimed,llims,ulims).
  125 +;=== Fixed parameters (fpd) and their values (fiv) are also set here.
  126 +;=== Refer to the DustEM and DustEMWrap User guides for an explanation
  127 +;=== of the physical meaning of dust model and plug-in parameters, and
  128 +;=== how to specify them.
  129 +
  130 +;=== Example is provided for J13
  131 +;=== If you choose a different dust model, you will need to adjust the
  132 +;=== parameter information accordingly.
  133 +
  134 +pd = [ $
  135 + '(*!dustem_params).grains(0).mdust_o_mh'$ ;CM20 -- power law size distribution
  136 + ,'(*!dustem_params).grains(1).mdust_o_mh'$ ;CM20 -- logN size distribution
  137 + ,'(*!dustem_params).grains(2).mdust_o_mh' $ ;aPyM5
  138 + ,'(*!dustem_params).grains(3).mdust_o_mh' ] ;aOlM5
  139 +
  140 +iv = [2.9e-3, 6.7e-4, 2.1e-3, 2.8e-3] ; starting guess for the fit
  141 +
  142 +Npar=n_elements(pd)
  143 +ulimed=replicate(0,Npar)
  144 +llimed=replicate(1,Npar)
  145 +llims=replicate(1.e-15,Npar)
  146 +
  147 +;=== Fixed parameters
  148 +fpd=[]
  149 +;==initial parameter values for fixed parameters
  150 +fiv=[]
  151 +
  152 +
  153 +if keyword_set(wait) then begin
  154 + message,'Finished setting dust model and plug-in parameters: '+use_model,/info
  155 + wait,wait
  156 +end
  157 +
  158 +;== INITIALISE DUSTEM
  159 +dustem_init,model=use_model,polarization=use_polarization
  160 +!dustem_nocatch=1
  161 +!dustem_verbose=use_verbose
  162 +IF keyword_set(noobj) THEN !dustem_noobj=1
  163 +!EXCEPT=2 ; for debugging
  164 +
  165 +;=== READ EXTINCTION DATA THAT WE WANT TO FIT
  166 +dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'
  167 +file=dir+'Mathis90Fitz99_DISM_NH20.xcat'
  168 +IF keyword_set(ext_file) THEN file=ext_file
  169 +ext=read_xcat(file,/silent)
  170 +
  171 +;=== HERE WE SET THE DATA SO THAT THE COMPUTE_ FUNCTIONS BELOW HAVE
  172 +;=== THE NECESSARY FIELDS AND INFORMATION
  173 +dustem_set_data,sed,sed,ext,ext
  174 +
  175 +;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS FOR THE FIT
  176 +;== AND ACTIVATE ANY PLUGINS
  177 +dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims,polarization=use_polarization
  178 +
  179 +if keyword_set(wait) then begin
  180 + message,'Finished initializing DustEMWrap, including plugins and fixed parameters',/info
  181 + wait,wait
  182 +end
  183 +
  184 +
  185 +;== INFORMATION TO RUN THE FIT
  186 +tol=1.e-10
  187 +
  188 +;=== INFORMATION TO MAKE THE PLOTS
  189 +;=== _x/_extinction means extinction plots, _m/_emissions means emission plots
  190 +xr_x=[0.01,30]
  191 +yr_x=[1.00E-10,10]
  192 +tit='Dust Optical Depth'
  193 +;ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
  194 +;xtit=textoidl('\lambda (\mum)')
  195 +
  196 +;=== RUN THE FIT
  197 +t1=systime(0,/sec)
  198 +res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $
  199 + ,/xlog,/ylog,xr_x=xr_x,yr_x=yr_x,xtit=xtit,ytit=ytit,title=tit $
  200 + ,legend_xpos=legend_xpos,legend_ypos=legend_ypos $
  201 + ,errors=errors,chi2=chi2,rchi2=rchi2,show_plot=show_plot)
  202 +t2=systime(0,/sec)
  203 +
  204 +if keyword_set(wait) then begin
  205 + message,'Finished running DustEMWrap, using Niters: '+strtrim(string(use_Nitermax),2),/info
  206 + message,'Time taken [sec]: '+sigfig(t2-t1,2,/sci),/info
  207 + wait,wait
  208 +end
  209 +
  210 +IF keyword_set(fits_save_and_restore) THEN BEGIN
  211 + message,'Writing out results structure: '+fits_save_and_restore,/info
  212 + dustem_write_fits_table,filename=fits_save_and_restore,help=help
  213 +;=== At this point, you could erase all dustem system variables, or exit idl... all the
  214 +;=== information needed to recover the results and remake the plots has been saved in the FITS table
  215 + ;stop
  216 + dustem_read_fits_table,filename=fits_save_and_restore,dustem_st=dustem_st
  217 + ;==== plot result taken from the saved fits table
  218 + res=*(*!dustem_fit).CURRENT_PARAM_VALUES
  219 + IF !dustem_noobj THEN BEGIN
  220 + dustemwrap_plot_noobj,res,dustem_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
  221 + ENDIF ELSE BEGIN
  222 + dustemwrap_plot,res,dustem_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
  223 + ENDELSE
  224 + IF keyword_set(wait) THEN BEGIN
  225 + message,'Saved the results as FITS in the file: '+fits_save_and_restore+', and made a plot using the data in this file',/info
  226 + wait,wait
  227 + ENDIF
  228 +ENDIF
  229 +
  230 +message,'Finished dustem_fit_ext_example',/info
  231 +
  232 +the_end:
  233 +
  234 +END
... ...
src/idl/dustem_fit_spectro_example.pro 0 → 100644
  Diff comments   View file @9d93919
... ... @@ -0,0 +1,241 @@
  1 +PRO dustem_fit_spectro_example,model=model $
  2 + ,sed_file=sed_file $
  3 + ,Nitermax=Nitermax $
  4 + ,fits_save_and_restore=fits_save_and_restore $
  5 + ,ucfactor=ucfactor $
  6 + ,wait=wait $
  7 + ,noobj=noobj $
  8 + ,help=help $
  9 + ,verbose=verbose
  10 +
  11 +;+
  12 +; NAME:
  13 +; dustem_fit_spectro_example
  14 +;
  15 +; PURPOSE: This routine is an example of how to fit an observational
  16 +; spectrometer data in emission (StokesI) with DustEM and
  17 +; DustEMWrap.
  18 +;
  19 +; See the DustEMWrap User Guide for more information.
  20 +;
  21 +; CATEGORY:
  22 +; DustEMWrap, Distributed, High-Level, User Example
  23 +;
  24 +; CALLING SEQUENCE:
  25 +; dustem_fit_spectro_example[,model=][sed_file=][ext_file=][,Nitermax=][,fits_save_and_restore=][,/help,/wait,/verbose,/noobj]
  26 +;
  27 +; INPUTS:
  28 +; None
  29 +;
  30 +; OPTIONAL INPUT PARAMETERS:
  31 +; None
  32 +;
  33 +; OUTPUTS:
  34 +; None
  35 +;
  36 +; OPTIONAL OUTPUT PARAMETERS:
  37 +; Plots, results structure in binary FITS table format
  38 +;
  39 +; ACCEPTED KEY-WORDS:
  40 +; model = specifies the interstellar dust mixture used by DustEM
  41 +; 'MC10' model from Compiegne et al 2010
  42 +; 'DBP90' model from Desert et al 1990
  43 +; 'DL01' model from Draine & Li 2001
  44 +; 'WD01_RV5P5B' model from Weingartner & Draine 2002 with Rv=5.5
  45 +; 'DL07' model from Draine & Li 2007
  46 +; 'J13' model from Jones et al 2013, as updated in
  47 +; Koehler et al 2014
  48 +; 'G17_ModelA' model A from Guillet et al (2018). Includes
  49 +; polarisation. See Tables 2 and 3 of that paper for details.
  50 +; 'G17_ModelB' model B from Guillet et al (2018)
  51 +; 'G17_ModelC' model C from Guillet et al (2018)
  52 +; 'G17_ModelD' model A from Guillet et al (2018)
  53 +; sed_file = string naming the path to text file in .xcat format that
  54 +; describes the observational SED. If not set, the file
  55 +; 'Data/EXAMPLE_OBSDATA/example_SED_2.xcat' is used.
  56 +; Nitermax = maximum number of fit iterations. Default is 5.
  57 +; ucfactor = Multiplicative factor for uncertainties on
  58 +; spectrometer data. Default is 1.
  59 +; fits_save_and_restore = if set, save the fit results in a binary
  60 +; FITS file. The code then restore this file and plots
  61 +; the results using the saved results information.
  62 +; help = if set, print this help
  63 +; wait = if set, wait this many seconds between each step of
  64 +; the code (for illustration purposes)
  65 +; verbose = if set, subroutines will run in verbose mode
  66 +; noobj = if set, runs with no object graphics
  67 +;
  68 +; COMMON BLOCKS:
  69 +; None
  70 +;
  71 +; SIDE EFFECTS:
  72 +; None
  73 +;
  74 +; RESTRICTIONS:
  75 +; The DustEM fortran code must be installed
  76 +; The DustEMWrap IDL code must be installed
  77 +;
  78 +; PROCEDURES AND SUBROUTINES USED:
  79 +;
  80 +; EXAMPLES
  81 +; dustem_fit_spectro_example
  82 +; dustem_fit_spectro_example,Nitermax=10,fits_save_and_restore='/tmp/mysavefile.fits'
  83 +; dustem_fit_spectro_example,model='DBP90'
  84 +;
  85 +; MODIFICATION HISTORY:
  86 +; Written by AH Oct 2022
  87 +; Evolution details on the DustEMWrap gitlab.
  88 +; See http://dustemwrap.irap.omp.eu/ for FAQ and help.
  89 +;-
  90 +
  91 +IF keyword_set(help) THEN BEGIN
  92 + doc_library,'dustem_fit_spectro_example'
  93 + goto,the_end
  94 +END
  95 +
  96 +IF keyword_set(model) THEN BEGIN
  97 + use_model=strupcase(model)
  98 +ENDIF ELSE BEGIN
  99 + use_model='WD01_RV5P5B'
  100 +ENDELSE
  101 +
  102 +known_mdls=['MC10','DBP90','DL01','WD01_RV5P5B','DL07','J13','G17_MODELA','G17_MODELB','G17_MODELC','G17_MODELD']
  103 +test_model = where(known_mdls eq use_model,ct)
  104 +if ct eq 0 then begin
  105 + message,'ISM dust model '+use_model+' unknown',/continue
  106 + message,'Known models are MC10,DBP90,DL01,WD01_RV5P5B,DL07,J13,G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD',/continue
  107 + stop
  108 +end
  109 +
  110 +use_polarization=0. ;default is no polarization in models
  111 +use_window=2 ; default graphics window number to use for plotting the results
  112 +use_verbose=0
  113 +if keyword_set(verbose) then use_verbose=1
  114 +use_ucfactor=1.
  115 +if keyword_set(ucfactor) then use_ucfactor=ucfactor
  116 +use_Nitermax=5 ; maximum number of iterations for the fit
  117 +IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax
  118 +
  119 +dustem_define_la_common
  120 +
  121 +;=== Set the (model-dependent) parameters that you want to fit
  122 +;=== Refer to the DustEM and DustEMWrap userguides for an explanation
  123 +; of the different grain types
  124 +
  125 +;=== Set the (model-dependent) parameters that you want to fit (pd),
  126 +;=== their initial values (iv)
  127 +;=== and whether they are bounded (ulimed,llimed,llims,ulims).
  128 +;=== Fixed parameters (fpd) and their values (fiv) are also set here.
  129 +;=== Refer to the DustEM and DustEMWrap User guides for an explanation
  130 +;=== of the physical meaning of dust model and plug-in parameters, and
  131 +;=== how to specify them.
  132 +
  133 +;=== Example is provided for WD01_RV5P5B
  134 +;=== If you choose a different dust model, you will need to adjust the
  135 +;=== parameter information accordingly.
  136 +
  137 +pd = [ '(*!dustem_params).G0', $ ; multiplicative factor of the radiation field
  138 + '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
  139 + '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
  140 + '(*!dustem_params).grains(2).mdust_o_mh',$ ;Gra
  141 + '(*!dustem_params).grains(3).mdust_o_mh',$ ;Gra
  142 + '(*!dustem_params).grains(4).mdust_o_mh'] ;aSil
  143 +
  144 +iv = [2., 5.4e-4, 5.4e-4, 1.8e-4, 1.65e-3, 7.e-3] ; starting guesses
  145 +
  146 +Npar=n_elements(pd)
  147 +ulimed=replicate(0,Npar)
  148 +llimed=replicate(1,Npar)
  149 +llims=replicate(1.e-15,Npar)
  150 +
  151 +;=== Fixed parameters
  152 +fpd=[]
  153 +;=== Initial parameter values for fixed parameters
  154 +fiv=[]
  155 +
  156 +
  157 +
  158 +if keyword_set(wait) then begin
  159 + message,'Finished setting dust model and plug-in parameters: '+use_model,/info
  160 + wait,wait
  161 +end
  162 +
  163 +;== INITIALISE DUSTEM
  164 +dustem_init,model=use_model,polarization=use_polarization
  165 +!dustem_nocatch=1
  166 +!dustem_verbose=use_verbose
  167 +IF keyword_set(noobj) THEN !dustem_noobj=1
  168 +!EXCEPT=2 ; for debugging
  169 +
  170 +
  171 +;=== READ EXAMPLE EMISSION DATA TO FIT
  172 +dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'
  173 +file=dir+'example_SED_2.xcat'
  174 +if keyword_set(sed_file) then file=sed_file
  175 +sed=read_xcat(file,/silent)
  176 +
  177 +;=== TWEAK THE UNCERTAINTIES ON THE SPECTROMETER DATA TO ADJUST THEIR
  178 +;=== RELATIVE WEIGHT IN THE FIT
  179 +fidx=where(sed.filter eq 'SPECTRUM',fct)
  180 +if fct gt 0 then sed[fidx].sigmaii = use_ucfactor*sed[fidx].sigmaii
  181 +
  182 +;== SET THE OBSERVATIONAL STRUCTURE
  183 +dustem_set_data, sed,sed
  184 +
  185 +;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS FOR THE FIT
  186 +;== AND ACTIVATE ANY PLUGINS
  187 +dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims,polarization=use_polarization
  188 +
  189 +if keyword_set(wait) then begin
  190 + message,'Finished initializing DustEMWrap, including plugins and fixed parameters',/info
  191 + wait,wait
  192 +end
  193 +
  194 +
  195 +;== INFORMATION TO RUN THE FIT
  196 +tol=1.e-10 ;fit tolerance
  197 +
  198 +;=== INFORMATION TO MAKE THE PLOT
  199 +;=== _x/_extinction means extinction plots, _m/_emissions means emission plots
  200 +xr_m = [1.,5e5]
  201 +yr_m = [5e-8,1.00e6]
  202 +tit_m='Spectral Energy Distribution'
  203 +
  204 +;=== RUN THE FIT
  205 +t1=systime(0,/sec)
  206 +res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $
  207 + ,/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 $
  208 + ,legend_xpos=legend_xpos,legend_ypos=legend_ypos $
  209 + ,errors=errors,chi2=chi2,rchi2=rchi2,show_plot=show_plot)
  210 +t2=systime(0,/sec)
  211 +
  212 +if keyword_set(wait) then begin
  213 + message,'Finished running DustEMWrap, using Niters: '+strtrim(string(use_Nitermax),2),/info
  214 + message,'Time taken [sec]: '+sigfig(t2-t1,2,/sci),/info
  215 + wait,wait
  216 +end
  217 +
  218 +
  219 +
  220 +IF keyword_set(fits_save_and_restore) THEN BEGIN
  221 + message,'Writing out results structure: '+fits_save_and_restore,/info
  222 + dustem_write_fits_table,filename=fits_save_and_restore,help=help
  223 + dustem_read_fits_table,filename=fits_save_and_restore,dustem_st=dustem_spectra_st
  224 + ;==== plot result taken from the saved fits table
  225 + res=*(*!dustem_fit).CURRENT_PARAM_VALUES
  226 + IF !dustem_noobj THEN BEGIN
  227 + dustemwrap_plot_noobj,res,dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,tit_m=tit_m+' (From Saved FITS file)'
  228 + ENDIF ELSE BEGIN
  229 + dustemwrap_plot,res,dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,tit_m=tit_m+' (From Saved FITS file)'
  230 + ENDELSE
  231 + IF keyword_set(wait) THEN BEGIN
  232 + message,'Saved the results as FITS in the file: '+fits_save_and_restore+', and made a plot using the data in this file',/info
  233 + wait,wait
  234 + ENDIF
  235 +ENDIF
  236 +
  237 +message,'Finished dustem_fit_spectro_example',/info
  238 +
  239 +the_end:
  240 +
  241 +END
... ...