Merge pull request #202 from karllark/add_astrodust

Add astrodust model

LICENSE

@@ -1,6 +1,6 @@
 BSD 3-Clause License
 
-Copyright (c) 2017, Karl Gordon
+Copyright (c) 2017-2023, Karl Gordon
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without

docs/dust_extinction/choose_model.rst

@@ -172,6 +172,8 @@ observationally (e.g., in the extreme UV below 912 A).
 +--------------+----------------+------------------+--------------+
 |   J13 MWRV31 |   0.00001 - 25 |    0.04 - 100000 |  MW R(V)=3.1 |
 +--------------+----------------+------------------+--------------+
+|  HD23 MWRV31 |  0.000033 - 10 |      0.1 - 30000 |  MW R(V)=3.1 |
++--------------+----------------+------------------+--------------+
 
 Shape Models
 ============

docs/dust_extinction/model_flavors.rst

@@ -351,7 +351,7 @@ Grain models
    import matplotlib.pyplot as plt
    import astropy.units as u
 
-   from dust_extinction.grain_models import DBP90, WD01, D03, ZDA04, C11, J13
+   from dust_extinction.grain_models import DBP90, WD01, D03, ZDA04, C11, J13, HD23
 
    fig, ax = plt.subplots()
 
@@ -363,12 +363,14 @@ Grain models
              WD01, WD01, WD01,
              D03, D03, D03,
              ZDA04,
-             C11, J13]
+             C11, J13,
+             HD23]
    modelnames = ["MWRV31",
                  "MWRV31", "MWRV40", "MWRV55",
                  "MWRV31", "MWRV40", "MWRV55",
                  "BARE-GR-S",
-                 "MWRV31", "MWRV31"]
+                 "MWRV31", "MWRV31",
+                 "MWRV31"]
 
    for cmodel, cname in zip(models, modelnames):
       ext_model = cmodel(cname)
@@ -400,7 +402,7 @@ Grain models
    from matplotlib.ticker import ScalarFormatter
    import astropy.units as u
 
-   from dust_extinction.grain_models import DBP90, WD01, D03, ZDA04, C11, J13
+   from dust_extinction.grain_models import DBP90, WD01, D03, ZDA04, C11, J13, HD23
 
    fig, ax = plt.subplots()
 
@@ -412,12 +414,14 @@ Grain models
              WD01, WD01, WD01,
              D03, D03, D03,
              ZDA04,
-             C11, J13]
+             C11, J13,
+             HD23]
    modelnames = ["MWRV31",
                  "MWRV31", "MWRV40", "MWRV55",
                  "MWRV31", "MWRV40", "MWRV55",
                  "BARE-GR-S",
-                 "MWRV31", "MWRV31"]
+                 "MWRV31", "MWRV31",
+                 "MWRV31"]
 
    for cmodel, cname in zip(models, modelnames):
       ext_model = cmodel(cname)

docs/dust_extinction/references.rst

@@ -52,8 +52,11 @@ G16: `Gordon et al. 2016, ApJ 826, 104
 G21: `Gordon et al. 2021, ApJ, 916, 33
 <https://ui.adsabs.harvard.edu/abs/2021ApJ...916...33G>`_
 
-G23: `Gordon et al. 2022, ApJ, in press
-<https://arxiv.org/abs/2304.01991>`_
+G23: `Gordon et al. 2022, ApJ, 950, 86
+<https://ui.adsabs.harvard.edu/abs/2023ApJ...950...86G/abstract>`_
+
+HD23: `Hensley & Draine 2023, ApJ, 948, 55
+<https://ui.adsabs.harvard.edu/abs/2023ApJ...948...55H/abstract>`_
 
 I05: `Indebetouw et al. 2005, ApJ, 619, 931
 <https://ui.adsabs.harvard.edu/abs/2005ApJ...619..931I>`_

dust_extinction/grain_models.py

