O94 model w/ tests and docs

docs/dust_extinction/model_flavors.rst

@@ -56,10 +56,11 @@ R(V) (+ other variables) dependent prediction models
    These models provide predictions of the shape of the dust extinction
    given input variable(s).  
 
-   These include CCM89 the original R(V) dependent model from 
-   Cardelli, Clayton, and Mathis (1989) and updated versions F99 
+   The R(V) dependent models include CCM89 the original such model
+   (Cardelli, Clayton, and Mathis 1989), the O94 model that updates the
+   optical portion of the CCM89 model (O'Donnell 1994), and F99 model
    (Fitzpatrick 1999).  These models are based on the average
-   behavior of extinction in the Milky Way.
+   behavior of extinction in the Milky Way as a function of R(V).
 
    In addition, the (R(V), f_A) two parameter relationship from 
    Gordon et al. (2016) is included.  This model is based on the average
@@ -93,6 +94,33 @@ R(V) (+ other variables) dependent prediction models
    plt.tight_layout()
    plt.show()
 
+.. plot::
+
+   import numpy as np
+   import matplotlib.pyplot as plt
+   import astropy.units as u
+
+   from dust_extinction.dust_extinction import O94
+
+   fig, ax = plt.subplots()
+
+   # generate the curves and plot them
+   x = np.arange(0.5,10.0,0.1)/u.micron
+
+   Rvs = ['2.0','3.0','4.0','5.0','6.0']
+   for cur_Rv in Rvs:
+      ext_model = O94(Rv=cur_Rv)
+      ax.plot(x,ext_model(x),label='R(V) = ' + str(cur_Rv))
+
+   ax.set_xlabel('$x$ [$\mu m^{-1}$]')
+   ax.set_ylabel('$A(x)/A(V)$')
+
+   ax.set_title('O94')
+
+   ax.legend(loc='best')
+   plt.tight_layout()
+   plt.show()
+
 .. plot::
 
    import numpy as np

dust_extinction/dust_extinction.py

@@ -11,13 +11,14 @@
 from astropy.modeling import (Model, Fittable1DModel,
                               Parameter, InputParameterError)
 
-__all__ = ['CCM89', 'FM90', 'P92', 'F99',
+__all__ = ['CCM89', 'FM90', 'P92', 'O94', 'F99',
            'G03_SMCBar', 'G03_LMCAvg', 'G03_LMC2',
            'G16']
 
 x_range_CCM89 = [0.3,10.0]
 x_range_FM90 = [1.0/0.32,1.0/0.0912]
 x_range_P92 = [1.0/1e3,1.0/1e-3]
+x_range_O94 = x_range_CCM89
 x_range_F99 = [0.3,10.0]
 x_range_G03 = [0.3,10.0]
 x_range_G16 = x_range_G03
@@ -927,6 +928,137 @@ def evaluate(self, in_x,
     # use numerical derivaties (need to add analytic)
     fit_deriv = None
 
+class O94(BaseExtRvModel):
+    """
+    O94 extinction model calculation
+
+    Parameters
+    ----------
+    Rv: float
+        R(V) = A(V)/E(B-V) = total-to-selective extinction
+
+    Raises
+    ------
+    InputParameterError
+       Input Rv values outside of defined range
+
+    Notes
+    -----
+    O94 Milky Way R(V) dependent extinction model
+
+    From O'Donnell (1994, ApJ, 422, 158)
+      Updates/improves the optical portion of the CCM89 model
+
+    Example showing O94 curves for a range of R(V) values.
+
+    .. plot::
+        :include-source:
+
+        import numpy as np
+        import matplotlib.pyplot as plt
+        import astropy.units as u
+
+        from dust_extinction.dust_extinction import O94
+
+        fig, ax = plt.subplots()
+
+        # generate the curves and plot them
+        x = np.arange(0.5,10.0,0.1)/u.micron
+
+        Rvs = ['2.0','3.0','4.0','5.0','6.0']
+        for cur_Rv in Rvs:
+           ext_model = O94(Rv=cur_Rv)
+           ax.plot(x,ext_model(x),label='R(V) = ' + str(cur_Rv))
+
+        ax.set_xlabel('$x$ [$\mu m^{-1}$]')
+        ax.set_ylabel('$A(x)/A(V)$')
+
+        ax.legend(loc='best')
+        plt.show()
+    """
+
+    Rv_range = [2.0,6.0]
+    x_range = x_range_O94
+
+    @staticmethod
+    def evaluate(in_x, Rv):
+        """
+        O94 function
+
+        Parameters
+        ----------
+        in_x: float
+           expects either x in units of wavelengths or frequency
+           or assumes wavelengths in wavenumbers [1/micron]
+
+           internally wavenumbers are used
+
+        Returns
+        -------
+        axav: np array (float)
+            A(x)/A(V) extinction curve [mag]
+
+        Raises
+        ------
+        ValueError
+           Input x values outside of defined range
+        """
+        # convert to wavenumbers (1/micron) if x input in units
+        # otherwise, assume x in appropriate wavenumber units
+        with u.add_enabled_equivalencies(u.spectral()):
+            x_quant = u.Quantity(in_x, 1.0/u.micron, dtype=np.float64)
+
+        # strip the quantity to avoid needing to add units to all the
+        #    polynomical coefficients
+        x = x_quant.value
+
+        # check that the wavenumbers are within the defined range
+        _test_valid_x_range(x, x_range_O94, 'O94')
+
+        # setup the a & b coefficient vectors
+        n_x = len(x)
+        a = np.zeros(n_x)
+        b = np.zeros(n_x)
+
+        # define the ranges
+        ir_indxs = np.where(np.logical_and(0.3 <= x,x < 1.1))
+        opt_indxs = np.where(np.logical_and(1.1 <= x,x < 3.3))
+        nuv_indxs = np.where(np.logical_and(3.3 <= x,x <= 8.0))
+        fnuv_indxs = np.where(np.logical_and(5.9 <= x,x <= 8))
+        fuv_indxs = np.where(np.logical_and(8 < x,x <= 10))
+
+        # Infrared
+        y = x[ir_indxs]**1.61
+        a[ir_indxs] = .574*y
+        b[ir_indxs] = -0.527*y
+
+        # NIR/optical
+        y = x[opt_indxs] - 1.82
+        a[opt_indxs] = np.polyval((-0.505, 1.647, -0.827, -1.718,
+                                   1.137, 0.701, -0.609, 0.104, 1), y)
+        b[opt_indxs] = np.polyval((3.347, -10.805, 5.491, 11.102,
+                                   -7.985, -3.989, 2.908, 1.952, 0), y)
+
+        # NUV
+        a[nuv_indxs] = 1.752-.316*x[nuv_indxs] \
+                       - 0.104/((x[nuv_indxs] - 4.67)**2 + .341)
+        b[nuv_indxs] = -3.09 + \
+                       1.825*x[nuv_indxs] \
+                       + 1.206/((x[nuv_indxs] - 4.62)**2 + .263)
+
+        # far-NUV
+        y = x[fnuv_indxs] - 5.9
+        a[fnuv_indxs] += -.04473*(y**2) - .009779*(y**3)
+        b[fnuv_indxs] += .2130*(y**2) + .1207*(y**3)
+
+        # FUV
+        y = x[fuv_indxs] - 8.0
+        a[fuv_indxs] = np.polyval((-.070, .137, -.628, -1.073), y)
+        b[fuv_indxs] = np.polyval((.374, -.42, 4.257, 13.67), y)
+
+        # return A(x)/A(V)
+        return a + b/Rv
+
 class F99(BaseExtRvModel):
     """
     F99 extinction model calculation

dust_extinction/tests/test_o94.py

@@ -0,0 +1,182 @@
+import numpy as np
+import pytest
+
+import astropy.units as u
+from astropy.modeling import InputParameterError
+
+from ..dust_extinction import O94
+
+@pytest.mark.parametrize("Rv_invalid", [-1.0,0.0,1.9,6.1,10.])
+def test_invalid_Rv_input(Rv_invalid):
+    with pytest.raises(InputParameterError) as exc:
+        tmodel = O94(Rv=Rv_invalid)
+    assert exc.value.args[0] == 'parameter Rv must be between 2.0 and 6.0'
+
+@pytest.mark.parametrize("x_invalid", [-1.0, 0.2, 10.1, 100.])
+def test_invalid_wavenumbers(x_invalid):
+    tmodel = O94(Rv=3.1)
+    with pytest.raises(ValueError) as exc:
+        tmodel(x_invalid)
+    assert exc.value.args[0] == 'Input x outside of range defined for O94' \
+                                + ' [' \
+                                + str(tmodel.x_range[0]) \
+                                +  ' <= x <= ' \
+                                + str(tmodel.x_range[1]) \
+                                + ', x has units 1/micron]'
+
+@pytest.mark.parametrize("x_invalid_wavenumber",
+                         [-1.0, 0.2, 10.1, 100.]/u.micron)
+def test_invalid_wavenumbers_imicron(x_invalid_wavenumber):
+    tmodel = O94(Rv=3.1)
+    with pytest.raises(ValueError) as exc:
+        tmodel(x_invalid_wavenumber)
+    assert exc.value.args[0] == 'Input x outside of range defined for O94' \
+                                + ' [' \
+                                + str(tmodel.x_range[0]) \
+                                +  ' <= x <= ' \
+                                + str(tmodel.x_range[1]) \
+                                + ', x has units 1/micron]'
+
+@pytest.mark.parametrize("x_invalid_micron",
+                         u.micron/[-1.0, 0.2, 10.1, 100.])
+def test_invalid_micron(x_invalid_micron):
+    tmodel = O94(Rv=3.1)
+    with pytest.raises(ValueError) as exc:
+        tmodel(x_invalid_micron)
+    assert exc.value.args[0] == 'Input x outside of range defined for O94' \
+                                + ' [' \
+                                + str(tmodel.x_range[0]) \
+                                +  ' <= x <= ' \
+                                + str(tmodel.x_range[1]) \
+                                + ', x has units 1/micron]'
+
+@pytest.mark.parametrize("x_invalid_angstrom",
+                         u.angstrom*1e4/[-1.0, 0.2, 10.1, 100.])
+def test_invalid_micron(x_invalid_angstrom):
+    tmodel = O94(Rv=3.1)
+    with pytest.raises(ValueError) as exc:
+        tmodel(x_invalid_angstrom)
+    assert exc.value.args[0] == 'Input x outside of range defined for O94' \
+                                + ' [' \
+                                + str(tmodel.x_range[0]) \
+                                +  ' <= x <= ' \
+                                + str(tmodel.x_range[1]) \
+                                + ', x has units 1/micron]'
+
+def test_axav_o94_rv31():
+    # values from Bastiaasen (1992) Table 6
+    x = np.array([2.939, 2.863, 2.778, 2.642, 2.476, 2.385,
+                  2.275, 2.224, 2.124, 2.000, 1.921, 1.849,
+                  1.785, 1.718, 1.637, 1.563, 1.497, 1.408,
+                  1.332, 1.270])
+    cor_vals = np.array([1.725, 1.651, 1.559, 1.431, 1.292, 1.206,
+                         1.100, 1.027, 0.907, 0.738, 0.606, 0.491,
+                         0.383, 0.301, 0.190, 0.098, -0.004, -0.128,
+                         -0.236, -0.327])
+
+    # initialize extinction model
+    tmodel = O94(Rv=3.1)
+
+    # get the model results and change to E(l-1.5)/E(2.2-1.5)
+    mod_vals = tmodel(x)
+    norm_vals = tmodel([1.5, 2.2])
+    mod_vals = (mod_vals - norm_vals[0])/(norm_vals[1] - norm_vals[0])
+
+    # test (table in paper has limited precision)
+    np.testing.assert_allclose(mod_vals, cor_vals, atol=6e-2)
+
+def get_axav_cor_vals(Rv):
+    # testing only NIR or UV wavenumbers (optical tested in previous test)
+    # O94 is the same as CCM89 for these wavelengths
+    x = np.array([ 10.  ,   9.  ,   8.  ,   7.  ,
+                   6.  ,   5.  ,   4.6 ,   4.  ,
+                   0.8 ,   0.63,
+                   0.46])
+
+    # add units
+    x = x/u.micron
+
+    # correct values
+    if Rv == 3.1:
+        cor_vals = np.array([ 5.23835484,  4.13406452,  3.33685933,  2.77962453,
+                              2.52195399,  2.84252644,  3.18598916,  2.31531711,
+                              0.28206957,  0.19200814,
+                              0.11572348])
+    elif Rv == 2.0:
+        cor_vals = np.array([ 9.407     ,  7.3065    ,  5.76223881,  4.60825807,
+                              4.01559036,  4.43845534,  4.93952892,  3.39275574,
+                              0.21678862,  0.14757062,
+                              0.08894094])
+    elif Rv == 3.0:
+        cor_vals = np.array([ 5.491     ,  4.32633333,  3.48385202,  2.8904508 ,
+                              2.6124774 ,  2.9392494 ,  3.2922643 ,  2.38061642,
+                              0.27811315,  0.18931496,
+                              0.11410029])
+    elif Rv == 4.0:
+        cor_vals = np.array([ 3.533     ,  2.83625   ,  2.34465863,  2.03154717,
+                              1.91092092,  2.18964643,  2.46863199,  1.87454675,
+                              0.30877542,  0.21018713,
+                              0.12667997])
+    elif Rv == 5.0:
+        cor_vals = np.array([ 2.3582    ,  1.9422    ,  1.66114259,  1.51620499,
+                              1.48998704,  1.73988465,  1.97445261,  1.57090496,
+                              0.32717278,  0.22271044,
+                              0.13422778])
+    elif Rv == 6.0:
+        cor_vals = np.array([ 1.575     ,  1.34616667,  1.20546523,  1.17264354,
+                              1.20936444,  1.44004346,  1.64499968,  1.36847709,
+                              0.33943769,  0.23105931,
+                              0.13925965])
+    else:
+        cor_vals = np.array([ 0.0 ])
+
+    return (x, cor_vals)
+
+
+@pytest.mark.parametrize("Rv", [2.0, 3.0, 3.1, 4.0, 5.0, 6.0])
+def test_extinction_O94_values(Rv):
+    # get the correct values
+    x, cor_vals = get_axav_cor_vals(Rv)
+
+    # initialize extinction model
+    tmodel = O94(Rv=Rv)
+
+    # test
+    np.testing.assert_allclose(tmodel(x), cor_vals)
+
+def test_extinguish_no_av_or_ebv():
+    tmodel = O94()
+    with pytest.raises(InputParameterError) as exc:
+        tmodel.extinguish([1.0])
+    assert exc.value.args[0] == 'neither Av or Ebv passed, one required'
+
+@pytest.mark.parametrize("Rv", [2.0, 3.0, 3.1, 4.0, 5.0, 6.0])
+def test_extinction_O94_extinguish_values_Av(Rv):
+    # get the correct values
+    x, cor_vals = get_axav_cor_vals(Rv)
+
+    # calculate the cor_vals in fractional units
+    Av = 1.0
+    cor_vals = np.power(10.0,-0.4*(cor_vals*Av))
+
+    # initialize extinction model
+    tmodel = O94(Rv=Rv)
+
+    # test
+    np.testing.assert_allclose(tmodel.extinguish(x, Av=Av), cor_vals)
+
+@pytest.mark.parametrize("Rv", [2.0, 3.0, 3.1, 4.0, 5.0, 6.0])
+def test_extinction_O94_extinguish_values_Ebv(Rv):
+    # get the correct values
+    x, cor_vals = get_axav_cor_vals(Rv)
+
+    # calculate the cor_vals in fractional units
+    Ebv = 1.0
+    Av = Ebv*Rv
+    cor_vals = np.power(10.0,-0.4*(cor_vals*Av))
+
+    # initialize extinction model
+    tmodel = O94(Rv=Rv)
+
+    # test
+    np.testing.assert_allclose(tmodel.extinguish(x, Ebv=Ebv), cor_vals)