From 8ccafdbfa9d314c9acde91372c210ee661367fdd Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Thu, 17 Feb 2022 16:13:27 +0000
Subject: [PATCH 01/36] Add feature to calculate the schechter params with de
 la Bella quenching model

---
 docs/galaxies.rst              |  1 +
 skypy/galaxies/stellar_mass.py | 52 ++++++++++++++++++++++++++++++++++
 2 files changed, 53 insertions(+)

diff --git a/docs/galaxies.rst b/docs/galaxies.rst
index bfddeeec6..229a94b04 100644
--- a/docs/galaxies.rst
+++ b/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
diff --git a/skypy/galaxies/stellar_mass.py b/skypy/galaxies/stellar_mass.py
index 2b1a33e86..9173dfe7f 100644
--- a/skypy/galaxies/stellar_mass.py
+++ b/skypy/galaxies/stellar_mass.py
@@ -8,6 +8,7 @@
 
 __all__ = [
     'schechter_smf_mass',
+    'schechter_smf_parameters',
 ]
 
 
@@ -74,3 +75,54 @@ 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_{*})`.
+
+    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 = (phi, 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}

From 90153d64b149c8b8aa5c144eaa4a4026ff27d0c2 Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Thu, 17 Feb 2022 16:34:00 +0000
Subject: [PATCH 02/36] Add feature to calculate the schechter params with de
 la Bella quenching model

---
 skypy/halo/quenching.py            | 106 +++++++++++++++++++++++++++++
 skypy/halo/tests/test_quenching.py |  30 ++++++++
 2 files changed, 136 insertions(+)
 create mode 100644 skypy/halo/quenching.py
 create mode 100644 skypy/halo/tests/test_quenching.py

diff --git a/skypy/halo/quenching.py b/skypy/halo/quenching.py
new file mode 100644
index 000000000..c3a332004
--- /dev/null
+++ b/skypy/halo/quenching.py
@@ -0,0 +1,106 @@
+"""Galaxy quenching.
+
+This module implements models for environment and mass
+quenching by dark matter halos.
+"""
+
+import numpy as np
+from scipy import special
+
+__all__ = [
+    'environment_quenched',
+    'mass_quenched',
+    ]
+
+
+def environment_quenched(nh, probability):
+    r'''Environment quenching function.
+    This function implements the model proposed by A.Amara where the
+    probability of a subhalo being quenched is a fixed
+    probability. The model is inspired on [1]_ and [2]_.
+
+    Parameters
+    ----------
+    nh: integer
+        Number of subhalos.
+    probability: float
+        Quenching probability.
+
+    Returns
+    -------
+    quenched: (nh,) array_like,  boolean
+        Boolean array indicating which subhalo's host galaxies are
+        (satellite) environment-quenched.
+
+    Examples
+    ---------
+
+    This example shows how many subhalos are environtment quenched (True)
+    and how many survive (False) from a list of 1000 halos:
+
+    >>> import numpy as np
+    >>> from skypy.halo.quenching import environment_quenched
+    >>> from collections import Counter
+    >>> quenched = environment_quenched(1000, 0.5)
+    >>> Counter(quenched)
+    Counter({...})
+
+    References
+    ----------
+    .. [1] Peng et al. 2010, doi 10.1088/0004-637X/721/1/193.
+    .. [2] Birrer et al. 2014, arXiv 1401.3162.
+
+    '''
+
+    return np.random.uniform(size=nh) < probability
+
+
+def mass_quenched(halo_mass, offset, width):
+    r'''Mass quenching function.
+    This function implements the model proposed by A.Amara where the
+    probability of a halo being quenched is related to the error function
+    of the logarithm of the halo's mass standardised by an offset and width
+    parameter.  The model is inspired on [1]_ and [2]_.
+
+    Parameters
+    ----------
+    halo_mass: (nh,) array_like
+        Array of halo masses in units of solar mass, :math:`M_{sun}`.
+    offset: float
+        Halo mass in :math:`M_{sun}` at which quenching probability is 50%.
+    width: float
+        Width of the error function.
+
+    Returns
+    -------
+    quenched: (nh,) array_like, boolean
+        Boolean array indicating which halo's host galaxies are mass-quenched.
+
+    Examples
+    ---------
+
+    This example shows how many halos are mass quenched (True)
+    and how many survive (False) from a list of 1000 halos:
+
+    >>> import numpy as np
+    >>> from astropy import units
+    >>> from skypy.halo.quenching import mass_quenched
+    >>> from collections import Counter
+    >>> offset, width = 1.0e12, 0.5
+    >>> halo_mass = np.random.lognormal(mean=np.log(offset), sigma=width,
+    ...                                 size=1000)
+    >>> quenched = mass_quenched(halo_mass, offset, width)
+    >>> Counter(quenched)
+    Counter({...})
+
+    References
+    ----------
+    .. [1] Peng et al. 2010, doi 10.1088/0004-637X/721/1/193.
+    .. [2] Birrer et al. 2014, arXiv 1401.3162.
+
+    '''
+
+    standardised_mass = np.log10(halo_mass / offset) / width
+    probability = 0.5 * (1.0 + special.erf(standardised_mass/np.sqrt(2)))
+    nh = len(halo_mass)
+    return np.random.uniform(size=nh) < probability
diff --git a/skypy/halo/tests/test_quenching.py b/skypy/halo/tests/test_quenching.py
new file mode 100644
index 000000000..a3f18b592
--- /dev/null
+++ b/skypy/halo/tests/test_quenching.py
@@ -0,0 +1,30 @@
+import numpy as np
+from scipy import stats
+from collections import Counter
+
+from skypy.halo.quenching import environment_quenched, mass_quenched
+
+
+def test_environment_quenched():
+    # Test the quenching process follows a binomial distribution
+    n, p = 1000, 0.7
+    quenched = environment_quenched(n, p)
+    number_quenched = Counter(quenched)[True]
+
+    p_value = stats.binom_test(number_quenched, n=n, p=p,
+                               alternative='two-sided')
+    assert p_value > 0.05
+
+
+def test_mass_quenched():
+    # Test the quenching process follows a binomial distribution if the
+    # logarithmic halo masses are symmetrical about the offset
+    n = 1000
+    offset, width = 1.0e12, 0.5
+    halo_mass = 10 ** np.random.uniform(11, 13, n)
+    quenched = mass_quenched(halo_mass, offset, width)
+    number_quenched = Counter(quenched)[True]
+
+    p_value = stats.binom_test(number_quenched, n=n, p=0.5,
+                               alternative='two-sided')
+    assert p_value > 0.05

From 1ada400195ac5bfc785f3a4aa7f212b76813747c Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Thu, 17 Feb 2022 16:40:46 +0000
Subject: [PATCH 03/36] Add placeholder for schechter_smf_parameters test

---
 skypy/galaxies/tests/test_stellar_mass.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py
index cb6679892..aba1461fa 100644
--- a/skypy/galaxies/tests/test_stellar_mass.py
+++ b/skypy/galaxies/tests/test_stellar_mass.py
@@ -60,3 +60,7 @@ 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()
+    # placeholder for test
\ No newline at end of file

From d69728a8ffee7fc1fd77e367ffeff943c6f8d772 Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Thu, 17 Feb 2022 16:42:03 +0000
Subject: [PATCH 04/36] Remove halo quenching file

---
 skypy/halo/quenching.py | 106 ----------------------------------------
 1 file changed, 106 deletions(-)
 delete mode 100644 skypy/halo/quenching.py

diff --git a/skypy/halo/quenching.py b/skypy/halo/quenching.py
deleted file mode 100644
index c3a332004..000000000
--- a/skypy/halo/quenching.py
+++ /dev/null
@@ -1,106 +0,0 @@
-"""Galaxy quenching.
-
-This module implements models for environment and mass
-quenching by dark matter halos.
-"""
-
-import numpy as np
-from scipy import special
-
-__all__ = [
-    'environment_quenched',
-    'mass_quenched',
-    ]
-
-
-def environment_quenched(nh, probability):
-    r'''Environment quenching function.
-    This function implements the model proposed by A.Amara where the
-    probability of a subhalo being quenched is a fixed
-    probability. The model is inspired on [1]_ and [2]_.
-
-    Parameters
-    ----------
-    nh: integer
-        Number of subhalos.
-    probability: float
-        Quenching probability.
-
-    Returns
-    -------
-    quenched: (nh,) array_like,  boolean
-        Boolean array indicating which subhalo's host galaxies are
-        (satellite) environment-quenched.
-
-    Examples
-    ---------
-
-    This example shows how many subhalos are environtment quenched (True)
-    and how many survive (False) from a list of 1000 halos:
-
-    >>> import numpy as np
-    >>> from skypy.halo.quenching import environment_quenched
-    >>> from collections import Counter
-    >>> quenched = environment_quenched(1000, 0.5)
-    >>> Counter(quenched)
-    Counter({...})
-
-    References
-    ----------
-    .. [1] Peng et al. 2010, doi 10.1088/0004-637X/721/1/193.
-    .. [2] Birrer et al. 2014, arXiv 1401.3162.
-
-    '''
-
-    return np.random.uniform(size=nh) < probability
-
-
-def mass_quenched(halo_mass, offset, width):
-    r'''Mass quenching function.
-    This function implements the model proposed by A.Amara where the
-    probability of a halo being quenched is related to the error function
-    of the logarithm of the halo's mass standardised by an offset and width
-    parameter.  The model is inspired on [1]_ and [2]_.
-
-    Parameters
-    ----------
-    halo_mass: (nh,) array_like
-        Array of halo masses in units of solar mass, :math:`M_{sun}`.
-    offset: float
-        Halo mass in :math:`M_{sun}` at which quenching probability is 50%.
-    width: float
-        Width of the error function.
-
-    Returns
-    -------
-    quenched: (nh,) array_like, boolean
-        Boolean array indicating which halo's host galaxies are mass-quenched.
-
-    Examples
-    ---------
-
-    This example shows how many halos are mass quenched (True)
-    and how many survive (False) from a list of 1000 halos:
-
-    >>> import numpy as np
-    >>> from astropy import units
-    >>> from skypy.halo.quenching import mass_quenched
-    >>> from collections import Counter
-    >>> offset, width = 1.0e12, 0.5
-    >>> halo_mass = np.random.lognormal(mean=np.log(offset), sigma=width,
-    ...                                 size=1000)
-    >>> quenched = mass_quenched(halo_mass, offset, width)
-    >>> Counter(quenched)
-    Counter({...})
-
-    References
-    ----------
-    .. [1] Peng et al. 2010, doi 10.1088/0004-637X/721/1/193.
-    .. [2] Birrer et al. 2014, arXiv 1401.3162.
-
-    '''
-
-    standardised_mass = np.log10(halo_mass / offset) / width
-    probability = 0.5 * (1.0 + special.erf(standardised_mass/np.sqrt(2)))
-    nh = len(halo_mass)
-    return np.random.uniform(size=nh) < probability

From 7e90c0a5512e85df363e90d3fc68e7c1de0547dd Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Thu, 17 Feb 2022 16:42:52 +0000
Subject: [PATCH 05/36] Fix pep8 errors for test

---
 skypy/galaxies/tests/test_stellar_mass.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py
index aba1461fa..61367df95 100644
--- a/skypy/galaxies/tests/test_stellar_mass.py
+++ b/skypy/galaxies/tests/test_stellar_mass.py
@@ -62,5 +62,6 @@ def calc_cdf(m):
     assert p_value >= 0.01
 
 
-def test_schechter_smf_parameters()
-    # placeholder for test
\ No newline at end of file
+def test_schechter_smf_parameters():
+    # placeholder for test
+    return 0
\ No newline at end of file

From 350310fc6d583a515b76d88dee9de99739197d08 Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Thu, 17 Feb 2022 16:43:43 +0000
Subject: [PATCH 06/36] Remove test halo quenching file

---
 skypy/halo/tests/test_quenching.py | 30 ------------------------------
 1 file changed, 30 deletions(-)
 delete mode 100644 skypy/halo/tests/test_quenching.py

diff --git a/skypy/halo/tests/test_quenching.py b/skypy/halo/tests/test_quenching.py
deleted file mode 100644
index a3f18b592..000000000
--- a/skypy/halo/tests/test_quenching.py
+++ /dev/null
@@ -1,30 +0,0 @@
-import numpy as np
-from scipy import stats
-from collections import Counter
-
-from skypy.halo.quenching import environment_quenched, mass_quenched
-
-
-def test_environment_quenched():
-    # Test the quenching process follows a binomial distribution
-    n, p = 1000, 0.7
-    quenched = environment_quenched(n, p)
-    number_quenched = Counter(quenched)[True]
-
-    p_value = stats.binom_test(number_quenched, n=n, p=p,
-                               alternative='two-sided')
-    assert p_value > 0.05
-
-
-def test_mass_quenched():
-    # Test the quenching process follows a binomial distribution if the
-    # logarithmic halo masses are symmetrical about the offset
-    n = 1000
-    offset, width = 1.0e12, 0.5
-    halo_mass = 10 ** np.random.uniform(11, 13, n)
-    quenched = mass_quenched(halo_mass, offset, width)
-    number_quenched = Counter(quenched)[True]
-
-    p_value = stats.binom_test(number_quenched, n=n, p=0.5,
-                               alternative='two-sided')
-    assert p_value > 0.05

From e9339dd83ae3c93d7da7be1b21a40d9aba4324b5 Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Thu, 17 Feb 2022 16:45:13 +0000
Subject: [PATCH 07/36] Fix no end of line error

---
 skypy/galaxies/tests/test_stellar_mass.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py
index 61367df95..34855d739 100644
--- a/skypy/galaxies/tests/test_stellar_mass.py
+++ b/skypy/galaxies/tests/test_stellar_mass.py
@@ -64,4 +64,4 @@ def calc_cdf(m):
 
 def test_schechter_smf_parameters():
     # placeholder for test
-    return 0
\ No newline at end of file
+    return 0

From 4ef115660634d4e49f34efd08e18e949ca3f62f0 Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Mon, 21 Feb 2022 09:57:43 +0000
Subject: [PATCH 08/36] Clarify docstring

---
 skypy/galaxies/stellar_mass.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/skypy/galaxies/stellar_mass.py b/skypy/galaxies/stellar_mass.py
index 9173dfe7f..0e81bd497 100644
--- a/skypy/galaxies/stellar_mass.py
+++ b/skypy/galaxies/stellar_mass.py
@@ -89,6 +89,8 @@ def schechter_smf_parameters(active_parameters, fsatellite, fenvironment):
     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.
 
     fsatellite: float, (n, ) array_like
         Fraction of active satellite galaxies between 0 and 1.

From 7670106a6ec4474b37d63205f2f00b33f3788c41 Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Mon, 21 Feb 2022 10:02:18 +0000
Subject: [PATCH 09/36] Clarify docstring for amplitude input

---
 skypy/galaxies/stellar_mass.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/skypy/galaxies/stellar_mass.py b/skypy/galaxies/stellar_mass.py
index 0e81bd497..5d94acc74 100644
--- a/skypy/galaxies/stellar_mass.py
+++ b/skypy/galaxies/stellar_mass.py
@@ -91,6 +91,8 @@ def schechter_smf_parameters(active_parameters, fsatellite, fenvironment):
         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.

From 6ce403b34d7a2ede465a8997b9bdb949cb070e05 Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Mon, 21 Feb 2022 11:10:15 +0000
Subject: [PATCH 10/36] Add test for the schechter parameters function

---
 skypy/galaxies/tests/test_stellar_mass.py | 57 ++++++++++++++++++++++-
 1 file changed, 55 insertions(+), 2 deletions(-)

diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py
index 34855d739..8698bb1b9 100644
--- a/skypy/galaxies/tests/test_stellar_mass.py
+++ b/skypy/galaxies/tests/test_stellar_mass.py
@@ -63,5 +63,58 @@ def calc_cdf(m):
 
 
 def test_schechter_smf_parameters():
-    # placeholder for test
-    return 0
+    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 satellite fraction
+    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])

From d97bb79be8a1885313fd81c42d2b9ff6bb7a09c2 Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Mon, 21 Feb 2022 11:13:24 +0000
Subject: [PATCH 11/36] Fix trailing whitespace

---
 skypy/galaxies/stellar_mass.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skypy/galaxies/stellar_mass.py b/skypy/galaxies/stellar_mass.py
index 5d94acc74..4c134911c 100644
--- a/skypy/galaxies/stellar_mass.py
+++ b/skypy/galaxies/stellar_mass.py
@@ -92,7 +92,7 @@ def schechter_smf_parameters(active_parameters, fsatellite, fenvironment):
         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 is scalar.
 
     fsatellite: float, (n, ) array_like
         Fraction of active satellite galaxies between 0 and 1.

From 829e6b46c351b951b63b8108d20ddceb1a147fbc Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Mon, 21 Feb 2022 11:21:52 +0000
Subject: [PATCH 12/36] update case II assertment

---
 skypy/galaxies/tests/test_stellar_mass.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py
index 8698bb1b9..9dc3c3963 100644
--- a/skypy/galaxies/tests/test_stellar_mass.py
+++ b/skypy/galaxies/tests/test_stellar_mass.py
@@ -116,5 +116,5 @@ def test_schechter_smf_parameters():
 
     # Case II: no satellite-quenched galaxies
     sp_rho0 = schechter_smf_parameters(blue, fsat_scalar, 0)
-    assert - sp_rho0['satellite_quenched'][0] == 0, \
+    assert sp_rho0['satellite_quenched'][0] == 0, \
         'The satellite-quenched is not zero {}'.format(sp_rho0['satellite_quenched'][0])

From 4a920b6e650c2d4c300aa425020ba5640c42b9e7 Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Mon, 21 Feb 2022 11:26:24 +0000
Subject: [PATCH 13/36] Fix comment

---
 skypy/galaxies/tests/test_stellar_mass.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py
index 9dc3c3963..ca27b26eb 100644
--- a/skypy/galaxies/tests/test_stellar_mass.py
+++ b/skypy/galaxies/tests/test_stellar_mass.py
@@ -90,7 +90,7 @@ def test_schechter_smf_parameters():
             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 satellite fraction
+    # 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)

From e504479065fccf34773ea1cfc694ac0f4ecd6100 Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Mon, 21 Feb 2022 11:29:40 +0000
Subject: [PATCH 14/36] Fix bug with mass-quenched galaxies

---
 skypy/galaxies/stellar_mass.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skypy/galaxies/stellar_mass.py b/skypy/galaxies/stellar_mass.py
index 4c134911c..8aeda4d50 100644
--- a/skypy/galaxies/stellar_mass.py
+++ b/skypy/galaxies/stellar_mass.py
@@ -125,7 +125,7 @@ def schechter_smf_parameters(active_parameters, fsatellite, fenvironment):
 
     centrals = (phic, alpha, mstar)
     satellites = (phis, alpha, mstar)
-    mass_quenched = (phi, alpha + 1, mstar)
+    mass_quenched = (phic + phis, alpha + 1, mstar)
     satellite_quenched = (- np.log(1 - fenvironment) * phis, alpha, mstar)
 
     return {'centrals': centrals, 'satellites': satellites,

From e85f24f1b2f0e2715cffe9f1e0886fc283f81290 Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Wed, 15 Jun 2022 16:31:00 +0100
Subject: [PATCH 15/36] Address major refactoring review

---
 skypy/galaxies/stellar_mass.py            | 147 +++++++++++++++-----
 skypy/galaxies/tests/test_stellar_mass.py | 162 ++++++++++++++--------
 2 files changed, 219 insertions(+), 90 deletions(-)

diff --git a/skypy/galaxies/stellar_mass.py b/skypy/galaxies/stellar_mass.py
index 8aeda4d50..f70416d6f 100644
--- a/skypy/galaxies/stellar_mass.py
+++ b/skypy/galaxies/stellar_mass.py
@@ -8,7 +8,10 @@
 
 __all__ = [
     'schechter_smf_mass',
-    'schechter_smf_parameters',
+    'schechter_smf_amplitude_centrals',
+    'schechter_smf_amplitude_satellites',
+    'schechter_smf_amplitude_mass_quenched',
+    'schechter_smf_amplitude_satellite_quenched',
 ]
 
 
@@ -77,39 +80,62 @@ def schechter_smf_mass(redshift, alpha, m_star, m_min, m_max, size=None,
     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]_.
+def schechter_smf_amplitude_centrals(phi_blue_total, fsatellite):
+    r'''Schechter amplitude of centrals.
+    This function returns the Schechter mass function amplitude
+    for active central population 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
+    phi_blue_total: float, (nt, ) array_like
+        Schechter mass function amplitude for the entire active
+        sample of galaxies, :math:`(\phi_b, \alpha_b, m_{*})`.
+
+    fsatellite: float, (nm, ) 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
+    -------
+    amplitude: float, (nt, nm) array_like
+        Amplitude of the Schechter mass function.
+
+     References
+    ----------
+    .. [1] de la Bella et al. 2021, Quenching and Galaxy Demographics,
+           arXiv 2112.11110.
+
+    '''
+
+    if np.ndim(phi_blue_total) == 1 and np.ndim(fsatellite) == 1:
+        phi_blue_total = phi_blue_total[:, np.newaxis]
+
+    sum_phics = (1 - fsatellite) * (1 - np.log(1 - fsatellite))
+
+    return (1 - fsatellite) * phi_blue_total / sum_phics
 
 
+def schechter_smf_amplitude_satellites(phi_centrals, fsatellite):
+    r'''Schechter amplitude of satellites.
+    This function returns the Schechter mass function amplitude
+    for active satellite population based on equation (15)
+    in de la Bella et al. 2021 [1]_.
+
+    Parameters
+    ----------
+    phi_centrals: float, (nt, nm) array_like
+        Schechter mass function amplitude of the central
+        active galaxies.
+
+    fsatellite: float, (nm, ) 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.
+
     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_{*})`.