@@ -4,11 +4,12 @@
 
 from astropy.table import Table
 from astropy.modeling import InputParameterError
+from astropy.io.fits import getdata
 
 from .helpers import _get_x_in_wavenumbers, _test_valid_x_range
 from .baseclasses import BaseExtModel
 
-__all__ = ["DBP90", "WD01", "D03", "ZDA04", "C11", "J13"]
+__all__ = ["DBP90", "WD01", "D03", "ZDA04", "C11", "J13", "HD23"]
 
 
 class GMBase(BaseExtModel):
@@ -114,7 +115,6 @@ class DBP90(GMBase):
     possnames = {"MWRV31": ("EXT_DBP90.RES.dat", 3.1)}
 
     def __init__(self, modelname="MWRV31", **kwargs):
-
         if modelname not in self.possnames.keys():
             raise InputParameterError("modelname not recognized")
         filename = self.possnames[modelname][0]
@@ -124,7 +124,10 @@ def __init__(self, modelname="MWRV31", **kwargs):
         data_path = pkg_resources.resource_filename("dust_extinction", "data/")
 
         a = Table.read(
-            data_path + filename, data_start=1, header_start=None, format="ascii.basic",
+            data_path + filename,
+            data_start=1,
+            header_start=None,
+            format="ascii.basic",
         )
 
         self.data_x = 1.0 / a["col1"].data
@@ -203,7 +206,6 @@ class WD01(GMBase):
     }
 
     def __init__(self, modelname="MWRV31", **kwargs):
-
         if modelname not in self.possnames.keys():
             raise InputParameterError("modelname not recognized")
         filename = self.possnames[modelname][0]
@@ -296,7 +298,6 @@ class D03(GMBase):
     }
 
     def __init__(self, modelname="MWRV31", **kwargs):
-
         if modelname not in self.possnames.keys():
             raise InputParameterError("modelname not recognized")
         filename = self.possnames[modelname][0]
@@ -385,7 +386,6 @@ class ZDA04(GMBase):
     possnames = {"BARE-GR-S": ("zubko2004_bare-gr-s_alam_av.dat", 3.1)}
 
     def __init__(self, modelname="BARE-GR-S", **kwargs):
-
         if modelname not in self.possnames.keys():
             raise InputParameterError("modelname not recognized")
         filename = self.possnames[modelname][0]
@@ -394,7 +394,10 @@ def __init__(self, modelname="BARE-GR-S", **kwargs):
         # get the tabulated information
         data_path = pkg_resources.resource_filename("dust_extinction", "data/")
 
-        a = Table.read(data_path + filename, format="ascii.basic",)
+        a = Table.read(
+            data_path + filename,
+            format="ascii.basic",
+        )
 
         self.data_x = 1.0 / a["lam[um]"].data
 
@@ -464,7 +467,6 @@ class C11(GMBase):
     possnames = {"MWRV31": ("EXT_C11.RES.dat", 3.1)}
 
     def __init__(self, modelname="MWRV31", **kwargs):
-
         if modelname not in self.possnames.keys():
             raise InputParameterError("modelname not recognized")
         filename = self.possnames[modelname][0]
@@ -474,7 +476,10 @@ def __init__(self, modelname="MWRV31", **kwargs):
         data_path = pkg_resources.resource_filename("dust_extinction", "data/")
 
         a = Table.read(
-            data_path + filename, data_start=1, header_start=None, format="ascii.basic",
+            data_path + filename,
+            data_start=1,
+            header_start=None,
+            format="ascii.basic",
         )
 
         self.data_x = 1.0 / a["col1"].data
@@ -545,7 +550,6 @@ class J13(GMBase):
     possnames = {"MWRV31": ("EXT_J13.RES.dat", 3.1)}
 
     def __init__(self, modelname="MWRV31", **kwargs):
-
         if modelname not in self.possnames.keys():
             raise InputParameterError("modelname not recognized")
         filename = self.possnames[modelname][0]
