Repository: skypyproject/skypy
Pull request #476 (2021-07-22 12:53:31)

ENH: compute kcorrect remaining stellar mass
## Description

Adds the `stellar_mass_remain()` function to the kcorrect template class which computes total surviving stellar mass

## Checklist
- [x] Follow the [Contributor Guidelines](https://github.com/skypyproject/skypy/blob/main/CONTRIBUTING.rst)
- [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


itrharrison: We are happy to merge this but I will write a test first!

Patch:
diff --git a/skypy/galaxies/spectrum.py b/skypy/galaxies/spectrum.py
index 3d0722e..f549276 100644
--- a/skypy/galaxies/spectrum.py
+++ b/skypy/galaxies/spectrum.py
@@ -143,6 +143,28 @@ class KCorrectTemplates(SpectrumTemplates):
         Mt = self.absolute_magnitudes(coefficients, filter)
         return np.power(10, 0.4*(Mt-magnitudes))
 
+    def stellar_mass_remain(self, coefficients, magnitudes, filter):
+        r'''Compute total surviving stellar mass from absolute magnitudes.
+
+        Parameters
+        ----------
+        coefficients : (ng, nt) array_like
+            Array of template coefficients.
+        magnitudes : (ng,) array_like
+            The magnitudes to match in the reference bandpass.
+        filter : str
+            A single reference bandpass filter specification for
+            `~speclite.filters.load_filters`.
+
+        Returns
+        -------
+        stellar_mass : (ng,) array_like
+            Total surviving stellar mass of each galaxy in template units.
+        '''
+        m = self.stellar_mass(coefficients, magnitudes, filter)
+        m *= np.dot(coefficients, self.mremain)
+        return m
+
     def metallicity(self, coefficients):
         r'''Galaxy metallicities from kcorrect templates.
 


Test patch:
diff --git a/skypy/galaxies/tests/test_spectrum.py b/skypy/galaxies/tests/test_spectrum.py
index 5aea0f4..f529530 100644
--- a/skypy/galaxies/tests/test_spectrum.py
+++ b/skypy/galaxies/tests/test_spectrum.py
@@ -158,6 +158,39 @@ def test_kcorrect_stellar_mass():
     np.testing.assert_allclose(stellar_mass, truth)
 
 
+@pytest.mark.skipif(not HAS_SPECLITE, reason='test requires speclite')
+def test_kcorrect_stellar_mass_remain():
+
+    from astropy import units
+    from skypy.galaxies.spectrum import kcorrect
+    from speclite.filters import FilterResponse
+
+    # Gaussian bandpass
+    filt_lam = np.logspace(3, 4, 1000) * units.AA
+    filt_mean = 5000 * units.AA
+    filt_width = 100 * units.AA
+    filt_tx = np.exp(-((filt_lam-filt_mean)/filt_width)**2)
+    filt_tx[[0, -1]] = 0
+    FilterResponse(wavelength=filt_lam, response=filt_tx,
+                   meta=dict(group_name='test', band_name='filt'))
+
+    # Using the identity matrix for the coefficients yields trivial test cases
+    coeff = np.eye(5)
+    Mt = kcorrect.absolute_magnitudes(coeff, 'test-filt')
+
+    # Using the absolute magnitudes of the templates as reference magnitudes
+    # should return one solar mass for each template.
+    stellar_mass = kcorrect.stellar_mass_remain(coeff, Mt, 'test-filt')
+    truth = kcorrect.mremain
+    np.testing.assert_allclose(stellar_mass, truth)
+
+    # Solution for given magnitudes without template mixing
+    Mb = np.array([10, 20, 30, 40, 50])
+    stellar_mass = kcorrect.stellar_mass_remain(coeff, Mb, 'test-filt')
+    truth = np.power(10, -0.4*(Mb-Mt)) * kcorrect.mremain
+    np.testing.assert_allclose(stellar_mass, truth)
+
+
 def test_kcorrect_metallicity():
 
     from skypy.galaxies.spectrum import kcorrect