Examples

Example 0: Running DustEM non-iteratively using DustEMWrap

In this example, we demonstrate how to run the DustEM fortran code non-iteratively from within IDL using DustEMWrap, and several ways of visualising the output of the fortran code. The example routine is called dustem_run_example.pro. 

To run the example, open an IDL session and type
IDLPrompt> dustem_run_example,'J13'
This will run DustEM using the J13 insterstellar dust model and produce the default set of plots that illustrates the model properties. In the plots, the different grain populations in the model are indicated with different line styles (if relevant). The overall model is indicated with a solid blue line, and typical observations of the Galactic diffuse ISM are indicated for comparison with red points. 

Looking at the code in dustem_run_example.pro, you will see that it is possible to modify the input properties of the dust model after reading the default model input files into an IDL structure. For example, we could increase the intensity of the dust-heating radiation field from G0 to 2*G0, and reduce the abundance of one of the amorphous silicate grain types in the J13 model before running DustEM to compute the emission and extinction. 

; read the default input files for the model into an IDL structure
st_model=dustem_read_all(dir_in)

st_model.G0=2.0 ; change the default value of the ISRF intensity used by DustEM
st_model.grains[3].mdust_O_mh=20.e-4 ; change the abundance of grain type 3 (aOlM5)

; write out the modified model properties to the input files used by the fortran code
dustem_write_all,st_model,!dustem_dat

; run the fortran using the new input parameters, and store the results into an IDL structure
 st_results=dustem_run() 

There are several ways to visualise the input/output of DustEM. Plots can be constructed directly from the variables in the st_model and st_results structures. In the routine dustem_run_example.pro, an example using st_model to plot the dust-heating radiation field is shown. 

The dustem_show_fortran procedure called by dustem_run_example currently allows the user to plot 

The default behaviour is show="all", which will generate all the plots that are appropriate for the model. A subset of plots can be generated by specifying the desired plots as an array of strings in the show keyword, e.g. show=["emis", "sdist", "alb", "align"]. The total intensity SED is always generated. We note that the dustem_show_fortran procedure will be deprecated in a future release of DustEMWrap, and will be replaced by dedicated functions that produce each of the different plot types.


Example 1: Fitting a broadband SED with DustEMWrap: Total Intensity

In this example, we demonstrate how DustEMWrap can be used to fit an SED with total intensity (Stokes I only) broadband measurements. The routine is called dustem_fit_intensity_example.pro. By default, this routine fits the data in the input file Data/EXAMPLE_OBSDATA/example_SED_1.xcat, which contains Stokes I and Stokes I uncertainty measurements for a range of measurements at infrared and millimetre wavelengths. 

As noted in the header of example_SED_1.xcat, the emission values are assumed to correspond to a fiducial column density N(H) = 1.e20 cm-2. Undefined or missing values are set to -32768. The data columns that are essential for this example are FILTER (the filter name), STOKESI (the Stokes I flux density measurement), and SIGMAII (the variance in Stokes I). Other columns (e.g. STOKESQ) are not used, although some columns may be completed for the user's convenience (e.g. WAVE).

To run the example, open an IDL session and type
IDLPrompt> dustem_fit_intensity_example

This will run Dustemwrap for Nitermax=5 iterations using the default DPB90 interstellar dust model, attempting to optimise the relative abundance of the three grain types and the intensity of the dust-heating radiation field. Looking at the code in dustem_fit_intensity_example.pro, the free parameters are defined in the parameter description vector (pd), and their initial values are set (iv), and lower-bounded to 0 (llimed, llims).  DustEMWrap initialises the dust grain information by reading the GRAIN_DBP90.DAT file, which is part of the DustEM fortran distribution. An important general rule of DustEMWrap is that dust model parameters that are not specified in either the pd or fpd vectors are held fixed at their default values, which are defined in the fortran GRAIN_XXX.dat files. 

While the DustEMWrap code runs, the results of each iteration of the DustEM fortran are printed to the terminal. A plot window is created, which shows the observational data from the input SED file (blue points), the expectation of the model for each band (including color corrections, red squares), the overall dust model spectrum (black line) and the spectra for each of the dust grain populations in the model (coloured lines). A lower sub-panel indicates the residuals of the fit. Two small additional windows display the current best-fitting values of the free parameters of the dust model, and of any plug-ins. (In this example, nothing is displayed in the plug-ins window). 

