Repository: skypyproject/skypy (python)

Issue:
ENH: Logistic Completeness function
## Description
This function returns the logistic completeness function:

                  1
    p(m) = -----------------
           1 + exp{-k(m-m0)}

that describes the probability `p(m)` that a galaxy of magnitude `m` is detected in the band.

See also the discussion around eq. (7) and fig. 2 in the ALHAMBRA paper by Lopez-Sanjuan et al. (2017).

## Inputs
 - apparent magnitudes
 - 95% completeness
 - 50% completeness

## Outputs
 - Logistic completeness function

## References
- [Lopez-Sanjuan et al. 2017](https://www.aanda.org/articles/aa/full_html/2017/03/aa29517-16/aa29517-16.html)


Pull request description:
ENH: Logistic completeness function
## Description
This PR adds a feature to calculate the logistic completeness function. Closes #519.

## Checklist
- [x] Follow the [Contributor Guidelines](https://github.com/skypyproject/skypy/blob/main/CONTRIBUTING.rst)
- [x] Write unit tests
- [x] Write documentation strings
- [x] Assign someone from your working team to review this pull request
- [x] Assign someone from the infrastructure team to review this pull request


Patch:
diff --git a/docs/utils/index.rst b/docs/utils/index.rst
index 3545bb6..090453e 100644
--- a/docs/utils/index.rst
+++ b/docs/utils/index.rst
@@ -144,6 +144,7 @@ create their own named `speclite.filters.FilterResponse`.
    mag_ab
    SpectrumTemplates
    magnitude_error_rykoff
+   logistic_completeness_function
 
 
 Random sampling (`skypy.utils.random`)
diff --git a/skypy/utils/photometry.py b/skypy/utils/photometry.py
index 1f79e2c..e596559 100644
--- a/skypy/utils/photometry.py
+++ b/skypy/utils/photometry.py
@@ -5,6 +5,7 @@
 
 from abc import ABCMeta, abstractmethod
 import numpy as np
+import scipy.special
 
 
 __all__ = [
@@ -14,6 +15,7 @@ __all__ = [
     'mag_ab',
     'SpectrumTemplates',
     'magnitude_error_rykoff',
+    'logistic_completeness_function',
 ]
 
 try:
@@ -392,3 +394,54 @@ def magnitude_error_rykoff(magnitude, magnitude_limit, magnitude_zp, a, b, error
     error = 2.5 / np.log(10) * np.sqrt((1 + flux_noise / flux) / (flux * t_eff))
 
     return np.minimum(error, error_limit)
+
+
+def logistic_completeness_function(magnitude, magnitude_95, magnitude_50):
+    r'''Logistic completeness function.
+
+    This function calculates the logistic completeness function (based on eq. (7) in
+    López-Sanjuan, C. et al. (2017) [1]_.
+
+    .. math::
+
+        p(m) = \frac{1}{1 + \exp[\kappa (m - m_{50})]}\;,
+
+    which describes the probability :math:`p(m)` that an object of magnitude :math:`m` is detected
+    in a specific band and with
+
+    .. math::
+
+        \kappa = \frac{\ln(\frac{1}{19})}{m_{95} - m_{50}}\;.
+
+    Here, :math:`m_{95}` and :math:`m_{50}` are the 95% and 50% completeness
+    magnitudes, respectively.
+
+    Parameters
+    ----------
+    magnitude : array_like
+        Magnitudes. Can be multidimensional for computing with multiple filter bands.
+    magnitude_95 : scalar or 1-D array_like
+        95% completeness magnitude.
+        If `magnitude_50` is 1-D array it has to be scalar or 1-D array of the same shape.
+    magnitude_50 : scalar or 1-D array_like
+        50% completeness magnitude.
+        If `magnitude_95` is 1-D array it has to be scalar or 1-D array of the same shape.
+
+    Returns
+    -------
+    probability : scalar or array_like
+        Probability of detecting an object with magnitude :math:`m`.
+        Returns array_like of the same shape as magnitude.
+        Exemption: If magnitude is scalar and `magnitude_95` or `magnitude_50`
+        is array_like of shape (nb, ) it returns array_like of shape (nb, ).
+
+    References
+    -----------
+    .. [1] López-Sanjuan, C. et al., `2017A&A…599A..62L`_
+    .. _2017A&A…599A..62L: https://ui.adsabs.harvard.edu/abs/2017A%26A...599A..62L
+
+    '''
+
+    kappa = np.log(1. / 19) / np.subtract(magnitude_95, magnitude_50)
+    arg = kappa * np.subtract(magnitude, magnitude_50)
+    return scipy.special.expit(-arg)


Test patch:
diff --git a/skypy/utils/tests/test_photometry.py b/skypy/utils/tests/test_photometry.py
index 560e2a8..c6ec04d 100644
--- a/skypy/utils/tests/test_photometry.py
+++ b/skypy/utils/tests/test_photometry.py
@@ -212,3 +212,21 @@ def test_magnitude_error_rykoff():
     error_limit = 1
     error = magnitude_error_rykoff(magnitude, magnitude_limit, magnitude_zp, a, b, error_limit)
     assert error == error_limit
+
+
+def test_logistic_completeness_function():
+    from skypy.utils.photometry import logistic_completeness_function
+
+    # Test that arguments broadcast correctly
+    m = np.full((2, 1, 1), 21)
+    m95 = np.full((3, 1), 22)
+    m50 = np.full(5, 23)
+    p = logistic_completeness_function(m, m95, m50)
+    assert p.shape == np.broadcast(m, m95, m50).shape
+
+    # Test result of completeness function for different given magnitudes
+    m95 = 24
+    m50 = 25
+    m = [np.finfo(np.float64).min, m95, m50, 2*m50-m95, np.finfo(np.float64).max]
+    p = logistic_completeness_function(m, m95, m50)
+    assert np.allclose(p, [1, 0.95, 0.5, 0.05, 0])