Repository: skypyproject/skypy
Pull request #242 (2020-07-22 09:20:46)

compute AB magnitude from spectrum and bandpass
## Description

This PR contains the functionality for computing AB magnitudes from spectra and filters.

## 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/spectrum.py`](https://github.com/skypyproject/skypy/blob/7fb2b050bf03387e1a9e8c3244d35ed80f03c225/skypy/galaxy/spectrum.py):

> [Line 235:80](https://github.com/skypyproject/skypy/blob/7fb2b050bf03387e1a9e8c3244d35ed80f03c225/skypy/galaxy/spectrum.py#L235): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (89 > 79 characters)
> [Line 239:80](https://github.com/skypyproject/skypy/blob/7fb2b050bf03387e1a9e8c3244d35ed80f03c225/skypy/galaxy/spectrum.py#L239): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (89 > 79 characters)


ntessore: @rrjbca This is a draft for now, need to discuss the mechanism for getting spectra and bandpasses, as well as assembling, storing, and retrieving the redshift tracks if we want to use that mechanism.

Patch:
diff --git a/skypy/galaxy/spectrum.py b/skypy/galaxy/spectrum.py
index 24c22bc..765bee4 100644
--- a/skypy/galaxy/spectrum.py
+++ b/skypy/galaxy/spectrum.py
@@ -8,7 +8,8 @@ from astropy.io.fits import getdata
 
 __all__ = [
     'dirichlet_coefficients',
-    'kcorrect_spectra'
+    'kcorrect_spectra',
+    'mag_ab',
 ]
 
 
@@ -195,3 +196,87 @@ def kcorrect_spectra(redshift, stellar_mass, coefficients):
                                     wavelength.reshape(1, len(wavelength)))
 
     return wavelength_observed, sed
+
+
+def mag_ab(spec_lam, spec_flux, band_lam, band_tx, redshift=None):
+    r'''Compute absolute AB magnitude from spectrum and bandpass.
+
+    This function takes an *emission* spectrum and an observation bandpass and
+    computes the AB magnitude for a source at 10pc (i.e. absolute magnitude).
+    The emission spectrum can optionally be redshifted. The definition of the
+    bandpass AB magnitude is taken from [1]_.
+
+    Both the spectrum and the bandpass must be given as functions of wavelength
+    in Angstrom. The spectrum must be given as flux in erg/s/cm2/A.
+
+    Parameters
+    ----------
+    spec_lam : (ns,) array_like
+        Vector of spectrum wavelengths in units of Angstrom.
+    spec_flux : (ns,) array_like
+        Vector of spectrum fluxes in units of erg/s/cm2/A.
+    band_lam : (nb,) array_like
+        Vector of bandpass wavelengths in units of Angstrom.
+    band_tx : (nb,) array_like
+        Vector of bandpass transmissions.
+    redshift : array_like, optional
+        Optional array of values for redshifting the source spectrum. Default
+        is a redshift of 0.0.
+
+    Returns
+    -------
+    mag_ab : array_like
+        The absolute AB magnitude. If redshifts are given, the output has the
+        same shape *and type* as the `redshift` argument.
+
+    References
+    ----------
+    .. [1] M. R. Blanton et al., 2003, AJ, 125, 2348
+    '''
+
+    assert np.ndim(spec_lam) == 1, 'spec_lam must be 1d array'
+    assert np.ndim(spec_flux) == 1, 'spec_flux must be 1d array'
+    assert len(spec_lam) == len(spec_flux), 'spec_lam and spec_flux must match'
+    assert np.all(np.less_equal(spec_lam[:-1], spec_lam[1:])), 'spec_lam must be ordered'
+    assert np.ndim(band_lam) == 1, 'band_lam must be 1d array'
+    assert np.ndim(band_tx) == 1, 'band_tx must be 1d array'
+    assert len(band_lam) == len(band_tx), 'band_lam and band_tx must match'
+    assert np.all(np.less_equal(band_lam[:-1], band_lam[1:])), 'band_lam must be ordered'
+
+    # redshift zero if not given
+    if redshift is None:
+        redshift = 0.
+
+    # allocate output array
+    mag_ab = np.empty_like(redshift)
+
+    # compute magnitude contribution from band normalisation [denominator of (2)]
+    m_band = -2.5*np.log10(np.trapz(band_tx/band_lam, band_lam))
+
+    # magnitude offset from band and AB definition [constant in (2)]
+    m_offs = -2.4079482426801846 - m_band
+
+    # compute flux integrand at emitted wavelengths
+    spec_intg = spec_lam*spec_flux
+
+    # go through redshifts ...
+    for i, z in np.ndenumerate(redshift):
+
+        # observed wavelength of spectrum
+        obs_lam = (1 + z)*spec_lam
+
+        # interpolate band to get transmission at observed wavelengths
+        obs_tx = np.interp(obs_lam, band_lam, band_tx, left=0, right=0)
+
+        # compute magnitude contribution from flux [numerator of (2)]
+        m_flux = -2.5*np.log10(np.trapz(spec_intg*obs_tx, obs_lam))
+
+        # combine AB magnitude [all of (2)]
+        mag_ab[i] = m_flux + m_offs
+
+    # simplify scalar output
+    if np.isscalar(redshift):
+        mag_ab = mag_ab.item()
+
+    # all done
+    return mag_ab


Test patch:
diff --git a/skypy/galaxy/tests/test_spectrum.py b/skypy/galaxy/tests/test_spectrum.py
index 64dfcaa..a964690 100644
--- a/skypy/galaxy/tests/test_spectrum.py
+++ b/skypy/galaxy/tests/test_spectrum.py
@@ -83,3 +83,41 @@ def test_kcorrect_spectra():
 
     assert np.allclose(lam_o, lam)
     assert np.allclose(sed, templates[0])
+
+
+def test_mag_ab_standard_source():
+
+    from skypy.galaxy.spectrum import mag_ab
+
+    # create a bandpass
+    bp_lam = np.logspace(0, 4, 1000)
+    bp_tx = np.exp(-((bp_lam - 1000)/(100))**2)
+
+    # test that the AB standard source has zero magnitude
+    lam = np.logspace(0, 4, 1000)
+    flam = 0.10884806248538730623/lam**2
+    m = mag_ab(lam, flam, bp_lam, bp_tx)
+
+    assert np.isclose(m, 0)
+
+
+def test_mag_ab_redshift_dependence():
+
+    from skypy.galaxy.spectrum import mag_ab
+
+    # make a wide tophat bandpass
+    th_lam = np.logspace(-10, 10, 3)
+    th_tx = np.ones(3)
+
+    # create a narrow gaussian source
+    lam = np.logspace(0, 3, 1000)
+    flam = np.exp(-((lam - 100)/10)**2)
+
+    # array of redshifts
+    z = np.linspace(0, 1, 11)
+
+    # compute the AB magnitude at different redshifts
+    m = mag_ab(lam, flam, th_lam, th_tx, z)
+
+    # compare with expected redshift dependence
+    np.testing.assert_allclose(m, m[0] - 2.5*np.log10(1 + z))