Add abstract base class for channels

changelog/98.feature.1.rst

@@ -0,0 +1,2 @@
+Added `sunkit_instruments.response.abstractions.AbstractChannel` to standardize an interface
+for computing wavelength and temperature response functions.

changelog/98.feature.2.rst

@@ -0,0 +1,3 @@
+Added `sunkit_instruments.response.SourceSpectra` to provide a container for
+spectra as a function of temperature and wavelength needed for computing temperature
+response functions.

docs/code_ref/index.rst

@@ -14,3 +14,4 @@ API Reference
    lyra
    rhessi
    suvi
+   response

docs/code_ref/response.rst

@@ -0,0 +1,7 @@
+========
+Response
+========
+
+.. automodapi:: sunkit_instruments.response
+
+.. automodapi:: sunkit_instruments.response.abstractions

pyproject.toml

@@ -16,7 +16,8 @@ authors = [
   { name = "The SunPy Community", email = "[email protected]" },
 ]
 dependencies = [
-  "sunpy[map,net,timeseries,visualization]>=6.0.0"
+  "sunpy[map,net,timeseries,visualization]>=6.0.0",
+  "xarray",
 ]
 dynamic = ["version"]
 

sunkit_instruments/response/__init__.py

@@ -0,0 +1,5 @@
+"""
+A subpackage for computing instrument responses
+"""
+
+from sunkit_instruments.response.thermal import *

sunkit_instruments/response/abstractions.py

@@ -0,0 +1,113 @@
+"""This module defines abstractions for computing instrument response."""
+import abc
+
+import astropy.units as u
+
+__all__ = ["AbstractChannel"]
+
+
+class AbstractChannel(abc.ABC):
+    """
+    An abstract base class for defining instrument channels.
+
+    For all methods and properties defined here, see the
+    topic guide on instrument response for more information.
+    """
+
+    @u.quantity_input
+    def wavelength_response(
+        self, obstime=None
+    ) -> u.cm**2 * u.DN * u.steradian / (u.photon * u.pixel):
+        """
+        Instrument response as a function of wavelength
+
+        The wavelength response is the effective area with
+        the conversion factors from photons to DN and steradians
+        to pixels.
+
+        Parameters
+        ----------
+        obstime: any format parsed by `~sunpy.time.parse_time`, optional
+            If specified, this is used to compute the time-dependent
+            instrument degradation.
+        """
+        area_eff = self.effective_area(obstime=obstime)
+        return (
+            area_eff
+            * self.energy_per_photon
+            * self.pixel_solid_angle
+            * self.camera_gain
+            / self.energy_per_electron
+        )
+
+    @u.quantity_input
+    def effective_area(self, obstime=None) -> u.cm**2:
+        """
+        Effective area as a function of wavelength.
+
+        The effective area is the geometrical collecting area
+        weighted by the mirror reflectance, filter transmittance,
+        quantum efficiency, and instrument degradation.
+
+        Parameters
+        ----------
+        obstime: any format parsed by `sunpy.time.parse_time`, optional
+            If specified, this is used to compute the time-dependent
+            instrument degradation.
+        """
+        return (
+            self.geometrical_area
+            * self.mirror_reflectance
+            * self.filter_transmittance
+            * self.effective_quantum_efficiency
+            * self.degradation(obstime=obstime)
+        )
+
+    @property
+    @u.quantity_input
+    def energy_per_photon(self) -> u.eV / u.photon:
+        return self.wavelength.to("eV", equivalencies=u.spectral()) / u.photon
+
+    @abc.abstractmethod
+    @u.quantity_input
+    def degradation(self, obstime=None): ...
+
+    @property
+    @abc.abstractmethod
+    @u.quantity_input
+    def geometrical_area(self) -> u.cm**2: ...
+
+    @property
+    @abc.abstractmethod
+    @u.quantity_input
+    def mirror_reflectance(self) -> u.dimensionless_unscaled: ...
+
+    @property
+    @abc.abstractmethod
+    @u.quantity_input
+    def filter_transmittance(self) -> u.dimensionless_unscaled: ...
+
+    @property
+    @abc.abstractmethod
+    @u.quantity_input
+    def effective_quantum_efficiency(self) -> u.dimensionless_unscaled: ...
+
+    @property
+    @abc.abstractmethod
+    @u.quantity_input
+    def camera_gain(self) -> u.DN / u.electron: ...
+
+    @property
+    @abc.abstractmethod
+    @u.quantity_input
+    def energy_per_electron(self) -> u.eV / u.electron: ...
+
+    @property
+    @abc.abstractmethod
+    @u.quantity_input
+    def pixel_solid_angle(self) -> u.steradian / u.pixel: ...
+
+    @property
+    @abc.abstractmethod
+    @u.quantity_input
+    def wavelength(self) -> u.Angstrom: ...

