From ba38b5a88afe210be6b9be16b780d9b616e56941 Mon Sep 17 00:00:00 2001
From: Alexandre Boucaud <3065310+aboucaud@users.noreply.github.com>
Date: Wed, 2 Mar 2022 23:42:21 +0100
Subject: [PATCH] Implement utility functions purely related to survey and
 filter numbers (#24)

* Tentative implementation of functions

* Simplify the computation for the conversion

* Add docstring explainations

* Cast value to integer

* Add tests
---
 .github/workflows/python-package.yml |  5 +-
 galcheat/utilities.py                | 77 ++++++++++++++++++++++++++++
 tests/test_utilities.py              | 33 ++++++++++++
 3 files changed, 114 insertions(+), 1 deletion(-)
 create mode 100644 galcheat/utilities.py
 create mode 100644 tests/test_utilities.py

diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
index f337364..4dfdcc8 100644
--- a/.github/workflows/python-package.yml
+++ b/.github/workflows/python-package.yml
@@ -27,7 +27,7 @@ jobs:
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
-        python -m pip install flake8
+        python -m pip install flake8 pytest pytest-cov
     - name: Lint with flake8
       run: |
         # stop the build if there are Python syntax errors or undefined names
@@ -43,3 +43,6 @@ jobs:
     - name: Test the console script
       run: |
         galcheat
+    - name: Test the utility functions
+      run: |
+        pytest --cov galcheat/ --cov-branch --cov-report term-missing -vv
diff --git a/galcheat/utilities.py b/galcheat/utilities.py
new file mode 100644
index 0000000..04cafe1
--- /dev/null
+++ b/galcheat/utilities.py
@@ -0,0 +1,77 @@
+import astropy.units as u
+
+from galcheat.helpers import get_filter, get_survey
+
+
+def mag2counts(magnitude, survey_name, filter_name):
+    """Convert source magnitude to counts for a given filter of a survey
+
+    To perform the computation, we use the filter zeropoint computed
+    with `speclite` under classical atmospheric conditions and at a
+    given airmass and we integrate over the survey lifetime using the
+    full filter exposure time.
+
+    Expect a rough estimate from this calculation since e.g. it does not
+    take into account the atmospheric extinction. Therefore the result
+    is casted to an integer.
+
+    Parameters
+    ----------
+    magnitude: float
+        magnitude of source
+    survey_name: str
+        Name of a given survey
+    filter_name: str
+        Name of the survey filter
+
+    Returns
+    -------
+    The corresponding flux in counts
+
+    References
+    ----------
+    The `WeakLensingDeblending` package
+    https://github.com/LSSTDESC/WeakLensingDeblending
+
+    """
+    if not isinstance(magnitude, u.Quantity):
+        magnitude *= u.mag(u.ct / u.s)
+    else:
+        magnitude = magnitude.value * u.mag(u.ct / u.s)
+
+    filter = get_filter(filter_name, survey_name)
+
+    flux = (magnitude - filter.zeropoint).to(u.ct / u.s)
+    counts = flux * filter.exposure_time
+
+    return counts.astype(int)
+
+
+def mean_sky_level(survey_name, filter_name):
+    """Computes the mean sky level for a given survey and a filter
+
+    This computation uses the sky brightness parameter from galcheat,
+    expressed as a magnitude per square arcminute, weights it by the
+    pixel area and converts it to counts.
+
+    Parameters
+    ----------
+    survey_name: str
+        Name of a given survey
+    filter_name: str
+        Name of the survey filter
+
+    Returns
+    -------
+    The corresponding mean sky level in counts
+
+    """
+    survey = get_survey(survey_name)
+    filter = get_filter(filter_name, survey_name)
+
+    sky_brightness_counts = mag2counts(filter.sky_brightness, survey_name, filter_name)
+    pixel_area = survey.pixel_scale.to_value(u.arcsec) ** 2
+
+    mean_sky_level = sky_brightness_counts * pixel_area
+
+    return mean_sky_level
diff --git a/tests/test_utilities.py b/tests/test_utilities.py
new file mode 100644
index 0000000..3d8f89b
--- /dev/null
+++ b/tests/test_utilities.py
@@ -0,0 +1,33 @@
+from galcheat.utilities import mag2counts, mean_sky_level
+
+BTK_COUNTS_MAG24 = {
+    "Rubin_u": 15321.782101179268,
+    "Rubin_g": 121397.91888074754,
+    "Rubin_r": 240956.73939657194,
+    "Rubin_i": 179448.18815675998,
+    "Rubin_z": 108953.51289042353,
+    "Rubin_y": 50727.240442255556,
+}
+
+BTK_MEAN_SKY_LEVEL = {
+    "Rubin_u": 1687.9876819744386,
+    "Rubin_g": 23241.878848667337,
+    "Rubin_r": 127057.1381640446,
+    "Rubin_i": 180301.3875959776,
+    "Rubin_z": 250784.8105213134,
+    "Rubin_y": 293292.6831817095,
+}
+
+
+def test_mag2counts():
+    survey = "Rubin"
+    for filt in "ugrizy":
+        counts = mag2counts(24, survey, filt).value
+        assert counts == int(BTK_COUNTS_MAG24[f"{survey}_{filt}"])
+
+
+def test_mean_sky_level():
+    survey = "Rubin"
+    for filt in "ugrizy":
+        sky_level = mean_sky_level(survey, filt).value
+        assert int(sky_level) == int(BTK_MEAN_SKY_LEVEL[f"{survey}_{filt}"])