Added optika.sensors.mean_charge_capture() function. (#79)

docs/refs.bib

@@ -177,3 +177,19 @@ @ARTICLE{Heymes2020
   keywords={Passivation;Silicon;Sensors;Photonics;Wavelength measurement;Absorption;Standards;CCD;CIS;EUV;quantum efficiency (QE);vacuum ultraviolet (VUV);X-ray},
   doi={10.1109/TNS.2020.3001622}
 }
+@INPROCEEDINGS{Stern2004,
+       author = {{Stern}, Robert A. and {Shing}, Lawrence and {Catura}, Paul R. and {Morrison}, Mons D. and {Duncan}, Dexter W. and {Lemen}, James R. and {Eaton}, Tim and {Pool}, Peter J. and {Steward}, Roy and {Walton}, Dave M. and {Smith}, Alan},
+        title = "{Characterization of the flight CCD detectors for the GOES N and O solar x-ray imagers}",
+    booktitle = {Telescopes and Instrumentation for Solar Astrophysics},
+         year = 2004,
+       editor = {{Fineschi}, Silvano and {Gummin}, Mark A.},
+       series = {Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series},
+       volume = {5171},
+        month = feb,
+        pages = {77-88},
+          doi = {10.1117/12.506346},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2004SPIE.5171...77S},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+

optika/sensors/__init__.py

@@ -5,6 +5,7 @@
 
 from ._materials import (
     charge_diffusion,
+    mean_charge_capture,
     energy_bandgap,
     energy_electron_hole,
     quantum_yield_ideal,
@@ -30,6 +31,7 @@
 
 __all__ = [
     "charge_diffusion",
+    "mean_charge_capture",
     "energy_bandgap",
     "energy_electron_hole",
     "quantum_yield_ideal",

optika/sensors/_materials/__init__.py

@@ -1,5 +1,6 @@
 from ._diffusion import (
     charge_diffusion,
+    mean_charge_capture,
 )
 from ._materials import (
     energy_bandgap,
@@ -22,6 +23,7 @@
 
 __all__ = [
     "charge_diffusion",
+    "mean_charge_capture",
     "energy_bandgap",
     "energy_electron_hole",
     "quantum_yield_ideal",

optika/sensors/_materials/_diffusion.py

@@ -1,10 +1,11 @@
 import astropy.units as u
 import numpy as np
-
+import scipy.special
 import named_arrays as na
 
 __all__ = [
     "charge_diffusion",
+    "mean_charge_capture",
 ]
 
 
@@ -127,4 +128,125 @@ def charge_diffusion(
 
     result = x_ff * np.sqrt(1 - depth_avg / x_ff)
 
-    return result
+    return np.nan_to_num(result)
+
+
+def mean_charge_capture(
+    width_diffusion: u.Quantity | na.AbstractScalar,
+    width_pixel: u.Quantity | na.AbstractScalar,
+) -> na.AbstractScalar:
+    r"""
+    A function to compute the mean charge capture :cite:p:`Stern2004`,
+    the fraction of charge from each photon event retained in the central pixel.
+
+    Parameters
+    ----------
+    width_diffusion
+        The standard deviation of the charge diffusion kernel.
+    width_pixel
+        The width of a pixel on the sensor.
+
+    Examples
+    --------
+
+    Plot the mean charge capture as a function of wavelength
+    and energy for the sensor parameters in :cite:t:`Heymes2020`.
+
+    .. jupyter-execute::
+
+        import matplotlib.pyplot as plt
+        import astropy.units as u
+        import astropy.visualization
+        import named_arrays as na
+        import optika
+
+        # Define a grid of wavelengths
+        wavelength = na.geomspace(1, 10000, axis="wavelength", num=1001) * u.AA
+
+        # Convert the grid to energies as well
+        energy = wavelength.to(u.eV, equivalencies=u.spectral())
+
+        # Load the optical properties of silicon
+        si = optika.chemicals.Chemical("Si")
+
+        # Retrieve the absorption coefficient of silicon
+        # for the given wavelenghts.
+        absorption = si.absorption(wavelength)
+
+        # Compute the charge diffusion
+        width_diffusion = optika.sensors.charge_diffusion(
+            absorption=absorption,
+            thickness_implant=40 * u.nm,
+            thickness_substrate=14 * u.um,
+            thickness_depletion=2.4 * u.um,
+        )
+
+        # Compute the mean charge capture
+        mcc = optika.sensors.mean_charge_capture(
+            width_diffusion=width_diffusion,
+            width_pixel=16 * u.um,
+        )
+
+        # Plot the mean charge capture as a function
+        # of wavelength and energy
+        with astropy.visualization.quantity_support():
+            fig, ax = plt.subplots()
+            ax2 = ax.twiny()
+            ax2.invert_xaxis()
+            na.plt.plot(
+                wavelength,
+                mcc,
+                ax=ax,
+            )
+            na.plt.plot(
+                energy,
+                mcc,
+                ax=ax2,
+                linestyle="None",
+            )
+            ax.set_xscale("log")
+            ax2.set_xscale("log")
+            ax.set_xlabel(f"wavelength ({ax.get_xlabel()})")
+            ax2.set_xlabel(f"energy ({ax2.get_xlabel()})")
+            ax.set_ylabel(f"mean charge capture")
+
+    Notes
+    -----
+    Naively, the mean charge capture (MCC) is the integral of the charge
+    diffusion kernel over the extent of a pixel.
+    However, since a photon can strike anywhere within the central pixel,
+    the charge diffusion kernel should be convolved with a rectangle function
+    the width of a pixel before integrating.
+    So, our definition for the MCC is
+
+    .. math::
+
+        P_\text{MCC} = \left\{ \frac{1}{d} \int_{-d/2}^{d/2} \left[ K(x') * \Pi \left( \frac{x'}{d} \right) \right](x) \, dx \right\}^2,
+
+    where :math:`K(x)` is the charge diffusion kernel,
+    :math:`\Pi(x)` is the `rectangle function <https://en.wikipedia.org/wiki/Rectangular_function>`_,
+    and :math:`d` is the width of a pixel.
+    If we assume that the charge diffusion kernel is a Gaussian with standard
+    deviation :math:`\sigma`,
+
+    .. math::
+
+        K(x) = \frac{1}{\sqrt{2\pi} \sigma} \exp \left( -\frac{x^2}{2 \sigma^2} \right),
+
+    then we can analytically solve for the MCC,
+
+    .. math::
+
+        P_\text{MCC} &= \left\{ \frac{1}{2d} \int_{-d/2}^{d/2} \left[ \text{erf} \left( \frac{d - 2x}{2 \sqrt{2} \sigma} \right) + \text{erf} \left( \frac{d + 2x}{2 \sqrt{2} \sigma} \right) \right] dx \right\}^2 \\
+                     &= \left\{ \sqrt{\frac{2}{\pi}} \frac{\sigma}{d} \left[ \exp \left( -\frac{d^2}{2 \sigma^2} \right) - 1 \right] + \text{erf} \left( \frac{d}{\sqrt{2} \sigma} \right) \right\}^2,
+
+    where :math:`\text{erf}(x)` is the `error function <https://en.wikipedia.org/wiki/Error_function>`_.
+    """
+    a = width_pixel / width_diffusion
+
+    t1 = np.sqrt(2 / np.pi) * (np.exp(-np.square(a) / 2) - 1) / a
+    t2 = scipy.special.erf(a / np.sqrt(2))
+
+    result = np.square(t1 + t2)
+
+    return result.to(u.dimensionless_unscaled)

optika/sensors/_materials/_diffusion_test.py

@@ -1,4 +1,5 @@
 import pytest
+import numpy as np
 import astropy.units as u
 import named_arrays as na
 import optika
@@ -42,3 +43,28 @@ def test_charge_diffusion(
     )
 
     assert result > 0 * u.um
+
+
+@pytest.mark.parametrize(
+    argnames="width_diffusion",
+    argvalues=[
+        10 * u.um,
+        na.linspace(1, 10, "width", 5) * u.um,
+    ],
+)
+@pytest.mark.parametrize(
+    argnames="width_pixel",
+    argvalues=[
+        15 * u.um,
+    ],
+)
+def test_mean_charge_capture(
+    width_diffusion: u.Quantity | na.AbstractScalar,
+    width_pixel: u.Quantity | na.AbstractScalar,
+):
+    result = optika.sensors.mean_charge_capture(
+        width_diffusion=width_diffusion,
+        width_pixel=width_pixel,
+    )
+    assert np.all(result > 0)
+    assert np.all(result < 1)