Modified optika.sensors.signal() to model charge diffusion using a Monte-Carlo simulation.

optika/rays/_ray_vectors.py

@@ -70,6 +70,22 @@ def unvignetted(self) -> na.AbstractScalar:
         and :obj:`False` indicates the ray has been blocked by an aperture.
         """
 
+    @property
+    def n(self) -> complex | na.AbstractScalar:
+        """
+        The complex index of refraction of the medium that the ray is
+        traveling in.
+
+        The real part is :attr:`index_refraction` and the imaginary part is the
+        extinction coefficient, :math:`k = \\alpha \\lambda / 4 \\pi`, computed
+        from the :attr:`attenuation` coefficient :math:`\\alpha` and the
+        :attr:`wavelength` :math:`\\lambda`.
+        """
+        return (
+            self.index_refraction
+            + self.attenuation * self.wavelength / (4 * np.pi) * 1j
+        )
+
     @property
     def type_abstract(self) -> type[na.AbstractArray]:
         return AbstractRayVectorArray

optika/rays/_tests/test_ray_vectors.py

@@ -80,6 +80,21 @@ def test_attenuation(self, array: optika.rays.AbstractRayVectorArray):
     def test_index_refraction(self, array: optika.rays.AbstractRayVectorArray):
         assert isinstance(na.as_named_array(array.index_refraction), na.AbstractScalar)
 
+    def test_n(self, array: optika.rays.AbstractRayVectorArray):
+        result = array.n
+        assert isinstance(na.as_named_array(result), na.AbstractScalar)
+
+        # The real part is the real index of refraction.
+        assert np.allclose(np.real(result), array.index_refraction)
+
+        # The imaginary part is the (non-negative) extinction coefficient,
+        # k = alpha * lambda / (4 pi).
+        extinction = (array.attenuation * array.wavelength / (4 * np.pi)).to(
+            u.dimensionless_unscaled
+        )
+        assert np.allclose(np.imag(result), extinction)
+        assert np.all(np.imag(result) >= 0)
+
     def test_unvignetted(self, array: optika.rays.AbstractRayVectorArray):
         assert isinstance(na.as_named_array(array.unvignetted), na.AbstractScalar)
 

optika/sensors/__init__.py

@@ -16,6 +16,7 @@
     fano_factor,
     fano_factor_inf,
     transmittance,
+    absorption_effective,
     absorbance,
     charge_collection_efficiency,
     quantum_efficiency_effective,
@@ -42,6 +43,7 @@
     "fano_factor",
     "fano_factor_inf",
     "transmittance",
+    "absorption_effective",
     "absorbance",
     "charge_collection_efficiency",
     "quantum_efficiency_effective",

optika/sensors/_sensors.py

@@ -81,70 +81,53 @@ def aperture(self):
             half_width=self.width_pixel * self.num_pixel / 2,
         )
 
-    def readout(
+    def collect(
         self,
         rays: optika.rays.RayVectorArray,
         wavelength: na.AbstractScalar,
-        timedelta: None | u.Quantity | na.AbstractScalar = None,
         axis: None | str | Sequence[str] = None,
         where: bool | na.AbstractScalar = True,
-        noise: bool = True,
-    ) -> na.FunctionArray[
-        na.SpectralPositionalVectorArray,
+    ) -> tuple[
+        na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar],
         na.AbstractScalar,
     ]:
         """
-        Given a set of rays incident on the sensor surface,
-        where each ray represents an expected number of photons per unit time,
-        simulate the number of electrons that would be measured by the sensor.
+        Bin a cloud of rays onto the pixel grid.
+
+        Returns the per-pixel photon image (the binned ray intensity) and the
+        flux-weighted mean cosine of the *refracted* angle inside the
+        light-sensitive region in each pixel: the two quantities :meth:`expose`
+        needs. Refracting each ray here (with its own ambient index of
+        refraction) and binning the result is what lets :meth:`expose` and the
+        material's ``signal`` model be shared with systems that have no rays,
+        without threading a separate ambient-index argument through them.
 
         Parameters
         ----------
         rays
-            A set of incident rays in local coordinates to measure.
+            A set of incident rays in local coordinates to bin.
         wavelength
             The edges of the wavelength bins to sample.
-        timedelta
-            The exposure time of the measurement.
-            If :obj:`None` (the default), the value in :attr:`timedelta_exposure`
-            will be used.
         axis
             The logical axes along which to collect photons.
         where
