Repository: skypyproject/skypy (python)

Issue:
ENH: Rykoff apparent magnitude uncertainty of galaxies
## Description
Add a function that returns the uncertainty of the galaxies' apparent magnitude in the given filter bands as described by equation (6) in Rykoff et al. (2015). Note: the normalisation constant k is absorbed in the included in parameter t_eff

## Inputs
 - Apparent magnitudes
 - 10 sigma apparent magnitude limits
 - Parameters a and b of the model
 - Uncertainty limit

## Outputs
 - Apparent magnitude uncertainty of galaxies

## References
- Rykoff et al. (2015), Assessing galaxy limiting magnitudes in large optical surveys, [arXiv 1509.00870](https://arxiv.org/abs/1509.00870)

Pull request:
ENH: Rykoff model of the magnitude uncertainty
## Description
This PR adds a feature to calculate the magnitude uncertainty following the model of Rykoff et al. (2015). Closes #517.
## 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 ac33e19..3545bb6 100644
--- a/docs/utils/index.rst
+++ b/docs/utils/index.rst
@@ -143,6 +143,7 @@ create their own named `speclite.filters.FilterResponse`.
    luminosity_in_band
    mag_ab
    SpectrumTemplates
+   magnitude_error_rykoff
 
 
 Random sampling (`skypy.utils.random`)
diff --git a/skypy/utils/photometry.py b/skypy/utils/photometry.py
index 8e978e3..1f79e2c 100644
--- a/skypy/utils/photometry.py
+++ b/skypy/utils/photometry.py
@@ -13,6 +13,7 @@ __all__ = [
     'luminosity_in_band',
     'mag_ab',
     'SpectrumTemplates',
+    'magnitude_error_rykoff',
 ]
 
 try:
@@ -304,3 +305,90 @@ def absolute_magnitude_from_luminosity(luminosity, zeropoint=None):
         zeropoint = -luminosity_in_band[zeropoint]
 
     return -2.5*np.log10(luminosity) - zeropoint
+
+
+def magnitude_error_rykoff(magnitude, magnitude_limit, magnitude_zp, a, b, error_limit=np.inf):
+    r"""Magnitude error acoording to the model from Rykoff et al. (2015).
+
+    Given an apparent magnitude calculate the magnitude error that is introduced
+    by the survey specifications and follows the model described in Rykoff et al. (2015).
+
+    Parameters
+    ----------
+    magnitude: array_like
+        Apparent magnitude. This and the other array_like parameters must
+        be broadcastable to the same shape.
+    magnitude_limit: array_like
+        :math:`10\sigma` limiting magnitude of the survey. This and the other
+        array_like parameters must be broadcastable to the same shape.
+    magnitude_zp: array_like
+        Zero-point magnitude of the survey. This and the other array_like parameters must
+        be broadcastable to the same shape.
+    a,b: array_like
+        Model parameters: a is the intercept and
+        b is the slope of the logarithmic effective time.
+        These and the other array_like parameters must be broadcastable to the same shape.
+    error_limit: float, optional
+        Upper limit of the returned error. If given, all values larger than this value
+        will be set to error_limit. Default is None.
+
+    Returns
+    -------
+    error: ndarray
+        The apparent magnitude error in the Rykoff et al. (2015) model. This is a scalar
+        if magnitude, magnitude_limit, magnitude_zp, a and b are scalars.
+
+    Notes
+    -----
+    Rykoff et al. (2015) (see [1]_) describe the error of the apparent magnitude :math:`m` as
+
+    .. math::
+
+        \sigma_m(m;m_{\mathrm{lim}}, t_{\mathrm{eff}}) &=
+            \sigma_m(F(m);F_{\mathrm{noise}}(m_{\mathrm{lim}}), t_{\mathrm{eff}}) \\
+        &= \frac{2.5}{\ln(10)} \left[ \frac{1}{Ft_{\mathrm{eff}}}
+            \left( 1 + \frac{F_{\mathrm{noise}}}{F} \right) \right]^{1/2} \;,
+
+    where
+
+    .. math::
+
+        F=10^{-0.4(m - m_{\mathrm{ZP}})}
+
+    is the source's flux,
+
+    .. math::
+
+        F_\mathrm{noise} = \frac{F_{\mathrm{lim}}^2 t_{\mathrm{eff}}}{10^2} - F_{\mathrm{lim}}
+
+    is the effective noise flux and :math:`t_\mathrm{eff}` is the effective exposure time
+    (we absorbed the normalisation constant :math:`k` in the definition of
+    :math:`t_\mathrm{eff}`).
+    Furthermore, :math:`m_\mathrm{ZP}` is the zero-point magnitude of the survey and
+    :math:`F_\mathrm{lim}` is the :math:`10\sigma` limiting flux.
+    Accordingly, :math:`m_\mathrm{lim}` is the :math:`10\sigma` limiting magnitud
+    associated with :math:`F_\mathrm{lim}`.
+
+    The effective exposure time is described by
+
+    .. math::
+
+        \ln{t_\mathrm{eff}} = a + b(m_\mathrm{lim} - 21)\;,
+
+    where :math:`a` and :math:`b` are free parameters.
+
+    Further note that the model was originally used for SDSS galaxy photometry.
+
+    References
+    ----------
+    .. [1] Rykoff E. S., Rozo E., Keisler R., 2015, eprint arXiv:1509.00870
+
+    """
+
+    flux = luminosity_from_absolute_magnitude(magnitude, -magnitude_zp)
+    flux_limit = luminosity_from_absolute_magnitude(magnitude_limit, -magnitude_zp)
+    t_eff = np.exp(a + b * np.subtract(magnitude_limit, 21.0))
+    flux_noise = np.square(flux_limit / 10) * t_eff - flux_limit
+    error = 2.5 / np.log(10) * np.sqrt((1 + flux_noise / flux) / (flux * t_eff))
+
+    return np.minimum(error, error_limit)


Test patch:
diff --git a/skypy/utils/tests/test_photometry.py b/skypy/utils/tests/test_photometry.py
index ad7e914..082949c 100644
--- a/skypy/utils/tests/test_photometry.py
+++ b/skypy/utils/tests/test_photometry.py
@@ -175,3 +175,40 @@ def test_speclite_not_installed():
     filter = 'bessell-B'
     with pytest.raises(ImportError):
         mag_ab(wavelength, spectrum, filter)
+
+
+def test_magnitude_error_rykoff():
+    from skypy.utils.photometry import magnitude_error_rykoff
+
+    # Test broadcasting to same shape given array for each parameter and
+    # test for correct result.
+    magnitude = np.full((2, 1, 1, 1, 1), 21)
+    magnitude_limit = np.full((3, 1, 1, 1), 21)
+    magnitude_zp = np.full((5, 1, 1), 21)
+    a = np.full((7, 1), np.log(200))
+    b = np.zeros(11)
+    error = magnitude_error_rykoff(magnitude, magnitude_limit, magnitude_zp, a, b)
+    # test result
+    assert np.allclose(error, 0.25 / np.log(10))
+    # test shape
+    assert error.shape == (2, 3, 5, 7, 11)
+
+    # second test for result
+    magnitude = 20
+    magnitude_limit = 22.5
+    magnitude_zp = 25
+    b = 2
+    a = np.log(10) - 1.5 * b
+    error = magnitude_error_rykoff(magnitude, magnitude_limit, magnitude_zp, a, b)
+    assert error == 2.5 / np.log(10) * np.sqrt(1 / 1000)
+
+    # test that error limit is returned if error is larger than error_limit
+    # The following set-up would give a value larger than 10
+    magnitude = 30
+    magnitude_limit = 25
+    magnitude_zp = 30
+    a = 0.5
+    b = 1.0
+    error_limit = 1
+    error = magnitude_error_rykoff(magnitude, magnitude_limit, magnitude_zp, a, b, error_limit)
+    assert error == error_limit