In this example, the fitting procedure stops because the maximum number of iterations is reached. The final solution is illustrated in the plot window, and the final parameter values are listed in the two parameter windows. 

Once you have run dustem_fit_intensity_example.pro with the default parameters to completion, you should explore how to modify the fixed and free parameters included in your fit, how to invoke a different interstellar dust model, include plugins, change the number of iterations, fit different observational data etc. 

Here are some examples of how to make minor adjustments to your dustem_fit_intensity_example run from the IDL command line:

To modify the parameters included in the fit and/or modify their initial values, you need to edit the dustem_fit_intensity_example.pro routine. The measurements in the example_SED_3.xcat file include data points at long wavelengths that seem incompatible with thermal dust emission alone. Let's run the DPB90 model again on the example_SED_3.xcat file again, but this time include a plug-in that adds synchrotron emission. For this, you need to edit dustem_fit_intensity_example.pro. 

The parameter description vector should be modified as follows:
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

And the initial values vector likewise requires input for the synchrotron parameters, e.g.
; initial values vector
; [G0,PAH0,VSG,BG,alpha_CR,Amp_syn]
iv =   [1.6, 2.2e-4, 5.7e-4, 3.4e-3, 2.7,0.01]

Having made this modifications, try running DustEMWrap again for three iterations on the input SED example_SED_3.xcat:
IDLPrompt> dustem_fit_intensity_example,sed_file='path_to_dustemwrap/Data/EXAMPLE_OBSDATA/example_SED_3.xcat',niter=3
This time, you should see the plugin emission component as a dashed coloured line in the plotting window. The additional windows now include the current best-fitting values of the free parameters of the synchrotron plugin. Once again, the fitting procedure stops because the maximum number of iterations is reached. 

Finally, let's run DustEMWrap including a 750K blackbody component and a radiation field with amplitude G0=0.5. In this case, we modify the  dustem_fit_intensity_example.pro routine to declare these as fixed parameters and set their values using the fpd and fiv vectors: 

fpd = [ $
     '(*!dustem_params).G0', $                           ;G0
    'dustem_plugin_continuum_1', $                    ;Temperature of the blackbody
     'dustem_plugin_continuum_2']                      ;Intensity at peak of the blackbody

fiv = [0.5, 750, 0.01]

The parameter description and initial values vectors should also be modified to remove G0  as a free parameter, i.e.
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',$         ;BG mass fraction
     'dustem_plugin_synchrotron_1', $                    ;Spectral index of CREs
     'dustem_plugin_synchrotron_2']                      ;Synchrotron amplitude at 10 mm

iv =   [2.2e-4, 5.7e-4, 3.4e-3, 2.7,0.01]

Having made this modifications, try running DustEMWrap again:
IDLPrompt> dustem_fit_intensity_example,sed_file='path_to_dustemwrap/Data/EXAMPLE_OBSDATA/example_SED_3.xcat',niters=10
This time, you should see both the synchrotron and blackbody plugin emission components as dashed coloured lines in the plotting window. The additional windows include the current best-fitting values of the free parameters of the dust model and synchrotron plugin, but also the values of the fixed parameters. 

DustEMWrap runs using other interstellar dust models and including different combinations of free/fixed parameters are included in the dustem_fit_intensity_example.pro routine. You can try them by uncommenting the relevant lines.

Example 2: Fitting a broadband SED with a modified blackbody model

In this example, we demonstrate how DustEMWrap can be used to fit an SED (Stokes I only) using a modified blackbody to model the emission in the far-infrared, rather than a physical interstellar dust model. Note that an interstellar dust model must still be invoked for DustEMWrap to work. In this example, we use the MC10 model and the input SED Data/EXAMPLE_OBSDATA/example_SED_4.xcat. We use MC10 to describe the observed emission at short wavelengths due to PAHs, but replace the three big grain components by a plugin that models a modified blackbody.

To run the example, open an IDL session and type
IDLPrompt> dustem_fit_intensity_mbb_example
This will run DustEMWrap for Nitermax=5 iterations, attempting to optimise the relative abundance of the PAH0 and PAH1 grain types in the MC10 model, the intensity of the dust-heating radiation field, and the three free parameters of the modified blackbody (temperature, peak intensity and spectral index).