sunkit_instruments/response/tests/test_channel.py

@@ -0,0 +1,92 @@
+import numpy as np
+import pytest
+
+import astropy.units as u
+
+from sunkit_instruments.response import SourceSpectra
+from sunkit_instruments.response.abstractions import AbstractChannel
+
+
+class TestChannel(AbstractChannel):
+    @property
+    @u.quantity_input
+    def wavelength(self) -> u.Angstrom:
+        return np.linspace(100, 200, 100) * u.AA
+
+    @u.quantity_input
+    def degradation(self, obstime=None) -> u.dimensionless_unscaled:
+        return 1.0
+
+    @property
+    @u.quantity_input
+    def geometrical_area(self) -> u.cm**2:
+        return 10 * u.cm**2
+
+    @property
+    @u.quantity_input
+    def mirror_reflectance(self) -> u.dimensionless_unscaled:
+        return np.exp(
+            -((self.wavelength - self.wavelength[0]) / self.wavelength[0]).decompose()
+        )
+
+    @property
+    @u.quantity_input
+    def filter_transmittance(self) -> u.dimensionless_unscaled:
+        return np.exp(
+            -((self.wavelength - 150 * u.AA) ** 2 / (1 * u.AA) ** 2).decompose()
+        )
+
+    @property
+    @u.quantity_input
+    def effective_quantum_efficiency(self) -> u.dimensionless_unscaled:
+        return np.ones(self.wavelength.shape)
+
+    @property
+    @u.quantity_input
+    def camera_gain(self) -> u.DN / u.electron:
+        return 2 * u.DN / u.electron
+
+    @property
+    @u.quantity_input
+    def energy_per_electron(self) -> u.eV / u.electron:
+        return 3.65 * u.eV / u.electron
+
+    @property
+    @u.quantity_input
+    def pixel_solid_angle(self) -> u.steradian / u.pixel:
+        return (1 * u.arcsec) ** 2 / u.pixel
+
+
+@pytest.fixture
+def fake_channel():
+    return TestChannel()
+
+
+def test_effective_area(fake_channel):
+    assert isinstance(fake_channel.effective_area(), u.Quantity)
+
+
+def test_wavelength_response(fake_channel):
+    assert isinstance(fake_channel.wavelength_response(), u.Quantity)
+
+
+@pytest.fixture
+def fake_spectra():
+    temperature = np.logspace(5, 8, 100) * u.K
+    density = 1e15 * u.cm ** (-3) * u.K / temperature
+    wavelength = np.linspace(50, 250, 1000) * u.AA
+    data = np.random.rand(*temperature.shape + wavelength.shape) * u.Unit(
+        "photon cm3 s-1 sr-1 Angstrom-1"
+    )
+    return SourceSpectra(temperature, wavelength, data, density=density)
+
+
+def test_spectra_repr(fake_spectra):
+    assert isinstance(fake_spectra.__repr__(), str)
+
+
+@pytest.mark.parametrize('obstime', [None, '2020-01-01'])
+def test_temperature_response(fake_channel, fake_spectra, obstime):
+    temp_response = fake_spectra.temperature_response(fake_channel, obstime=obstime)
+    assert isinstance(temp_response, u.Quantity)
+    assert temp_response.shape == fake_spectra.temperature.shape

sunkit_instruments/response/thermal.py

