Repository: skypyproject/skypy
Pull request #208 (2020-07-03 09:13:48)

absolute magnitude standard band zeropoints
## Description

Adds an optional `band` parameter to functions working with absolute magnitudes and luminosities. If a standard filter band is given, luminosities are returned in units of solar luminosity for that band.

## Checklist
- [x] Follow the [Contributor Guidelines](https://github.com/skypyproject/skypy/blob/master/CONTRIBUTING.md)
- [x] Write unit tests
- [x] Write documentation strings
- [x] Assign someone from your working team to review this pull request
- [x] Assign someone from the infrastructure team to review this pull request

skypybot: Hello @ntessore! Thanks for opening this PR. We checked the lines you've touched for [PEP 8](https://www.python.org/dev/peps/pep-0008) issues, and found:





* In the file [`skypy/galaxy/tests/test_luminosity.py`](https://github.com/skypyproject/skypy/blob/4749ae9f83b359df801f9c3fee84eaf09c1f3160/skypy/galaxy/tests/test_luminosity.py):

> [Line 29:17](https://github.com/skypyproject/skypy/blob/4749ae9f83b359df801f9c3fee84eaf09c1f3160/skypy/galaxy/tests/test_luminosity.py#L29): [E127](https://duckduckgo.com/?q=pep8%20E127) continuation line over-indented for visual indent





Patch:
diff --git a/skypy/galaxy/luminosity.py b/skypy/galaxy/luminosity.py
index 3c0307a..74eb8d5 100644
--- a/skypy/galaxy/luminosity.py
+++ b/skypy/galaxy/luminosity.py
@@ -44,7 +44,8 @@ def herbel_luminosities(redshift, alpha, a_m, b_m, size=None,
     Returns
     -------
     luminosity : array_like
-        Drawn luminosities from the Schechter luminosity function.
+        Drawn luminosities from the Schechter luminosity function. Luminosities
+        are B-band luminosities in units of solar luminosity.
 
     Notes
     -----
@@ -98,20 +99,21 @@ def herbel_luminosities(redshift, alpha, a_m, b_m, size=None,
     >>> luminosities = lum.herbel_luminosities(z, -1.3, -0.9408582,
     ...                                         -20.40492365)
 
-
-
     """
 
     if size is None and np.shape(redshift):
         size = np.shape(redshift)
 
-    luminosity_star = _calculate_luminosity_star(redshift, a_m, b_m)
+    # this is the AB zeropoint of the B-band in units of solar luminosity
+    zeropoint = -4.73
+
+    luminosity_star = _calculate_luminosity_star(redshift, a_m, b_m, zeropoint)
 
     x_sample = schechter(alpha, x_min, x_max, resolution=resolution, size=size)
 
     return luminosity_star * x_sample
 
 
-def _calculate_luminosity_star(redshift, a_m, b_m):
+def _calculate_luminosity_star(redshift, a_m, b_m, zeropoint=None):
     absolute_magnitude_star = a_m * redshift + b_m
-    return astro.luminosity_from_absolute_magnitude(absolute_magnitude_star)
+    return astro.luminosity_from_absolute_magnitude(absolute_magnitude_star, zeropoint)
diff --git a/skypy/galaxy/redshift.py b/skypy/galaxy/redshift.py
index 8da23c6..245a501 100644
--- a/skypy/galaxy/redshift.py
+++ b/skypy/galaxy/redshift.py
@@ -9,7 +9,6 @@ import scipy.special as sc
 import scipy.integrate
 
 import skypy.utils.special as special
-import skypy.utils.astronomy as astro
 
 
 __all__ = [
@@ -159,11 +158,9 @@ def herbel_redshift(alpha, a_phi, b_phi, a_m, b_m, cosmology, low=0.0,
 
     """
 
-    luminosity_min = astro.luminosity_from_absolute_magnitude(
-        absolute_magnitude_max)
     redshift = np.linspace(low, high, resolution)
     pdf = herbel_pdf(redshift, alpha, a_phi, b_phi, a_m, b_m, cosmology,
-                     luminosity_min)
+                     absolute_magnitude_max)
     cdf = scipy.integrate.cumtrapz(pdf, redshift, initial=0)
     cdf = cdf / cdf[-1]
     u = np.random.uniform(size=size)
@@ -173,7 +170,7 @@ def herbel_redshift(alpha, a_phi, b_phi, a_m, b_m, cosmology, low=0.0,
 
 
 def herbel_pdf(redshift, alpha, a_phi, b_phi, a_m, b_m, cosmology,
-               luminosity_min):
+               absolute_magnitude_max):
     r"""Calculates the redshift pdf of the Schechter luminosity function
     according to the model of Herbel et al. [1]_ equation (3.6).
 
@@ -226,7 +223,6 @@ def herbel_pdf(redshift, alpha, a_phi, b_phi, a_m, b_m, cosmology,
     Examples
     --------
     >>> from skypy.galaxy.redshift import herbel_pdf
-    >>> import skypy.utils.astronomy as astro
     >>> from astropy.cosmology import FlatLambdaCDM
     >>> import numpy as np
 
@@ -236,15 +232,15 @@ def herbel_pdf(redshift, alpha, a_phi, b_phi, a_m, b_m, cosmology,
 
     >>> cosmology = FlatLambdaCDM(H0=70, Om0=0.3, Tcmb0=2.725)
     >>> redshift = np.linspace(0, 2, 100)
-    >>> luminosity_min = astro.luminosity_from_absolute_magnitude(-16.0)
+    >>> mag_lim = -16.0
     >>> redshift = herbel_pdf(redshift=redshift, alpha=-1.3,
     ...                     a_phi=-0.10268436,a_m=-0.9408582, b_phi=0.00370253,
     ...                     b_m=-20.40492365, cosmology=cosmology,
-    ...                     luminosity_min=luminosity_min)
+    ...                     absolute_magnitude_max=mag_lim)
     """
+    abs_mag = a_m*redshift + b_m
     dv = cosmology.differential_comoving_volume(redshift).value
-    x = luminosity_min \
-        * 1./astro.luminosity_from_absolute_magnitude(a_m*redshift + b_m)
+    x = 10.**(-0.4*(absolute_magnitude_max - abs_mag))
     lg = sc.gammaln(alpha+1)
     gx = np.fabs(special.gammaincc(alpha+1, x))
     pdf = dv * b_phi * np.exp(a_phi * redshift + lg) * gx
diff --git a/skypy/utils/astronomy.py b/skypy/utils/astronomy.py
index 1cddff2..a863d52 100644
--- a/skypy/utils/astronomy.py
+++ b/skypy/utils/astronomy.py
@@ -18,33 +18,45 @@ Utility functions
 import numpy as np
 
 
-def luminosity_from_absolute_magnitude(absolute_magnitude):
+def luminosity_from_absolute_magnitude(absolute_magnitude, zeropoint=None):
     """ Converts absolute magnitudes into luminosities
 
     Parameters
     ----------
     absolute_magnitude : array_like
-                    Input absolute magnitudes
+        Input absolute magnitudes
+    zeropoint : float, optional
+        Zeropoint for the conversion.
+
     Returns
     -------
     ndarray, or float if input is scalar
     Luminosity values.
     """
 
-    return np.power(10, -0.4*absolute_magnitude)
+    if zeropoint is None:
+        zeropoint = 0.
+
+    return 10.**(-0.4*np.add(absolute_magnitude, zeropoint))
 
 
-def absolute_magnitude_from_luminosity(luminosity):
+def absolute_magnitude_from_luminosity(luminosity, zeropoint=None):
     """ Converts luminosities into absolute magnitudes
 
     Parameters
     ----------
     luminosity : array_like
-            Input luminosity
+        Input luminosity
+    zeropoint : float, optional
+        Zeropoint for the conversion.
 
     Returns
     -------
     ndarray, or float if input is scalar
-    Absolute magnitude values
+    Absolute magnitude values.
     """
-    return -np.log(luminosity)/(0.4 * np.log(10))
+
+    if zeropoint is None:
+        zeropoint = 0.
+
+    return -2.5*np.log10(luminosity) - zeropoint


Test patch:
diff --git a/skypy/galaxy/tests/test_luminosity.py b/skypy/galaxy/tests/test_luminosity.py
index 66711c1..3205283 100644
--- a/skypy/galaxy/tests/test_luminosity.py
+++ b/skypy/galaxy/tests/test_luminosity.py
@@ -29,14 +29,15 @@ def test_herbel_luminosities():
     # Test that sampling corresponds to sampling from the right pdf.
     # For this, we sample an array of luminosities for redshift z = 1.0 and we
     # compare it to the corresponding cdf.
+    # The B-band zeropoint for the model is -4.73.
 
     def calc_pdf(luminosity, z, alpha, a_m, b_m, absolute_magnitude_max=-16.0,
                  absolute_magnitude_min=-28.0):
         luminosity_min = astro.luminosity_from_absolute_magnitude(
-            absolute_magnitude_max)
+            absolute_magnitude_max, -4.73)
         luminosity_max = astro.luminosity_from_absolute_magnitude(
-            absolute_magnitude_min)
-        lum_star = lum._calculate_luminosity_star(z, a_m, b_m)
+            absolute_magnitude_min, -4.73)
+        lum_star = lum._calculate_luminosity_star(z, a_m, b_m, -4.73)
         c = np.fabs(special.gammaincc(alpha+1, luminosity_min/lum_star))
         d = np.fabs(special.gammaincc(alpha+1, luminosity_max/lum_star))
         return 1./lum_star*np.power(luminosity/lum_star, alpha)*np.exp(
@@ -59,8 +60,9 @@ def test_herbel_luminosities():
 
 def test_exponential_distribution():
     # When alpha=0, L*=1 and x_min~0 we get a truncated exponential
+    # the b_m value offsets the B-band zeropoint of -4.73
     q_max = 1e2
-    sample = lum.herbel_luminosities(0, 0, 0, 0, size=1000,
+    sample = lum.herbel_luminosities(0, 0, 0, 4.73, size=1000,
                                      x_min=1e-10, x_max=q_max,
                                      resolution=1000)
     d, p_value = scipy.stats.kstest(sample, 'truncexpon', args=(q_max,))
diff --git a/skypy/galaxy/tests/test_redshift.py b/skypy/galaxy/tests/test_redshift.py
index 894650d..2ac184d 100644
--- a/skypy/galaxy/tests/test_redshift.py
+++ b/skypy/galaxy/tests/test_redshift.py
@@ -39,7 +39,7 @@ def test_herbel_pdf():
     pdf = herbel_pdf(np.array([0.01, 0.5, 1, 2]),
                      -0.5, -0.70596888,
                      0.0035097, -0.70798041,
-                     -20.37196157, cosmology, np.power(10, -0.4 * -16.0))
+                     -20.37196157, cosmology, -16.0)
     result = np.array(
         [4.09063927e+04, 4.45083420e+07, 7.26629445e+07, 5.40766813e+07])
     np.testing.assert_allclose(pdf, result)
@@ -80,7 +80,7 @@ def test_herbel_redshift_sampling():
                          b_phi=0.00370253,
                          b_m=-20.40492365,
                          cosmology=cosmology,
-                         luminosity_min=2511886.4315095823)
+                         absolute_magnitude_max=-16)
         cdf = scipy.integrate.cumtrapz(pdf, z, initial=0)
         cdf = cdf / cdf[-1]
         return cdf