Looking at the code of dustem_fit_intensity_mbb_example.pro, you will see that parameter description vector pd includes the PAH abundances as free parameters, and the intensity of the interstellar radiation field. The big grain abundances specified in the fixed parameter description vector fpd and their values set to 1.e-12 (i.e. a very small number). 

pd = ['(*!dustem_params).G0', $                           ;G0
      '(*!dustem_params).grains(0).mdust_o_mh',$          ;PAH0 mass fraction
      '(*!dustem_params).grains(1).mdust_o_mh',$          ;PAH1 mass fraction
      'dustem_plugin_mbbdy_1', $                          ;Amplitude of mBB
       'dustem_plugin_mbbdy_2', $                          ;T of mBB
       'dustem_plugin_mbbdy_3']                            ;Emissivity index of mBB

iv =   [1.8, 7.8e-4, 7.8e-4, 4.e-5, 22., 1.8]

fpd = [(*!dustem_params).grains(2).mdust_o_mh',$          ;amCBEx
      '(*!dustem_params).grains(3).mdust_o_mh',$          ;amCBEx
      '(*!dustem_params).grains(4).mdust_o_mh']           ;aSilx

fiv =   [1.e-12,1.e-12,1.e-12]

In this example, the fitting procedure stops because the maximum number of iterations is reached. The final solution is illustrated in the plot window, and the final parameter values are listed in the two parameter windows.  

Example 3: Fitting an SED with combined broadband and spectrometer data

In this example, we demonstrate how DustEMWrap can be used to fit an SED (Stokes I only) that has a combination of spectrometer and broadband measurements. By default, this routine fits the data in the input file Data/EXAMPLE_OBSDATA/example_SED_2.xcat, which contains measurements from the FIRAS and ISO spectrometers, as well as DIRBE and WMAP broadband measurements.

As noted in the header of example_SED_2.xcat, the emission values are assumed to correspond to a fiducial column density N(H) = 1.e20 cm-2. Undefined or missing values are set to -32768. Unlike the other example .xcat files that we have seen, the filter name for the spectrometer data is always SPECTRUM, since color corrections are negligible for very narrow bandpasses (like spectrometer measurements).

To run the example, open an IDL session and type
IDLPrompt> dustem_fit_spectro_example
This will run DustEMWrap for Nitermax=5 iterations, attempting to optimise the relative abundance of the different grain types in the WD01_RV5P5B model and the intensity of the dust-heating radiation field.

The art of fitting a combination of broadband and spectrometer SED is to define appropriate uncertainties in the input SED such that the relative weight given to the different data types is appropriate. Modifying the uncertainties can be done by editing the example_SED_2.xcat file directly (here the values in the SIGMAII column) or from within the dustem_fit_spectro_example.pro routine after the SED file has been read into the sed IDL variable i.e.

;=== READ EXAMPLE EMISSION DATA TO FIT
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)

;=== TWEAK THE UNCERTAINTIES ON THE SPECTROMETER DATA TO ADJUST THE RELATIVE WEIGHT IN THE FIT
fidx=where(sed.filter eq 'SPECTRUM',fct)
if fct gt 0 then sed[fidx].sigmaii = use_ucfactor*sed[fidx].sigmaii

By default use_ucfactor = 1.0, but you can give less weight to the spectrometer data by increasing the associated uncertainties (i.e. by setting use_ucfactor to a value greater than 1). In our example, setting use_ucfactor = 2.0 leads to a fit with a slightly higher abundance of PAH0 and aSil grain types, and a slightly lower abundance of Gra grains with a power law size distribution (compared to use_ucfactor = 1.0).

Example 4: Fitting extinction curve data with DustEMWrap

In this example, we demonstrate how DustEMWrap can be used to fit an extinction curve. By default, this routine fits the data in the input file Data/EXAMPLE_OBSDATA/Mathis90Fitz99_DISM_NH20.xcat, which corresponds to a representative extinction curve for the diffuse ISM, as described by Mathis 1990 and Fitzpatrick & Massa 1999. 

As noted in the header of Mathis90Fitz99_DISM_NH20.xcat, the extincation values are assumed to correspond to a fiducial column density N(H) = 1.e20 cm-2. Undefined or missing values are set to -32768. Unlike the other example .xcat files that we have seen, the instrument name for extinction data is always EXTINCTION, and the filter name is  SPECTRUM, since the curve is continuous.