+    amplitude: float, (nt, nm) array_like
+        Amplitude of the Schechter mass function.
 
     References
     ----------
@@ -117,16 +143,69 @@ def schechter_smf_parameters(active_parameters, fsatellite, fenvironment):
            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))
+    if np.ndim(phi_centrals) == 2 and np.ndim(fsatellite) == 1:
+        phi_centrals = phi_centrals[:, np.newaxis]
+
+    return phi_centrals * np.log(1 / (1 - fsatellite))
+
+
+def schechter_smf_amplitude_mass_quenched(phi_centrals, phi_satellites):
+    r'''Schechter amplitude of satellites.
+    This function returns the Schechter mass function amplitude
+    for passive mass-quenched population based on equation (15)
+    in de la Bella et al. 2021 [1]_.
+
+    Parameters
+    ----------
+    phi_centrals: float, (nt, nm) array_like
+        Schechter mass function amplitude of the central
+        active galaxies.
 
-    centrals = (phic, alpha, mstar)
-    satellites = (phis, alpha, mstar)
-    mass_quenched = (phic + phis, alpha + 1, mstar)
-    satellite_quenched = (- np.log(1 - fenvironment) * phis, alpha, mstar)
+    phi_satellites: float, (nt, nm) array_like
+        Schechter mass function amplitude of the satellite
+        active galaxies.
+
+    Returns
+    -------
+    amplitude: float, (nt, nm) array_like
+        Amplitude of the Schechter mass function.
+
+    References
+    ----------
+    .. [1] de la Bella et al. 2021, Quenching and Galaxy Demographics,
+           arXiv 2112.11110.
+
+    '''
+
+    return phi_centrals + phi_satellites
+
+
+def schechter_smf_amplitude_satellite_quenched(phi_satellites, fenvironment):
+    r'''Schechter amplitude of satellites.
+    This function returns the Schechter mass function amplitude
+    for active central population based on equation (15)
+    in de la Bella et al. 2021 [1]_.
+
+    Parameters
+    ----------
+    phi_satellites: float, (nt, nm) array_like
+        Schechter mass function amplitude of the satellite
+        active galaxies.
+
+    fenvironment: float
+        Fraction of satellite-quenched galaxies between 0 and 1.
+
+    Returns
+    -------
+    amplitude: float, (nt, nm) array_like
+        Amplitude of the Schechter mass function.
+
+    References
+    ----------
+    .. [1] de la Bella et al. 2021, Quenching and Galaxy Demographics,
+           arXiv 2112.11110.
+
+    '''
 
