PSF model

docs/changes/2643.feature.rst

@@ -0,0 +1,2 @@
+Adds psf models to ``ctapipe.instrument.optics`` with the parent class ``PSFModel`` and a psf model based on pure coma aberration called ``ComaPSFModel``
+- The function ``PSFModel.pdf`` gives the value of the PSF in a given location

docs/references.bib

@@ -137,6 +137,17 @@ @article{bright-star-catalog
   url={https://heasarc.gsfc.nasa.gov/W3Browse/star-catalog/bsc5p.html},
 }
 
+
+@article{startracker,
+	author={Abe, K. and others},
+  title={Star tracking for pointing determination of Imaging Atmospheric Cherenkov Telescopes - Application to the Large-Sized Telescope of the Cherenkov Telescope Array},
+	doi={10.1051/0004-6361/202347128},
+	journal={A\&A},
+	year={2023},
+	volume={679},
+	pages={A90},
+}
+
 @manual{ctao-software-standards,
   title={Software Programming Standards},
   number={CTA-STD-OSO-000000-0001},

src/ctapipe/instrument/__init__.py

@@ -1,7 +1,14 @@
 from .atmosphere import get_atmosphere_profile_functions
 from .camera import CameraDescription, CameraGeometry, CameraReadout, PixelShape
 from .guess import guess_telescope
-from .optics import FocalLengthKind, OpticsDescription, ReflectorShape, SizeType
+from .optics import (
+    ComaPSFModel,
+    FocalLengthKind,
+    OpticsDescription,
+    PSFModel,
+    ReflectorShape,
+    SizeType,
+)
 from .subarray import SubarrayDescription, UnknownTelescopeID
 from .telescope import TelescopeDescription
 from .trigger import SoftwareTrigger
@@ -23,4 +30,6 @@
     "SizeType",
     "SoftwareTrigger",
     "FromNameWarning",
+    "PSFModel",
+    "ComaPSFModel",
 ]

src/ctapipe/instrument/optics.py

@@ -3,21 +3,29 @@
 """
 
 import logging
+from abc import abstractmethod
 from enum import Enum, auto, unique
 
 import astropy.units as u
 import numpy as np
 from astropy.table import QTable
+from scipy.stats import laplace, laplace_asymmetric
+from traitlets import validate
 
 from ..compat import StrEnum
+from ..core import TelescopeComponent
+from ..core.traits import Float, List, TraitError
 from ..utils import get_table_dataset
+from ..utils.quantities import all_to_value
 from .warnings import warn_from_name
 
 logger = logging.getLogger(__name__)
 
 __all__ = [
     "OpticsDescription",
     "FocalLengthKind",
+    "PSFModel",
+    "ComaPSFModel",
 ]
 
 
@@ -263,3 +271,174 @@
 
     def __str__(self):
         return self.name
+
+
+class PSFModel(TelescopeComponent):
+    """
+    Base component to describe image distortion due to the optics of the different cameras.
+    """
+
+    @u.quantity_input(x=u.m, y=u.m, x0=u.m, y0=u.m)
+    @abstractmethod
+    def pdf(self, x, y, x0, y0) -> np.ndarray:
+        """
+        Calculates the value of the psf at a given location
+
+        Parameters
+        ----------
+        x : u.Quantity[length]
+            x-coordinate of the point on the focal plane where the psf is evaluated
+        y : u.Quantity[length]
+            y-coordinate of the point on the focal plane where the psf is evaluated
+        x0 : u.Quantity[length]
+            x-coordinate of the point source on the focal plane
+        y0 : u.Quantity[length]
+            y-coordinate of the point source on the focal plane
+        Returns
+        ----------
+        psf : np.ndarray
+            value of the PSF at the specified location with the specified position of the point source
+        """
+
+        pass
+
+
+def _cartesian_to_polar(x, y):
+    r = np.sqrt(x**2 + y**2)
+    phi = np.arctan2(y, x)
+    return r, phi
+
+
+class ComaPSFModel(PSFModel):
+    r"""
+    PSF model describing pure coma aberrations PSF effect.
+
+    The PSF is described by a product of an asymmetric Laplacian for the radial part and a symmetric Laplacian in the polar direction.
+    Explicitly, the radial part is given by
+
+    .. math:: f_{R}(r, K) = \begin{cases}\frac{1}{S_{R}(K+K^{-1})}e^{-K\frac{r-r_0}{S_{R}}}, r\ge r_0\\ \frac{1}{S_{R}(K+K^{-1})}e^{\frac{r-r_0}{KS_{R}}}, r < r_0\end{cases}
+
+    and the polar part is given by
+
+    .. math:: f_{\Phi}(\phi) = \frac{1}{2S_\phi}e^{-|\frac{\phi-\phi_0}{S_\phi}|}
+
+    The parameters :math:`K`, :math:`S_{R}`, and :math:`S_{\phi}` are functions of the distance :math:`r` to the optical axis.
+    Their detailed description is provided in the attributes section.
+
+    Attributes
+    ----------
+    asymmetry_params : list
+        Describes the dependency of the PSF on the distance to the center of the camera.
+        Used to calculate a PDF asymmetry parameter K of the asymmetric radial Laplacian
+        of the PSF as a function of the distance r to the optical axis.
+
+        .. math:: K(r) = 1 - c_0 \tanh(c_1 r) - c_2 r
+
+    radial_scale_params : list
+        Describes the dependency of the radial scale on the distance to the center of the camera.
+        Used to calculate width Sr of the asymmetric radial Laplacian in the PSF as a function of the distance :math:`r` to the optical axis.
+
+        .. math:: S_{R}(r) = b_1 + b_2\,r + b_3\,r^2 + b_4\,r^3
+
+    phi_scale_params : list
+        Describes the dependency of the polar angle (:math:`\phi`) scale on the distance to the center of the camera.
+        Used to calculate the width Sf of the polar Laplacian in the PSF as a function of the distance :math:`r` to the optical axis.
+
+        .. math:: S_{\phi}(r) = a_1\,\exp{(-a_2\,r)}+\frac{a_3}{a_3+r}
+
+    Parameters
+    ----------
+    subarray : ctapipe.instrument.SubarrayDescription
+        Description of the subarray.
+
+    References
+    ----------
+    For reference, see :cite:p:`startracker`
+    """
+
+    asymmetry_params = List(
+        help=(
+            "Describes the dependency of the PSF on the distance "
+            "to the center of the camera. Used to calculate a PDF "
+            "asymmetry parameter :math:`K` of the asymmetric radial Laplacian "
+            "of the PSF as a function of the distance r to the optical axis"
+        )
+    ).tag(config=True)
+
+    radial_scale_params = List(
+        help=(
+            "Describes the dependency of the radial scale on the "
+            "distance to the center of the camera. Used to calculate "
+            "width :math:`S_R` of the asymmetric radial Laplacian in the PSF "
+            "as a function of the distance r to the optical axis"
+        )
+    ).tag(config=True)
+
+    phi_scale_params = List(
+        help=(
+            "Describes the dependency of the polar scale on the "
+            "distance to the center of the camera. Used to calculate "
+            "the width :math:`S_\phi` of the polar Laplacian in the PSF "
+            "as a function of the distance r to the optical axis"
+        )
+    ).tag(config=True)
+
+    pixel_width = Float(
+        default_value=0.05,
+        help="Width of a pixel of the camera in meters",
+    ).tag(config=True)
+
+    def _k(self, r):
+        c1, c2, c3 = self.asymmetry_params
+        return 1 - c1 * np.tanh(c2 * r) - c3 * r
+
+    def _s_r(self, r):
+        return np.polyval(self.radial_scale_params[::-1], r)
+
+    def _s_phi(self, r):
+        a1, a2, a3 = self.phi_scale_params
+        return a1 * np.exp(-a2 * r) + a3 / (a3 + r)
+
+    @u.quantity_input(x=u.m, y=u.m, x0=u.m, y0=u.m)
+    def pdf(self, x, y, x0, y0) -> np.ndarray:
+        x, y, x0, y0 = all_to_value(x, y, x0, y0, unit=u.m)
+        r, phi = _cartesian_to_polar(x, y)
+        r0, phi0 = _cartesian_to_polar(x0, y0)
+
+        k = self._k(r0)
+        s_r = self._s_r(r0)
+        s_phi = self._s_phi(r0)
+
+        radial_pdf = laplace_asymmetric.pdf(r, k, r0, s_r)
+        polar_pdf = laplace.pdf(phi, phi0, s_phi)
+
+        # Phi is not defined at the center
+        at_center = np.isclose(r0, 0, atol=self.pixel_width)
+        polar_pdf = np.where(at_center, 1 / (2 * s_phi), polar_pdf)
+        # Polar PDF is valid under approximation that the polar axis is orthogonal to the radial axis
+        # Thus, we limit the PDF to a chord of 6 pixels or covering ~30deg around the radial axis, whichever is smaller
+        chord_length = min(6 * self.pixel_width, 0.5 * r0)
+        if r0 != 0:
+            dphi = np.arcsin(chord_length / (2 * r0))
+            polar_pdf[phi < phi0 - dphi] = 0
+            polar_pdf[phi > phi0 + dphi] = 0
+
+        return radial_pdf * polar_pdf
+
+    @validate("asymmetry_params")
+    def _check_asymmetry_params(self, proposal):
+        if not len(proposal["value"]) == 3:
+            raise TraitError("asymmetry_params needs to have length 3")
+        return proposal["value"]
+
+    @validate("radial_scale_params")
+    def _check_radial_scale_params(self, proposal):
+        if not len(proposal["value"]) == 4:
+            raise TraitError("radial_scale_params needs to have length 4")
+        return proposal["value"]
+
+    @validate("phi_scale_params")
+    def _check_phi_scale_params(self, proposal):
+        if not len(proposal["value"]) == 3:
+            raise TraitError("phi_scale_params needs to have length 3")
+        return proposal["value"]

src/ctapipe/instrument/tests/test_psf_model.py

@@ -0,0 +1,62 @@
+"""
+This module contains the ctapipe.image.psf_model unit tests
+"""
+
+import astropy.units as u
+import numpy as np
+import pytest
+
+from ctapipe.core.traits import TraitError
+from ctapipe.instrument.optics import PSFModel
+
+
+@pytest.fixture(scope="session")
+def coma_psf(example_subarray):
+    psf = PSFModel.from_name(
+        "ComaPSFModel",
+        subarray=example_subarray,
+        asymmetry_params=[0.5, 10, 0.15],
+        radial_scale_params=[0.015, -0.1, 0.06, 0.03],
+        phi_scale_params=[0.25, 7.5, 0.02],
+    )
+    return psf
+
+
+def test_psf(example_subarray):
+    with pytest.raises(
+        TraitError,
+        match="phi_scale_params needs to have length 3",
+    ):
+        PSFModel.from_name(
+            "ComaPSFModel",
+            subarray=example_subarray,
+            asymmetry_params=[0.0, 0.0, 0.0],
+            radial_scale_params=[0.0, 0.0, 0.0, 0.0],
+            phi_scale_params=[0.0],
+        )
+    with pytest.raises(
+        TraitError,
+        match="radial_scale_params needs to have length 4",
+    ):
+        PSFModel.from_name(
+            "ComaPSFModel",
+            subarray=example_subarray,
+            asymmetry_params=[0.0, 0.0, 0.0],
+            radial_scale_params=[0.0, 0.0, 0.0],
+            phi_scale_params=[0.0, 0.0, 0.0],
+        )
+    with pytest.raises(
+        TraitError,
+        match="asymmetry_params needs to have length 3",
+    ):
+        PSFModel.from_name(
+            "ComaPSFModel",
+            subarray=example_subarray,
+            asymmetry_params=[0.0, 0.0, 0.0, 0.0],
+            radial_scale_params=[0.0, 0.0, 0.0, 0.0],
+            phi_scale_params=[0.0, 0.0, 0.0],
+        )
+
+
+def test_asymptotic_behavior(coma_psf):
+    assert np.isclose(coma_psf.pdf(*([10.0, 0.0, 1.0, 0.0] * u.m)), 0.0)