ENH: Time evolution Schechter amplitude

docs/galaxies.rst

@@ -89,6 +89,7 @@ The following models are found in the `skypy.galaxies.stellar_mass` package.
   :nosignatures:
 
   schechter_smf_mass
+  schechter_smf_phi_active_redshift
 
 
 Velocity dispersion

skypy/galaxies/stellar_mass.py

@@ -2,13 +2,15 @@
 """
 
 import numpy as np
+from scipy.integrate import trapz
 
 from ..utils.random import schechter
 from ..utils import dependent_argument
 
 
 __all__ = [
     'schechter_smf_mass',
+    'schechter_smf_phi_active_redshift',
 ]
 
 
@@ -78,3 +80,67 @@ def schechter_smf_mass(redshift, alpha, m_star, m_min, m_max, size=None,
     m = schechter(alpha, x_min, x_max, resolution, size=size, scale=m_star)
 
     return m
+
+
+def schechter_smf_phi_active_redshift(redshift,
+                                      phi_active_today,
+                                      probability_satellites,
+                                      cosmology,
+                                      redshift_initial=10,
+                                      resolution=100):
+    r'''Redshift-dependent Schechter amplitude of active galaxies.
+    This function returns the time-dependent Schechter mass function amplitude
+    for the total active population based on equation (17)
+    in de la Bella et al. 2021 [1]_.
+
+    Parameters
+    ----------
+    redshift:
+        Values of redshift at which to evaluate the amplitude.
+
+    phi_active_today: array_like
+        Schechter mass function amplitude for the entire active
+        sample of galaxies in the past when density and mass was very low.
+
+    probability_satellites: func
+        Probability of active galaxies becoming satellites
+        as a function of redshift.
+
+    cosmology: astropy.cosmology.Cosmology
+        Cosmology object providing methods for the evolution history
+        of omega_matter and omega_lambda with redshift.
+
+    redshift_initial: float
+        Value of redshift in the past when density and mass was very low.
+        Default is 10.
+
+    resolution: float
+        Resolution of the integral. Default is 100.
+
+    Returns
+    -------
+    amplitude: array_like
+        Amplitude of the Schechter mass function.
+
+     References
+    ----------
+    .. [1] de la Bella et al. 2021, Quenching and Galaxy Demographics,
+           arXiv 2112.11110.
+
+    '''
+
+    def integrand(z):
+        return probability_satellites(z) / cosmology.H(z).value / (1 + z)
+
+    # Calculate the amplitude in the past
+    redshift_today = np.linspace(0, redshift_initial, resolution)
+    integrand_today = integrand(redshift_today)
+    integral_today = trapz(integrand_today)
+    B = phi_active_today * np.exp(- integral_today) / (1 - integral_today)
+
+    # Calculate the amplitude for the given redshift
+    redshift_array = np.linspace(redshift, redshift_initial, resolution)
+    integrand_redshift = integrand(redshift_array)
+    integral = trapz(integrand_redshift)
+
+    return B * (1 - integral) * np.exp(integral)

skypy/galaxies/tests/test_stellar_mass.py

@@ -4,6 +4,8 @@
 from scipy.special import gammaln
 import pytest
 from astropy.modeling.models import Exponential1D
+from hypothesis import given
+from hypothesis.strategies import integers
 
 from skypy.galaxies import stellar_mass
 from skypy.utils import special
@@ -79,3 +81,30 @@ def calc_cdf(m):
                                              size=1000, resolution=100)
     p_value = scipy.stats.kstest(sample, calc_cdf)[1]
     assert p_value >= 0.01
+
+
+@given(integers(1, 100))
+def test_schechter_smf_phi_active_redshift(redshift):
+    # Amplitude today and cosmology
+    from astropy.cosmology import Planck15
+    phi_today = 10**-2.423
+
+    # Simple probability of satellites
+    def prob_sat(z, a=0.4):
+        return a * np.exp(- z)
+
+    # Test that today the output of the active amplitude equals the input
+    output_today = stellar_mass.schechter_smf_phi_active_redshift(0, phi_today, prob_sat, Planck15)
+    assert output_today == pytest.approx(phi_today)
+
+    # Test that at any redshift for a null probability of creating satellite galaxies
+    # the output of the active amplitude equals the input
+    def prob_null(z):
+        return prob_sat(z, 0)
+
+    output_null = stellar_mass.schechter_smf_phi_active_redshift(redshift,
+                                                                 phi_today,
+                                                                 prob_null,
+                                                                 Planck15
+                                                                 )
+    assert output_null == pytest.approx(phi_today)