-    return {'centrals': centrals, 'satellites': satellites,
-            'mass_quenched': mass_quenched, 'satellite_quenched': satellite_quenched}
+    return - np.log(1 - fenvironment) * phi_satellites
diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py
index ca27b26eb..5268b0545 100644
--- a/skypy/galaxies/tests/test_stellar_mass.py
+++ b/skypy/galaxies/tests/test_stellar_mass.py
@@ -62,59 +62,109 @@ def calc_cdf(m):
     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])
+def test_schechter_smf_amplitude_centrals():
+    # Scalar inputs
+    phiblue_scalar = 10**-2.423
+    fsat_scalar = 0.4
+
+    # Array inputs
+    phiblue_array = np.array([10**-2.423, 10**-2.422])
+    fsat_array = np.array([0.40, 0.41, 0.42])
+
+    # Test for scalar output
+    phic_scalar = stellar_mass.schechter_smf_amplitude_centrals(phiblue_scalar, fsat_scalar)
+    assert np.isscalar(phic_scalar)
+
+    # Test for 1 dim output
+    phic_1d_phib = stellar_mass.schechter_smf_amplitude_centrals(phiblue_array, fsat_scalar)
+    phic_1d_fsat = stellar_mass.schechter_smf_amplitude_centrals(phiblue_scalar, fsat_array)
+    assert phic_1d_phib.shape == phiblue_array.shape
+    assert phic_1d_fsat.shape == fsat_array.shape
+
+    # Test for 2 dim output
+    phic_2d = stellar_mass.schechter_smf_amplitude_centrals(phiblue_array, fsat_array)
+    assert phic_2d.shape == (len(phiblue_array), len(fsat_array))
+
+
+def test_schechter_smf_amplitude_satellites():
+    # Scalar inputs
+    phic_scalar = 10**-2.423
+    fsat_scalar = 0.4
+
+    # Array inputs
+    phic_array = np.array([10**-2.423, 10**-2.422])
+    fsat_array = np.array([0.40, 0.41, 0.42])
+
+    # Test for scalar output
+    phis_scalar = stellar_mass.schechter_smf_amplitude_satellites(phic_scalar, fsat_scalar)
+    assert np.isscalar(phis_scalar)
+
+    # Test for 1 dim output
+    phis_1d_phib = stellar_mass.schechter_smf_amplitude_satellites(phic_array, fsat_scalar)
+    phis_1d_fsat = stellar_mass.schechter_smf_amplitude_satellites(phic_scalar, fsat_array)
+    assert phis_1d_phib.shape == phic_array.shape
+    assert phis_1d_fsat.shape == fsat_array.shape
+
+    # Test for 2 dim output
+    phis_2d = stellar_mass.schechter_smf_amplitude_satellites(phic_array, fsat_array)
+    assert phis_2d.shape == (len(phic_array), len(fsat_array))
+
+    # Corner case: no satellite galaxies
+    fsat_null = 0
+    phis_null_scalar = stellar_mass.schechter_smf_amplitude_satellites(phic_scalar, fsat_null)
+    phis_null_array = stellar_mass.schechter_smf_amplitude_satellites(phic_array, fsat_null)
+    assert phis_null_scalar == 0
+    assert np.all(phis_null_array == np.zeros(len(phic_array)))
+
+
+def test_schechter_smf_amplitude_mass_quenched():
+    # Scalar inputs
+    phic_scalar = 10**-2.5
+    phis_scalar = 10**-2.4
+
+    # 1D Array inputs
+    phic_1d = np.array([10**-2.50, 10**-2.51])
+    phis_1d = np.array([10**-2.40, 10**-2.41])
+
+    # 2D Array inputs
+    phic_2d = np.array([[10**-2.50, 10**-2.51], [10**-2.50, 10**-2.51]])
+    phis_2d = np.array([[10**-2.40, 10**-2.41], [10**-2.40, 10**-2.41]])
+
+    # Test for scalar output
+    phimq_scalar = stellar_mass.schechter_smf_amplitude_mass_quenched(phic_scalar, phis_scalar)
+    assert np.isscalar(phis_scalar)
+    assert phimq_scalar == phic_scalar + phis_scalar
+
+    # Test for 1 dim output
+    phimq_1d = stellar_mass.schechter_smf_amplitude_mass_quenched(phic_1d, phis_1d)
+    assert phimq_1d.shape == phic_1d.shape == phis_1d.shape
+    assert np.all(phimq_1d == phic_1d + phis_1d)
+
+    # Test for 2 dim output
+    phimq_2d = stellar_mass.schechter_smf_amplitude_mass_quenched(phic_2d, phis_2d)
+    assert phimq_2d.shape == phic_2d.shape == phis_2d.shape
+    assert np.all(phimq_2d == phic_2d + phis_2d)
+
+
+def test_schechter_smf_amplitude_satellite_quenched():
+    # Scalar inputs
+    phis_scalar = 10**-2.423
+    frho_scalar = 0.2
+
+    # Array inputs
+    phis_array = np.array([10**-2.423, 10**-2.422])
+
+    # Test for scalar output
+    phisq_scalar = stellar_mass.schechter_smf_amplitude_satellite_quenched(phis_scalar, frho_scalar)
+    assert np.isscalar(phisq_scalar)
+
+    # Test for 1 dim output
+    phisq_1d_phib = stellar_mass.schechter_smf_amplitude_satellite_quenched(phis_array, frho_scalar)
+    assert phisq_1d_phib.shape == phis_array.shape
+
+    # Corner case no satellite-quenched galaxies
+    frho_null = 0
+    phisq_null_scalar = stellar_mass.schechter_smf_amplitude_satellites(phis_scalar, frho_null)
+    phisq_null_array = stellar_mass.schechter_smf_amplitude_satellites(phis_array, frho_null)
+    assert phisq_null_scalar == 0
+    assert np.all(phisq_null_array == np.zeros(len(phis_array)))