@@ -555,7 +559,10 @@ def __init__(self, modelname="MWRV31", **kwargs):
         data_path = pkg_resources.resource_filename("dust_extinction", "data/")
 
         a = Table.read(
-            data_path + filename, data_start=1, header_start=None, format="ascii.basic",
+            data_path + filename,
+            data_start=1,
+            header_start=None,
+            format="ascii.basic",
         )
 
         self.data_x = 1.0 / a["col1"].data
@@ -569,3 +576,81 @@ def __init__(self, modelname="MWRV31", **kwargs):
         self.data_tolerance = 1e-6
 
         super().__init__(**kwargs)
+
+
+class HD23(GMBase):
+    r"""
+    Hensley & Draine (2023) Grain Model
+
+    Parameters
+    ----------
+    None
+
+    Raises
+    ------
+    None
+
+    Notes
+    -----
+    From Hensley & Draine (2023, ApJ, 948, 55).  File from
+    https://dataverse.harvard.edu/dataverse/astrodust
+
+    Example showing the possible curves
+
+    .. plot::
+        :include-source:
+
+        import numpy as np
+        import matplotlib.pyplot as plt
+        import astropy.units as u
+
+        from dust_extinction.grain_models import HD23
+
+        fig, ax = plt.subplots()
+
+        ext_model = HD23()
+
+        lam = np.logspace(np.log10(1.0/ext_model.x_range[1]),
+                          np.log10(1.0/ext_model.x_range[0]),
+                          num=1000)
+        x = (1.0 / lam) / u.micron
+
+        # define the extinction model
+        ax.plot(lam,ext_model(x),label=ext_model.__class__.__name__)
+
+        ax.set_xlabel(r'$\lambda$ [$\mu m$]')
+        ax.set_ylabel(r'$A(x)/A(V)$')
+
+        ax.set_xscale('log')
+        ax.set_yscale('log')
+
+        ax.legend(loc='best')
+        plt.show()
+    """
+
+    x_range = [1.0 / 3e4, 1.0 / 0.1]
+
+    possnames = {"MWRV31": ("astrodust+PAH_MW_RV3.1.fits", 3.1)}
+
+    def __init__(self, modelname="MWRV31", **kwargs):
+        if modelname not in self.possnames.keys():
+            raise InputParameterError("modelname not recognized")
+        filename = self.possnames[modelname][0]
+        self.Rv = self.possnames[modelname][1]
+
+        # get the tabulated information
+        data_path = pkg_resources.resource_filename("dust_extinction", "data/")
+
+        a = getdata(data_path + filename, 2)
+
+        self.data_x = 1.0 / a[:, 0]
+
+        # normalized by wavelength closest to V band
+        sindxs = np.argsort(np.absolute(self.data_x - 1.0 / 0.55))
+
+        self.data_axav = a[:, 3] / a[sindxs[0], 3]
+
+        # accuracy of the data based on calculations
+        self.data_tolerance = 1e-6
+
+        super().__init__(**kwargs)

dust_extinction/tests/helpers.py

@@ -29,7 +29,7 @@
     G21_MWAvg,
     D22_MWAvg,
 )
-from dust_extinction.grain_models import DBP90, WD01, D03, ZDA04, C11, J13
+from dust_extinction.grain_models import DBP90, WD01, D03, ZDA04, C11, J13, HD23
 
 param_ave_models_Rv = [CCM89, O94, F99, F04, VCG04, GCC09, M14, F19, D22, G23]
 param_ave_models_Rv_fA = [G16]
@@ -50,7 +50,7 @@
     G21_MWAvg,
     D22_MWAvg,
 ]
-grain_models = [DBP90, WD01, D03, ZDA04, C11, J13]
+grain_models = [DBP90, WD01, D03, ZDA04, C11, J13, HD23]
 
 all_models = (
     param_ave_models_Rv