To run the example, open an IDL session and type
IDLPrompt> dustem_fit_ext_example
This will run DustEMWrap for Nitermax=5 iterations, attempting to optimise the relative abundance of the different grain types in the J13 model.

In this example, the fitting procedure stops because the maximum number of iterations is reached. The final solution is illustrated in the plot window, and the final parameter values are listed in the dust parameter window.  

Example 5: Generating and fitting a Stokes IQU SED 

In this example, we demonstrate how DustEMWrap can be used to generate an SED using an interstellar dust model, and to fit the resulting SED.  Here we specifically work with a Stokes IQU SED that includes polarization information, but a similar approach can be used to generate a Stokes I (total intensity only) SED.

In general, we note that the SED fitting procedure implemented in DustEMWrap works exclusively with Stokes parameters, consistent with observational data. Interstellar dust models, however, typically express their predictions in terms of a polarisation fraction. When generating polarized SEDs, the user must specify a dust polarization angle, which depends on the hypothetical magnetic field geometry and is not intrinsic to the interstellar dust model. Comparing polarised interstellar dust models with observational data therefore requires converting between the Stokes IQU and (p=P/I,psi) conventions for describing the polarisation of the emission.  

In this example, we use the G17_MODELD model and the synchrotron plugin to generate an SED. To run the example, open an IDL session and type
IDLPrompt> dustem_make_polarization_sed_example
This will generate an SED that corresponds to the G17_MODELD model with the parameters set to the values in the pd_true vector. A polarized synchrotron component is added to the SED using a plugin, and the dust polarization angle is set. The synchrotron and dust polarization angles are independent. Looking at the code, you will see

pd = [  '(*!dustem_params).G0', $                   ;G0
 '(*!dustem_params).grains(0).mdust_o_mh',$  ;PAH0_MC10 mass fraction
        '(*!dustem_params).grains(1).mdust_o_mh',$  ;amCBE_0.3333x mass fraction
        '(*!dustem_params).grains(2).mdust_o_mh', $ ;aSil2001BE6pctG_0.4x mass fraction
        'dustem_plugin_modify_dust_pol_2',  $     ;Dust polarization angle
        'dustem_plugin_synchrotron_1', $   ;Synchrotron spectral index
        'dustem_plugin_synchrotron_2', $   ;Synchrotron amplitude at 10 mm
        'dustem_plugin_synchrotron_3', $    ;Synchrotron polarization fraction
        'dustem_plugin_synchrotron_4'  $    ;Synchrotron polarization angle
     ]                   

pd_true=[1., $ ;G0
        7.8000E-04, $ ;PAH0_MC10 mass fraction
         7.8000E-04, $ ;amCBE_0.3333x mass fraction
        7.8000E-04, $ ;aSil2001BE6pctG_0.4x mass fraction
        10., $ ;Dust polarization angle
        3., $ ;Synchrotron spectral index
        1.e-2, $ ;Synchrotron amplitude at 10 mm
         0.3, $ ;Synchrotron polarization fraction
        45.] ;Synchrotron polarization angle

To construct an observational SED, DustEMWrap further requires the set of instrument filters (filters, which can be passed as a keyword), and the name of the ouput file for the SED (outfile). Based on the input dust model and plugin parameters, and the chosen list of filters, the IQU SEDs are computed using dustem_compute_sed (for Stokes I) and dustem_compute_stokes (for Stokes QU). The SEDs are stored in the sed IDL structure. In our example, uncertainties (i.e. variances and co-variances in IQU) are manually specified, e.g. sed.sigmaQQ = abs(sed.StokesQ*0.000001). 

Having computed IQU, the corresponding polarised intensity, polarisation fraction and polarisation angle columns of the sed structure are completed using the dustem_fill_sed_dependent_columns function (for user convenience only, since DustEMWrap only uses Stokes IQU during the fitting).

We illustrate how to fit a Stokes IQU SED using the example routine dustem_fit_polarization_example.pro. The code is very similar to Example 1 (dustem_fit_intensity_example.pro), except the graphical output display of DustEMWrap is different since additional panels displaying Stokes QU, P, p=P/I and psi as a function of wavelength are shown in addition to Stokes I.