From 133be5a5bcb2ae4005a185f89d2c245be403e11b Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Wed, 15 Jun 2022 16:33:28 +0100
Subject: [PATCH 16/36] Fix docs by adding new functions to the galaxy rst file

---
 docs/galaxies.rst | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/docs/galaxies.rst b/docs/galaxies.rst
index 229a94b04..c31a36a38 100644
--- a/docs/galaxies.rst
+++ b/docs/galaxies.rst
@@ -89,7 +89,10 @@ The following models are found in the `skypy.galaxies.stellar_mass` package.
   :nosignatures:
 
   schechter_smf_mass
-  schechter_smf_parameters
+  schechter_smf_amplitude_centrals
+  schechter_smf_amplitude_satellites
+  schechter_smf_amplitude_mass_quenched
+  schechter_smf_amplitude_satellite_quenched
 
 
 Reference/API

From 03fd46764e05828da586ba302494b77851832d88 Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Wed, 15 Jun 2022 16:55:34 +0100
Subject: [PATCH 17/36] Fix broadcast issue

---
 skypy/galaxies/stellar_mass.py            |  3 ---
 skypy/galaxies/tests/test_stellar_mass.py | 25 +++++++++++++----------
 2 files changed, 14 insertions(+), 14 deletions(-)

diff --git a/skypy/galaxies/stellar_mass.py b/skypy/galaxies/stellar_mass.py
index f70416d6f..e2dd2bdc3 100644
--- a/skypy/galaxies/stellar_mass.py
+++ b/skypy/galaxies/stellar_mass.py
@@ -144,9 +144,6 @@ def schechter_smf_amplitude_satellites(phi_centrals, fsatellite):
 
     '''
 
-    if np.ndim(phi_centrals) == 2 and np.ndim(fsatellite) == 1:
-        phi_centrals = phi_centrals[:, np.newaxis]
-
     return phi_centrals * np.log(1 / (1 - fsatellite))
 
 
diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py
index 5268b0545..1bab25aac 100644
--- a/skypy/galaxies/tests/test_stellar_mass.py
+++ b/skypy/galaxies/tests/test_stellar_mass.py
@@ -91,30 +91,33 @@ def test_schechter_smf_amplitude_satellites():
     phic_scalar = 10**-2.423
     fsat_scalar = 0.4
 
-    # Array inputs
-    phic_array = np.array([10**-2.423, 10**-2.422])
-    fsat_array = np.array([0.40, 0.41, 0.42])
+    # 1D Array inputs
+    phic_1d = np.array([10**-2.423, 10**-2.422])
+    fsat_1d = np.array([0.40, 0.41, 0.42])
+
+    # 2D Array inputs
+    phic_2d = np.array([[10**-2.50, 10**-2.51, 10**-2.51], [10**-2.50, 10**-2.51, 10**-2.51]])
 
     # Test for scalar output
     phis_scalar = stellar_mass.schechter_smf_amplitude_satellites(phic_scalar, fsat_scalar)
     assert np.isscalar(phis_scalar)
 
     # Test for 1 dim output
-    phis_1d_phib = stellar_mass.schechter_smf_amplitude_satellites(phic_array, fsat_scalar)
-    phis_1d_fsat = stellar_mass.schechter_smf_amplitude_satellites(phic_scalar, fsat_array)
-    assert phis_1d_phib.shape == phic_array.shape
-    assert phis_1d_fsat.shape == fsat_array.shape
+    phis_1d_phib = stellar_mass.schechter_smf_amplitude_satellites(phic_1d, fsat_scalar)
+    phis_1d_fsat = stellar_mass.schechter_smf_amplitude_satellites(phic_scalar, fsat_1d)
+    assert phis_1d_phib.shape == phic_1d.shape
+    assert phis_1d_fsat.shape == fsat_1d.shape
 
     # Test for 2 dim output
-    phis_2d = stellar_mass.schechter_smf_amplitude_satellites(phic_array, fsat_array)
-    assert phis_2d.shape == (len(phic_array), len(fsat_array))
+    phis_2d = stellar_mass.schechter_smf_amplitude_satellites(phic_2d, fsat_1d)
+    assert phis_2d.shape == (len(phic_2d), len(fsat_1d))
 
     # Corner case: no satellite galaxies
     fsat_null = 0
     phis_null_scalar = stellar_mass.schechter_smf_amplitude_satellites(phic_scalar, fsat_null)
-    phis_null_array = stellar_mass.schechter_smf_amplitude_satellites(phic_array, fsat_null)
+    phis_null_1d = stellar_mass.schechter_smf_amplitude_satellites(phic_1d, fsat_null)
     assert phis_null_scalar == 0
-    assert np.all(phis_null_array == np.zeros(len(phic_array)))
+    assert np.all(phis_null_1d == np.zeros(len(phic_1d)))
 
 
 def test_schechter_smf_amplitude_mass_quenched():

From 9a2d528c807ea453bf3a316b7f117c656f07db89 Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Fri, 17 Jun 2022 15:48:04 +0100
Subject: [PATCH 18/36] Add more tests for centrals

---
 skypy/galaxies/tests/test_stellar_mass.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py
index 1bab25aac..73a60dea9 100644
--- a/skypy/galaxies/tests/test_stellar_mass.py
+++ b/skypy/galaxies/tests/test_stellar_mass.py
@@ -85,6 +85,12 @@ def test_schechter_smf_amplitude_centrals():
     phic_2d = stellar_mass.schechter_smf_amplitude_centrals(phiblue_array, fsat_array)
     assert phic_2d.shape == (len(phiblue_array), len(fsat_array))
 
+    # Special case
+    fsat_special = 1 - np.exp(-1)
+    phic_special = stellar_mass.schechter_smf_amplitude_centrals(phiblue_scalar, fsat_special)
+    expectation = 0.5 * phiblue_scalar
+    assert phic_special == expectation
+
 
 def test_schechter_smf_amplitude_satellites():
     # Scalar inputs

From 9e8db46357021b75e1d3c7b1274bdccefdaac69a Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Fri, 17 Jun 2022 15:54:28 +0100
Subject: [PATCH 19/36] Delete redundant variable

---
 skypy/galaxies/tests/test_stellar_mass.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py
index 73a60dea9..a79d78045 100644
--- a/skypy/galaxies/tests/test_stellar_mass.py
+++ b/skypy/galaxies/tests/test_stellar_mass.py
@@ -88,8 +88,7 @@ def test_schechter_smf_amplitude_centrals():
     # Special case
     fsat_special = 1 - np.exp(-1)
     phic_special = stellar_mass.schechter_smf_amplitude_centrals(phiblue_scalar, fsat_special)
-    expectation = 0.5 * phiblue_scalar
-    assert phic_special == expectation
+    assert phic_special == 0.5 * phiblue_scalar
 
 
 def test_schechter_smf_amplitude_satellites():

From 63d92d8d335ff756044066e1dcd2b520bd745a92 Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Fri, 17 Jun 2022 15:55:39 +0100
Subject: [PATCH 20/36] Add extra test for mass quenched function

---
 skypy/galaxies/tests/test_stellar_mass.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py
index a79d78045..89de2f3f2 100644
--- a/skypy/galaxies/tests/test_stellar_mass.py
+++ b/skypy/galaxies/tests/test_stellar_mass.py
@@ -153,6 +153,12 @@ def test_schechter_smf_amplitude_mass_quenched():
     assert phimq_2d.shape == phic_2d.shape == phis_2d.shape
     assert np.all(phimq_2d == phic_2d + phis_2d)
 
+    # Special case
+    phic_null = 0
+    phis_null = 0
+    phimq_null = stellar_mass.schechter_smf_amplitude_mass_quenched(phic_null, phis_null)
+    assert phimq_null == 0
+
 
 def test_schechter_smf_amplitude_satellite_quenched():
     # Scalar inputs

From 84055187945d09953dc36f1e830dda4e1d51e0cc Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Tue, 5 Jul 2022 14:42:58 +0100
Subject: [PATCH 21/36] Change function names from amplitude to phi

---
 skypy/galaxies/stellar_mass.py | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/skypy/galaxies/stellar_mass.py b/skypy/galaxies/stellar_mass.py
index e2dd2bdc3..4829e64f3 100644
--- a/skypy/galaxies/stellar_mass.py
+++ b/skypy/galaxies/stellar_mass.py
@@ -8,10 +8,10 @@
 
 __all__ = [
     'schechter_smf_mass',
-    'schechter_smf_amplitude_centrals',
-    'schechter_smf_amplitude_satellites',
-    'schechter_smf_amplitude_mass_quenched',
-    'schechter_smf_amplitude_satellite_quenched',
+    'schechter_smf_phi_centrals',
+    'schechter_smf_phi_satellites',
+    'schechter_smf_phi_mass_quenched',
+    'schechter_smf_phi_satellite_quenched',
 ]
 
 
@@ -80,7 +80,7 @@ def schechter_smf_mass(redshift, alpha, m_star, m_min, m_max, size=None,
     return m
 
 
-def schechter_smf_amplitude_centrals(phi_blue_total, fsatellite):
+def schechter_smf_phi_centrals(phi_blue_total, fsatellite):
     r'''Schechter amplitude of centrals.
     This function returns the Schechter mass function amplitude
     for active central population based on equation (15)
