Merge pull request #172 from sot/new-scaling

Add routine for modern dark scaling

chandra_aca/dark_model.py

@@ -1,21 +1,27 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 """
-Routines related to the canonical Chandra ACA dark current model.
+Routines related to the Chandra ACA dark current models.
 
-The model is based on smoothed twice-broken power-law fits of
-dark current histograms from Jan-2007 though Aug-2017.  This analysis
-was done entirely with dark current maps scaled to -14 C.
+The canonical model for ACA dark current is a based on smoothed
+twice-broken power-law fits of dark current histograms from Jan-2007
+though Aug-2017.  This analysis was done entirely with dark current maps scaled to -14 C.
 
 See: /proj/sot/ska/analysis/dark_current_model/dark_model.ipynb
 and other files in that directory.
 
 Alternatively:
 http://nbviewer.ipython.org/url/asc.harvard.edu/mta/ASPECT/analysis/dark_current_model/dark_model.ipynb
+
+The dark_temp_scale_img method in this file uses a more recent 2023 approach with per-pixel scaling.
+
+https://nbviewer.org/url/asc.harvard.edu/mta/ASPECT/analysis/dark_current_model/dark_model-2023.ipynb
+
 """
 
 import warnings
 
 import numpy as np
+import numpy.typing as npt
 from Chandra.Time import DateTime
 
 # Define a common fixed binning of dark current distribution
@@ -340,3 +346,48 @@ def synthetic_dark_image(date, t_ccd_ref=None):
         dark *= dark_temp_scale(-14, t_ccd_ref)
 
     return dark
+
+
+def dark_temp_scale_img(img: float | npt.ArrayLike, t_ccd: float, t_ccd_ref: float):
+    """
+    Get dark current taken at ``t_ccd`` scaled to reference temperature ``t_ccd_ref``
+
+    This scales the dark current based on an exponential scaling factor that depends on
+    the dark current value of each pixel.  This is a more accurate way to scale dark
+    current images than using a global scaling factor as in dark_temp_scale().
+
+    See the reference notebook at:
+    https://nbviewer.org/url/asc.harvard.edu/mta/ASPECT/analysis/dark_current_model/dark_model-2023.ipynb
+
+    Parameters
+    ----------
+    img : float, ArrayLike
+        Dark current image or value in e-/sec
+    t_ccd : float
+        CCD temperature (degC) of the input image
+    t_ccd_ref : float
+        CCD temperature (degC) of the scaled output image
+
+    Returns
+    -------
+    float, np.ndarray
+        Dark current image scaled to ``t_ccd_ref``
+    """
+
+    # Confirm t_ccd and t_ref are just floats
+    t_ccd = float(t_ccd)
+    t_ccd_ref = float(t_ccd_ref)
+
+    # this comes from the simple fit to DC averages, with fixed T_CCD=265.15 (-8.0 C)
+    a, b, c, d, e = [
+        -4.88802057e00,
+        -1.66791619e-04,
+        -2.22596103e-01,
+        -2.45720364e-03,
+        1.90718453e-01,
+    ]
+    t = 265.15
+    y = np.log(np.clip(img, 20, 1e4)) - e * t
+    scale = a + b * t + c * y + d * y**2
+
+    return img * np.exp(scale * (t_ccd_ref - t_ccd))

chandra_aca/tests/test_dark_model.py

@@ -1,10 +1,15 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 import os
+import pickle
+from pathlib import Path
 
 import numpy as np
 import pytest
+from mica.archive.aca_dark import get_dark_cal_props
 from mica.common import MICA_ARCHIVE
 
+from chandra_aca.dark_model import dark_temp_scale_img
+
 from ..dark_model import dark_temp_scale, get_warm_fracs, synthetic_dark_image
 
 HAS_MICA = os.path.exists(MICA_ARCHIVE)
@@ -72,3 +77,37 @@ def test_get_warm_fracs_2017185():
     exp = np.array([207845, 79207, 635, 86, 26])
 
     assert np.allclose(wps, exp, rtol=0.001, atol=2)
+
+
+@pytest.mark.parametrize(
+    "img, t_ccd, t_ref, expected",
+    [
+        (np.array([100, 1000, 500]), -10.0, -5.0, np.array([171.47, 1668.62, 852.79])),
+        (np.array([100, 1000, 500]), -10.0, -15.0, np.array([58.31, 599.29, 293.15])),
+        ([[100, 1000], [500, 600]], -15.0, -15.0, [[100, 1000], [500, 600]]),
+        ([200, 2000, 600], -6.0, 2.0, np.array([478.16, 4299.02, 1399.40])),
+        (2000, -6, 2, 4299.02),
+        ([2000, np.nan], -6, 2, [4299.02, np.nan]),
+        (np.nan, -6, 2, np.nan),
+    ],
+)
+def test_dark_temp_scale_img(img, t_ccd, t_ref, expected):
+    scaled_img = dark_temp_scale_img(img, t_ccd, t_ref)
+    assert np.allclose(scaled_img, expected, atol=0.1, rtol=0, equal_nan=True)
+    if np.shape(img):
+        assert isinstance(scaled_img, np.ndarray)
+    else:
+        assert isinstance(scaled_img, float)
+
+
+@pytest.mark.skipif("not HAS_MICA")
+def test_dark_temp_scale_img_real_dc():
+    test_data = {}
+    with open((Path(__file__).parent / "data" / "dark_scaled_img.pkl"), "rb") as f:
+        test_data.update(pickle.load(f))
+
+    dc = get_dark_cal_props("2024:001", include_image=True)
+    # just use 100 square pixels instead of 1024x1024
+    img = dc["image"][100:200, 100:200]
+    scaled_img = dark_temp_scale_img(img, dc["t_ccd"], -3.0)
+    assert np.allclose(scaled_img, test_data["scaled_img"], atol=0.1, rtol=0)