To run the example, open an IDL session and type
IDLPrompt> dustem_fit_polarization_example

By default, this will fit the SED in the Data/EXAMPLE_OBSDATA/example_SED_3.xcat file with the G17_MODELD dust model. (This input SED file was generated using the same input parameters for the dust model and synchrotron plugin as described above.) In the example routine, the dust grain abundances, the dust polarisation angle and some of the synchrotron characteristics are set as free parameters, while the intensity of the radiation field and the spectral index of the cosmic rays are held fixed. Our starting guess for the free parameters is close to the values that we used to generate the SED, i.e.
iv = [7.2000E-04, $ ;PAH0_MC10 mass fraction
         7.2000E-04, $ ;amCBE_0.3333x mass fraction
        8.5000E-04, $ ;aSil2001BE6pctG_0.4x mass fraction
        -10., $ ;Dust polarization angle
        1.5e-2, $ ;Synchrotron amplitude at 10 mm
         0.25, $ ;Synchrotron polarization fraction
        35.] ;Synchrotron polarization angle

We constrain the possible values of the synchrotron emission, and place a lower-bound of zero on the dust grain abundances, i.e.
ulimed=[0,0,0,0,1,1,1]
llimed=[1,1,1,0,1,1,1]
ulims=[0,0,0,0,0.02,0.5,60.]
llims=[0,0,0,0,0.005,0.05,30.]

Example 6: Generating and fitting a Stokes IQU extinction curve

In this example, we demonstrate how DustEMWrap can be used to generate an extinction curve using an interstellar dust model, and to fit the resulting curve. This example is very similar to Example 5, but here we work with extinction measurements.  Our extinction curve in this example includes polarization information, but a similar approach can be used to generate Stokes I extinction data. 

In this example, we use the G17_MODELD model to generate an extinction curve. To run the example, open an IDL session and type
IDLPrompt> dustem_make_polarization_ext_example

This will generate an extinction curve that corresponds to the G17_MODELD model with the parameters set to the values in the pd_true vector (including a dust polarization angle for which we use a plugin). Looking at the code, you will see

pd = [
'(*!dustem_params).grains(0).mdust_o_mh',$  ;PAH0_MC10 mass fraction
        '(*!dustem_params).grains(1).mdust_o_mh',$  ;amCBE_0.3333x mass fraction
        '(*!dustem_params).grains(2).mdust_o_mh', $ ;aSil2001BE6pctG_0.4x mass fraction
        'dustem_plugin_modify_dust_polx_2'  $     ;Dust polarization angle in extinction
    ]                   

true_vals=[
        7.8000E-04, $ ;PAH0_MC10 mass fraction
         7.8000E-04, $ ;amCBE_0.3333x mass fraction
        7.8000E-04, $ ;aSil2001BE6pctG_0.4x mass fraction
        10. $ ; Dust polarization angle
  ]

To construct an observational extinction curve, DustEMWrap further requires the wavelength range (wave) over which the extinction curve is defined, the number of measurements in the extinction curve (Next), and the name of the ouput file for the extinction curve  (outfile).  By default, the extinction curve is defined on a logarithmically spaced wavelength vector with 100 elements between 0.01 and 50 microns.

Based on the input dust model and plugin parameters, the IQU extinction curves are computed using dustem_compute_ext (for Stokes I) and dustem_compute_stokext (for Stokes QU). The extinction curves are stored in the ext IDL structure. In our example, uncertainties (i.e. variances and co-variances in IQU) are manually specified, e.g. ext.sigmaextqq = abs(ext.EXT_Q*0.000001). 

Having computed IQU, the corresponding polarisation fraction and polarisation angle columns of the ext structure may be completed using the dustem_fill_ext_dependent_columns function (for user convenience only, since DustEMWrap only uses Stokes IQU during the fitting).

We illustrate how to fit a Stokes IQU extinction curve using the example routine dustem_fit_ext_pol_example.pro. The code is very similar to Example 4 (dustem_fit_ext_example.pro), except the graphical output display of DustEMWrap is different since additional panels displaying polarisation information are included. 

To run the example, open an IDL session and type
IDLPrompt> dustem_fit_ext_pol_example

By default, this will fit an extinction curve using the G17_MODELD dust model. In the example routine, we are interested in the dust grain abundances and the dust polarisation angle parameters (pd). 