@@ -116,7 +116,7 @@ def schechter_smf_amplitude_centrals(phi_blue_total, fsatellite):
     return (1 - fsatellite) * phi_blue_total / sum_phics
 
 
-def schechter_smf_amplitude_satellites(phi_centrals, fsatellite):
+def schechter_smf_phi_satellites(phi_centrals, fsatellite):
     r'''Schechter amplitude of satellites.
     This function returns the Schechter mass function amplitude
     for active satellite population based on equation (15)
@@ -147,7 +147,7 @@ def schechter_smf_amplitude_satellites(phi_centrals, fsatellite):
     return phi_centrals * np.log(1 / (1 - fsatellite))
 
 
-def schechter_smf_amplitude_mass_quenched(phi_centrals, phi_satellites):
+def schechter_smf_phi_mass_quenched(phi_centrals, phi_satellites):
     r'''Schechter amplitude of satellites.
     This function returns the Schechter mass function amplitude
     for passive mass-quenched population based on equation (15)
@@ -178,7 +178,7 @@ def schechter_smf_amplitude_mass_quenched(phi_centrals, phi_satellites):
     return phi_centrals + phi_satellites
 
 
-def schechter_smf_amplitude_satellite_quenched(phi_satellites, fenvironment):
+def schechter_smf_phi_satellite_quenched(phi_satellites, fenvironment):
     r'''Schechter amplitude of satellites.
     This function returns the Schechter mass function amplitude
     for active central population based on equation (15)

From 3792b61520c2d3bad7372ca4cb4238bae8fa72f9 Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Tue, 5 Jul 2022 14:45:28 +0100
Subject: [PATCH 22/36] Modify docstrings and types

---
 skypy/galaxies/stellar_mass.py | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/skypy/galaxies/stellar_mass.py b/skypy/galaxies/stellar_mass.py
index 4829e64f3..04634d0cd 100644
--- a/skypy/galaxies/stellar_mass.py
+++ b/skypy/galaxies/stellar_mass.py
@@ -88,17 +88,17 @@ def schechter_smf_phi_centrals(phi_blue_total, fsatellite):
 
     Parameters
     ----------
-    phi_blue_total: float, (nt, ) array_like
+    phi_blue_total: array_like
         Schechter mass function amplitude for the entire active
         sample of galaxies, :math:`(\phi_b, \alpha_b, m_{*})`.
 
-    fsatellite: float, (nm, ) array_like
+    fsatellite: 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.
 
     Returns
     -------
-    amplitude: float, (nt, nm) array_like
+    amplitude: array_like
         Amplitude of the Schechter mass function.
 
      References
@@ -124,7 +124,7 @@ def schechter_smf_phi_satellites(phi_centrals, fsatellite):
 
     Parameters
     ----------
-    phi_centrals: float, (nt, nm) array_like
+    phi_centrals: array_like
         Schechter mass function amplitude of the central
         active galaxies.
 
@@ -134,7 +134,7 @@ def schechter_smf_phi_satellites(phi_centrals, fsatellite):
 
     Returns
     -------
-    amplitude: float, (nt, nm) array_like
+    amplitude: array_like
         Amplitude of the Schechter mass function.
 
     References
@@ -155,17 +155,17 @@ def schechter_smf_phi_mass_quenched(phi_centrals, phi_satellites):
 
     Parameters
     ----------
-    phi_centrals: float, (nt, nm) array_like
+    phi_centrals: array_like
         Schechter mass function amplitude of the central
         active galaxies.
 
-    phi_satellites: float, (nt, nm) array_like
+    phi_satellites: array_like
         Schechter mass function amplitude of the satellite
         active galaxies.
 
     Returns
     -------
-    amplitude: float, (nt, nm) array_like
+    amplitude: array_like
         Amplitude of the Schechter mass function.
 
     References
@@ -186,7 +186,7 @@ def schechter_smf_phi_satellite_quenched(phi_satellites, fenvironment):
 
     Parameters
     ----------
-    phi_satellites: float, (nt, nm) array_like
+    phi_satellites: array_like
         Schechter mass function amplitude of the satellite
         active galaxies.
 
@@ -195,7 +195,7 @@ def schechter_smf_phi_satellite_quenched(phi_satellites, fenvironment):
 
     Returns
     -------
-    amplitude: float, (nt, nm) array_like
+    amplitude: array_like
         Amplitude of the Schechter mass function.
 
     References

From 6c6f7a4cd9da5e643013780ef280a2345d2ef383 Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Tue, 5 Jul 2022 14:48:31 +0100
Subject: [PATCH 23/36] Fix function to allow for numpy broadcastability

---
 skypy/galaxies/stellar_mass.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/skypy/galaxies/stellar_mass.py b/skypy/galaxies/stellar_mass.py
index 04634d0cd..c418bc536 100644
--- a/skypy/galaxies/stellar_mass.py
+++ b/skypy/galaxies/stellar_mass.py
@@ -108,9 +108,6 @@ def schechter_smf_phi_centrals(phi_blue_total, fsatellite):
 
     '''
 
-    if np.ndim(phi_blue_total) == 1 and np.ndim(fsatellite) == 1:
-        phi_blue_total = phi_blue_total[:, np.newaxis]
-
     sum_phics = (1 - fsatellite) * (1 - np.log(1 - fsatellite))
 
     return (1 - fsatellite) * phi_blue_total / sum_phics

From 4711b6336672cff6507f7ebd263822e0419eb5bd Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Tue, 5 Jul 2022 14:54:21 +0100
Subject: [PATCH 24/36] Update function names for tests

---
 skypy/galaxies/tests/test_stellar_mass.py | 46 +++++++++++------------
 1 file changed, 23 insertions(+), 23 deletions(-)

diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py
index 89de2f3f2..c0b4cbdb5 100644
--- a/skypy/galaxies/tests/test_stellar_mass.py
+++ b/skypy/galaxies/tests/test_stellar_mass.py
@@ -62,7 +62,7 @@ def calc_cdf(m):
     assert p_value >= 0.01
 
 
-def test_schechter_smf_amplitude_centrals():
+def test_schechter_smf_phi_centrals():
     # Scalar inputs
     phiblue_scalar = 10**-2.423
     fsat_scalar = 0.4
@@ -72,26 +72,26 @@ def test_schechter_smf_amplitude_centrals():
     fsat_array = np.array([0.40, 0.41, 0.42])
 
     # Test for scalar output
-    phic_scalar = stellar_mass.schechter_smf_amplitude_centrals(phiblue_scalar, fsat_scalar)
+    phic_scalar = stellar_mass.schechter_smf_phi_centrals(phiblue_scalar, fsat_scalar)
     assert np.isscalar(phic_scalar)
 
     # Test for 1 dim output
-    phic_1d_phib = stellar_mass.schechter_smf_amplitude_centrals(phiblue_array, fsat_scalar)
-    phic_1d_fsat = stellar_mass.schechter_smf_amplitude_centrals(phiblue_scalar, fsat_array)
+    phic_1d_phib = stellar_mass.schechter_smf_phi_centrals(phiblue_array, fsat_scalar)
+    phic_1d_fsat = stellar_mass.schechter_smf_phi_centrals(phiblue_scalar, fsat_array)
     assert phic_1d_phib.shape == phiblue_array.shape
     assert phic_1d_fsat.shape == fsat_array.shape
 
     # Test for 2 dim output
-    phic_2d = stellar_mass.schechter_smf_amplitude_centrals(phiblue_array, fsat_array)
+    phic_2d = stellar_mass.schechter_smf_phi_centrals(phiblue_array, fsat_array)
     assert phic_2d.shape == (len(phiblue_array), len(fsat_array))
 
     # Special case
     fsat_special = 1 - np.exp(-1)
-    phic_special = stellar_mass.schechter_smf_amplitude_centrals(phiblue_scalar, fsat_special)
+    phic_special = stellar_mass.schechter_smf_phi_centrals(phiblue_scalar, fsat_special)
     assert phic_special == 0.5 * phiblue_scalar
 
 
-def test_schechter_smf_amplitude_satellites():
+def test_schechter_smf_phi_satellites():
     # Scalar inputs
     phic_scalar = 10**-2.423
     fsat_scalar = 0.4
@@ -104,28 +104,28 @@ def test_schechter_smf_amplitude_satellites():
     phic_2d = np.array([[10**-2.50, 10**-2.51, 10**-2.51], [10**-2.50, 10**-2.51, 10**-2.51]])
 
     # Test for scalar output
