Merge b060b98 into 52bd24d

desihub · Nov 29, 2018 · a51ff75 · a51ff75
2 parents 52bd24d + b060b98
commit a51ff75
Show file tree

Hide file tree

Showing 3 changed files with 144 additions and 20 deletions.
diff --git a/doc/changes.rst b/doc/changes.rst
@@ -10,10 +10,13 @@ Change Log
 1.9.15 (unreleased)
 -------------------
 
+* Add :func:`desiutil.dust.ext_odonnell` and :func:`desiutil.dust.ext_ccm`
+  originally from desispec (PR `#128`_).
 * Update permissions set by :command:`fix_permissions.sh` (PR `#126`_).
 * Set read-only permissions on all Module files, and unlock them as needed (PR `#125`_).
 * Draw ecliptic in all-sky plots (PR `#124`_).
 
+.. _`#128`: https://github.com/desihub/desiutil/pull/128
 .. _`#126`: https://github.com/desihub/desiutil/pull/126
 .. _`#125`: https://github.com/desihub/desiutil/pull/125
 .. _`#124`: https://github.com/desihub/desiutil/pull/124

diff --git a/py/desiutil/dust.py b/py/desiutil/dust.py
@@ -1,7 +1,123 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 # -*- coding: utf-8 -*-
 
-# Some code was copied on Nov/20/2016 from https://github.com/kbarbary/sfdmap/ commit: bacdbbd
+"""
+=============
+desiutil.dust
+=============
+
+Get :math:`E(B-V)` values from the `Schlegel, Finkbeiner & Davis (1998; SFD98)`_ dust map.
+
+.. _`Schlegel, Finkbeiner & Davis (1998; SFD98)`: http://adsabs.harvard.edu/abs/1998ApJ...500..525S.
+"""
+
+import os
+import numpy as np
+from astropy.io.fits import getdata
+from astropy.coordinates import SkyCoord
+from astropy import units as u
+
+from .log import get_logger
+log = get_logger()
+
+def ext_odonnell(wave,Rv=3.1) :
+    """Return extinction curve from Odonnell (1994), defined in the wavelength
+    range [3030,9091] Angstroms.  Outside this range, use CCM (1989).
+
+    Args:
+        wave : 1D array of vacuum wavelength [Angstroms]
+        Rv   : Value of R_V (scalar); default is 3.1
+
+     Returns:
+        1D array of A(lambda)/A(V)
+    """
+
+    # direct python translation of idlutils/pro/dust/ext_odonnell.pro
+
+    A  = np.zeros(wave.shape)
+    xx = 10000. / wave
+
+    optical_waves = (xx>=1.1)&(xx<=3.3)
+    other_waves   = (xx<1.1)|(xx>3.3)
+
+    if np.sum(optical_waves)>0 :
+        yy = xx[optical_waves]-1.82
+        afac = 1.0 + 0.104*yy - 0.609*yy**2 + 0.701*yy**3 + 1.137*yy**4 - 1.718*yy**5 - 0.827*yy**6 + 1.647*yy**7 - 0.505*yy**8
+        bfac = 1.952*yy + 2.908*yy**2 - 3.989*yy**3 - 7.985*yy**4 + 11.102*yy**5 + 5.491*yy**6 - 10.805*yy**7 + 3.347*yy**8
+        A[optical_waves] = afac + bfac / Rv
+    if np.sum(other_waves)>0 :
+        A[other_waves]   =  ext_ccm(wave[other_waves],Rv=Rv)
+
+    return A
+
+def ext_ccm(wave,Rv=3.1) :
+    """Return extinction curve from CCM (1989), defined in the wavelength
+    range [1250,33333] Angstroms.
+
+    Args:
+        wave : 1D array of vacuum wavelength [Angstroms]
+        Rv   : Value of R_V (scalar); default is 3.1
+
+     Returns:
+        1D array of A(lambda)/A(V)
+    """
+
+    # direct python translation of idlutils/pro/dust/ext_ccm.pro
+    # numeric values checked with other implementation
+
+    A  = np.zeros(wave.shape)
+    xx = 10000. / wave
+
+    # Limits for CCM fitting function
+    qLO  = (xx > 8.0)                # No data, lambda < 1250 Ang
+    qUV  = (xx > 3.3)&(xx <= 8.0)    # UV + FUV
+    qOPT = (xx > 1.1)&(xx <= 3.3)    #  Optical/NIR
+    qIR  = (xx > 0.3)&(xx <= 1.1)    # IR
+    qHI  = (xx <= 0.3)               # No data, lambda > 33,333 Ang
+
+    # For lambda < 1250 Ang, arbitrarily return Alam=5
+    if np.sum(qLO) > 0 :
+        A[qLO] = 5.0
+
+    if np.sum(qUV) > 0 :
+        xt = xx[qUV]
+        afac = 1.752 - 0.316*xt - 0.104 / ( (xt-4.67)**2 + 0.341 )
+        bfac = -3.090 + 1.825*xt + 1.206 / ( (xt-4.62)**2 + 0.263 )
+
+        qq = (xt >= 5.9)&(xt <= 8.0)
+        if np.sum(qq)> 0 :
+            Fa = -0.04473*(xt[qq]-5.9)**2 - 0.009779*(xt[qq]-5.9)**3
+            Fb = 0.2130*(xt[qq]-5.9)**2 + 0.1207*(xt[qq]-5.9)**3
+            afac[qq] += Fa
+            bfac[qq] += Fb
+
+        A[qUV] = afac + bfac / Rv
+
+    if np.sum(qOPT) > 0 :
+        yy = xx[qOPT] - 1.82
+        afac = 1.0 + 0.17699*yy - 0.50447*yy**2 - 0.02427*yy**3 + 0.72085*yy**4 + 0.01979*yy**5 - 0.77530*yy**6 + 0.32999*yy**7
+        bfac = 1.41338*yy + 2.28305*yy**2 + 1.07233*yy**3 - 5.38434*yy**4 - 0.62251*yy**5 + 5.30260*yy**6 - 2.09002*yy**7
+        A[qOPT] = afac + bfac / Rv
+
+    if np.sum(qIR) > 0 :
+        yy = xx[qIR]**1.61
+        afac = 0.574*yy
+        bfac = -0.527*yy
+        A[qIR] = afac + bfac / Rv
+
+    # For lambda > 33,333 Ang, arbitrarily extrapolate the IR curve
+    if np.sum(qHI) > 0 :
+        yy = xx[qHI]**1.61
+        afac = 0.574*yy
+        bfac = -0.527*yy
+        A[qHI] = afac + bfac / Rv
+
+    return A
+
+#-------------------------------------------------------------------------
+# The SFDMap and _Hemisphere classes and the _bilinear_interpolate and ebv
+# functions below were copied on Nov/20/2016 from
+# https://github.com/kbarbary/sfdmap/ commit: bacdbbd
 # which was originally Licensed under an MIT "Expat" license:
 
 # Copyright (c) 2016 Kyle Barbary
