Added iris.sg.SpectrographObservation.radiance property which conve…

…rts from DN to radiometric units. (#21)

iris/sg/__init__.py

@@ -3,11 +3,17 @@
 """
 
 from . import background
-from ._effective_area import effective_area
+from ._effective_area import (
+    gain,
+    width_slit,
+    effective_area,
+)
 from ._spectrograph import SpectrographObservation
 
 __all__ = [
     "background",
+    "gain",
+    "width_slit",
     "effective_area",
     "SpectrographObservation",
 ]

iris/sg/_effective_area.py

@@ -4,9 +4,23 @@
 import iris
 
 __all__ = [
+    "gain",
+    "width_slit",
     "effective_area",
 ]
 
+gain = 12 * u.photon / u.DN
+"""
+The conversion factor between counts and photons measured by the sensor
+:cite:p:`Wulser2018`.
+"""
+
+
+width_slit = 1 / 3 * u.arcsec
+"""
+The angular subtent of the spectrographic slit.
+"""
+
 
 def effective_area(wavelength: u.Quantity | na.AbstractScalar) -> na.AbstractScalar:
     """

iris/sg/_spectrograph.py

@@ -1,7 +1,9 @@
+from typing_extensions import Self
 import dataclasses
 import pathlib
 import numpy as np
 import astropy.units as u
+import astropy.constants
 import astropy.time
 import astropy.wcs
 import astropy.io.fits
@@ -432,3 +434,41 @@ def empty(
             axis_detector_x=axis_detector_x,
             axis_detector_y=axis_detector_y,
         )
+
+    @property
+    def radiance(self) -> Self:
+        """
+        Convert to radiometric units using :func:`iris.sg.effective_area`.
+        """
+
+        wavelength = self.inputs.wavelength
+
+        lower = {self.axis_wavelength: slice(None, ~0)}
+        upper = {self.axis_wavelength: slice(+1, None)}
+        wavelength = (wavelength[lower] + wavelength[upper]) / 2
+
+        energy = astropy.constants.h * astropy.constants.c / wavelength / u.ph
+
+        gain = iris.sg.gain
+
+        area_eff = iris.sg.effective_area(wavelength)
+
+        pix_xy = np.diff(self.inputs.position, axis=self.axis_detector_y).length
+        pix_xy = pix_xy.mean(self.axis_detector_x)
+
+        pix_lambda = np.diff(self.inputs.wavelength, axis=self.axis_wavelength)
+
+        t_exp = self.timedelta
+
+        w_slit = iris.sg.width_slit
+
+        factor = energy * gain / (area_eff * pix_xy * pix_lambda * t_exp * w_slit)
+
+        outputs = self.outputs * factor
+
+        outputs = outputs.to(u.erg / (u.cm**2 * u.sr * u.nm * u.s))
+
+        return dataclasses.replace(
+            self,
+            outputs=outputs,
+        )

iris/sg/_spectrograph_test.py

@@ -1,4 +1,6 @@
 import pytest
+import numpy as np
+import astropy.units as u
 import astropy.time
 import iris
 
@@ -25,3 +27,9 @@ def test_axis_detector_x(self, array: iris.sg.SpectrographObservation):
 
     def test_axis_detector_y(self, array: iris.sg.SpectrographObservation):
         assert isinstance(array.axis_detector_y, str)
+
+    def test_radiance(self, array: iris.sg.SpectrographObservation):
+        result = array.radiance
+        assert isinstance(result, iris.sg.SpectrographObservation)
+        assert np.all(result.inputs == array.inputs)
+        assert np.nansum(result.outputs) > 0 * u.erg / (u.cm**2 * u.sr * u.s * u.nm)