@@ -0,0 +1,152 @@
+"""
+Classes for computing the temperature response
+"""
+import xarray
+
+import astropy.units as u
+
+__all__ = ["SourceSpectra", "get_temperature_response"]
+
+
+def get_temperature_response(channel, spectra, obstime=None):
+    """
+    Calculate the temperature response function for a given instrument channel
+    and input spectra
+
+    Parameters
+    ----------
+    channel: `~sunkit_instruments.`
+    spectra: `~sunkit_instruments.response.SourceSpectra`
+    obstime: any format parsed by `sunpy.time.parse_time`, optional
+    """
+    return spectra.temperature, spectra.temperature_response(channel, obstime=obstime)
+
+
+class SourceSpectra:
+    """
+    Source spectra as a function of temperature and wavelength.
+
+    The source spectra describes how a plasma emits (under the optically-thin
+    assumption) as a function of both temperature and wavelength. This source
+    spectra is typically computed using a database like CHIANTI by summing the
+    emission spectra of many ions as well as the continuum emission. For more
+    information, see the `topic guide on instrument response <>`_.
+
+    Parameters
+    ----------
+    temperature: `~astropy.units.Quantity`
+        1D array describing the variation along the temperature axis
+    wavelength: `~astropy.units.Quantity`
+        1D array describing the variation along the wavelength axis
+    spectra: `~astropy.units.Quantity`
+        Source spectra as a 2D array. The first axis should correspond to temperature and the
+        second axis should correspond to wavelength.
+    density: `~astropy.units.Quantity`, optional
+        1D array describing the variation in density along the temperature axis. It is assumed
+        that temperature and density are dependent.
+    meta: `dict`, optional
+        Any optional metadata to attach to the spectra, e.g. abundance model, CHIANTI version
+    """
+
+    @u.quantity_input
+    def __init__(
+        self,
+        temperature: u.K,
+        wavelength: u.Angstrom,
+        spectra: u.photon * u.cm**3 / (u.s * u.Angstrom * u.steradian),
+        density: u.cm ** (-3) = None,
+        meta=None,
+    ):
+        self.meta = meta
+        coords = {
+            "temperature": xarray.Variable(
+                "temperature",
+                temperature.value,
+                attrs={"unit": temperature.unit.to_string()},
+            ),
+            "wavelength": xarray.Variable(
+                "wavelength",
+                wavelength.value,
+                attrs={"unit": wavelength.unit.to_string()},
+            ),
+        }
+        if density is not None:
+            coords["density"] = xarray.Variable(
+                "temperature", density.value, attrs={"unit": density.unit.to_string()}
+            )
+        self._da = xarray.DataArray(
+            spectra.data,
+            dims=["temperature", "wavelength"],
+            coords=coords,
+            attrs={"unit": spectra.unit.to_string(), **self.meta},
+        )
+
+    def __repr__(self):
+        return self._da.__repr__()
+
+    def __str__(self):
+        return self._da.__str__()
+
+    def _repr_html_(self):
+        return self._da._repr_html_()
+
+    @property
+    def meta(self):
+        return self._meta
+
+    @meta.setter
+    def meta(self, x):
+        if x is None:
+            self._meta = {}
+        elif isinstance(x, dict):
+            self._mata = x
+        else:
+            raise TypeError(f'Unsupported metadata type {type(x)}')
+
+    @property
+    @u.quantity_input
+    def temperature(self) -> u.K:
+        return u.Quantity(self._da.temperature.data, self._da.temperature.attrs["unit"])
+
+    @property
+    @u.quantity_input
+    def wavelength(self) -> u.Angstrom:
+        return u.Quantity(self._da.wavelength.data, self._da.wavelength.attrs["unit"])
+
+    @property
+    @u.quantity_input
+    def data(self) -> u.photon * u.cm**3 / (u.s * u.Angstrom * u.steradian):
+        return u.Quantity(self._da.data, self._da.attrs["unit"])
+
+    @u.quantity_input
+    def temperature_response(
+        self, channel, obstime=None
+    ) -> u.cm**5 * u.DN / (u.pixel * u.s):
+        """
+        Temperature response function for a given instrument channel.
+
+        The temperature response function describes the sensitivity of an imaging
+        instrument as a function of temperature. The temperature response is
+        calculated by integrating the source spectra over the wavelength dimension,
+        weighted by the wavelength response of the instrument.
+
+        Parameters
+        ----------
+        channel: `~sunkit_instruments.response.abstractions.AbstractChannel`
+            The relevant instrument channel object used to compute the wavelength
+            response function.
+        obstime: any format parsed by `sunpy.time.parse_time`, optional
+            A time of a particular observation. This is used to calculated any
+            time-dependent instrument degradation.
+        """
+        spec_interp = self._da.interp(
+            wavelength=channel.wavelength.to_value(self.wavelength.unit),
+            kwargs={'bounds_error':False, 'fill_value':0.0},
+        )
+        wave_response = channel.wavelength_response(obstime=obstime)
+        spec_interp_weighted = spec_interp * wave_response.to_value()
+        temp_response = spec_interp_weighted.integrate(coord='wavelength')
+        final_unit = (u.Unit(spec_interp.unit)
+                      * u.Unit(wave_response.unit)
+                      * u.Unit(spec_interp_weighted.wavelength.unit))
+        return u.Quantity(temp_response.data, final_unit)