Merge 16ee6dd into 843cf82

.rtd-environment.yml

@@ -8,5 +8,6 @@ dependencies:
     - Cython
     - matplotlib
     - numpy
+    - scipy
 #   - pip:
 #       - dependency_from_pip

.travis.yml

@@ -38,7 +38,8 @@ env:
         # For this package-template, we include examples of Cython modules,
         # so Cython is required for testing. If your package does not include
         # Cython code, you can set CONDA_DEPENDENCIES=''
-        - CONDA_DEPENDENCIES='Cython'
+        - CONDA_DEPENDENCIES='Cython scipy'
+        #- CONDA_ALL_DEPENDENCIES='Cython scipy'
 
         # Conda packages for affiliated packages are hosted in channel
         # "astropy" while builds for astropy LTS with recent numpy versions

dust_extinction/dust_extinction.py

@@ -9,14 +9,16 @@
 import warnings
 
 import numpy as np
+from scipy import interpolate
 
 import astropy.units as u
 from astropy.modeling import Model, Parameter, InputParameterError
 
-__all__ = ['BaseExtRvModel', 'CCM89', 'FM90']
+__all__ = ['BaseExtRvModel', 'CCM89', 'FM90', 'F99']
 
 x_range_CCM89 = [0.3,10.0]
 x_range_FM90 = [1.0/0.32,1.0/0.0912]
