Merge pull request #25 from karllark/fix_f99_bug

Fix 1 coefficient in f99 and update tests

.travis.yml

@@ -79,10 +79,10 @@ matrix:
         - os: linux
           env: SETUP_CMD='build_docs -w'
 
-        # Now try Astropy dev and LTS vesions with the latest 3.x and 2.7.
-        - os: linux
-          env: PYTHON_VERSION=2.7 ASTROPY_VERSION=development
-               EVENT_TYPE='pull_request push cron'
+        # Now try Astropy dev and LTS vesions with the latest 3.x.
+        #- os: linux
+        #  env: PYTHON_VERSION=2.7 ASTROPY_VERSION=development
+        #       EVENT_TYPE='pull_request push cron'
         - os: linux
           env: ASTROPY_VERSION=development
                EVENT_TYPE='pull_request push cron'
@@ -117,8 +117,6 @@ matrix:
           env: MAIN_CMD='pycodestyle dust_extinction --count' SETUP_CMD=''
 
     allow_failures:
-        - os: linux
-          env: PYTHON_VERSION=3.3 NUMPY_VERSION=1.9
         # Do a PEP8 test with pycodestyle
         # (allow to fail unless your code completely compliant)
         - os: linux

dust_extinction/dust_extinction.py

@@ -559,34 +559,34 @@ def fit_deriv(in_x, C1, C2, C3, C4, xo, gamma):
         Derivatives of the FM90 function with respect to the parameters
         """
         x = in_x
-        
+
         # useful quantitites
         x2 = x**2
         xo2 = xo**2
         g2 = gamma**2
         x2mxo2_2 = (x2 - xo2)**2
         denom = (x2mxo2_2 - x2*g2)**2
-        
+
         # derivatives
         d_C1 = np.full((len(x)),1.)
         d_C2 = x
 
         d_C3 = (x2/(x2mxo2_2 + x2*g2))
 
         d_xo = (4.*C2*x2*xo*(x2 - xo2))/denom
-        
+
         d_gamma = (2.*C2*(x2**2)*gamma)/denom
-        
+
         d_C4 = np.zeros((len(x)))
         fuv_indxs = np.where(x >= 5.9)
         if len(fuv_indxs) > 0:
             y = x[fuv_indxs] - 5.9
             d_C4[fuv_indxs] = (0.5392*(y**2) + 0.05644*(y**3))
 
         return [d_C1, d_C2, d_C3, d_C4, d_xo, d_gamma]
-    
+
     #fit_deriv = None
-        
+
 class F99(BaseExtRvModel):
     """
     F99 extinction model calculation
@@ -687,7 +687,7 @@ def evaluate(self, in_x, Rv):
         # determine optical/IR values at spline points
         opt_axebv_y = np.array([-0.426 + 1.0044*Rv,
                                 -0.050 + 1.0016*Rv,
-                                0.701 + 1.016*Rv,
+                                0.701 + 1.0016*Rv,
                                 1.208 + 1.0032*Rv - 0.00033*(Rv**2)])
         nir_axebv_y = np.array([0.265,0.829])*Rv/3.1
         optnir_axebv_y = np.concatenate([nir_axebv_y,opt_axebv_y])
@@ -767,6 +767,8 @@ class G03_SMCBar(BaseExtAve):
                              4.243, 4.472, 4.776, 5.000,
                              5.272, 5.575, 5.795, 6.074,
                              6.297, 6.436, 6.992])
+    # accuracy of the observed data based on published table
+    obsdata_tolerance = 6e-2
 
     def evaluate(self, in_x):
         """
@@ -877,6 +879,8 @@ class G03_LMCAvg(BaseExtAve):
                              2.607, 2.668, 2.787, 2.874,
                              2.983, 3.118, 3.231, 3.374,
                              3.366])
+    # accuracy of the observed data based on published table
+    obsdata_tolerance = 6e-2
 
     def evaluate(self, in_x):
         """
@@ -991,6 +995,8 @@ class G03_LMC2(BaseExtAve):
                              3.060, 3.110, 3.299, 3.408,
                              3.515, 3.670, 3.862, 3.937,
                              4.055])
