Merge pull request #118 from karllark/update_calc_av

update A(V) calc to use latest - Decleir et al. (2022)

measure_extinction/extdata.py

@@ -502,17 +502,18 @@ def calc_EBV(self):
         else:
             self.columns["EBV"] = (self.exts["BAND"][bindx], self.uncs["BAND"][bindx])
 
-    def calc_AV(self, akav=0.112):
+    def calc_AV(self, akav=0.105):
         """
         Calculate A(V) from the observed extinction curve:
             - fit a powerlaw to the SpeX extinction curve, if available
             - otherwise: extrapolate the K-band extinction
 
         Parameters
         ----------
-        akav : float  [default = 0.112]
+        akav : float
            Value of A(K)/A(V)
-           default is from Rieke & Lebofsky (1985)
+           default is from Decleir et al. (2022)
+           for Rieke & Lebofsky (1985) use akav=0.112
            van de Hulst No. 15 curve has A(K)/A(V) = 0.0885
 
         Returns
@@ -539,23 +540,33 @@ def calc_AV(self, akav=0.112):
                 avunc = np.absolute(av * (self.uncs["BAND"][kindx] / ekv))
                 self.columns["AV"] = (av, avunc)
 
-    def calc_AV_JHK(self):
+    def calc_AV_JHK(
+        self,
+        ref_waves=np.array([1.25, 1.6, 2.2]) * u.micron,
+        ref_alav=[0.269, 0.163, 0.105],
+    ):
         """
         Calculate A(V) from the observed extinction curve:
             - extrapolate from J, H, & K photometry
-            - assumes functional form from Rieke, Rieke, & Paul (1989)
+
+        Parameters
+        ----------
+        ref_waves : floats
+            wavelengths for reference values (default = JHK)
+        ref_alav : floats
+            A(lambda)/A(V) values for reference
+            default is for JHK from Decleir et al. (2022)
+            for Rieke, Rieke, & Paul (1989) use ref_alav=[0.2815534, 0.17475728, 0.11197411],
 
         Returns
         -------
         Updates self.columns["AV"]
         """
         # J, H, K
-        rrp89_waves = np.array([1.25, 1.6, 2.2]) * u.micron
-        rrp89_alav = [0.2815534, 0.17475728, 0.11197411]
 
         avs = []
         avs_unc = []
-        for cwave, calav in zip(rrp89_waves, rrp89_alav):
+        for cwave, calav in zip(ref_waves, ref_alav):
             dwaves = np.absolute(self.waves["BAND"] - cwave)
             kindx = dwaves.argmin()
             if dwaves[kindx] < 0.1 * u.micron:

measure_extinction/tests/test_EbvAvRv.py

@@ -25,8 +25,8 @@ def test_calc_av():
     text.calc_AV()
     av = text.columns["AV"]
     assert isinstance(av, tuple)
-    assert np.isclose(av[0], 0.56306, atol=1e-3)
-    assert np.isclose(av[1], 0.02252, atol=1e-3)
+    assert np.isclose(av[0], 0.55866, atol=1e-3)
+    assert np.isclose(av[1], 0.02235, atol=1e-3)
 
 
 def test_calc_rv_fromelv():
@@ -38,8 +38,8 @@ def test_calc_rv_fromelv():
     text.calc_RV()
     rv = text.columns["RV"]
     assert isinstance(rv, tuple)
-    assert np.isclose(rv[0], 3.37837, atol=1e-3)
-    assert np.isclose(rv[1], 0.207647, atol=1e-3)
+    assert np.isclose(rv[0], 3.351955, atol=1e-3)
+    assert np.isclose(rv[1], 0.206023, atol=1e-3)
 
 
 def test_calc_rv_fromebvav():