+x_range_F99 = [1.0/6.0,8.7]
 
 class BaseExtRvModel(Model):
     """
@@ -108,10 +110,6 @@ class CCM89(BaseExtRvModel):
 
     From Cardelli, Clayton, and Mathis (1989, ApJ, 345, 245)
 
-    FM07 should be used instead as it is based on 10x more observations 
-    and a better treatment of the optical/NIR photometry based portion 
-    of the curves.
-
     Example showing CCM89 curves for a range of R(V) values.
 
     .. plot::
@@ -274,7 +272,7 @@ class FM90(Model):
         fig, ax = plt.subplots()
 
         # generate the curves and plot them
-        x = np.arange(3.2,10.0,0.1)/u.micron
+        x = np.arange(3.8,8.6,0.1)/u.micron
     
         ext_model = FM90()
         ax.plot(x,ext_model(x),label='total')
@@ -369,3 +367,175 @@ def evaluate(in_x, C1, C2, C3, C4, xo, gamma):
         # return E(x-V)/E(B-V)
         return exvebv
 
+class F99(BaseExtRvModel):
+    """
+    F99 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
+    -----
+    F99 Milky Way R(V) dependent extinction model
+
+    From Fitzpatrick (1999, PASP, 111, 63)
+
+    Updated for the C1 vs C2 correlation in
+       Fitzpatrick & Massa (2007, ApJ, 663, 320)
+
+    Example showing F99 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 F99
+
+        fig, ax = plt.subplots()
+
+        # generate the curves and plot them
+        x = np.arange(0.2,8.7,0.1)/u.micron
+    
+        Rvs = ['2.0','3.0','4.0','5.0','6.0']
+        for cur_Rv in Rvs:
+           ext_model = F99(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_F99
+
+    @staticmethod
+    def evaluate(in_x, Rv):
+        """
+        F99 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
+        if np.logical_or(np.any(x < x_range_F99[0]),
+                         np.any(x > x_range_F99[1])):
+            raise ValueError('Input x outside of range defined for F99' \
+                             + ' ['
+                             + str(x_range_F99[0])
+                             +  ' <= x <= '
+                             + str(x_range_F99[1])
+                             + ', x has units 1/micron]')
+
+        # ensure Rv is a single element, not numpy array
+        Rv = Rv[0]
+
+        # initialize extinction curve storage
+        axebv = np.zeros(len(x))
+
+        # **Ultraviolet** portion
+        #   calculated using the FM90 parameterization
+
+        # constant terms
+        C3 = 3.23
+        C4 = 0.41
+        xo = 4.596
+        gamma = 0.99
+
+        # terms depending on Rv
+        C2 = -0.824 + 4.717/Rv
+        # updated for FM07 correlation between C1 and C2
+        C1 = 2.030 - 3.007*C2
+
+        # x value above which FM90 parametrization used
+        x_cutval_uv = 10000.0/2700.0
+
+        # required UV points for spline interpolation
+        x_splineval_uv = 10000.0/np.array([2700.0,2600.0])
+
+        # UV points in input x
+        indxs_uv = np.where(x >= x_cutval_uv)
+
+        # add in required spline points, otherwise just spline points
+        if len(indxs_uv) > 0:
+            xuv = np.concatenate([x_splineval_uv, x[indxs_uv]])
+        else:
+            xuv = x_splineval_uv
+
+        # FM90 model and values
+        fm90_model = FM90(C1=C1, C2=C2, C3=C3, C4=C4, xo=xo, gamma=gamma)
+        # evaluate model and get results in A(x)/E(B-V)
+        axebv_fm90 = fm90_model(xuv) + Rv
+
+        # save spline points
+        y_splineval_uv = axebv_fm90[0:2]
+
+        # ingore the spline points
+        if len(indxs_uv) > 0:
+            axebv[indxs_uv] = axebv_fm90[2:]
+
+        # **Optical Portion**
+        #   using cubic spline anchored in UV, optical, and IR
+
+        # optical/NIR points in input x
+        indxs_opir = np.where(x < x_cutval_uv)
+
+        # spline points
+        x_splineval_optir = 10000./np.array([26500.0,12200.0,6000.0,
+                                             5470.0,4670.0,4110.0])
+        # add in zero extinction at infinite wavelength
+        x_splineval_optir = np.insert(x_splineval_optir, 0, 0.0)
+
+        # determine optical/IR values at splint points
+        y_splineval_opt   = np.array([-0.426 + 1.0044*Rv,
+                                      -0.050 + 1.0016*Rv,
+                                      0.701 + 1.016*Rv,
+                                      1.208 + 1.0032*Rv - 0.00033*(Rv**2)])
+        y_splineval_ir    = np.array([0.0,0.265,0.829])*Rv/3.1 
+        y_splineval_optir = np.concatenate([y_splineval_ir,y_splineval_opt])
+
+        if len(indxs_opir) > 0:
+            spline_x = np.concatenate([x_splineval_optir, x_splineval_uv])
+            spline_y = np.concatenate([y_splineval_optir, y_splineval_uv])
+            spline_rep = interpolate.splrep(spline_x, spline_y)
+            axebv[indxs_opir] = interpolate.splev(x[indxs_opir],
+                                                  spline_rep, der=0)
+
+        # return A(x)/A(V)
+        return axebv/Rv
+

dust_extinction/tests/test_f99.py

@@ -0,0 +1,82 @@
+import numpy as np
+import pytest
+
+import astropy.units as u
+from astropy.modeling import InputParameterError
+
+from ..dust_extinction import F99
+
+x_bad = [-1.0, 0.1, 9.0, 100.]
+@pytest.mark.parametrize("x_invalid", x_bad)
+def test_invalid_wavenumbers(x_invalid):
+    tmodel = F99()
+    with pytest.raises(ValueError) as exc:
+        tmodel(x_invalid)
+    assert exc.value.args[0] == 'Input x outside of range defined for F99' \
+                                + ' [' \
+                                + str(tmodel.x_range[0]) \
+                                +  ' <= x <= ' \
+                                + str(tmodel.x_range[1]) \
+                                + ', x has units 1/micron]'
+
+@pytest.mark.parametrize("x_invalid_wavenumber", x_bad/u.micron)
+def test_invalid_wavenumbers_imicron(x_invalid_wavenumber):
+    tmodel = F99()
+    with pytest.raises(ValueError) as exc:
+        tmodel(x_invalid_wavenumber)
+    assert exc.value.args[0] == 'Input x outside of range defined for F99' \
+                                + ' [' \
+                                + str(tmodel.x_range[0]) \
+                                +  ' <= x <= ' \
+                                + str(tmodel.x_range[1]) \
+                                + ', x has units 1/micron]'
+
+@pytest.mark.parametrize("x_invalid_micron", u.micron/x_bad)
+def test_invalid_micron(x_invalid_micron):
+    tmodel = F99()
+    with pytest.raises(ValueError) as exc:
+        tmodel(x_invalid_micron)
+    assert exc.value.args[0] == 'Input x outside of range defined for F99' \
+                                + ' [' \
+                                + 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/x_bad)
+def test_invalid_micron(x_invalid_angstrom):
+    tmodel = F99()
+    with pytest.raises(ValueError) as exc:
+        tmodel(x_invalid_angstrom)
+    assert exc.value.args[0] == 'Input x outside of range defined for F99' \
+                                + ' [' \
+                                + str(tmodel.x_range[0]) \
+                                +  ' <= x <= ' \
+                                + str(tmodel.x_range[1]) \
+                                + ', x has units 1/micron]'
+
+def get_axav_cor_vals():
+    # testing wavenumbers
+    x = np.array([0.5, 1.0, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0])
+
+    # add units
+    x = x/u.micron
+
+    # correct values
+    cor_vals = np.array([0.124997,  0.377073,  1.130636,  1.419644,
+                         1.630377,  1.888546,  2.275900,  3.014577,
+                         2.762256,  2.475272,  2.711508,  3.197144])
+
+    return (x, cor_vals)
+
+
+def test_extinction_F99_values():
+    # get the correct values
+    x, cor_vals = get_axav_cor_vals()
+
+    # initialize extinction model    
+    tmodel = F99()
+
+    # test
+    np.testing.assert_allclose(tmodel(x), cor_vals, rtol=1e-05)
+

setup.cfg

@@ -45,7 +45,7 @@ license = BSD
 url = http://astropy.org/
 edit_on_github = False
 github_project = astropy/astropy
-install_requires = astropy
+install_requires = astropy scipy
 # version should be PEP440 compatible (http://www.python.org/dev/peps/pep-0440)
 version = 0.1.dev