Factor out generalised template magnitude functions

skypyproject · rrjbca · Jan 19, 2021 · Dec 1, 2020 · Dec 1, 2020 · Dec 1, 2020
commit 578c06047436218b0b7bee5ac6899556fa945c81
diff --git a/skypy/galaxy/spectrum.py b/skypy/galaxy/spectrum.py
@@ -13,6 +13,8 @@
     'dirichlet_coefficients',
     'load_spectral_data',
     'ab_maggies_redshift',
+    'template_absolute_magnitudes',
+    'template_apparent_magnitudes',
     'kcorrect_absolute_magnitudes',
     'kcorrect_apparent_magnitudes',
     'kcorrect_stellar_mass',
@@ -219,6 +221,26 @@ def ab_maggies_redshift(wavelength, spectrum, filters, redshift, interpolate=100
     return m * np.reshape((1 + redshift), np.shape(redshift) + (1,)*(nd_s+nd_f))
 
 
+def template_absolute_magnitudes(lambda_, spec, coefficients, filters, stellar_mass=None):
+
+    flux = ab_maggies_redshift(lambda_, spec, filters, 0)
+    flux = np.sum((coefficients.T * flux.T).T, axis=1)
+    mass_modulus = 0 if stellar_mass is None else -2.5 * np.log10(stellar_mass)
+
+    return (-2.5*np.log10(flux).T + mass_modulus).T
+
+
+def template_apparent_magnitudes(lambda_, spec, coefficients, redshift, filters, cosmology,
+                                 *, stellar_mass=None, resolution=1000):
+
+    flux = ab_maggies_redshift(lambda_, spec, filters, redshift, resolution)
+    flux = np.sum((coefficients.T * flux.T).T, axis=1)
+    distance_modulus = cosmology.distmod(redshift).value
+    mass_modulus = 0 if stellar_mass is None else -2.5 * np.log10(stellar_mass)
+
+    return (-2.5*np.log10(flux).T + distance_modulus + mass_modulus).T
+
+
 def kcorrect_absolute_magnitudes(coefficients, filters, stellar_mass=None):
     '''Galaxy AB absolute magnitudes from kcorrect template spectra.
 
@@ -253,11 +275,7 @@ def kcorrect_absolute_magnitudes(coefficients, filters, stellar_mass=None):
         spec = hdul[1].data * units.Unit('erg s-1 cm-2 angstrom-1')
         lambda_ = hdul[11].data * units.Unit('angstrom')
 
-    flux = ab_maggies_redshift(lambda_, spec, filters, 0)
-    flux = np.sum((coefficients.T * flux.T).T, axis=1)
-    mass_modulus = 0 if stellar_mass is None else -2.5 * np.log10(stellar_mass)
-
-    return (-2.5*np.log10(flux).T + mass_modulus).T
+    return template_absolute_magnitudes(lambda_, spec, coefficients, filters, stellar_mass)
 
 
 def kcorrect_apparent_magnitudes(coefficients, redshift, filters, cosmology,
@@ -304,12 +322,8 @@ def kcorrect_apparent_magnitudes(coefficients, redshift, filters, cosmology,
         spec = hdul[1].data * units.Unit('erg s-1 cm-2 angstrom-1')
         lambda_ = hdul[11].data * units.Unit('angstrom')
 
-    flux = ab_maggies_redshift(lambda_, spec, filters, redshift, resolution)
-    flux = np.sum((coefficients.T * flux.T).T, axis=1)
-    distance_modulus = cosmology.distmod(redshift).value
-    mass_modulus = 0 if stellar_mass is None else -2.5 * np.log10(stellar_mass)
-
-    return (-2.5*np.log10(flux).T + distance_modulus + mass_modulus).T
+    return template_apparent_magnitudes(lambda_, spec, coefficients, redshift, filters,
+                                        cosmology, stellar_mass, resolution)
 
 
 def kcorrect_stellar_mass(coefficients, magnitudes, filter):

diff --git a/skypy/galaxy/tests/test_spectrum.py b/skypy/galaxy/tests/test_spectrum.py
@@ -156,47 +156,47 @@ def test_ab_maggies_redshift_multi():
 def test_template_spectra():
 
     from astropy import units
-    from skypy.galaxy.spectrum import mag_ab, magnitudes_from_templates
+    from skypy.galaxy.spectrum import ab_maggies_redshift, template_absolute_magnitudes, template_apparent_magnitudes
     from astropy.cosmology import Planck15
-    from specutils import Spectrum1D
     from speclite.filters import FilterResponse
 
     # 3 Flat Templates
     lam = np.logspace(-1, 4, 1000)*units.AA
     A = np.array([[2], [3], [4]])
     flam = A * 0.10885464149979998*units.Unit('erg s-1 cm-2 AA')/lam**2
-    spec = Spectrum1D(spectral_axis=lam, flux=flam)
 
     # Gaussian bandpass
     filt_lam = np.logspace(0, 4, 1000)*units.AA
     filt_tx = np.exp(-((filt_lam - 1000*units.AA)/(100*units.AA))**2)
     filt_tx[[0, -1]] = 0
-    filt = FilterResponse(wavelength=filt_lam, response=filt_tx,
-                          meta=dict(group_name='test', band_name='filt'))
+    FilterResponse(wavelength=filt_lam, response=filt_tx,
+                   meta=dict(group_name='test', band_name='filt'))
 
     # Each test galaxy is exactly one of the templates
     coefficients = np.eye(3)
-    mt = magnitudes_from_templates(coefficients, spec, 'test-filt')
-    m = mag_ab(lam, flam, 'test-filt')
-    np.testing.assert_allclose(mt, m)
 
-    # Test distance modulus
+    # Absolute Magnitudes
+    mt = ab_maggies_redshift(lam, flam, 'test-filt', 0)
+    m = template_absolute_magnitudes(lam, flam, coefficients, 'test-filt')
+    np.testing.assert_allclose(-2.5*np.log10(mt), m)
+
+    # Apparent Magnitudes
     redshift = np.array([1, 2, 3])
     dm = Planck15.distmod(redshift).value
-    mt = magnitudes_from_templates(coefficients, spec, 'test-filt', distance_modulus=dm)
-    np.testing.assert_allclose(mt, m + dm)
+    mt = template_apparent_magnitudes(lam, flam, coefficients, redshift, 'test-filt', Planck15)
+    np.testing.assert_allclose(mt, m - 2.5*np.log10(1+redshift) + dm)
 
     # Test stellar mass
     sm = np.array([1, 2, 3])
-    mt = magnitudes_from_templates(coefficients, spec, 'test-filt', stellar_mass=sm)
+    mt = template_absolute_magnitudes(lam, flam, coefficients, 'test-filt', stellar_mass=sm)
     np.testing.assert_allclose(mt, m - 2.5*np.log10(sm))
 
     # Redshift interpolation test; linear interpolation sufficient over a small
     # redshift range at low relative tolerance
-    z = np.linspace(0, 0.1, 3)
-    m_true = magnitudes_from_templates(coefficients, spec, 'test-filt', redshift=z, resolution=4)
-    m_interp = magnitudes_from_templates(coefficients, spec, 'test-filt', redshift=z, resolution=2)
-    np.testing.assert_allclose(m_true, m_interp, rtol=1e-2)
+    z = np.linspace(0.1, 0.2, 3)
+    m_true = template_apparent_magnitudes(lam, flam, coefficients, z, 'test-filt', Planck15, resolution=4)
+    m_interp = template_apparent_magnitudes(lam, flam, coefficients, z, 'test-filt', Planck15, resolution=2)
+    np.testing.assert_allclose(m_true, m_interp, rtol=1e-7)
     assert not np.all(m_true == m_interp)