-    phis_scalar = stellar_mass.schechter_smf_amplitude_satellites(phic_scalar, fsat_scalar)
+    phis_scalar = stellar_mass.schechter_smf_phi_satellites(phic_scalar, fsat_scalar)
     assert np.isscalar(phis_scalar)
 
     # Test for 1 dim output
-    phis_1d_phib = stellar_mass.schechter_smf_amplitude_satellites(phic_1d, fsat_scalar)
-    phis_1d_fsat = stellar_mass.schechter_smf_amplitude_satellites(phic_scalar, fsat_1d)
+    phis_1d_phib = stellar_mass.schechter_smf_phi_satellites(phic_1d, fsat_scalar)
+    phis_1d_fsat = stellar_mass.schechter_smf_phi_satellites(phic_scalar, fsat_1d)
     assert phis_1d_phib.shape == phic_1d.shape
     assert phis_1d_fsat.shape == fsat_1d.shape
 
     # Test for 2 dim output
-    phis_2d = stellar_mass.schechter_smf_amplitude_satellites(phic_2d, fsat_1d)
+    phis_2d = stellar_mass.schechter_smf_phi_satellites(phic_2d, fsat_1d)
     assert phis_2d.shape == (len(phic_2d), len(fsat_1d))
 
     # Corner case: no satellite galaxies
     fsat_null = 0
-    phis_null_scalar = stellar_mass.schechter_smf_amplitude_satellites(phic_scalar, fsat_null)
-    phis_null_1d = stellar_mass.schechter_smf_amplitude_satellites(phic_1d, fsat_null)
+    phis_null_scalar = stellar_mass.schechter_smf_phi_satellites(phic_scalar, fsat_null)
+    phis_null_1d = stellar_mass.schechter_smf_phi_satellites(phic_1d, fsat_null)
     assert phis_null_scalar == 0
     assert np.all(phis_null_1d == np.zeros(len(phic_1d)))
 
 
-def test_schechter_smf_amplitude_mass_quenched():
+def test_schechter_smf_phi_mass_quenched():
     # Scalar inputs
     phic_scalar = 10**-2.5
     phis_scalar = 10**-2.4
@@ -139,28 +139,28 @@ def test_schechter_smf_amplitude_mass_quenched():
     phis_2d = np.array([[10**-2.40, 10**-2.41], [10**-2.40, 10**-2.41]])
 
     # Test for scalar output
-    phimq_scalar = stellar_mass.schechter_smf_amplitude_mass_quenched(phic_scalar, phis_scalar)
+    phimq_scalar = stellar_mass.schechter_smf_phi_mass_quenched(phic_scalar, phis_scalar)
     assert np.isscalar(phis_scalar)
     assert phimq_scalar == phic_scalar + phis_scalar
 
     # Test for 1 dim output
-    phimq_1d = stellar_mass.schechter_smf_amplitude_mass_quenched(phic_1d, phis_1d)
+    phimq_1d = stellar_mass.schechter_smf_phi_mass_quenched(phic_1d, phis_1d)
     assert phimq_1d.shape == phic_1d.shape == phis_1d.shape
     assert np.all(phimq_1d == phic_1d + phis_1d)
 
     # Test for 2 dim output
-    phimq_2d = stellar_mass.schechter_smf_amplitude_mass_quenched(phic_2d, phis_2d)
+    phimq_2d = stellar_mass.schechter_smf_phi_mass_quenched(phic_2d, phis_2d)
     assert phimq_2d.shape == phic_2d.shape == phis_2d.shape
     assert np.all(phimq_2d == phic_2d + phis_2d)
 
     # Special case
     phic_null = 0
     phis_null = 0
-    phimq_null = stellar_mass.schechter_smf_amplitude_mass_quenched(phic_null, phis_null)
+    phimq_null = stellar_mass.schechter_smf_phi_mass_quenched(phic_null, phis_null)
     assert phimq_null == 0
 
 
-def test_schechter_smf_amplitude_satellite_quenched():
+def test_schechter_smf_phi_satellite_quenched():
     # Scalar inputs
     phis_scalar = 10**-2.423
     frho_scalar = 0.2
@@ -169,16 +169,16 @@ def test_schechter_smf_amplitude_satellite_quenched():
     phis_array = np.array([10**-2.423, 10**-2.422])
 
     # Test for scalar output
-    phisq_scalar = stellar_mass.schechter_smf_amplitude_satellite_quenched(phis_scalar, frho_scalar)
+    phisq_scalar = stellar_mass.schechter_smf_phi_satellite_quenched(phis_scalar, frho_scalar)
     assert np.isscalar(phisq_scalar)
 
     # Test for 1 dim output
-    phisq_1d_phib = stellar_mass.schechter_smf_amplitude_satellite_quenched(phis_array, frho_scalar)
+    phisq_1d_phib = stellar_mass.schechter_smf_phi_satellite_quenched(phis_array, frho_scalar)
     assert phisq_1d_phib.shape == phis_array.shape
 
     # Corner case no satellite-quenched galaxies
     frho_null = 0
-    phisq_null_scalar = stellar_mass.schechter_smf_amplitude_satellites(phis_scalar, frho_null)
-    phisq_null_array = stellar_mass.schechter_smf_amplitude_satellites(phis_array, frho_null)
+    phisq_null_scalar = stellar_mass.schechter_smf_phi_satellites(phis_scalar, frho_null)
+    phisq_null_array = stellar_mass.schechter_smf_phi_satellites(phis_array, frho_null)
     assert phisq_null_scalar == 0
     assert np.all(phisq_null_array == np.zeros(len(phis_array)))

From 2358e197bde6c83da5ce099e04517b27f2c44a43 Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Tue, 5 Jul 2022 15:06:38 +0100
Subject: [PATCH 25/36] Modify satellite functions by defining common
 functionality

---
 skypy/galaxies/stellar_mass.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/skypy/galaxies/stellar_mass.py b/skypy/galaxies/stellar_mass.py
index c418bc536..b4b28e830 100644
--- a/skypy/galaxies/stellar_mass.py
+++ b/skypy/galaxies/stellar_mass.py
@@ -141,7 +141,7 @@ def schechter_smf_phi_satellites(phi_centrals, fsatellite):
 
     '''
 
-    return phi_centrals * np.log(1 / (1 - fsatellite))
+    return _satellite_computation(phi_centrals, fsatellite)
 
 
 def schechter_smf_phi_mass_quenched(phi_centrals, phi_satellites):
@@ -202,4 +202,8 @@ def schechter_smf_phi_satellite_quenched(phi_satellites, fenvironment):
 
     '''
 
-    return - np.log(1 - fenvironment) * phi_satellites
+    return _satellite_computation(phi_satellites, fenvironment)
+
+
+def _satellite_computation(amplitude, fraction):
+    return - np.log(1 - fraction) * amplitude

From cda6c68b226c0742bfcf3ea3877a6108610fd20e Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Tue, 5 Jul 2022 15:24:58 +0100
Subject: [PATCH 26/36] Merge the satellite related tests into a single one

---
 skypy/galaxies/tests/test_stellar_mass.py | 67 +++++++----------------
 1 file changed, 20 insertions(+), 47 deletions(-)

diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py
index c0b4cbdb5..3c4b420e8 100644
--- a/skypy/galaxies/tests/test_stellar_mass.py
+++ b/skypy/galaxies/tests/test_stellar_mass.py
@@ -91,40 +91,6 @@ def test_schechter_smf_phi_centrals():
     assert phic_special == 0.5 * phiblue_scalar
 
 
-def test_schechter_smf_phi_satellites():
-    # Scalar inputs
-    phic_scalar = 10**-2.423
-    fsat_scalar = 0.4
-
-    # 1D Array inputs
-    phic_1d = np.array([10**-2.423, 10**-2.422])
-    fsat_1d = np.array([0.40, 0.41, 0.42])
-
-    # 2D Array inputs
-    phic_2d = np.array([[10**-2.50, 10**-2.51, 10**-2.51], [10**-2.50, 10**-2.51, 10**-2.51]])
-
-    # Test for scalar output
-    phis_scalar = stellar_mass.schechter_smf_phi_satellites(phic_scalar, fsat_scalar)
-    assert np.isscalar(phis_scalar)
-
-    # Test for 1 dim output
-    phis_1d_phib = stellar_mass.schechter_smf_phi_satellites(phic_1d, fsat_scalar)
-    phis_1d_fsat = stellar_mass.schechter_smf_phi_satellites(phic_scalar, fsat_1d)
-    assert phis_1d_phib.shape == phic_1d.shape
-    assert phis_1d_fsat.shape == fsat_1d.shape
-
-    # Test for 2 dim output
-    phis_2d = stellar_mass.schechter_smf_phi_satellites(phic_2d, fsat_1d)
-    assert phis_2d.shape == (len(phic_2d), len(fsat_1d))
-
-    # Corner case: no satellite galaxies
-    fsat_null = 0
-    phis_null_scalar = stellar_mass.schechter_smf_phi_satellites(phic_scalar, fsat_null)
-    phis_null_1d = stellar_mass.schechter_smf_phi_satellites(phic_1d, fsat_null)
-    assert phis_null_scalar == 0
-    assert np.all(phis_null_1d == np.zeros(len(phic_1d)))
-
-
 def test_schechter_smf_phi_mass_quenched():
     # Scalar inputs
     phic_scalar = 10**-2.5
@@ -160,25 +126,32 @@ def test_schechter_smf_phi_mass_quenched():
     assert phimq_null == 0
 
 
