Merge pull request #1 from eblur/rogantini_fek

Add Rogantini (2018) Fe K olivine cross-section to absorption model

README.md

@@ -8,7 +8,52 @@ To download, use git to create an ismdust folder:
 
         git clone https://github.com/eblur/ismdust.git ismdust
 
-## Setup for ISMdust
+## About the model parameters
+
+Each model applies the value `exp(-tau md)` to the flux, where tau is the path-integrated
+absorption, scattering, or extinction optical depth and `md` is the
+dust mass column in units of 10^-4 g cm^-2.
+
+To estimate the appropriate dust mass column for a given ISM column (NH),
+one must choose a dust-to-gas mass ratio. For for the Milky Way ISM,
+this value is typically ~0.01.
+
+For example, NH = 10^22 cm^-2 corresponds to NH * m_p (proton mass) * 0.01 =
+ 1.67e-4 g cm^-2. One can then set, for example, the `md_sil`
+parameter in ISMdust to 1.67.
+
+**ISMdust model parameters**
+
+    sil_md : dust mass column for Silicate in units of 1.e-4
+      default = 0.6
+
+    gra_md : dust mass column for Graphite in units of 1.e-4
+      default = 0.4
+
+    redshift : redshift of the obscuring dust (XSPEC only)
+      default = 0.0
+
+The default mixture of 60% / 40% silicate / graphite is described in
+[Corrales et al. (2016)](http://arxiv.org/abs/1602.01100).
+Please cite this paper if you use this model.
+
+**OlivineAbs model parameters**
+
+    olv_md : dust mass column for Olivine grains in units of 1.e-4
+      default = 1.0
+
+    redshift : redshift of the obscuring dust (XSPEC only)
+      default = 0.0
+
+A simple model for olivine absorption that uses the silicate absorption model
+from ISMdust with the [Rogantini et al. (2018)](http://adsabs.harvard.edu/abs/2018A%26A...609A..22R)
+cross-section for the Fe K edge incorporated.
+Please cite both Corrales et al. (2016) and Rogantini et al. (2018)
+if you use this model.
+
+## XSPEC installation
+
+### ISMdust setup
 
 Enter the ismdust directory and start XSPEC
 
@@ -21,14 +66,37 @@ Finally, set an environment variable to point to location of ismdust, e.g.
 
 Now, when you want to load the ISMDUST model in XSPEC, load the correct library. The following example loads ismdust and then applies the extinction model to a power law component:
 
-    XSPEC12> load libismdust.dylib
+    XSPEC12> lmod ismdust $::env(ISMDUSTROOT)
     XSPEC12> mo ismdust*pow
-
+
+Note: You can omit the `$::env(ISMDUSTROOT)` portion if you have the `LMODDIR`
+environment variable set to the location where ISMdust is installed.
+
 Try the test file to make sure it's working.
 
     XSPEC12> @test.xcm
 
-## Setup for ismdust fits with ISIS (Interactive Spectral Interpretation System) models
+### OlivineAbs setup
+
+Follow the instructions above to set your `ISMDUSTROOT` environment variable
+(and `LMODDIR`, if you choose).
+
+Now install the local model
+
+    XSPEC12> initpackage olivineabs lmodel_olivineabs.dat .
+
+To load the model in XSPEC:
+
+    XSPEC12> lmod olivineabs $::env(ISMDUSTROOT)
+    XSPEC12> mo olivineabs*pow
+
+Try the test file to make sure it's working.
+
+    XSPEC12> @test_olivine.xcm
+
+## ISIS (Interactive Spectral Interpretation System) installation
+
+### ISMdust setup
 
 Add a line to your .isisrc file
 
@@ -44,22 +112,39 @@ When you want to invoke the model, use the require function in ISIS to load ismd
 
 To set up the model extinction model with a power law continuum, for example, do:
 
+    isis> require("ismdust");
     isis> fit_fun("ismdust(1, powerlaw(1))");
 
-## Setup for Fe-L edge fits with ISIS (Interactive Spectral Interpretation System) models
+**For absorption component only**
 
-Add a line to your .isisrc file
+    isis> require("ismdust");
+    isis> fit_fun("ismdust_abs(1, powerlaw(1))");
 
-    add_to_isis_load_path("/path/to/ismdust/ismdust_isis");
+**For scattering component only**
 
-Set an environment variable to point to the location of the Fe-L edge templates:
+    isis> require("ismdust");
+    isis> fit_fun("ismdust_sca(1, powerlaw(1))");
 
-    export FEPATH=/path/to/ismdust/ismdust_isis/
+See `ismdust_isis/test_ismdust.sl`
 
-When you want to invoke the model, use the require function in ISIS to load the model.
+### OlivineAbs setup
 
-    isis> require("FeLedge");
+Use the same set up instructions as above.
+
+To run the model with a power law continuum, for example, do:
+
+    isis> require("olivineabs");
+    isis> fit_fun("olivineabs(1, powerlaw(1))");
+
+See also `ismdust_isis/test_olivine.sl`
+
+### Fe-L edge fits with ISIS
+
+Set an environment variable to point to the location of the Fe-L edge templates:
+
+    export FEPATH=/path/to/ismdust/ismdust_isis/
 
 To invoke the model:
 
+    isis> require("FeLedge");
     isis> fit_fun("FeLedge(1, powerlaw(1))");

ismdust_isis/olivineabs.sl

@@ -0,0 +1,48 @@
+#! /usr/bin/env isis
+
+%% olivineabs.sl -- Reads in table model data olivine absorption and applies it.
+%%    This is a user defined model for fitting with ISIS.
+%%    If using this model, please cite Corrales et al. (2016)
+%%    http://adsabs.harvard.edu/abs/2016MNRAS.458.1345C
+%%    and Rogantini et al. (2018)
+%%    http://adsabs.harvard.edu/abs/2018A%26A...609A..22R
+%%
+%% 2018.06.23 - [email protected]
+%%---------------------------------------------------------------------
+
+_traceback = 1;
+
+variable ISMPATH = getenv("ISMDUSTROOT");
+variable OLIV_FILE = ISMPATH + "/edge_files/olivine_abs.fits";
+
+private define read_xsect( filename )
+{
+    variable wavel, xabs;
+    (wavel, xabs) = fits_read_col(filename, "angstroms", "abs");
+
+    return struct{wavel=wavel[[::-1]], tau_abs=xabs[[::-1]]};
+}
+
+static variable X_lo, X_hi, Value;
+static variable olv_ext = read_xsect(OLIV_FILE);
+
+define olivineabs_fit(lo, hi, par, fun)
+{
+    variable olv_md  = par[0];
+    variable Angs = 0.5*(lo+hi);
+    variable tau  = olv_md * interpol(Angs,olv_ext.wavel,olv_ext.tau_abs);
+
+    return fun * exp(-tau);
+}
+
+define olivineabs_defaults(i)
+{
+    switch(i)
+    { case 0: % olv_md
+    return ( 1.0, 0, 0, 10000 );
+    }
+}
+
+add_slang_function ("olivineabs", ["olv_md"]);
+set_function_category("olivineabs", ISIS_FUN_OPERATOR);
+set_param_default_hook("olivineabs", "olivineabs_defaults");

ismdust_isis/test_ismdust.sl

@@ -0,0 +1,21 @@
+
+variable x1, x2, yext, yabs, ysca;
+
+require("ismdust");
+(x1, x2) = linear_grid(1.0, 40.0, 10000);
+
+fit_fun("ismdust(1, constant(1))");
+yext = eval_fun(x1, x2);
+
+fit_fun("ismdust_abs(1, constant(1))");
+yabs = eval_fun(x1, x2);
+
+fit_fun("ismdust_sca(1, constant(1))");
+ysca = eval_fun(x1, x2);
+
+% plots the cross-sections (tau, unitless)
+
+xlog; ylog;
+hplot(x1, x2, -log(yext), 1);
+ohplot(x1, x2, -log(yabs), 2);
+ohplot(x1, x2, -log(ysca), 4);

ismdust_isis/test_olivine.sl

@@ -0,0 +1,10 @@
+
+variable x1, x2, y;
+
+require("olivineabs");
+(x1, x2) = linear_grid(1.0, 2.5, 1000);
+fit_fun("olivineabs(1, powerlaw(1))");
+set_par(3, 100.0);
+y = eval_fun(x1, x2);
+xlog; ylog;
+hplot(x1, x2, y);

lmodel_olivineabs.dat

@@ -0,0 +1,3 @@
+olivineabs       2   0.01    1.e6      olivineabs     mul  0
+moliv    10^-4    1.0    0.    0.    1E4  1E5    1E-3
+redshift  " "    0.0    -1.0  0.0  10.   10.    -0.01