pd = [ $
     '(*!dustem_params).grains(0).mdust_o_mh',$         ;PAH0 mass fraction
    '(*!dustem_params).grains(1).mdust_o_mh',$         ;amCBE mass fraction
    '(*!dustem_params).grains(2).mdust_o_mh',$          ;aSil
    'dustem_plugin_modify_dust_polx_2' ]

If an extinction file is not specified (via the ext_file= keyword), then the example generates the data using the parameter values

true_vals = [7.1e-4, 1.3e-3, 5.6e-3, 35.] ; these are the values that will be used to generate data if no ext_file is specified

and the dustem_compute_ext (for Stokes I) and dustem_compute_stokext (for Stokes QU) routines. 

Our starting guess for the free parameters is close to the values that we used to generate the SED, i.e.
iv = true_vals+[5.00E-4,-7.00E-4,8.00E-4,10] ; this is the vector of starting guesses for the fit

Example 7: Fitting extinction and emission data simultaneously

In this example, we demonstrate how to fit observations of emission and extinction with an interstellar dust model using DustEMWrap. This example is very similar to the preceding ones, but here we work with Stokes IQU extinction and Stokes IQU emission input data all at the same time. The example describes observational data the includes polarization information, but a similar approach can be used if only Stokes I measurements in emission and/or extinction are present (see the routine dustem_fit_sed_ext_stokesi_example.pro for such a Stokes I only example).

In this example, we use the G17_MODELD model to generate and fit the observational data . To run the example, open an IDL session and type
IDLPrompt> dustem_fit_sed_ext_pol_example

This will generate an extinction curve and an SED that corresponds to the G17_MODELD model with the parameters set to the values in the true_vals vector (including independent dust polarization angles for emission and extinction, for which we use plugins). Looking at the code, you will see

pd = [ '*!dustem_params).grains(0).mdust_o_mh',$  ;PAH0_MC10 mass fraction
        '(*!dustem_params).grains(1).mdust_o_mh',$  ;amCBE_0.3333x mass fraction
        '(*!dustem_params).grains(2).mdust_o_mh', $ ;aSil2001BE6pctG_0.4x mass fraction
        'dustem_plugin_modify_dust_pol_2'  $     ;Dust polarization angle in emission
         'dustem_plugin_modify_dust_polx_2'  $     ;Dust polarization angle in extinction   ]                   

true_vals = [5.4e-4, 6.4e-4, 4.8e-4, 48., 52] ; these are the true values that will be used to generate the data

Based on the input dust model and plugin parameters, the IQU curves are computed using dustem_compute_sed (for emission Stokes I) , dustem_compute_stokes (for emission Stokes QU), dustem_compute_ext (for extinction Stokes I) and dustem_compute_stokext (for extinction Stokes QU). The SEDs are stored in the sed IDL structure, and the extinction curves are stored in the ext IDL structure. In our example, we have set the uncertainties and covariances to a constant small value (1.e-15). 

Our starting guess for the free parameters is close to the values that we used to generate the SED, i.e.
iv = true_vals+[2.00E-4,-2.00E-4,2.00E-4, 10, -10] ; this is the vector of starting guesses for the fit

As for Example 6, this example is written so that a user can input their own SED and/or extinction curve data using text files in.xcat format. These are specified with the ext_file= and  sed_file= keywords. For example, 

IDLPrompt> dustem_fit_sed_ext_pol_example,ext_file='myextinction.xcat',sed_file='mysed.xcat'

Example 8: Using the modify_ISRF plugin

In the next two examples, we demonstrate how to use the distributed plugins that modify the dust-heating ISRF. There are two plugins of this type : dustem_plugin_modify_isrf.pro and dustem_plugin_stellar_population.pro. The first reads an ISRF from a user-provided text file, while the second constructs an ISRF internally based on a user-specified population of local main-sequence stars. 

In the current release of DustEMWrap, the amplitude of an ISRF provided via a text-file can be included as a fit parameter when using dustem_plugin_modify_isrf.pro, i.e. the spectral shape of the user-provided ISRF is not free to vary. The other plugin, dustem_plugin_stellar_population.pro is more powerful since both the shape and the amplitude of the ISRF are free to vary according to the user-specified stellar parameters (effective temperature, stellar radius, distance and number) that are included as fit parameters.