@@ -24,25 +140,6 @@
 # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 # SOFTWARE.
 
-"""
-=============
-desiutil.dust
-=============
-
-Get :math:`E(B-V)` values from the `Schlegel, Finkbeiner & Davis (1998; SFD98)`_ dust map.
-
-.. _`Schlegel, Finkbeiner & Davis (1998; SFD98)`: http://adsabs.harvard.edu/abs/1998ApJ...500..525S.
-"""
-
-import os
-import numpy as np
-from astropy.io.fits import getdata
-from astropy.coordinates import SkyCoord
-from astropy import units as u
-
-from .log import get_logger
-log = get_logger()
-
 def _bilinear_interpolate(data, y, x):
     """Map a two-dimensional integer pixel-array at float coordinates.
 

diff --git a/py/desiutil/test/test_dust.py b/py/desiutil/test/test_dust.py
@@ -101,6 +101,30 @@ def test_class(self):
         # ADM assert that the tests worked 
         self.assertTrue(testcnt == 2)
 
+    def test_extinction(self):
+        """Test ext_odonnel and ext_ccm functions"""
+        wave = np.arange(2000,10001,100)
+        ext_odl_31 = dust.ext_odonnell(wave, Rv=3.1)
+        ext_odl_33 = dust.ext_odonnell(wave, Rv=3.3)
+        ext_ccm_31 = dust.ext_ccm(wave, Rv=3.1)
+        ext_ccm_33 = dust.ext_ccm(wave, Rv=3.3)
+
+        #- Sanity check on ranges
+        self.assertTrue(np.all(0.4<ext_odl_31) and np.all(ext_odl_31<4.0))
+        self.assertTrue(np.all(0.4<ext_odl_33) and np.all(ext_odl_33<4.0))
+        self.assertTrue(np.all(0.4<ext_ccm_31) and np.all(ext_ccm_31<4.0))
+        self.assertTrue(np.all(0.4<ext_ccm_33) and np.all(ext_ccm_33<4.0))
+
+        #- Changing Rv should change answer
+        self.assertTrue(np.all(ext_odl_31 != ext_odl_33))
+        self.assertTrue(np.all(ext_ccm_31 != ext_ccm_33))
+
+        #- Odonnell == CCM for some but not all wavelengths
+        self.assertTrue(np.any(ext_odl_31 == ext_ccm_31))
+        self.assertTrue(np.any(ext_odl_31 != ext_ccm_31))
+        self.assertTrue(np.any(ext_odl_33 == ext_ccm_33))
+        self.assertTrue(np.any(ext_odl_33 != ext_ccm_33))
+
 def test_suite():
     """Allows testing of only this module with the command::