-            A boolean mask used to indicate which photons should be considered
-            when calculating the signal measured by the sensor.
-        noise
-            Whether to add noise to the result
+            A boolean mask used to indicate which rays should be considered.
         """
-        if timedelta is None:
-            timedelta = self.timedelta_exposure
-
         where = where & rays.unvignetted
 
-        sgn = np.sign(rays.intensity)
-
-        rays = dataclasses.replace(
-            rays,
-            intensity=np.abs(rays.intensity) * timedelta,
-        )
-
         normal = self.sag.normal(rays.position)
 
-        rays = self.material.signal(
-            rays=rays,
+        # Cosine of the refracted angle inside the light-sensitive region,
+        # folding in each ray's ambient index of refraction. This is generally
+        # complex, so its real and imaginary parts are binned separately below.
+        direction = self.material.direction_refracted(
+            wavelength=rays.wavelength,
+            direction=rays.direction,
+            n=rays.n,
             normal=normal,
-            noise=noise,
-        )
-
-        rays = dataclasses.replace(
-            rays,
-            intensity=sgn * rays.intensity,
         )
 
-        rays = self.material.charge_diffusion(
-            rays=rays,
-            normal=normal,
-        )
+        flux = rays.intensity * where
 
         bins = na.SpectralPositionalVectorArray(
             wavelength=wavelength,
@@ -155,15 +138,149 @@ def readout(
                 num=self.num_pixel + 1,
             ),
         )
+        a = na.SpectralPositionalVectorArray(
+            wavelength=rays.wavelength,
+            position=rays.position.xy,
+        )
 
-        return na.histogram(
-            a=na.SpectralPositionalVectorArray(
-                wavelength=rays.wavelength,
-                position=rays.position.xy,
-            ),
-            bins=bins,
+        image = na.histogram(a, bins=bins, axis=axis, weights=flux)
+        moment_real = na.histogram(
+            a, bins=bins, axis=axis, weights=flux * np.real(direction)
+        )
+        moment_imag = na.histogram(
+            a, bins=bins, axis=axis, weights=flux * np.imag(direction)
+        )
+
+        # flux-weighted mean refracted cosine; empty pixels carry no flux, so
+        # assume normal incidence (cosine of 1).
+        nonempty = image.outputs > 0
+        with np.errstate(invalid="ignore", divide="ignore"):
+            direction_real = np.where(nonempty, moment_real.outputs / image.outputs, 1)
+            direction_imag = np.where(nonempty, moment_imag.outputs / image.outputs, 0)
+        direction = direction_real + direction_imag * 1j
+
+        return image, direction
+
+    def expose(
+        self,
+        image: na.FunctionArray[
+            na.SpectralPositionalVectorArray,
+            na.AbstractScalar,
+        ],
+        direction: float | na.AbstractScalar = 1,
+        axis_wavelength: None | str = None,
+        timedelta: None | u.Quantity | na.AbstractScalar = None,
+        noise: bool = True,
+    ) -> na.FunctionArray[
+        na.SpectralPositionalVectorArray,
+        na.AbstractScalar,
+    ]:
+        """
+        Convert a per-pixel photon image into the electrons measured by the sensor.
+
+        This is the detector-physics step shared by every optical system: it
+        operates on a pixel grid (a photon image plus a refracted-cosine map),
+        so it can be driven by :meth:`collect` for a ray-traced system, or by
+        any model that produces a per-pixel photon image directly.
+
+        The photon flux is multiplied by the exposure time and converted to
+        electrons using
+        :meth:`~optika.sensors.materials.AbstractSensorMaterial.signal`, which
+        applies the quantum efficiency, noise, and charge-diffusion models.
+
+        Parameters
+        ----------
+        image
+            The expected photon flux incident on each pixel, as a function of
+            wavelength and pixel position.
+            The wavelength inputs (``image.inputs.wavelength``) must be the
+            bin *edges*, not the centers.
+        direction
+            The cosine of the refracted angle inside the light-sensitive region
+            in each pixel, as produced by :meth:`collect`.
+        axis_wavelength
+            The logical axis of `image` corresponding to changing wavelength.
+            If :obj:`None` (the default), ``image.inputs.wavelength`` must have
+            only one logical axis.
+        timedelta
+            The exposure time of the measurement.
+            If :obj:`None` (the default), the value in :attr:`timedelta_exposure`
+            will be used.
+        noise
+            Whether to add shot noise and intrinsic sensor noise to the result.
+        """
+        if axis_wavelength is None:
+            shape_wavelength = na.shape(image.inputs.wavelength)
+            if len(shape_wavelength) != 1:  # pragma: nocover
+                raise ValueError(
+                    f"if `axis_wavelength` is `None`, `image.inputs.wavelength` "
+                    f"must have exactly one logical axis, got {shape_wavelength}."
+                )
+            (axis_wavelength,) = shape_wavelength
+
+        if timedelta is None:
+            timedelta = self.timedelta_exposure
+
+        photons = image.outputs * timedelta
+
+        electrons = self.material.signal(
+            photons=photons,
+            wavelength=image.inputs.wavelength.cell_centers(axis_wavelength),
+            direction=direction,
+            width_pixel=self.width_pixel,
+            axis_xy=(self.axis_pixel.x, self.axis_pixel.y),
+            noise=noise,
+        )
+
+        return dataclasses.replace(image, outputs=electrons)
+
+    def measure(
+        self,
+        rays: optika.rays.RayVectorArray,
+        wavelength: na.AbstractScalar,
+        axis: None | str | Sequence[str] = None,
+        where: bool | na.AbstractScalar = True,
+        timedelta: None | u.Quantity | na.AbstractScalar = None,
+        noise: bool = True,
+    ) -> na.FunctionArray[
+        na.SpectralPositionalVectorArray,
+        na.AbstractScalar,
+    ]:
+        """
+        Bin a set of rays onto the pixel grid and convert them to the electrons
+        measured by the sensor.
+
+        This composes :meth:`collect` (gather rays into the pixel grid) with
+        :meth:`expose` (apply the detector physics).
+
+        Parameters
+        ----------
+        rays
+            A set of incident rays in local coordinates to measure.
+        wavelength
+            The edges of the wavelength bins to sample.
+        axis
+            The logical axes along which to collect photons.
+        where
+            A boolean mask used to indicate which rays should be considered.
+        timedelta
+            The exposure time of the measurement.
+            If :obj:`None` (the default), the value in :attr:`timedelta_exposure`
+            will be used.
+        noise
+            Whether to add shot noise and intrinsic sensor noise to the result.
+        """
+        image, direction = self.collect(
+            rays=rays,
+            wavelength=wavelength,
             axis=axis,
-            weights=rays.intensity * where,
+            where=where,
+        )
+        return self.expose(
+            image,
+            direction=direction,
+            timedelta=timedelta,
+            noise=noise,
         )
 
 

optika/sensors/_sensors_test.py

@@ -26,6 +26,7 @@ def test_timedelta_exposure(self, a: optika.sensors.AbstractImagingSensor):
                 intensity=100 * u.photon / u.s,
                 wavelength=500 * u.nm,
                 position=na.Cartesian3dVectorArray() * u.mm,
+                direction=na.Cartesian3dVectorArray(0, 0, 1),
             ),
             optika.rays.RayVectorArray(
                 intensity=na.random.poisson(100, shape_random=dict(t=11)) * u.erg / u.s,
@@ -35,6 +36,7 @@ def test_timedelta_exposure(self, a: optika.sensors.AbstractImagingSensor):
                     y=na.random.uniform(-1, 1, shape_random=dict(t=11)) * u.mm,
                     z=0 * u.mm,
                 ),
+                direction=na.Cartesian3dVectorArray(0, 0, 1),
             ),
         ],
     )
@@ -44,13 +46,13 @@ def test_timedelta_exposure(self, a: optika.sensors.AbstractImagingSensor):
             na.linspace(500, 600, axis="wavelength", num=11) * u.nm,
         ],
     )
-    def test_readout(
+    def test_measure(
         self,
         a: optika.sensors.AbstractImagingSensor,
         rays: optika.rays.RayVectorArray,
         wavelength: u.Quantity | na.AbstractScalar,
     ):
-        result = a.readout(rays, wavelength)
+        result = a.measure(rays, wavelength)
         assert isinstance(result, na.FunctionArray)
         assert isinstance(result.inputs, na.SpectralPositionalVectorArray)
         assert isinstance(result.outputs, na.AbstractScalar)

optika/sensors/materials/_e2v_ccd203/_e2v_ccd203.py

@@ -88,10 +88,8 @@ def e2v_ccd203(
 
         # Evaluate the fitted QE using the given wavelengths
         eqe_fit = material.quantum_efficiency_effective(
-            rays=optika.rays.RayVectorArray(
-                wavelength=wavelength_fit,
-                direction=na.Cartesian3dVectorArray(0, 0, 1),
-            ),
+            wavelength=wavelength_fit,
+            direction=na.Cartesian3dVectorArray(0, 0, 1),
             normal=na.Cartesian3dVectorArray(0, 0, -1),
         )
 
@@ -154,11 +152,7 @@ def e2v_ccd203(
         # Compute the width of the charge diffusion kernel
         # for each wavelength.
         width = material.width_charge_diffusion(
-            rays=optika.rays.RayVectorArray(
-                wavelength=wavelength_fit,
-                direction=na.Cartesian3dVectorArray(0, 0, 1),
-            ),
-            normal=na.Cartesian3dVectorArray(0, 0, -1),
+            wavelength=wavelength_fit,
         )
 
         # Plot the results