Let's start by looking at dustem_plugin_modify_isrf.pro. To run the example, open an IDL session and type
IDLPrompt> dustem_myisrf_example

In this example, we fit the SED described by the input file Data/EXAMPLE_OBSDATA/example_SED_1.xcat using the DBP90 dust model and a synchrotron component (by invoking dustem_plugin_synchrotron). We provide our own ISRF by including dustem_plugin_modify_isrf in the pd vector, and setting the system variable !dustem_isrf_file=ptr_new(isrf_file), where isrf_file is the name of a text file describing the ISRF that we have created using the routine dustem_create_rfield.pro. It is not necessary to use dustem_create_rfield.pro. to generate an ISRF, but the input ISRF text file must conform to the same format as the output of this routine, i.e. column 1 listing the wavelength (in microns) covering 0.01 to 1e5 microns in 200 steps and column 2 listing the corresponding ISRF intensity in erg/s/cm2/Hz.   

The free parameters in the fit for this example are the abundances of the PAH component in the DBP90 dust model, the amplitude of the synchrotron at 1cm, and the amplitude of the ISRF that we have provided as input. The fit converges to a solution where the amplitude of the ISRF is ~0.001. This makes sense since the ISRF that we generated is ~1000 times stronger than the standard Mathis field, while the observational data in example_SED_1.xcat corresponds to G0~1. 

When using the plugins that modify the ISRF, it is important to remember the general rule of DustEMWrap that parameters that are not specified in either the pd or fpd vectors are held fixed at their default values. In our example, we have not specified (*!dustem_params).G0 (the amplitude of the standard Mathis radiation field) in either pd or fpd, so it is held fixed at its default value G0=1 during the DustEM run. To suppress the Mathis component of the ISRF altogether, we should specify fpd=['(*!dustem_params).G0'] and fix its value to a very small number, e.g. fiv=[1.e-10].

Example 9: Using the stellar population ISRF plugin

In this example, we demonstrate how to use dustem_plugin_stellar_population.pro.  Here we use the DBP90 model to generate and fit the observational data. We do not try to fit any parameters of the dust model, however, but instead try to determine ISRF parameters that are most consistent with the data.

To run the example, open an IDL session and type
IDLPrompt> dustem_stellarpopisrf_example

This will generate an SED that corresponds to the DBP90 model with the default grain properties. We activate the stellar population plugin for a single O7V star, i.e.

pd = ['dustem_plugin_stellar_population_O7V3']
true_vals =   [20.]                ; true distance to the star

fpd=['dustem_plugin_stellar_population_O7V4' , $ ;number of O7V stars (FIXED)
'(*!dustem_params).G0'] ; Mathis field (FIXED TO ~ZERO)
fiv=[1.,1.e-12]

We bracket the allowed values of the distance to the star between 15 and 25pc using ulims[0]=25. and llims[0]=15.

Recall that parameters that are not specified in either the pd or fpd vectors are held fixed at their default values. For the amplitude of the standard Mathis radiation field(*!dustem_params).G0 , the default value G0=1. This is why we fix G0 to ~zero,

Based on the input DBP90 dust model and plugin parameters, the SED is computed using dustem_compute_sed. In our example, we set the Stokes I relative uncertainties to 1%, sed.SigmaII = sed.StokesI*0.01. Our starting guess for the free parameters is close to the values that we used to generate the SED, i.e.
iv = true_vals+[4.]               ; distance we use as an initial guess

The fit converges to a solution with the right distance in a few iterations although the uncertainty on the distance is quite large since our data does not provide very strong constraints. 

The dustem_plugin_modify_isrf and dustem_plugin_stellar_population plugins can be combined. The dustem_stellarpopisrf_example.pro code includes a second part where we include a Habing field as well as the radiation from a single O7V star. Looking at the code, we see

pd = [ 'dustem_plugin_stellar_population_O7V3' , $          ;distance to O7V star
       'dustem_plugin_modify_isrf_1']                        

As before, we fix the number of stars, the amplitude of the Mathis field and limit the allowed values of the distance to the star. We use dustem_create_rfield.pro. to generate a Habing ISRF within the code itself, i.e.

my_isrf_file='./myisrf_habing.dat'
myisrf=dustem_create_rfield([0],isrf=2,fname=my_isrf_file,x=mywaves)
isrf_file=my_isrf_file