+    # accuracy of the observed data based on published table
+    obsdata_tolerance = 6e-2
 
     def evaluate(self, in_x):
         """

dust_extinction/tests/test_f99.py

@@ -57,41 +57,50 @@ def test_invalid_micron(x_invalid_angstrom):
 
 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])
+    #  from previous version of code
+    #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])
+
+    # correct values (not quite right)
+    #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])
+
+    # from Fitzpatrick (1999) Table 3
+    x = np.array([0.377, 0.820, 1.667, 1.828, 2.141, 2.433,
+                  3.704, 3.846])
+
+    cor_vals = np.array([0.265, 0.829, 2.688, 3.055, 3.806, 4.315,
+                         6.265, 6.591])
+    tolerance = 2e-3
+
+    # convert from A(x)/E(B-V) to A(x)/A(V)
+    cor_vals /= 3.1
 
     # 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)
+    return (x, cor_vals, tolerance)
 
 
 def test_extinction_F99_values():
     # get the correct values
-    x, cor_vals = get_axav_cor_vals()
+    x, cor_vals, tolerance = get_axav_cor_vals()
 
     # initialize extinction model
     tmodel = F99()
 
     # test
-    np.testing.assert_allclose(tmodel(x), cor_vals, rtol=1e-05)
+    np.testing.assert_allclose(tmodel(x), cor_vals, rtol=tolerance)
 
 
-test_vals = zip([0.5, 1.0, 2.0, 2.5, 3.0, 3.5, 4.0,
-                 4.5, 5.0, 6.0, 7.0, 8.0],
-                [0.124997,  0.377073,  1.130636,  1.419644,
-                 1.630377,  1.888546,  2.275900,  3.014577,
-                 2.762256,  2.475272,  2.711508,  3.197144])
+x_vals, axav_vals, tolerance = get_axav_cor_vals()
+test_vals = zip(x_vals, axav_vals, np.full(len(x_vals),tolerance))
 @pytest.mark.parametrize("test_vals", test_vals)
 def test_extinction_F99_single_values(test_vals):
-    x, cor_val = test_vals
+    x, cor_val, tolerance = test_vals
 
     # initialize extinction model
     tmodel = F99()
 
     # test
-    np.testing.assert_allclose(tmodel(x), cor_val, rtol=1e-05)
+    np.testing.assert_allclose(tmodel(x), cor_val, rtol=tolerance)

dust_extinction/tests/test_fm90.py

@@ -91,13 +91,13 @@ def test_FM90_fitting():
     y = (g03_model.obsdata_axav - 1.0)*g03_model.Rv
     # only fit the UV portion (FM90 only valid in UV)
     gindxs, = np.where(x > 3.125)
-    
+
     fm90_init = FM90()
-    fit = LevMarLSQFitter()    
+    fit = LevMarLSQFitter()
     g03_fit = fit(fm90_init, x[gindxs], y[gindxs])
     fit_vals = [g03_fit.C1.value, g03_fit.C2.value, g03_fit.C3.value,
                 g03_fit.C4.value, g03_fit.xo.value, g03_fit.gamma.value]
-    
+
     good_vals = np.array([-0.958016797002, 1.0109751831, 2.96430606652,
                           0.313137860902, 4.59996300532, 0.99000982258])
 

dust_extinction/tests/test_g03.py

@@ -70,13 +70,15 @@ def test_extinction_G03_values(tmodel):
     #  not to numerical precision as we are using the FM90 fits
     #  and spline functions and the correct values are the data
     np.testing.assert_allclose(tmodel(tmodel.obsdata_x),
-                               tmodel.obsdata_axav, rtol=6e-02)
+                               tmodel.obsdata_axav,
+                               rtol=tmodel.obsdata_tolerance)
 
 @pytest.mark.parametrize("tmodel", models)
 def test_extinction_G03_single_values(tmodel):
     # test
     for x, cor_val in zip(tmodel.obsdata_x, tmodel.obsdata_axav):
-        np.testing.assert_allclose(tmodel(x), cor_val, rtol=6e-02)
+        np.testing.assert_allclose(tmodel(x), cor_val,
+                                   rtol=tmodel.obsdata_tolerance)
 
 @pytest.mark.parametrize("tmodel", models)
 def test_extinction_G03_extinguish_values_Av(tmodel):

dust_extinction/tests/test_g16.py

@@ -5,6 +5,8 @@
 from astropy.modeling import InputParameterError
 
 from ..dust_extinction import G16
+from ..dust_extinction import G03_SMCBar
+from .test_f99 import get_axav_cor_vals as get_axav_cor_vals_fA_1
 
 @pytest.mark.parametrize("RvA_invalid", [-1.0,0.0,1.9,6.1,10.])
 def test_invalid_RvA_input(RvA_invalid):
@@ -67,43 +69,38 @@ def test_invalid_micron(x_invalid_angstrom):
                                 + 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
+def test_extinction_G16_fA_1_values():
+    # get the correct values
+    x, cor_vals, tolerance = get_axav_cor_vals_fA_1()
 
-    # 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])
+    # initialize extinction model
+    tmodel = G16(RvA=3.1, fA=1.0)
 
-    return (x, cor_vals)
+    # test
+    np.testing.assert_allclose(tmodel(x), cor_vals, rtol=tolerance)
 
+def test_extinction_G16_fA_0_values():
+    # initialize the model
+    tmodel = G16(fA=0.0)
 
-def test_extinction_G16_values():
     # get the correct values
-    x, cor_vals = get_axav_cor_vals()
-
-    # initialize extinction model
-    tmodel = G16()
+    gmodel = G03_SMCBar()
+    x = gmodel.obsdata_x
+    cor_vals = gmodel.obsdata_axav
+    tolerance = gmodel.obsdata_tolerance
 
     # test
-    np.testing.assert_allclose(tmodel(x), cor_vals, rtol=1e-05)
+    np.testing.assert_allclose(tmodel(x), cor_vals, rtol=tolerance)
 
-
-test_vals = zip([0.5, 1.0, 2.0, 2.5, 3.0, 3.5, 4.0,
-                 4.5, 5.0, 6.0, 7.0, 8.0],
-                [0.124997,  0.377073,  1.130636,  1.419644,
-                 1.630377,  1.888546,  2.275900,  3.014577,
-                 2.762256,  2.475272,  2.711508,  3.197144])
+x_vals, axav_vals, tolerance = get_axav_cor_vals_fA_1()
+test_vals = zip(x_vals, axav_vals, np.full(len(x_vals),tolerance))
 @pytest.mark.parametrize("test_vals", test_vals)
-def test_extinction_G16_single_values(test_vals):
-    x, cor_val = test_vals
+def test_extinction_G16_fA_1_single_values(test_vals):
+    x, cor_val, tolerance = test_vals
 
     # initialize extinction model
-    tmodel = G16()
+    tmodel = G16(RvA=3.1, fA=1.0)
 
     # test
-    np.testing.assert_allclose(tmodel(x), cor_val, rtol=1e-05)
+    np.testing.assert_allclose(tmodel(x), cor_val, rtol=tolerance)