-def test_schechter_smf_phi_satellite_quenched():
+SATELLITE_FUNCTIONS = [
+    stellar_mass.schechter_smf_phi_satellites,
+    stellar_mass.schechter_smf_phi_satellite_quenched,
+]
+
+
+@pytest.mark.parametrize('satellite_function', SATELLITE_FUNCTIONS)
+def test_schechter_smf_phi_satellites_common(satellite_function):
     # Scalar inputs
     phis_scalar = 10**-2.423
-    frho_scalar = 0.2
+    fraction_scalar = 0.2
 
     # Array inputs
     phis_array = np.array([10**-2.423, 10**-2.422])
 
     # Test for scalar output
-    phisq_scalar = stellar_mass.schechter_smf_phi_satellite_quenched(phis_scalar, frho_scalar)
-    assert np.isscalar(phisq_scalar)
+    phis_sat_scalar = satellite_function(phis_scalar, fraction_scalar)
+    assert np.isscalar(phis_sat_scalar)
 
     # Test for 1 dim output
-    phisq_1d_phib = stellar_mass.schechter_smf_phi_satellite_quenched(phis_array, frho_scalar)
-    assert phisq_1d_phib.shape == phis_array.shape
-
-    # Corner case no satellite-quenched galaxies
-    frho_null = 0
-    phisq_null_scalar = stellar_mass.schechter_smf_phi_satellites(phis_scalar, frho_null)
-    phisq_null_array = stellar_mass.schechter_smf_phi_satellites(phis_array, frho_null)
-    assert phisq_null_scalar == 0
-    assert np.all(phisq_null_array == np.zeros(len(phis_array)))
+    phis_sat_1d_phib = satellite_function(phis_array, fraction_scalar)
+    assert phis_sat_1d_phib.shape == phis_array.shape
+
+    # Corner case no satellite galaxies
+    fraction_null = 0
+    phis_sat_null_scalar = satellite_function(phis_scalar, fraction_null)
+    phis_sat_null_array = satellite_function(phis_array, fraction_null)
+    assert phis_sat_null_scalar == 0
+    assert np.all(phis_sat_null_array == np.zeros(len(phis_array)))

From dbfb875d7bb54a1c83abc8814355c78ae4192bf8 Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Tue, 5 Jul 2022 15:29:47 +0100
Subject: [PATCH 27/36] Add broadcastable example test

---
 skypy/galaxies/tests/test_stellar_mass.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py
index 3c4b420e8..4ebce8bf1 100644
--- a/skypy/galaxies/tests/test_stellar_mass.py
+++ b/skypy/galaxies/tests/test_stellar_mass.py
@@ -82,7 +82,7 @@ def test_schechter_smf_phi_centrals():
     assert phic_1d_fsat.shape == fsat_array.shape
 
     # Test for 2 dim output
-    phic_2d = stellar_mass.schechter_smf_phi_centrals(phiblue_array, fsat_array)
+    phic_2d = stellar_mass.schechter_smf_phi_centrals(phiblue_array[:, np.newaxis], fsat_array)
     assert phic_2d.shape == (len(phiblue_array), len(fsat_array))
 
     # Special case

From 4a19bd58642c9a7c53a28f35ab2416e65761b3eb Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Tue, 5 Jul 2022 16:08:58 +0100
Subject: [PATCH 28/36] Modify test for mass quenched galaxies by usng the
 hypothesis package

---
 skypy/galaxies/tests/test_stellar_mass.py | 37 +++++++----------------
 1 file changed, 11 insertions(+), 26 deletions(-)

diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py
index 4ebce8bf1..0280ebf02 100644
--- a/skypy/galaxies/tests/test_stellar_mass.py
+++ b/skypy/galaxies/tests/test_stellar_mass.py
@@ -3,6 +3,8 @@
 import scipy.integrate
 from scipy.special import gammaln
 import pytest
+from hypothesis import given
+from hypothesis.strategies import integers
 
 from skypy.galaxies import stellar_mass
 from skypy.utils import special
@@ -91,40 +93,23 @@ def test_schechter_smf_phi_centrals():
     assert phic_special == 0.5 * phiblue_scalar
 
 
-def test_schechter_smf_phi_mass_quenched():
-    # Scalar inputs
-    phic_scalar = 10**-2.5
-    phis_scalar = 10**-2.4
-
-    # 1D Array inputs
-    phic_1d = np.array([10**-2.50, 10**-2.51])
-    phis_1d = np.array([10**-2.40, 10**-2.41])
+@given(integers(), integers())
+def test_schechter_smf_phi_mass_quenched(phic, phis):
 
-    # 2D Array inputs
-    phic_2d = np.array([[10**-2.50, 10**-2.51], [10**-2.50, 10**-2.51]])
-    phis_2d = np.array([[10**-2.40, 10**-2.41], [10**-2.40, 10**-2.41]])
+    # Array inputs
+    phic_1d = np.array([phic, phic])
+    phis_1d = np.array([phis, phis])
 
     # Test for scalar output
-    phimq_scalar = stellar_mass.schechter_smf_phi_mass_quenched(phic_scalar, phis_scalar)
-    assert np.isscalar(phis_scalar)
-    assert phimq_scalar == phic_scalar + phis_scalar
+    phimq_scalar = stellar_mass.schechter_smf_phi_mass_quenched(phic, phis)
+    assert np.isscalar(phimq_scalar)
+    assert phimq_scalar == phic + phis
 
-    # Test for 1 dim output
+    # Test for array output
     phimq_1d = stellar_mass.schechter_smf_phi_mass_quenched(phic_1d, phis_1d)
     assert phimq_1d.shape == phic_1d.shape == phis_1d.shape
     assert np.all(phimq_1d == phic_1d + phis_1d)
 