The ISRF could also be provided as a text file via the isrf_file= keyword. Within a few iterations, the fit converges to a solution with a stronger Habing ISRF amplitude (~3) and more distant O-star (~23pc)  than we used to generate the model spectrum. This indicates that our data does not provide very strong constraints to separate the two ISRF contributions.

Example 10: Using the SED extractor helper tool

In this example dustem_extract_sed_example.pro, we demonstrate how to use a helper tool (dustem_sed_extractor.pro) to generate SEDs (Stokes IQU, or StokesI only) that conform to the input requirements of DustEMWrap. 

The SED extractor provides three different options for generating a SED, specified using the method keyword: 

The regions of interest are specified using pixel indices. The RoI where the SED will be extracted is a required input to dustem_sed_extractor (index_on), along with the data structures containing the input maps (maps) and the corresponding list of DustEMWrap filters (filters). The OFF region (only needed when using method=MEAN_ONOFF) is specified using the index_off keyword.

In this example, we generate some fake Stokes I and uncertainty data using the function dustem_planck_function to describe a blackbody at the wavelengths of the specified filters, and manually define a mask to extract the SED using pixel indices, i.e.

;; Define the region of interest (ON) and OFF region
mask=maps[*,*,0,0]
mask[5:35,5:35]=1 ; OFF will be a guard band of 5 pixel around the ON defined below
mask[10:30,10:30]=mask[10:30,10:30]+1 ; ON source

...

;; Loop through the maps corresponding to different filters
;; ON pixels are assigned the emission of a BB at 40K
;; OFF pixels are assigned the emission of a BB at 20K

FOR i=0L,Ndata_sets-1 DO BEGIN
 wave=dustem_filter2wav(filters[i])
 this_map0=maps[*,*,indI,i]
 this_map=this_map0
 this_map[index_on]=dustem_planck_function(temp_on,wave)
 this_map[index_off]=dustem_planck_function(temp_off,wave)
 maps[*,*,indI,i]=this_map
;; assign a variance value to each valid pixel
  maps[*,*,indII,i]=la_power(la_mul(maps[*,*,indI,i],0.3),2.)
ENDFOR

For real analysis, these lines of the dustem_extract_sed_example code would typically be replaced by reading FITS files that correspond to the requested filters and the masks. Note that the dustem_sed_extractor tool does not perform any smoothing or reprojection, i.e. it expects that all input data already at the desired resolution and on the same astrometric grid. 

To run the example, open an IDL session and type
IDLPrompt> dustem_extract_sed_example

In the example, there are two calls to the SED extractor. The first

sed=dustem_sed_extractor(maps,index_on,filters,maps_order=maps_order,/total_intensity_only)

returns a DustEMWrap SED structure using the default MEAN method for pixels inside the RoI specified by  index_on. The wavelengths of the SED correspond to the filters.  In this example, we defined the filters and index_on at the start of the code, i.e.

filters=['IRAS3','IRAS4','PACS1','PACS2','PACS3','PILOT1','SPIRE1','SPIRE2','SPIRE3','HFI3','HFI4','HFI5','HFI6']

mask[5:35,5:35]=1 ; OFF is a guard band of 5 pixel around the ON defined below
mask[10:30,10:30]=mask[10:30,10:30]+1 ; ON source
index_on=where(mask EQ 2,count_on)
index_off=where(mask EQ 1,count_off)

The second SED is extracted using the MEAN_ONOFF method for exactly the same region and filter set, and the OFF region that surrounds the ON, i.e.

sed2=dustem_sed_extractor(maps,index_on,filters,maps_order=maps_order,/total_intensity_only,index_off=index_off,method='MEAN_ONOFF')

The example first generates a simple plot. The first SED is shown with red data points, and the second SED is shown with blue data points. 

The example then initialises fits the sed using DustEM with the blackbody continuum plugin, with an initial guess of 15K for the blackbody temperature and 1e6 MJy/sr for the peak blackbody emission. Note that DustEM requires invoking a physical dust model to run (in this example, we use DBP90) but we disable the physical dust model from participating in the fit by fixing all the grain abundances to ~zero.

Image: NIRCAM/JWST, NASA, ESA, CSA, STScI; https://webbtelescope.org/news/first-images