GCC09 basic code - still need docs and tests

dust_extinction/dust_extinction.py

@@ -13,13 +13,14 @@
 
 __all__ = ['CCM89', 'FM90', 'P92', 'F99',
            'G03_SMCBar', 'G03_LMCAvg', 'G03_LMC2',
-           'G16']
+           'GCC09', '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_F99 = [0.3,10.0]
 x_range_G03 = [0.3,10.0]
+x_range_GCC09 = [3.3,11.0]
 x_range_G16 = x_range_G03
 
 def _test_valid_x_range(x, x_range, outname):
@@ -340,7 +341,7 @@ def fA(self, value):
                                       + str(self.fA_range[0])
                                       + " and "
                                       + str(self.fA_range[1]))
-        
+
 class CCM89(BaseExtRvModel):
     """
     CCM89 extinction model calculation
@@ -1388,6 +1389,132 @@ def evaluate(self, in_x):
                                  self.x_range, 'G03')
 
 
+class GCC09(BaseExtRvModel):
+    """
+    GCC09 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
+    -----
+    GCC09 Milky Way R(V) dependent extinction model
+
+    From Gordon, Cartledge, & Clayton (2009, ApJ, 705, 1320)
+       and errata 2014, ApJ, 781, 128
+
+    Example showing GCC09 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 GCC09
+
+        fig, ax = plt.subplots()
+
+        # generate the curves and plot them
+        x = np.arange(3.3,11.0,0.1)/u.micron
+
+        Rvs = ['2.0','3.0','4.0','5.0','6.0']
+        for cur_Rv in Rvs:
+           ext_model = GCC09(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_GCC09
+
+    @staticmethod
+    def evaluate(in_x, Rv):
+        """
+        GCC09 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_GCC09, 'GCC09')
+
+        # 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 <= 11.0))
+        fuv_indxs = np.where(np.logical_and(5.9 <= x,x <= 11.0))
+
+        # 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((.32999, -.7753, .01979, .72085, -.02427,
+        #                           -.50447, .17699, 1), y)
+        #b[opt_indxs] = np.polyval((-2.09002, 5.3026, -.62251, -5.38434,
+        #                           1.07233, 2.28305, 1.41338, 0), y)
+
+        # NUV
+        a[nuv_indxs] = 1.894-.373*x[nuv_indxs] \
+                       - 0.0101/((x[nuv_indxs] - 4.57)**2 + .0384)
+        b[nuv_indxs] = -3.490 + \
+                       2.057*x[nuv_indxs] \
+                       + 0.706/((x[nuv_indxs] - 4.59)**2 + .169)
+
+        # far-NUV
+        y = x[fuv_indxs] - 5.9
+        a[fuv_indxs] += -.110*(y**2) - .00100*(y**3)
+        b[fuv_indxs] += .531*(y**2) + .0544*(y**3)
+
+        # return A(x)/A(V)
+        return a + b/Rv
+
+
 class G16(BaseExtRvAfAModel):
     """
     G16 extinction model calculation