ENH: galaxy-quenching model

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_parameters
 
 
 Reference/API

skypy/galaxies/stellar_mass.py

@@ -8,6 +8,7 @@
 
 __all__ = [
     'schechter_smf_mass',
+    'schechter_smf_parameters',
 ]
 
 
@@ -74,3 +75,58 @@ 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_parameters(active_parameters, fsatellite, fenvironment):
+    r'''Schechter parameters.
+    This function returns the Schechter mass function parameters
+    for active galaxies (centrals and satellites)
+    and passive galaxies (mass- and satellite-quenched)
+    based on equation (15) in de la Bella et al. 2021 [1]_.
+
+    Parameters
+    ----------
+    active_parameters: tuple
+        Schechter mass function parameters for the entire active
+        sample of galaxies: :math:`(\phi_b, \alpha_b, m_{*})`.
+        Note: in [1], :math:`\alpha_b` and :math:`m_{*}` are
+        constants. The elements of the tuple should be scalars.
+        Exception: :math:`\phi_b` can be array only when
+        fsatellite is scalar.
+
+    fsatellite: float, (n, ) array_like
+        Fraction of active satellite galaxies between 0 and 1.
+        It could be a float or an array, depending on the model you choose.
+
+    fenvironment: float
+        Fraction of satellite-quenched galaxies between 0 and 1.
+
+
+    Returns
+    -------
+    parameters: dic
+        It returns a  dictionary with the parameters of the
+        Schechter mass function. The dictionary keywords:
+        `centrals`, `satellites`, `mass_quenched` and
+        `satellite_quenched`. The values correspond to a
+        tuple :math:`(\phi, \alpha, m_{*})`.
+
+    References
+    ----------
+    .. [1] de la Bella et al. 2021, Quenching and Galaxy Demographics,
+           arXiv 2112.11110.
+
+    '''
+    phi, alpha, mstar = active_parameters
+
+    sum_phics = (1 - fsatellite) * (1 - np.log(1 - fsatellite))
+    phic = (1 - fsatellite) * phi / sum_phics
+    phis = phic * np.log(1 / (1 - fsatellite))
+
+    centrals = (phic, alpha, mstar)
+    satellites = (phis, alpha, mstar)
+    mass_quenched = (phic + phis, alpha + 1, mstar)
+    satellite_quenched = (- np.log(1 - fenvironment) * phis, alpha, mstar)
+
+    return {'centrals': centrals, 'satellites': satellites,
+            'mass_quenched': mass_quenched, 'satellite_quenched': satellite_quenched}

skypy/galaxies/tests/test_stellar_mass.py

@@ -60,3 +60,61 @@ def calc_cdf(m):
                                              size=1000, resolution=100)
     p_value = scipy.stats.kstest(sample, calc_cdf)[1]
     assert p_value >= 0.01
+
+
+def test_schechter_smf_parameters():
+    from skypy.galaxies.stellar_mass import schechter_smf_parameters
+    # Check scalar inputs
+    blue = (10**-2.423, 10**10.60, -1.21)
+    fsat_scalar, frho = 0.4, 0.2
+    sp_scalar = schechter_smf_parameters(blue, fsat_scalar, frho)
+
+    assert type(sp_scalar) == dict, 'active_parameters is not a dict: {}'.format(type(sp_scalar))
+    assert len(sp_scalar) == 4, 'The length of the tuple is not four: {}'.format(len(sp_scalar))
+    for p in sp_scalar:
+        for k in range(3):
+            assert type(sp_scalar[p][k]) is not np.ndarray, \
+                '{} tuple is not a scalar {}'.format(p, type(sp_scalar[p][k]))
+
+    # Check array input for the satellite fraction
+    fsat_array = np.array([0.4, 0.5])
+    sp_array = schechter_smf_parameters(blue, fsat_array, frho)
+    assert type(sp_array) is dict, 'active_parameters is not a dict: {}'.format(type(sp_array))
+    assert len(sp_array) == 4, 'The length of the tuple is not four: {}'.format(len(sp_array))
+    for p in sp_array:
+        assert type(sp_array[p][0]) == np.ndarray, \
+            '{} amplitude is not an array {}'.format(p, type(sp_array[p][0]))
+        assert len(sp_array[p][0]) == len(fsat_array), \
+            '{} amplitude does not match input length {}'.format(p, len(fsat_array))
+        for k in range(1, 3):
+            assert type(sp_array[p][k]) != np.ndarray, \
+                '{} slope or mstar not a scalar {}'.format(p, type(sp_array[p][k]))
+
+    # Check array input for the amplitude
+    amplitude_array = np.array([10**-2.4, 10**-2.3])
+    blue_array = (amplitude_array, 10**10.60, -1.21)
+    sp_array = schechter_smf_parameters(blue_array, fsat_scalar, frho)
+    assert type(sp_array) is dict, 'active_parameters is not a dict: {}'.format(type(sp_array))
+    assert len(sp_array) == 4, 'The length of the tuple is not four: {}'.format(len(sp_array))
+    for p in sp_array:
+        assert type(sp_array[p][0]) is np.ndarray, \
+            '{} amplitude is not an array {}'.format(p, type(sp_array[p][0]))
+        assert len(sp_array[p][0]) is len(fsat_array), \
+            '{} amplitude does not match input length {}'.format(p, len(amplitude_array))
+        for k in range(1, 3):
+            assert type(sp_array[p][k]) is not np.ndarray, \
+                '{} slope or mstar not a scalar {}'.format(p, type(sp_array[p][k]))
+
+    # Corner cases
+    # Case I: no satellite galaxies
+    sp_sat0 = schechter_smf_parameters(blue, 0, frho)
+    assert sp_sat0['centrals'][0] == sp_sat0['mass_quenched'][0] == blue[0], \
+        'The amplitude of centrals and mass-quenched are not equal to the \
+         amplitude of the blue sample'
+    assert sp_sat0['satellites'][0] == sp_sat0['satellite_quenched'][0] == 0, \
+        'The amplitude of satellites and satellite-quenched are not zero'
+
+    # Case II: no satellite-quenched galaxies
+    sp_rho0 = schechter_smf_parameters(blue, fsat_scalar, 0)
+    assert sp_rho0['satellite_quenched'][0] == 0, \
+        'The satellite-quenched is not zero {}'.format(sp_rho0['satellite_quenched'][0])