-    # Test for 2 dim output
-    phimq_2d = stellar_mass.schechter_smf_phi_mass_quenched(phic_2d, phis_2d)
-    assert phimq_2d.shape == phic_2d.shape == phis_2d.shape
-    assert np.all(phimq_2d == phic_2d + phis_2d)
-
-    # Special case
-    phic_null = 0
-    phis_null = 0
-    phimq_null = stellar_mass.schechter_smf_phi_mass_quenched(phic_null, phis_null)
-    assert phimq_null == 0
-
 
 SATELLITE_FUNCTIONS = [
     stellar_mass.schechter_smf_phi_satellites,

From 54621710fc03dcf6a4663481c33af6c692f3b705 Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Tue, 5 Jul 2022 16:17:52 +0100
Subject: [PATCH 29/36] Update docs file

---
 docs/galaxies.rst | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/docs/galaxies.rst b/docs/galaxies.rst
index c31a36a38..9d29885cb 100644
--- a/docs/galaxies.rst
+++ b/docs/galaxies.rst
@@ -89,10 +89,10 @@ The following models are found in the `skypy.galaxies.stellar_mass` package.
   :nosignatures:
 
   schechter_smf_mass
-  schechter_smf_amplitude_centrals
-  schechter_smf_amplitude_satellites
-  schechter_smf_amplitude_mass_quenched
-  schechter_smf_amplitude_satellite_quenched
+  schechter_smf_phi_centrals
+  schechter_smf_phi_satellites
+  schechter_smf_phi_mass_quenched
+  schechter_smf_phi_satellite_quenched
 
 
 Reference/API

From d0b5aef27fb8e05b1fff7622061b4ea1c0d4a11b Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Tue, 5 Jul 2022 16:21:45 +0100
Subject: [PATCH 30/36] Fix title of functions

---
 skypy/galaxies/stellar_mass.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/skypy/galaxies/stellar_mass.py b/skypy/galaxies/stellar_mass.py
index b4b28e830..744e47970 100644
--- a/skypy/galaxies/stellar_mass.py
+++ b/skypy/galaxies/stellar_mass.py
@@ -145,7 +145,7 @@ def schechter_smf_phi_satellites(phi_centrals, fsatellite):
 
 
 def schechter_smf_phi_mass_quenched(phi_centrals, phi_satellites):
-    r'''Schechter amplitude of satellites.
+    r'''Schechter amplitude of mass-quenched galaxies.
     This function returns the Schechter mass function amplitude
     for passive mass-quenched population based on equation (15)
     in de la Bella et al. 2021 [1]_.
@@ -176,7 +176,7 @@ def schechter_smf_phi_mass_quenched(phi_centrals, phi_satellites):
 
 
 def schechter_smf_phi_satellite_quenched(phi_satellites, fenvironment):
-    r'''Schechter amplitude of satellites.
+    r'''Schechter amplitude of satellite-quenched galaxies.
     This function returns the Schechter mass function amplitude
     for active central population based on equation (15)
     in de la Bella et al. 2021 [1]_.

From 7f70ffde4055af2f8a1b7660a2d551f571a9bb37 Mon Sep 17 00:00:00 2001
From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com>
Date: Tue, 5 Jul 2022 16:25:18 +0100
Subject: [PATCH 31/36] Fix title of functions for active galaxies

---
 skypy/galaxies/stellar_mass.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/skypy/galaxies/stellar_mass.py b/skypy/galaxies/stellar_mass.py
index 744e47970..433b3b67a 100644
--- a/skypy/galaxies/stellar_mass.py
+++ b/skypy/galaxies/stellar_mass.py
@@ -81,7 +81,7 @@ def schechter_smf_mass(redshift, alpha, m_star, m_min, m_max, size=None,
 
 
 def schechter_smf_phi_centrals(phi_blue_total, fsatellite):
-    r'''Schechter amplitude of centrals.
+    r'''Schechter amplitude of central galaxies.
     This function returns the Schechter mass function amplitude
     for active central population based on equation (15)
     in de la Bella et al. 2021 [1]_.
@@ -114,7 +114,7 @@ def schechter_smf_phi_centrals(phi_blue_total, fsatellite):
 
 
 def schechter_smf_phi_satellites(phi_centrals, fsatellite):
-    r'''Schechter amplitude of satellites.
+    r'''Schechter amplitude of satellite galaxies.
     This function returns the Schechter mass function amplitude
     for active satellite population based on equation (15)
     in de la Bella et al. 2021 [1]_.

From 8ea593e01f8861584c6fba60edf4d0fff28a3937 Mon Sep 17 00:00:00 2001
From: Ian Harrison <itrharrison@gmail.com>
Date: Thu, 18 Aug 2022 17:31:35 +0100
Subject: [PATCH 32/36] Update test_skypy.py

E275 fixes
---
 skypy/pipeline/tests/test_skypy.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/skypy/pipeline/tests/test_skypy.py b/skypy/pipeline/tests/test_skypy.py
index f84a5ea55..988b01a64 100644
--- a/skypy/pipeline/tests/test_skypy.py
+++ b/skypy/pipeline/tests/test_skypy.py
@@ -58,7 +58,7 @@ def test_logging(capsys, tmp_path):
     output_filename = str(tmp_path / 'logging.fits')
     skypy.main([config_filename, output_filename])
     out, err = capsys.readouterr()
-    assert(not err)
+    assert (not err)
 
     # Run again with increased verbosity and capture log. Force an exception by
     # not using the "--overwrite" flag when the output file already exists.
@@ -78,19 +78,19 @@ def test_logging(capsys, tmp_path):
     # Check all jobs appear in the log
     for job in list(config) + list(tables) + columns:
         log_string = f"[INFO] skypy.pipeline: Generating {job}"
-        assert(log_string in err)
+        assert (log_string in err)
 
     # Check all functions appear in the log
     for f in functions:
         log_string = f"[INFO] skypy.pipeline: Calling {f.function.__name__}"
-        assert(log_string in err)
+        assert (log_string in err)
 
     # Check cosmology appears in the log
     if cosmology:
-        assert("[INFO] skypy.pipeline: Setting cosmology" in err)
+        assert ("[INFO] skypy.pipeline: Setting cosmology" in err)
 
     # Check writing output file is in the log
-    assert(f"[INFO] skypy: Writing {output_filename}" in err)
+    assert (f"[INFO] skypy: Writing {output_filename}" in err)
 
     # Check error for existing output file is in the log
     try:
@@ -100,10 +100,10 @@ def test_logging(capsys, tmp_path):
     except ImportError:
         # Fallback on old error message from astropy v4.x
         error_string = f"[ERROR] skypy: File {output_filename!r} already exists."
-    assert(error_string in err)
+    assert (error_string in err)
 
     # Run again with decreased verbosity and check the log is empty
     with pytest.raises(SystemExit):
         skypy.main([config_filename, output_filename, '-qq'])
     out, err = capsys.readouterr()
-    assert(not err)
+    assert (not err)

From 31be5062c167d8fc432c2b7cb79d44c053fa0858 Mon Sep 17 00:00:00 2001
From: Ian Harrison <itrharrison@gmail.com>
Date: Thu, 18 Aug 2022 17:32:16 +0100
Subject: [PATCH 33/36] Update skypy.py

E275 fixes (2)
---
 skypy/pipeline/scripts/skypy.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skypy/pipeline/scripts/skypy.py b/skypy/pipeline/scripts/skypy.py
index 62319de34..abc35364d 100644
--- a/skypy/pipeline/scripts/skypy.py
+++ b/skypy/pipeline/scripts/skypy.py
@@ -48,4 +48,4 @@ def main(args=None):
         logger.exception(e)
         raise SystemExit(2) from e
 
-    return(0)
+    return (0)

From f07ad416665a53d31c469addfea4039e509a06bb Mon Sep 17 00:00:00 2001
From: Ian Harrison <itrharrison@gmail.com>
Date: Thu, 18 Aug 2022 17:32:48 +0100
Subject: [PATCH 34/36] Update photometry.py

E275 fixes (3)
---
 skypy/utils/photometry.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skypy/utils/photometry.py b/skypy/utils/photometry.py
index b6df107b1..3262a1b42 100644
--- a/skypy/utils/photometry.py
+++ b/skypy/utils/photometry.py
@@ -131,7 +131,7 @@ def mag_ab(wavelength, spectrum, filters, *, redshift=None, coefficients=None,
         u = u.reshape(u.shape + (1,)*(nd_s+nd_f))
         m = np.ascontiguousarray(m[n])
         m += u*dm[n]
-        del(dm, n, u)
+        del (dm, n, u)
 
     # combine spectra if asked to
     if coefficients is not None:

From c09748cfacd168119ed3a030965c77ecb0f599db Mon Sep 17 00:00:00 2001
From: Ian Harrison <itrharrison@gmail.com>
Date: Thu, 25 Aug 2022 10:44:21 +0100
Subject: [PATCH 35/36] revert flake8 failures

---
 skypy/pipeline/scripts/skypy.py    |  2 +-
 skypy/pipeline/tests/test_skypy.py | 14 +++++++-------
 skypy/utils/photometry.py          |  2 +-
 3 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/skypy/pipeline/scripts/skypy.py b/skypy/pipeline/scripts/skypy.py
index abc35364d..62319de34 100644
--- a/skypy/pipeline/scripts/skypy.py
+++ b/skypy/pipeline/scripts/skypy.py
@@ -48,4 +48,4 @@ def main(args=None):
         logger.exception(e)
         raise SystemExit(2) from e
 
-    return (0)
+    return(0)
diff --git a/skypy/pipeline/tests/test_skypy.py b/skypy/pipeline/tests/test_skypy.py
index 988b01a64..f84a5ea55 100644
--- a/skypy/pipeline/tests/test_skypy.py
+++ b/skypy/pipeline/tests/test_skypy.py
@@ -58,7 +58,7 @@ def test_logging(capsys, tmp_path):
     output_filename = str(tmp_path / 'logging.fits')
     skypy.main([config_filename, output_filename])
     out, err = capsys.readouterr()
-    assert (not err)
+    assert(not err)
 
     # Run again with increased verbosity and capture log. Force an exception by
     # not using the "--overwrite" flag when the output file already exists.
@@ -78,19 +78,19 @@ def test_logging(capsys, tmp_path):
     # Check all jobs appear in the log
     for job in list(config) + list(tables) + columns:
         log_string = f"[INFO] skypy.pipeline: Generating {job}"
-        assert (log_string in err)
+        assert(log_string in err)
 
     # Check all functions appear in the log
     for f in functions:
         log_string = f"[INFO] skypy.pipeline: Calling {f.function.__name__}"
-        assert (log_string in err)
+        assert(log_string in err)
 
     # Check cosmology appears in the log
     if cosmology:
-        assert ("[INFO] skypy.pipeline: Setting cosmology" in err)
+        assert("[INFO] skypy.pipeline: Setting cosmology" in err)
 
     # Check writing output file is in the log
-    assert (f"[INFO] skypy: Writing {output_filename}" in err)
+    assert(f"[INFO] skypy: Writing {output_filename}" in err)
 
     # Check error for existing output file is in the log
     try:
@@ -100,10 +100,10 @@ def test_logging(capsys, tmp_path):
     except ImportError:
         # Fallback on old error message from astropy v4.x
         error_string = f"[ERROR] skypy: File {output_filename!r} already exists."
-    assert (error_string in err)
+    assert(error_string in err)
 
     # Run again with decreased verbosity and check the log is empty
     with pytest.raises(SystemExit):
         skypy.main([config_filename, output_filename, '-qq'])
     out, err = capsys.readouterr()
-    assert (not err)
+    assert(not err)
diff --git a/skypy/utils/photometry.py b/skypy/utils/photometry.py
index 3262a1b42..b6df107b1 100644
--- a/skypy/utils/photometry.py
+++ b/skypy/utils/photometry.py
@@ -131,7 +131,7 @@ def mag_ab(wavelength, spectrum, filters, *, redshift=None, coefficients=None,
         u = u.reshape(u.shape + (1,)*(nd_s+nd_f))
         m = np.ascontiguousarray(m[n])
         m += u*dm[n]
-        del (dm, n, u)
+        del(dm, n, u)
 
     # combine spectra if asked to
     if coefficients is not None:

From a47b0eff2204e12e2bc73cf3910f5300c2d5fadc Mon Sep 17 00:00:00 2001
From: Lucia Fonseca de la Bella <lfdelabella@Lucias-MacBook-Pro.local>
Date: Mon, 16 Jan 2023 15:51:34 +0000
Subject: [PATCH 36/36] Add corner cases

---
 skypy/galaxies/tests/test_stellar_mass.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py
index c064fdbd3..dc5e898ba 100644
--- a/skypy/galaxies/tests/test_stellar_mass.py
+++ b/skypy/galaxies/tests/test_stellar_mass.py
@@ -130,6 +130,12 @@ def test_schechter_smf_phi_mass_quenched(phic, phis):
     assert phimq_1d.shape == phic_1d.shape == phis_1d.shape
     assert np.all(phimq_1d == phic_1d + phis_1d)
 
+    # Corner cases
+    phi0 = 1
+    assert phi0 == 0.5 * stellar_mass.schechter_smf_phi_mass_quenched(phi0, phi0)
+    assert 0.0 == stellar_mass.schechter_smf_phi_mass_quenched(phi0, -phi0)
+    assert phi0 == stellar_mass.schechter_smf_phi_mass_quenched(phi0, 0.0)
+
 
 SATELLITE_FUNCTIONS = [
     stellar_mass.schechter_smf_phi_satellites,