From 5efed43c8bd5441ef0a6b47bb710332a8d5d3694 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Mon, 1 Jun 2026 19:11:53 -0600 Subject: [PATCH 01/19] Modified `optika.sensors.signal()` to model charge diffusion using a Monte-Carlo simulation. --- optika/sensors/materials/_materials.py | 20 +- .../_ramanathan_2020/_ramanathan_2020.py | 195 ++++++++++++------ optika/sensors/materials/_stern_1994.py | 1 + 3 files changed, 147 insertions(+), 69 deletions(-) diff --git a/optika/sensors/materials/_materials.py b/optika/sensors/materials/_materials.py index 38606333..a854daab 100644 --- a/optika/sensors/materials/_materials.py +++ b/optika/sensors/materials/_materials.py @@ -13,6 +13,7 @@ _thickness_oxide, _thickness_implant, _thickness_substrate, + _width_pixel, _cce_backsurface, ) from ._ramanathan_2020 import ( @@ -945,15 +946,18 @@ def electrons_measured_approx( return result -_absorbance = absorbance +_transmittance = transmittance def signal( photons_expected: u.Quantity | na.AbstractScalar, wavelength: u.Quantity | na.ScalarArray, - absorbance: None | float | na.AbstractScalar = None, + transmittance: None | float | na.AbstractScalar = None, absorption: None | u.Quantity | na.AbstractScalar = None, thickness_implant: u.Quantity | na.AbstractScalar = _thickness_implant, + thickness_depletion: u.Quantity | na.AbstractScalar = _thickness_substrate, + thickness_substrate: None | na.AbstractScalar = _thickness_substrate, + width_pixel: u.Quantity | na.AbstractScalar = _width_pixel, cce_backsurface: u.Quantity | na.AbstractScalar = _cce_backsurface, temperature: u.Quantity | na.ScalarArray = 300 * u.K, method: Literal["exact", "approx", "expected"] = "exact", @@ -974,9 +978,9 @@ def signal( The `expected` number of photons incident on the detector surface. wavelength The vacuum wavelength of the absorbed photons. - absorbance - The fraction of incident energy absorbed by the light-sensitive layer - of the detector computed using the average of :func:`absorbance`. + transmittance + The fraction of incident energy transmitted to the light-sensitive layer + of the detector. If :obj:`None` (the default), the result of :func:`absorbance` called with default values will be used. absorption @@ -1048,8 +1052,8 @@ def signal( ax.set_ylabel(f"variance-to-mean ratio ({electrons.unit:latex_inline})"); """ - if absorbance is None: - absorbance = _absorbance(wavelength).average + if transmittance is None: + transmittance = _transmittance(wavelength).average if absorption is None: absorption = optika.chemicals.Chemical("Si").absorption(wavelength) @@ -1064,7 +1068,7 @@ def signal( thickness_implant=thickness_implant, cce_backsurface=cce_backsurface, ) - return iqy * absorbance * cce * photons_expected.to(u.ph) + return iqy * transmittance * cce * photons_expected.to(u.ph) photons_absorbed_expected = absorbance * photons_expected.to(u.ph) photons_absorbed = na.random.poisson( diff --git a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py index af58fa3c..75a316ca 100644 --- a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py +++ b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py @@ -8,6 +8,8 @@ import optika from .._stern_1994 import ( _thickness_implant, + _thickness_substrate, + _width_pixel, _cce_backsurface, ) @@ -464,12 +466,16 @@ def probability_of_n_pairs( def electrons_measured( - photons_absorbed: u.Quantity | na.AbstractScalar, + photons_transmitted: u.Quantity | na.AbstractScalar, wavelength: u.Quantity | na.ScalarArray, absorption: None | u.Quantity | na.AbstractScalar = None, thickness_implant: u.Quantity | na.AbstractScalar = _thickness_implant, + thickness_depletion: u.Quantity | na.AbstractScalar = _thickness_substrate, + thickness_substrate: u.Quantity | na.AbstractScalar = _thickness_substrate, + width_pixel: u.Quantity | na.ScalarArray = _width_pixel, cce_backsurface: u.Quantity | na.AbstractScalar = _cce_backsurface, temperature: u.Quantity | na.ScalarArray = 300 * u.K, + axis_xy: None | tuple[str, str] = None, shape_random: None | dict[str, int] = None, ) -> na.AbstractScalar: r""" @@ -482,20 +488,34 @@ def electrons_measured( Parameters ---------- - photons_absorbed - The number of photons absorbed by the light-sensitive layer of the sensor. + photons_transmitted + The number of photons transmitted into the light-sensitive region + of the sensor. wavelength The vacuum wavelength of the absorbed photons. absorption The absorption coefficient of silicon at the given wavelength. thickness_implant The thickness of the implant layer, where partial-charge collection occurs. + thickness_depletion + The thickness of the depletion region, the region with significant electric + field. + This is sent to `thickness_substrate` by default, which is equivalent to + there being no field-free region in which charge diffusion occurs. + thickness_substrate + The thickness of the entire light-sensitive region of the device. + width_pixel + The size of a single pixel on the sensor. cce_backsurface The differential charge collection efficiency on the back surface of the sensor. temperature The temperature of the silicon detector. Default is room temperature. + axis_xy + The two logical axes corresponding to the pixel grid of the sensor + along which electrons will diffuse. + If :obj:`None` (the default), there is no charge diffusion. shape_random Additional shape used to specify the number of samples to draw. @@ -518,7 +538,7 @@ def electrons_measured( # Define the expected number of photons # for each experiment - photons_absorbed = (100 * u.photon).astype(int) + photons_transmitted = (100 * u.photon).astype(int) # Define the wavelength at which to sample the distribution wavelength = 5.9 * u.keV @@ -526,7 +546,7 @@ def electrons_measured( # Compute the actual number of electrons measured for each experiment electrons = optika.sensors.electrons_measured( - photons_absorbed=photons_absorbed, + photons_transmitted=photons_transmitted, wavelength=wavelength, shape_random=dict(experiment=num_experiments), ) @@ -565,18 +585,34 @@ def electrons_measured( shape_random = dict() shape = na.broadcast_shapes( - na.shape(photons_absorbed), + na.shape(photons_transmitted), na.shape(wavelength), na.shape(absorption), na.shape(thickness_implant), + na.shape(thickness_depletion), + na.shape(thickness_substrate), + na.shape(width_pixel), na.shape(cce_backsurface), na.shape(temperature), shape_random, ) - photons_absorbed = na.broadcast_to(photons_absorbed, shape) + if axis_xy is not None: + axis_x, axis_y = axis_xy + shape[axis_x] = shape.pop(axis_x) + shape[axis_y] = shape.pop(axis_y) + else: + axis_x = "__dummy_x__" + axis_y = "__dummy_y__" + shape[axis_x] = 1 + shape[axis_y] = 1 + + photons_transmitted = na.broadcast_to(photons_transmitted, shape) absorption = na.broadcast_to(absorption, shape) thickness_implant = na.broadcast_to(thickness_implant, shape) + thickness_depletion = na.broadcast_to(thickness_depletion, shape) + thickness_substrate = na.broadcast_to(thickness_substrate, shape) + width_pixel = na.broadcast_to(width_pixel, shape) cce_backsurface = na.broadcast_to(cce_backsurface, shape) if not isinstance(cce_backsurface, u.Quantity): @@ -600,10 +636,13 @@ def electrons_measured( fano_inf = fano_inf.broadcast_to(shape) result = _electrons_measured_quantity( - photons_absorbed=photons_absorbed.ndarray, + photons_transmitted=photons_transmitted.ndarray, wavelength=wavelength.ndarray, absorption=absorption.ndarray, thickness_implant=thickness_implant.ndarray, + thickness_depletion=thickness_depletion.ndarray, + thickness_substrate=thickness_substrate.ndarray, + width_pixel=width_pixel, cce_backsurface=cce_backsurface.ndarray, p_n=p_n.ndarray, n=n.ndarray, @@ -618,10 +657,13 @@ def electrons_measured( def _electrons_measured_quantity( - photons_absorbed: u.Quantity, + photons_transmitted: u.Quantity, wavelength: u.Quantity, absorption: u.Quantity, thickness_implant: u.Quantity, + thickness_depletion: u.Quantity, + thickness_substrate: u.Quantity, + width_pixel: u.Quantity, cce_backsurface: u.Quantity, p_n: np.ndarray, n: np.ndarray, @@ -630,32 +672,44 @@ def _electrons_measured_quantity( ) -> u.Quantity: shape = np.broadcast_shapes( - photons_absorbed.shape, + photons_transmitted.shape, absorption.shape, thickness_implant.shape, + thickness_depletion.shape, + thickness_substrate.shape, cce_backsurface.shape, + width_pixel.shape, ) + num_x = shape[~1] + num_y = shape[~0] + unit_length = u.mm - photons_absorbed = photons_absorbed.to_value(u.photon) + photons_transmitted = photons_transmitted.to_value(u.photon) energy = wavelength.to(u.eV, equivalencies=u.spectral()) absorption = absorption.to_value(1 / unit_length) thickness_implant = thickness_implant.to_value(unit_length) + thickness_depletion = thickness_depletion.to_value(unit_length) + thickness_substrate = thickness_substrate.to_value(unit_length) + width_pixel = width_pixel.to_value(unit_length) cce_backsurface = cce_backsurface.to_value(u.dimensionless_unscaled) energy_pair_inf = energy_pair_inf.to_value(u.eV) fano_inf = fano_inf.to_value(u.electron / u.photon) result = _electrons_measured_numba( - photons_absorbed=photons_absorbed.reshape(-1), - energy=energy.reshape(-1), - absorption=absorption.reshape(-1), - thickness_implant=thickness_implant.reshape(-1), - cce_backsurface=cce_backsurface.reshape(-1), - p_n=p_n.reshape(-1, p_n.shape[~0]), - n=n.reshape(-1, n.shape[~0]), - energy_pair_inf=energy_pair_inf.reshape(-1), - fano_inf=fano_inf.reshape(-1), + photons_transmitted=photons_transmitted.reshape(-1, num_x, num_y), + energy=energy.reshape(-1, num_x, num_y), + absorption=absorption.reshape(-1, num_x, num_y), + thickness_implant=thickness_implant.reshape(-1, num_x, num_y), + thickness_depletion=thickness_depletion.reshape(-1, num_x, num_y), + thickness_substrate=thickness_substrate.reshape(-1, num_x, num_y), + width_pixel=width_pixel.reshape(-1, num_x, num_y), + cce_backsurface=cce_backsurface.reshape(-1, num_x, num_y), + p_n=p_n.reshape(-1, num_x, num_y, p_n.shape[~0]), + n=n.reshape(-1, num_x, num_y, n.shape[~0]), + energy_pair_inf=energy_pair_inf.reshape(-1, num_x, num_y), + fano_inf=fano_inf.reshape(-1, num_x, num_y), ) result = result.reshape(shape) @@ -667,73 +721,92 @@ def _electrons_measured_quantity( @numba.njit(fastmath=True, parallel=True) def _electrons_measured_numba( - photons_absorbed: np.ndarray, + photons_transmitted: np.ndarray, energy: np.ndarray, absorption: np.ndarray, thickness_implant: np.ndarray, + thickness_depletion: np.ndarray, + thickness_substrate: np.ndarray, + width_pixel: np.ndarray, cce_backsurface: np.ndarray, p_n: np.ndarray, n: np.ndarray, energy_pair_inf: np.ndarray, fano_inf: np.ndarray, -) -> np.ndarray: # pragma: nocover +) -> np.ndarray: + + num_i, num_x, num_y, num_n = p_n.shape - num_i, num_n = p_n.shape - - result = np.empty(num_i) + result = np.zeros((num_i, num_x, num_y)) for i in numba.prange(num_i): + for x in range(num_x): + for y in range(num_y): + num_photon = int(photons_transmitted[i, x, y]) + energy_i = energy[i, x, y] + a = absorption[i, x, y] + W = thickness_implant[i, x, y] + h_0 = cce_backsurface[i, x, y] + cmf_i = np.cumsum(p_n[i, x, y]) + n_i = n[i, x, y] + energy_pair_inf_i = energy_pair_inf[i, x, y] + fano_inf_i = fano_inf[i, x, y] + z_substrate = thickness_substrate[i, x, y] + z_ff = z_substrate - thickness_depletion[i, x, y] + + d = 1 / a - num_photon = int(photons_absorbed[i]) - energy_i = energy[i] - a = absorption[i] - W = thickness_implant[i] - h_0 = cce_backsurface[i] - cmf_i = np.cumsum(p_n[i]) - n_i = n[i] - energy_pair_inf_i = energy_pair_inf[i] - fano_inf_i = fano_inf[i] + mean_inf = energy_i / energy_pair_inf_i + std_inf = math.sqrt(fano_inf_i * mean_inf) - d = 1 / a + low_energy = energy_i <= 50 - num_electron_i = 0 + for j in range(num_photon): + if low_energy: + x_ij = random.uniform(0, 1) - mean_inf = energy_i / energy_pair_inf_i - std_inf = math.sqrt(fano_inf_i * mean_inf) + for k_ij in range(num_n): + if cmf_i[k_ij] > x_ij: + break - low_energy = energy_i <= 50 + n_ij = n_i[k_ij] - for j in range(num_photon): + else: + n_ij = random.normalvariate( + mu=mean_inf, + sigma=std_inf, + ) + n_ij = round(n_ij) - if low_energy: + y_ij = random.uniform(0, 1) + z_ij = -d * math.log(1 - y_ij) - x_ij = random.uniform(0, 1) + if z_ij < W: + h_ij = h_0 + (1 - h_0) * z_ij / W + else: + h_ij = 1 - for k_ij in range(num_n): - if cmf_i[k_ij] > x_ij: - break + m_ij = np.random.binomial(n=n_ij, p=h_ij) - n_ij = n_i[k_ij] + u = random.uniform(-0.5, 0.5) + v = random.uniform(-0.5, 0.5) - else: - n_ij = random.normalvariate( - mu=mean_inf, - sigma=std_inf, - ) - n_ij = round(n_ij) + for e in range(m_ij): + if z_ij < z_ff: + w = z_ff * math.sqrt(1 - z_ij / z_ff) / width_pixel - y_ij = random.uniform(0, 1) - z_ij = -d * math.log(1 - y_ij) + p = random.gauss(u, w) + q = random.gauss(v, w) - if z_ij < W: - h_ij = h_0 + (1 - h_0) * z_ij / W - else: - h_ij = 1 + p = round(p) + q = round(q) - m_ij = np.random.binomial(n=n_ij, p=h_ij) + elif z_ij > z_substrate: + continue - num_electron_i = num_electron_i + m_ij + else: + p = q = 0 - result[i] = num_electron_i + result[i, (x + p) % num_x, (y + q) % num_y] += 1 return result diff --git a/optika/sensors/materials/_stern_1994.py b/optika/sensors/materials/_stern_1994.py index 99f66151..ee336e19 100644 --- a/optika/sensors/materials/_stern_1994.py +++ b/optika/sensors/materials/_stern_1994.py @@ -3,4 +3,5 @@ _thickness_oxide = 50 * u.AA _thickness_substrate = 7 * u.um _thickness_implant = 2317 * u.AA +_width_pixel = 27 * u.um _cce_backsurface = 0.21 From 8ce94faac7174439a8edc3cce7cd8347ab11fc03 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 2 Jun 2026 15:29:58 -0600 Subject: [PATCH 02/19] wip --- optika/sensors/materials/_materials.py | 202 +++++++++++------- .../_ramanathan_2020/_ramanathan_2020.py | 9 +- 2 files changed, 134 insertions(+), 77 deletions(-) diff --git a/optika/sensors/materials/_materials.py b/optika/sensors/materials/_materials.py index a854daab..716ac7fd 100644 --- a/optika/sensors/materials/_materials.py +++ b/optika/sensors/materials/_materials.py @@ -960,7 +960,8 @@ def signal( width_pixel: u.Quantity | na.AbstractScalar = _width_pixel, cce_backsurface: u.Quantity | na.AbstractScalar = _cce_backsurface, temperature: u.Quantity | na.ScalarArray = 300 * u.K, - method: Literal["exact", "approx", "expected"] = "exact", + method: Literal["monte-carlo", "expected"] = "monte-carlo", + axis_xy: None | tuple[str, str] = None, shape_random: None | dict[str, int] = None, ) -> na.AbstractScalar: r""" @@ -991,6 +992,16 @@ def signal( thickness_implant The thickness of the implant layer. Default is the value given in :cite:t:`Stern1994`. + thickness_depletion + The thickness of the depletion region, the region with significant electric + field. + If :obj:`None` (the default), this is set to the same value as + `thickness_substrate`. + thickness_substrate + The thickness of the entire light-sensitive region of the device. + The default + width_pixel + The size of a single pixel on the sensor. cce_backsurface The differential charge collection efficiency on the back surface of the sensor. @@ -1005,6 +1016,10 @@ def signal( of photons is high. The `expected` method does not add any noise to the signal and just returns the expected number of electrons. + axis_xy + The two logical axes corresponding to the pixel grid of the sensor + along which electrons will diffuse. + If :obj:`None` (the default), there is no charge diffusion. shape_random Additional shape used to specify the number of samples to draw. @@ -1063,33 +1078,36 @@ def signal( wavelength=wavelength, temperature=temperature, ) + absorbance = transmittance * np.exp(-absorption * thickness_substrate) cce = charge_collection_efficiency( absorption=absorption, thickness_implant=thickness_implant, cce_backsurface=cce_backsurface, ) - return iqy * transmittance * cce * photons_expected.to(u.ph) + return iqy * absorbance * cce * photons_expected.to(u.ph) - photons_absorbed_expected = absorbance * photons_expected.to(u.ph) - photons_absorbed = na.random.poisson( - lam=photons_absorbed_expected, - shape_random=shape_random, - ).astype(int) + elif method == "monte-carlo": - kwargs = dict( - photons_absorbed=photons_absorbed, - wavelength=wavelength, - absorption=absorption, - thickness_implant=thickness_implant, - cce_backsurface=cce_backsurface, - temperature=temperature, - shape_random=shape_random, - ) + photons_expected = transmittance * photons_expected.to(u.ph) + photons = na.random.poisson( + lam=photons_expected, + shape_random=shape_random, + ).astype(int) + + return electrons_measured( + photons_transmitted=photons, + wavelength=wavelength, + absorption=absorption, + thickness_implant=thickness_implant, + thickness_depletion=thickness_depletion, + thickness_substrate=thickness_substrate, + width_pixel=width_pixel, + cce_backsurface=cce_backsurface, + temperature=temperature, + axis_xy=axis_xy, + shape_random=shape_random, + ) - if method == "exact": - return electrons_measured(**kwargs) - elif method == "approx": - return electrons_measured_approx(**kwargs) else: # pragma: nocover raise ValueError(f"Unrecognized method: {method}") @@ -1262,22 +1280,31 @@ class AbstractSensorMaterial( @abc.abstractmethod def signal( self, - rays: optika.rays.RayVectorArray, + photons: u.Quantity | na.AbstractScalar, + wavelength: u.Quantity | na.AbstractScalar, + direction: na.AbstractCartesian3dVectorArray, normal: na.AbstractCartesian3dVectorArray, + axis_xy: None | tuple[str, str] = None, noise: bool = True, - ) -> optika.rays.RayVectorArray: + ) -> na.AbstractScalar: """ Given a set of absorbed rays, compute the number of electrons measured by the sensor using :func:`signal`. Parameters ---------- - rays - The rays absorbed by the light-sensitive silicon layer. - The :attr:`optika.rays.RayVectorArray.intensity` field should - either be in units of photons or energy. + photons + The number of photons incident on each pixel. + wavelength + An assumed grid of wavelengths for the incident photons. + direction + An assumed propagation direction for the incident photons. normal The vector perpendicular to the surface of the sensor. + axis_xy + The two logical axes corresponding to the pixel grid of the sensor. + If provided, charge diffusion will occur along these two axes. + If :obj:`None` (the default), no diffusion is performed. noise Whether to add noise to the result. """ @@ -1298,7 +1325,7 @@ def photons_incident( Parameters ---------- electrons - The number of electrons measured by the sensor. + The number of electrons measured by each pixel. wavelength An assumed grid of wavelengths for the incident photons. direction @@ -1307,24 +1334,6 @@ def photons_incident( The vector perpendicular to the surface of the sensor. """ - @abc.abstractmethod - def charge_diffusion( - self, - rays: optika.rays.RayVectorArray, - normal: na.AbstractCartesian3dVectorArray, - ) -> optika.rays.RayVectorArray: - """ - Given a set of incident rays, compute how much the position of each ray - is perturbed due to charge diffusion within the sensor. - - Parameters - ---------- - rays - The rays incident on the sensor surface. - normal - The vector perpendicular to the surface of the sensor. - """ - @dataclasses.dataclass(eq=False, repr=False) class IdealSensorMaterial( @@ -1338,15 +1347,18 @@ class IdealSensorMaterial( def signal( self, - rays: optika.rays.RayVectorArray, + photons: u.Quantity | na.AbstractScalar, + wavelength: u.Quantity | na.AbstractScalar, + direction: na.AbstractCartesian3dVectorArray, normal: na.AbstractCartesian3dVectorArray, - noise: bool = False, + axis_xy: None | tuple[str, str] = None, + noise: bool = True, ) -> optika.rays.RayVectorArray: - intensity = rays.intensity - if not intensity.unit.is_equivalent(u.photon): + + if not photons.unit.is_equivalent(u.photon): h = astropy.constants.h c = astropy.constants.c - intensity = intensity / (h * c / rays.wavelength) * u.photon + intensity = photons / (h * c / wavelength) * u.photon if noise: intensity = na.random.poisson(intensity.to(u.ph)).astype(int) @@ -1354,9 +1366,7 @@ def signal( electrons = intensity * u.electron / u.photon electrons = electrons.to(u.electron) - result = dataclasses.replace(rays, intensity=electrons) - - return result + return electrons def photons_incident( self, @@ -1367,13 +1377,6 @@ def photons_incident( ) -> na.AbstractScalar: return electrons * u.photon / u.electron - def charge_diffusion( - self, - rays: optika.rays.RayVectorArray, - normal: na.AbstractCartesian3dVectorArray, - ) -> optika.rays.RayVectorArray: - return rays - @dataclasses.dataclass(eq=False, repr=False) class AbstractSiliconSensorMaterial( @@ -1503,8 +1506,7 @@ def depletion(self) -> AbstractDepletionModel: def width_charge_diffusion( self, - rays: optika.rays.RayVectorArray, - normal: na.AbstractCartesian3dVectorArray, + wavelength: u.Quantity | na.AbstractScalar, ) -> na.AbstractScalar: """ The standard deviation of the charge diffusion kernel for this sensor. @@ -1512,21 +1514,62 @@ def width_charge_diffusion( Parameters ---------- - rays - The rays incident on the sensor surface. - normal - The vector perpendicular to the surface of the sensor. + wavelength + The wavelength of the incident light in vacuum. """ return optika.sensors.charge_diffusion( - self._chemical.absorption(rays.wavelength), + self._chemical.absorption(wavelength), thickness_substrate=self.thickness_substrate, thickness_depletion=self.depletion.thickness, ) + def transmittance( + self, + wavelength: u.Quantity | na.AbstractScalar, + direction: None | na.AbstractCartesian3dVectorArray = None, + index_refraction: float | na.AbstractScalar = 1, + normal: None | na.AbstractCartesian3dVectorArray = None, + ) -> optika.vectors.PolarizationVectorArray: + """ + Compute the fraction of energy transmitted to the light-sensitive region + of the sensor. + + Parameters + ---------- + wavelength + The wavelength of the incident light in vacuum. + direction + The propagation direction of the incident light in the ambient medium. + If :obj:`None` (default), normal incidence (:math:`\hat{z}`) is assumed. + index_refraction + The complex index of refraction of the ambient medium. + normal + The vector perpendicular to the surface of the CCD sensor. + """ + if direction is None: + direction = na.Cartesian3dVectorArray(0, 0, 1) + + if normal is None: + normal = na.Cartesian3dVectorArray(0, 0, 1) + + return transmittance( + wavelength=wavelength, + direction=-direction @ normal, + n=index_refraction, + thickness_oxide=self.thickness_oxide, + thickness_substrate=self.thickness_substrate, + chemical_oxide=self._chemical_oxide, + chemical_substrate=self._chemical, + roughness_oxide=self.roughness_oxide, + roughness_substrate=self.roughness_substrate, + ) + def absorbance( self, - rays: optika.rays.AbstractRayVectorArray, - normal: na.AbstractCartesian3dVectorArray, + wavelength: u.Quantity | na.AbstractScalar, + direction: None | na.AbstractCartesian3dVectorArray = None, + index_refraction: float | na.AbstractScalar = 1, + normal: None | na.AbstractCartesian3dVectorArray = None, ) -> optika.vectors.PolarizationVectorArray: """ Compute the fraction of energy absorbed by the light-sensitive region @@ -1534,15 +1577,26 @@ def absorbance( Parameters ---------- - rays - The light rays incident on the CCD surface. + wavelength + The wavelength of the incident light in vacuum. + direction + The propagation direction of the incident light in the ambient medium. + If :obj:`None` (default), normal incidence (:math:`\hat{z}`) is assumed. + n + The complex index of refraction of the ambient medium. normal The vector perpendicular to the surface of the CCD sensor. """ + if direction is None: + direction = na.Cartesian3dVectorArray(0, 0, 1) + + if normal is None: + normal = na.Cartesian3dVectorArray(0, 0, 1) + return absorbance( - wavelength=rays.wavelength, - direction=-rays.direction @ normal, - n=rays.index_refraction, + wavelength=wavelength, + direction=-direction @ normal, + n=index_refraction, thickness_oxide=self.thickness_oxide, thickness_substrate=self.thickness_substrate, chemical_oxide=self._chemical_oxide, diff --git a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py index 75a316ca..795cda98 100644 --- a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py +++ b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py @@ -470,7 +470,7 @@ def electrons_measured( wavelength: u.Quantity | na.ScalarArray, absorption: None | u.Quantity | na.AbstractScalar = None, thickness_implant: u.Quantity | na.AbstractScalar = _thickness_implant, - thickness_depletion: u.Quantity | na.AbstractScalar = _thickness_substrate, + thickness_depletion: None | u.Quantity | na.AbstractScalar = None, thickness_substrate: u.Quantity | na.AbstractScalar = _thickness_substrate, width_pixel: u.Quantity | na.ScalarArray = _width_pixel, cce_backsurface: u.Quantity | na.AbstractScalar = _cce_backsurface, @@ -500,8 +500,8 @@ def electrons_measured( thickness_depletion The thickness of the depletion region, the region with significant electric field. - This is sent to `thickness_substrate` by default, which is equivalent to - there being no field-free region in which charge diffusion occurs. + If :obj:`None` (the default), this is set to the same value as + `thickness_substrate`. thickness_substrate The thickness of the entire light-sensitive region of the device. width_pixel @@ -581,6 +581,9 @@ def electrons_measured( if absorption is None: absorption = optika.chemicals.Chemical("Si").absorption(wavelength) + if thickness_depletion is None: + thickness_depletion = thickness_substrate + if shape_random is None: shape_random = dict() From feb22beb9e92d529328f208e459d39e0041fb9ed Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Sat, 6 Jun 2026 18:52:03 -0600 Subject: [PATCH 03/19] wip --- optika/sensors/__init__.py | 2 + optika/sensors/_sensors.py | 178 +++-- optika/sensors/_sensors_test.py | 2 + optika/sensors/materials/_materials.py | 611 ++++++++++++------ optika/sensors/materials/_materials_test.py | 241 +++---- .../_ramanathan_2020/_ramanathan_2020.py | 66 +- .../_ramanathan_2020/_ramanathan_2020_test.py | 16 +- 7 files changed, 725 insertions(+), 391 deletions(-) diff --git a/optika/sensors/__init__.py b/optika/sensors/__init__.py index 42f8d154..b4139dd7 100644 --- a/optika/sensors/__init__.py +++ b/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", diff --git a/optika/sensors/_sensors.py b/optika/sensors/_sensors.py index 58adf079..8d266e28 100644 --- a/optika/sensors/_sensors.py +++ b/optika/sensors/_sensors.py @@ -81,70 +81,41 @@ def aperture(self): half_width=self.width_pixel * self.num_pixel / 2, ) - def readout( + def bin_rays( 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 incidence angle in each pixel: the two + quantities :meth:`expose` needs. Performing the ray-to-cosine projection + here is what lets :meth:`expose` be shared with systems that have no rays. 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) + direction = -(rays.direction @ normal) - rays = self.material.signal( - rays=rays, - 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 +126,124 @@ 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 = na.histogram(a, bins=bins, axis=axis, weights=flux * direction) + + # flux-weighted mean cosine; empty pixels carry no flux, assume normal incidence + nonempty = image.outputs > 0 + direction = np.where(nonempty, moment.outputs / image.outputs, 1) + + return image, direction + + def expose( + self, + image: na.FunctionArray[ + na.SpectralPositionalVectorArray, + na.AbstractScalar, + ], + direction: float | na.AbstractScalar = 1, + 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 an incidence-cosine map), + so it can be driven by :meth:`bin_rays` 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. + direction + The cosine of the incidence angle in each pixel. + 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 timedelta is None: + timedelta = self.timedelta_exposure + + photons = image.outputs * timedelta + sgn = np.sign(photons) + + electrons = self.material.signal( + photons=np.abs(photons), + wavelength=image.inputs.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=sgn * electrons) + + def readout( + 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:`bin_rays` (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.bin_rays( + rays=rays, + wavelength=wavelength, axis=axis, - weights=rays.intensity * where, + where=where, + ) + return self.expose( + image, + direction=direction, + timedelta=timedelta, + noise=noise, ) diff --git a/optika/sensors/_sensors_test.py b/optika/sensors/_sensors_test.py index 4256e0be..a5c2995c 100644 --- a/optika/sensors/_sensors_test.py +++ b/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), ), ], ) diff --git a/optika/sensors/materials/_materials.py b/optika/sensors/materials/_materials.py index 716ac7fd..f937c393 100644 --- a/optika/sensors/materials/_materials.py +++ b/optika/sensors/materials/_materials.py @@ -34,6 +34,7 @@ "quantum_yield_ideal", "fano_factor", "fano_factor_inf", + "absorption_effective", "absorbance", "charge_collection_efficiency", "quantum_efficiency_effective", @@ -50,10 +51,51 @@ ] +def absorption_effective( + wavelength: u.Quantity | na.AbstractScalar, + n_substrate: complex | na.AbstractScalar, + direction_substrate: complex | na.AbstractScalar, +) -> na.AbstractScalar: + r""" + The absorption coefficient per unit perpendicular depth for a wave + refracted into the (absorbing) light-sensitive region at oblique incidence. + + This is the rigorous generalization of :math:`4 \pi \, \text{Im}(n) / \lambda` + to oblique incidence: it uses the complex refracted cosine + ``direction_substrate`` so that the geometric :math:`\sec\theta` path-length + factor is no longer needed downstream. Reduces to the normal-incidence + coefficient when ``direction_substrate`` is 1. + + Parameters + ---------- + wavelength + The wavelength of the incident light in vacuum. + n_substrate + The complex index of refraction of the light-sensitive region. + direction_substrate + The (generally complex) cosine of the refracted angle inside the + light-sensitive region, e.g. from :func:`optika.materials.snells_law_scalar`. + + Notes + ----- + For a wave refracted from a transparent ambient medium into an absorbing + medium, the planes of constant amplitude remain parallel to the interface, + so the intensity attenuates along the surface normal with coefficient + + .. math:: + + \alpha_\perp = \frac{4 \pi}{\lambda} \, \text{Im}(n_2 \cos\theta_2), + + where :math:`n_2` is the complex index of the light-sensitive region and + :math:`\cos\theta_2` is its (complex) refracted cosine. + """ + return 4 * np.pi * np.imag(n_substrate * direction_substrate) / wavelength + + def transmittance( wavelength: u.Quantity | na.AbstractScalar, direction: float | na.AbstractScalar = 1, - n: float | na.AbstractScalar = 1, + n: complex | na.AbstractScalar = 1, thickness_oxide: u.Quantity | na.AbstractScalar = _thickness_oxide, thickness_substrate: u.Quantity | na.AbstractScalar = _thickness_substrate, chemical_oxide: str | optika.chemicals.AbstractChemical = "SiO2", @@ -70,11 +112,10 @@ def transmittance( wavelength The wavelength of the incident light in vacuum. direction - The component of the incident light's propagation direction antiparallel - to the surface normal of the sensor. + The cosine of the incidence angle. Default is normal incidence. n - The index of refraction in the ambient medium. + The complex index of refraction in the ambient medium. thickness_oxide The thickness of the oxide layer on the illuminated surface of the sensor. Default is the value given in :cite:t:`Stern1994`. @@ -158,13 +199,14 @@ def transmittance( def absorbance( wavelength: u.Quantity | na.AbstractScalar, direction: float | na.AbstractScalar = 1, - n: float | na.AbstractScalar = 1, + n: complex | na.AbstractScalar = 1, thickness_oxide: u.Quantity | na.AbstractScalar = _thickness_oxide, thickness_substrate: u.Quantity | na.AbstractScalar = _thickness_substrate, chemical_oxide: str | optika.chemicals.AbstractChemical = "SiO2", chemical_substrate: str | optika.chemicals.AbstractChemical = "Si", roughness_oxide: u.Quantity | na.AbstractScalar = 0 * u.nm, roughness_substrate: u.Quantity | na.AbstractScalar = 0 * u.nm, + method: Literal["exact", "Beer-Lambert"] = "Beer-Lambert", ) -> optika.vectors.PolarizationVectorArray: """ The fraction of incident energy absorbed by the light-sensitive @@ -175,11 +217,10 @@ def absorbance( wavelength The wavelength of the incident light in vacuum. direction - The component of the incident light's propagation direction antiparallel - to the surface normal of the sensor. + The cosine of the incidence angle. Default is normal incidence. n - The index of refraction in the ambient medium. + The complex index of refraction in the ambient medium. thickness_oxide The thickness of the oxide layer on the illuminated surface of the sensor. Default is the value given in :cite:t:`Stern1994`. @@ -196,6 +237,13 @@ def absorbance( The RMS roughness the oxide layer surface. roughness_substrate The RMS roughness of the substrate surface. + method + The method to use to compute the absorbance. + If ``exact``, this method allows thin-film interference effects + inside the light-sensitive region. + If ``Beer-Lambert``, this method assumes no interference effects. + These methods only differ in the infrared, where the wavelength is + commensurate with the thickness of the light-sensitive region. Examples -------- @@ -219,8 +267,14 @@ def absorbance( ) # Compute the absorbance vs wavelength - absorbance = optika.sensors.absorbance( + absorbance_exact = optika.sensors.absorbance( + wavelength=wavelength, + method="exact", + ) + + absorbance_beer = optika.sensors.absorbance( wavelength=wavelength, + method="Beer-Lambert", ) # Plot the average absorbance vs. wavelength @@ -233,9 +287,15 @@ def absorbance( ); na.plt.plot( wavelength, - absorbance.average, + absorbance_exact.average, + ax=ax, + label="exact absorbance", + ); + na.plt.plot( + wavelength, + absorbance_beer.average, ax=ax, - label="absorbance", + label="Beer-Lambert absorbance", ); ax.set_xscale("log"); ax.set_xlabel(f"wavelength ({wavelength.unit:latex_inline})"); @@ -248,28 +308,65 @@ def absorbance( if not isinstance(chemical_substrate, optika.chemicals.AbstractChemical): chemical_substrate = optika.chemicals.Chemical(chemical_substrate) - result = optika.materials.layer_absorbance( - index=1, - wavelength=wavelength, - direction=direction, - n=n, - layers=[ - optika.materials.Layer( - chemical=chemical_oxide, - thickness=thickness_oxide, - interface=optika.materials.profiles.ErfInterfaceProfile( - width=roughness_oxide, + if method == "exact": + + result = optika.materials.layer_absorbance( + index=1, + wavelength=wavelength, + direction=direction, + n=n, + layers=[ + optika.materials.Layer( + chemical=chemical_oxide, + thickness=thickness_oxide, + interface=optika.materials.profiles.ErfInterfaceProfile( + width=roughness_oxide, + ), ), - ), - optika.materials.Layer( - chemical=chemical_substrate, - thickness=thickness_substrate, - interface=optika.materials.profiles.ErfInterfaceProfile( - width=roughness_substrate, + optika.materials.Layer( + chemical=chemical_substrate, + thickness=thickness_substrate, + interface=optika.materials.profiles.ErfInterfaceProfile( + width=roughness_substrate, + ), ), - ), - ], - ) + ], + ) + + elif method == "Beer-Lambert": + + _transmittance = transmittance( + wavelength=wavelength, + direction=direction, + n=n, + thickness_oxide=thickness_oxide, + thickness_substrate=thickness_substrate, + chemical_oxide=chemical_oxide, + chemical_substrate=chemical_substrate, + roughness_oxide=roughness_oxide, + roughness_substrate=roughness_substrate, + ) + + n_substrate = chemical_substrate.n(wavelength) + + direction_substrate = optika.materials.snells_law_scalar( + cos_incidence=direction, + index_refraction=n, + index_refraction_new=n_substrate, + ) + + absorption = absorption_effective( + wavelength=wavelength, + n_substrate=n_substrate, + direction_substrate=direction_substrate, + ) + + _transmittance_total = _transmittance * np.exp(-absorption * thickness_substrate) + + result = _transmittance - _transmittance_total + + else: + raise ValueError(f"Method {method} not implemented") return np.real(result) @@ -278,7 +375,6 @@ def charge_collection_efficiency( absorption: u.Quantity | na.AbstractScalar, thickness_implant: u.Quantity | na.AbstractScalar = _thickness_implant, cce_backsurface: u.Quantity | na.AbstractScalar = _cce_backsurface, - cos_incidence: float | na.AbstractScalar = 1, ) -> na.AbstractScalar: r""" Compute the average charge collection efficiency using the differential @@ -287,8 +383,11 @@ def charge_collection_efficiency( Parameters ---------- absorption - The absorption coefficient of the light-sensitive material for the - wavelength of interest. + The absorption coefficient of the light-sensitive material per unit + perpendicular depth. + For oblique incidence, supply the effective coefficient from + :func:`absorption_effective`, which folds in the refracted angle, so no + separate angle argument is needed. thickness_implant The thickness of the implant layer, the layer where recombination can occur. @@ -297,9 +396,6 @@ def charge_collection_efficiency( The differential charge collection efficiency on the back surface of the sensor. Default is the value given in :cite:t:`Stern1994`. - cos_incidence - The cosine of the angle of the incident light's propagation direction - inside the substrate with the surface normal Examples -------- @@ -382,27 +478,14 @@ def charge_collection_efficiency( with the addition of an :math:`e^{-\alpha W}` term which represents photons absorbed inside the epitaxial layer but outside the implant layer. - Equation :eq:`cce` is only valid for normally-incident light. - We can generalize it to obliquely-incident light by making the substitution - - .. math:: - :label: x-oblique - - x \rightarrow \frac{x}{\cos \theta} - - where :math:`\theta` is the angle between the propagation direction - inside the silicon substrate and the normal vector. - - Substituting :eq:`x-oblique` into Equation :eq:`cce` and solving yields - - .. math:: - :label: eqe-oblique - - \text{CCE}(\lambda, \theta) = - \eta_0 - + \left( \frac{1 - \eta_0}{\alpha W \sec \theta} \right) (1 - e^{-\alpha W \sec \theta}) + In Equation :eq:`cce`, :math:`\alpha` and :math:`z` are measured along the + surface normal. Oblique incidence is therefore handled entirely by the + absorption coefficient: passing the perpendicular-depth coefficient + :math:`\alpha_\perp` from :func:`absorption_effective` (which folds in the + refracted angle, and reduces to :math:`\alpha` at normal incidence) makes + Equation :eq:`cce` valid at any incidence angle without further modification. """ - z0 = absorption * thickness_implant / np.real(cos_incidence) + z0 = absorption * thickness_implant exp_z0 = np.exp(-z0) term_1 = cce_backsurface @@ -413,8 +496,8 @@ def charge_collection_efficiency( def quantum_efficiency_effective( wavelength: u.Quantity | na.AbstractScalar, - direction: None | na.AbstractCartesian3dVectorArray = None, - n: float | na.AbstractScalar = 1, + direction: float | na.AbstractScalar = 1, + n: complex | na.AbstractScalar = 1, thickness_oxide: u.Quantity | na.AbstractScalar = _thickness_oxide, thickness_implant: u.Quantity | na.AbstractScalar = _thickness_implant, thickness_substrate: u.Quantity | na.AbstractScalar = _thickness_substrate, @@ -423,7 +506,6 @@ def quantum_efficiency_effective( chemical_substrate: str | optika.chemicals.AbstractChemical = "Si", roughness_oxide: u.Quantity | na.AbstractScalar = 0 * u.nm, roughness_substrate: u.Quantity | na.AbstractScalar = 0 * u.nm, - normal: None | na.AbstractCartesian3dVectorArray = None, ) -> na.AbstractScalar: r""" Calculate the effective quantum efficiency of a backilluminated detector. @@ -433,8 +515,8 @@ def quantum_efficiency_effective( wavelength The wavelength of the incident light in vacuum. direction - The propagation direction of the incident light in the ambient medium. - If :obj:`None` (default), normal incidence (:math:`\hat{z}`) is assumed. + The cosine of the incidence angle. + Default is normal incidence. n The complex index of refraction of the ambient medium. thickness_oxide @@ -460,9 +542,6 @@ def quantum_efficiency_effective( The RMS roughness the oxide layer surface. roughness_substrate The RMS roughness of the substrate surface. - normal - The vector perpendicular to the surface of the sensor. - If :obj:`None`, then the normal is assumed to be :math:`-\hat{z}` Examples -------- @@ -565,8 +644,6 @@ def quantum_efficiency_effective( and :math:`\text{CCE}(\lambda)` is the charge collection efficiency (computed by :func:`charge_collection_efficiency`). """ - if direction is None: - direction = na.Cartesian3dVectorArray(0, 0, 1) if not isinstance(chemical_oxide, optika.chemicals.AbstractChemical): chemical_oxide = optika.chemicals.Chemical(chemical_oxide) @@ -574,11 +651,6 @@ def quantum_efficiency_effective( if not isinstance(chemical_substrate, optika.chemicals.AbstractChemical): chemical_substrate = optika.chemicals.Chemical(chemical_substrate) - if normal is None: - normal = na.Cartesian3dVectorArray(0, 0, -1) - - direction = -direction @ normal - absorbance_substrate = absorbance( wavelength=wavelength, direction=direction, @@ -590,10 +662,8 @@ def quantum_efficiency_effective( roughness_oxide=roughness_oxide, roughness_substrate=roughness_substrate, ) - + n_substrate = chemical_substrate.n(wavelength) - wavenumber_substrate = np.imag(n_substrate) - absorption_substrate = 4 * np.pi * wavenumber_substrate / wavelength direction_substrate = optika.materials.snells_law_scalar( cos_incidence=direction, @@ -601,11 +671,16 @@ def quantum_efficiency_effective( index_refraction_new=n_substrate, ) + absorption_substrate = absorption_effective( + wavelength=wavelength, + n_substrate=n_substrate, + direction_substrate=direction_substrate, + ) + cce = charge_collection_efficiency( absorption=absorption_substrate, thickness_implant=thickness_implant, cce_backsurface=cce_backsurface, - cos_incidence=direction_substrate, ) result = absorbance_substrate.average * cce @@ -952,12 +1027,16 @@ def electrons_measured_approx( def signal( photons_expected: u.Quantity | na.AbstractScalar, wavelength: u.Quantity | na.ScalarArray, + direction: float | na.AbstractScalar = 1, + n: complex | na.AbstractScalar = 1, + n_substrate: None | complex | na.AbstractScalar = None, transmittance: None | float | na.AbstractScalar = None, - absorption: None | u.Quantity | na.AbstractScalar = None, thickness_implant: u.Quantity | na.AbstractScalar = _thickness_implant, thickness_depletion: u.Quantity | na.AbstractScalar = _thickness_substrate, thickness_substrate: None | na.AbstractScalar = _thickness_substrate, - width_pixel: u.Quantity | na.AbstractScalar = _width_pixel, + width_pixel: ( + u.Quantity | na.AbstractScalar | na.AbstractCartesian2dVectorArray + ) = _width_pixel, cce_backsurface: u.Quantity | na.AbstractScalar = _cce_backsurface, temperature: u.Quantity | na.ScalarArray = 300 * u.K, method: Literal["monte-carlo", "expected"] = "monte-carlo", @@ -979,16 +1058,19 @@ def signal( The `expected` number of photons incident on the detector surface. wavelength The vacuum wavelength of the absorbed photons. + direction + The cosine of the incidence angle. + n + The complex index of refraction of the ambient medium. + n_substrate + The complex index of refraction of the light-sensitive material + If :obj:`None` (the default), the result of + :meth:`optika.chemicals.Chemical.n` for silicon will be used. transmittance The fraction of incident energy transmitted to the light-sensitive layer of the detector. - If :obj:`None` (the default), the result of :func:`absorbance` + If :obj:`None` (the default), the result of :func:`transmittance` called with default values will be used. - absorption - The absorption coefficient of the light-sensitive material for the - wavelength of interest. - If :obj:`None` (the default), the result of - :meth:`optika.chemicals.Chemical.absorption` for silicon will be used. thickness_implant The thickness of the implant layer. Default is the value given in :cite:t:`Stern1994`. @@ -1068,17 +1150,37 @@ def signal( """ if transmittance is None: - transmittance = _transmittance(wavelength).average + transmittance = _transmittance( + wavelength=wavelength, + direction=direction, + n=n, + ).average - if absorption is None: - absorption = optika.chemicals.Chemical("Si").absorption(wavelength) + if n_substrate is None: + n_substrate = optika.chemicals.Chemical("Si").n(wavelength) + + direction_substrate = optika.materials.snells_law_scalar( + cos_incidence=direction, + index_refraction=n, + index_refraction_new=n_substrate, + ) + + absorption = absorption_effective( + wavelength=wavelength, + n_substrate=n_substrate, + direction_substrate=direction_substrate, + ) if method == "expected": iqy = quantum_yield_ideal( wavelength=wavelength, temperature=temperature, ) - absorbance = transmittance * np.exp(-absorption * thickness_substrate) + + transmittance_total = transmittance * np.exp(-absorption * thickness_substrate) + + absorbance = transmittance - transmittance_total + cce = charge_collection_efficiency( absorption=absorption, thickness_implant=thickness_implant, @@ -1114,7 +1216,9 @@ def signal( def vmr_signal( wavelength: u.Quantity | na.ScalarArray, - absorption: None | u.Quantity | na.AbstractScalar = None, + direction: float | na.AbstractScalar = 1, + n: complex | na.AbstractScalar = 1, + n_substrate: None | complex | na.AbstractScalar = None, thickness_implant: u.Quantity | na.AbstractScalar = _thickness_implant, cce_backsurface: u.Quantity | na.AbstractScalar = _cce_backsurface, temperature: u.Quantity | na.ScalarArray = 300 * u.K, @@ -1130,11 +1234,14 @@ def vmr_signal( ---------- wavelength The vacuum wavelength of the absorbed photons. - absorption - The absorption coefficient of the light-sensitive material for the - wavelength of interest. + direction + The cosine of the incidence angle. + n + The complex index of refraction of the ambient medium. + n_substrate + The complex index of refraction of the light-sensitive material. If :obj:`None` (the default), the result of - :meth:`optika.chemicals.Chemical.absorption` for silicon will be used. + :meth:`optika.chemicals.Chemical.n` for silicon will be used. thickness_implant The thickness of the implant layer. Default is the value given in :cite:t:`Stern1994`. @@ -1227,10 +1334,22 @@ def vmr_signal( and :math:`\eta_0` is the CCE at the back surface. """ - if absorption is None: - absorption = optika.chemicals.Chemical("Si").absorption(wavelength) + if n_substrate is None: + n_substrate = optika.chemicals.Chemical("Si").n(wavelength) - n = quantum_yield_ideal( + direction_substrate = optika.materials.snells_law_scalar( + cos_incidence=direction, + index_refraction=n, + index_refraction_new=n_substrate, + ) + + absorption = absorption_effective( + wavelength=wavelength, + n_substrate=n_substrate, + direction_substrate=direction_substrate, + ) + + iqy = quantum_yield_ideal( wavelength=wavelength, temperature=temperature, ) @@ -1246,7 +1365,7 @@ def vmr_signal( result = 0 if shot: - F_shot = n * cce + F_shot = iqy * cce result = result + F_shot @@ -1257,12 +1376,13 @@ def vmr_signal( if pcc: n0 = cce_backsurface - aW = (absorption * thickness_implant).to(u.dimensionless_unscaled).value + aW = absorption * thickness_implant + aW = aW.to(u.dimensionless_unscaled).value F_cce = 2 * np.exp(-aW) * np.square((n0 - 1) / aW) * (np.sinh(aW) - aW) / cce unit = u.electron / u.photon - F_pcc = 1 * unit - cce * unit - F_cce * unit + n * F_cce + F * F_cce + F_pcc = 1 * unit - cce * unit - F_cce * unit + iqy * F_cce + F * F_cce result = result + F_pcc @@ -1282,14 +1402,16 @@ def signal( self, photons: u.Quantity | na.AbstractScalar, wavelength: u.Quantity | na.AbstractScalar, - direction: na.AbstractCartesian3dVectorArray, - normal: na.AbstractCartesian3dVectorArray, + direction: float | na.AbstractScalar = 1, + width_pixel: ( + u.Quantity | na.AbstractScalar | na.AbstractCartesian2dVectorArray + ) = 0 * u.um, axis_xy: None | tuple[str, str] = None, noise: bool = True, ) -> na.AbstractScalar: """ - Given a set of absorbed rays, compute the number of electrons - measured by the sensor using :func:`signal`. + Given the photons incident on each pixel, compute the number of + electrons measured by the sensor using :func:`signal`. Parameters ---------- @@ -1298,9 +1420,9 @@ def signal( wavelength An assumed grid of wavelengths for the incident photons. direction - An assumed propagation direction for the incident photons. - normal - The vector perpendicular to the surface of the sensor. + The cosine of the incidence angle. + width_pixel + The physical size of each pixel, used by the charge-diffusion model. axis_xy The two logical axes corresponding to the pixel grid of the sensor. If provided, charge diffusion will occur along these two axes. @@ -1349,21 +1471,23 @@ def signal( self, photons: u.Quantity | na.AbstractScalar, wavelength: u.Quantity | na.AbstractScalar, - direction: na.AbstractCartesian3dVectorArray, - normal: na.AbstractCartesian3dVectorArray, + direction: float | na.AbstractScalar = 1, + width_pixel: ( + u.Quantity | na.AbstractScalar | na.AbstractCartesian2dVectorArray + ) = 0 * u.um, axis_xy: None | tuple[str, str] = None, noise: bool = True, - ) -> optika.rays.RayVectorArray: + ) -> na.AbstractScalar: if not photons.unit.is_equivalent(u.photon): h = astropy.constants.h c = astropy.constants.c - intensity = photons / (h * c / wavelength) * u.photon + photons = photons / (h * c / wavelength) * u.photon if noise: - intensity = na.random.poisson(intensity.to(u.ph)).astype(int) + photons = na.random.poisson(photons.to(u.ph)).astype(int) - electrons = intensity * u.electron / u.photon + electrons = photons * u.electron / u.photon electrons = electrons.to(u.electron) return electrons @@ -1527,10 +1651,10 @@ def transmittance( self, wavelength: u.Quantity | na.AbstractScalar, direction: None | na.AbstractCartesian3dVectorArray = None, - index_refraction: float | na.AbstractScalar = 1, + n: complex | na.AbstractScalar = 1, normal: None | na.AbstractCartesian3dVectorArray = None, ) -> optika.vectors.PolarizationVectorArray: - """ + r""" Compute the fraction of energy transmitted to the light-sensitive region of the sensor. @@ -1541,7 +1665,7 @@ def transmittance( direction The propagation direction of the incident light in the ambient medium. If :obj:`None` (default), normal incidence (:math:`\hat{z}`) is assumed. - index_refraction + n The complex index of refraction of the ambient medium. normal The vector perpendicular to the surface of the CCD sensor. @@ -1550,12 +1674,12 @@ def transmittance( direction = na.Cartesian3dVectorArray(0, 0, 1) if normal is None: - normal = na.Cartesian3dVectorArray(0, 0, 1) + normal = na.Cartesian3dVectorArray(0, 0, -1) return transmittance( wavelength=wavelength, direction=-direction @ normal, - n=index_refraction, + n=n, thickness_oxide=self.thickness_oxide, thickness_substrate=self.thickness_substrate, chemical_oxide=self._chemical_oxide, @@ -1568,10 +1692,10 @@ def absorbance( self, wavelength: u.Quantity | na.AbstractScalar, direction: None | na.AbstractCartesian3dVectorArray = None, - index_refraction: float | na.AbstractScalar = 1, + n: complex | na.AbstractScalar = 1, normal: None | na.AbstractCartesian3dVectorArray = None, ) -> optika.vectors.PolarizationVectorArray: - """ + r""" Compute the fraction of energy absorbed by the light-sensitive region of the sensor. @@ -1591,12 +1715,12 @@ def absorbance( direction = na.Cartesian3dVectorArray(0, 0, 1) if normal is None: - normal = na.Cartesian3dVectorArray(0, 0, 1) + normal = na.Cartesian3dVectorArray(0, 0, -1) return absorbance( wavelength=wavelength, direction=-direction @ normal, - n=index_refraction, + n=n, thickness_oxide=self.thickness_oxide, thickness_substrate=self.thickness_substrate, chemical_oxide=self._chemical_oxide, @@ -1607,49 +1731,88 @@ def absorbance( def charge_collection_efficiency( self, - rays: optika.rays.AbstractRayVectorArray, - normal: na.AbstractCartesian3dVectorArray, + wavelength: u.Quantity | na.AbstractScalar, + direction: None | na.AbstractCartesian3dVectorArray = None, + n: complex | na.AbstractScalar = 1, + normal: None | na.AbstractCartesian3dVectorArray = None, ) -> na.AbstractScalar: - """ + r""" Compute the charge collection efficiency of this CCD sensor material using :func:`charge_collection_efficiency`. Parameters ---------- - rays - The light rays incident on the CCD surface. + wavelength + The wavelength of the incident light in vacuum. + direction + The propagation direction of the incident light in the ambient medium. + If :obj:`None` (default), normal incidence (:math:`\hat{z}`) is assumed. + n + The complex index of refraction of the ambient medium. normal The vector perpendicular to the surface of the CCD sensor. """ + + if direction is None: + direction = na.Cartesian3dVectorArray(0, 0, 1) + + if normal is None: + normal = na.Cartesian3dVectorArray(0, 0, -1) + + direction = -direction @ normal + + n_substrate = self._chemical.n(wavelength) + + direction_substrate = optika.materials.snells_law_scalar( + cos_incidence=direction, + index_refraction=n, + index_refraction_new=n_substrate, + ) + return charge_collection_efficiency( - absorption=self._chemical.absorption(rays.wavelength), + absorption=absorption_effective( + wavelength=wavelength, + n_substrate=n_substrate, + direction_substrate=direction_substrate, + ), thickness_implant=self.thickness_implant, cce_backsurface=self.cce_backsurface, - cos_incidence=-rays.direction @ normal, ) def quantum_efficiency_effective( self, - rays: optika.rays.AbstractRayVectorArray, - normal: na.AbstractCartesian3dVectorArray, + wavelength: u.Quantity | na.AbstractScalar, + direction: None | na.AbstractCartesian3dVectorArray = None, + n: complex | na.AbstractScalar = 1, + normal: None | na.AbstractCartesian3dVectorArray = None, ) -> na.AbstractScalar: - """ + r""" Compute the effective quantum efficiency of this CCD material using :func:`quantum_efficiency_effective`. Parameters ---------- - rays - The light rays incident on the CCD surface + wavelength + The wavelength of the incident light in vacuum. + direction + The propagation direction of the incident light in the ambient medium. + If :obj:`None` (default), normal incidence (:math:`\hat{z}`) is assumed. + n + The complex index of refraction of the ambient medium. normal The vector perpendicular to the surface of the CCD. """ - k = rays.attenuation * rays.wavelength / (4 * np.pi) - n = rays.index_refraction + k * 1j + if direction is None: + direction = na.Cartesian3dVectorArray(0, 0, 1) + + if normal is None: + normal = na.Cartesian3dVectorArray(0, 0, -1) + + direction = -direction @ normal return quantum_efficiency_effective( - wavelength=rays.wavelength, - direction=rays.direction, + wavelength=wavelength, + direction=direction, n=n, thickness_oxide=self.thickness_oxide, thickness_implant=self.thickness_implant, @@ -1659,117 +1822,168 @@ def quantum_efficiency_effective( chemical_substrate=self._chemical, roughness_oxide=self.roughness_oxide, roughness_substrate=self.roughness_substrate, - normal=normal, ) def quantum_efficiency( self, - rays: optika.rays.AbstractRayVectorArray, - normal: na.AbstractCartesian3dVectorArray, + wavelength: u.Quantity | na.AbstractScalar, + direction: None | na.AbstractCartesian3dVectorArray = None, + n: complex | na.AbstractScalar = 1, + normal: None | na.AbstractCartesian3dVectorArray = None, ) -> na.AbstractScalar: - """ + r""" Compute the quantum efficiency of this CCD material using :meth:`quantum_efficiency_effective` and :meth:`quantum_yield_ideal`. Parameters ---------- - rays - The light rays incident on the CCD surface + wavelength + The wavelength of the incident light in vacuum. + direction + The propagation direction of the incident light in the ambient medium. + If :obj:`None` (default), normal incidence (:math:`\hat{z}`) is assumed. + n + The complex index of refraction of the ambient medium. normal The vector perpendicular to the surface of the CCD. """ - iqy = self.quantum_yield_ideal(rays.wavelength) - eqe = self.quantum_efficiency_effective(rays, normal) + iqy = self.quantum_yield_ideal(wavelength) + eqe = self.quantum_efficiency_effective( + wavelength=wavelength, + direction=direction, + n=n, + normal=normal, + ) return iqy * eqe def probability_measurement( self, - rays: optika.rays.AbstractRayVectorArray, - normal: na.AbstractCartesian3dVectorArray, + wavelength: u.Quantity | na.AbstractScalar, + direction: None | na.AbstractCartesian3dVectorArray = None, + normal: None | na.AbstractCartesian3dVectorArray = None, ) -> na.AbstractScalar: - """ + r""" Compute the probability of measuring an absorbed photon for this sensor using :func:`probability_measurement`. Parameters ---------- - rays - The light rays incident on the CCD surface + wavelength + The wavelength of the incident light in vacuum. + direction + The propagation direction of the incident light in the ambient medium. + If :obj:`None` (default), normal incidence (:math:`\hat{z}`) is assumed. normal The vector perpendicular to the surface of the CCD. """ return probability_measurement( - iqy=self.quantum_yield_ideal(rays.wavelength), - cce=self.charge_collection_efficiency(rays, normal), + iqy=self.quantum_yield_ideal(wavelength), + cce=self.charge_collection_efficiency( + wavelength=wavelength, + direction=direction, + normal=normal, + ), ) def electrons_measured( self, - rays: optika.rays.RayVectorArray, - normal: na.AbstractCartesian3dVectorArray, - ) -> optika.rays.RayVectorArray: + photons_transmitted: na.AbstractScalar, + wavelength: u.Quantity | na.AbstractScalar, + direction: None | na.AbstractCartesian3dVectorArray = None, + n: complex | na.AbstractScalar = 1, + normal: None | na.AbstractCartesian3dVectorArray = None, + width_pixel: ( + u.Quantity | na.AbstractScalar | na.AbstractCartesian2dVectorArray + ) = 0 * u.um, + axis_xy: None | tuple[str, str] = None, + ) -> na.AbstractScalar: """ Randomly sample the number of measured electrons given the number of absorbed photons using :func:`electrons_measured`. """ - intensity = rays.intensity - wavelength = rays.wavelength + if direction is None: + direction = na.Cartesian3dVectorArray(0, 0, 1) + + if normal is None: + normal = na.Cartesian3dVectorArray(0, 0, -1) + + direction = -direction @ normal - electrons = electrons_measured( - photons_absorbed=intensity, + n_substrate = self._chemical.n(wavelength) + + direction_substrate = optika.materials.snells_law_scalar( + cos_incidence=direction, + index_refraction=n, + index_refraction_new=n_substrate, + ) + + return electrons_measured( + photons_transmitted=photons_transmitted, wavelength=wavelength, - absorption=self._chemical.absorption(wavelength), + absorption=absorption_effective( + wavelength=wavelength, + n_substrate=n_substrate, + direction_substrate=direction_substrate, + ), thickness_implant=self.thickness_implant, + thickness_depletion=self.depletion.thickness, + thickness_substrate=self.thickness_substrate, + width_pixel=width_pixel, cce_backsurface=self.cce_backsurface, + axis_xy=axis_xy, temperature=self.temperature, ) - result = dataclasses.replace(rays, intensity=electrons) - - return result - def signal( self, - rays: optika.rays.RayVectorArray, - normal: na.AbstractCartesian3dVectorArray, + photons: u.Quantity | na.AbstractScalar, + wavelength: u.Quantity | na.AbstractScalar, + direction: float | na.AbstractScalar = 1, + n: complex | na.AbstractScalar = 1, + width_pixel: ( + u.Quantity | na.AbstractScalar | na.AbstractCartesian2dVectorArray + ) = 0 * u.um, + axis_xy: None | tuple[str, str] = None, noise: bool = True, - ) -> optika.rays.RayVectorArray: - intensity = rays.intensity - wavelength = rays.wavelength + ) -> na.AbstractScalar: - if not intensity.unit.is_equivalent(u.photon): + if not photons.unit.is_equivalent(u.photon): h = astropy.constants.h c = astropy.constants.c - intensity = intensity / (h * c / wavelength) * u.photon + photons = photons / (h * c / wavelength) * u.photon if noise: - method = "exact" + method = "monte-carlo" else: method = "expected" - electrons = signal( - photons_expected=intensity, + return signal( + photons_expected=photons, wavelength=wavelength, - absorbance=1, - absorption=self._chemical.absorption(wavelength), + direction=direction, + n=n, + n_substrate=self._chemical.n(wavelength), + transmittance=1, thickness_implant=self.thickness_implant, + thickness_depletion=self.depletion.thickness, + thickness_substrate=self.thickness_substrate, + width_pixel=width_pixel, cce_backsurface=self.cce_backsurface, temperature=self.temperature, method=method, + axis_xy=axis_xy, ) - result = dataclasses.replace(rays, intensity=electrons) - - return result def photons_incident( self, electrons: u.Quantity | na.AbstractScalar, wavelength: u.Quantity | na.AbstractScalar, - direction: na.AbstractCartesian3dVectorArray, - normal: na.AbstractCartesian3dVectorArray, + direction: None | na.AbstractCartesian3dVectorArray = None, + n: complex | na.AbstractScalar = 1, + normal: None | na.AbstractCartesian3dVectorArray = None, ) -> na.AbstractScalar: """ Compute the expected number of incident photons for a given number @@ -1783,47 +1997,30 @@ def photons_incident( The assumed wavelength of the incident photons. direction The assumed direction of the incident photons. + n + The complex index of refraction of the ambient medium. normal The vector perpendicular to the surface of the sensor. """ qe = self.quantum_efficiency( - rays=optika.rays.RayVectorArray( - wavelength=wavelength, - direction=direction, - ), + wavelength=wavelength, + direction=direction, + n=n, normal=normal, ) return electrons / qe - def charge_diffusion( - self, - rays: optika.rays.RayVectorArray, - normal: na.AbstractCartesian3dVectorArray, - ) -> optika.rays.RayVectorArray: - width = self.width_charge_diffusion(rays, normal) - - position = dataclasses.replace( - rays.position, - x=na.random.normal(rays.position.x, width), - y=na.random.normal(rays.position.y, width), - ) - - rays = dataclasses.replace( - rays, - position=position, - ) - - return rays - def efficiency( self, rays: optika.rays.AbstractRayVectorArray, normal: na.AbstractCartesian3dVectorArray, ) -> na.ScalarLike: - result = self.absorbance( - rays=rays, + result = self.transmittance( + wavelength=rays.wavelength, + direction=rays.direction, + n=rays.n, normal=normal, ) diff --git a/optika/sensors/materials/_materials_test.py b/optika/sensors/materials/_materials_test.py index 1f662c2c..ba4154ff 100644 --- a/optika/sensors/materials/_materials_test.py +++ b/optika/sensors/materials/_materials_test.py @@ -127,23 +127,15 @@ def test_absorbance( 1, ], ) -@pytest.mark.parametrize( - argnames="cos_incidence", - argvalues=[ - 1, - ], -) def test_charge_collection_efficiency( absorption: u.Quantity | na.AbstractScalar, thickness_implant: u.Quantity | na.AbstractScalar, cce_backsurface: u.Quantity | na.AbstractScalar, - cos_incidence: float | na.AbstractScalar, ): result = optika.sensors.charge_collection_efficiency( absorption=absorption, thickness_implant=thickness_implant, cce_backsurface=cce_backsurface, - cos_incidence=cos_incidence, ) assert np.all(result >= 0) @@ -160,8 +152,8 @@ def test_charge_collection_efficiency( @pytest.mark.parametrize( argnames="direction", argvalues=[ - None, - na.Cartesian3dVectorArray(0, 0, 1), + 1, + 0.5, ], ) @pytest.mark.parametrize( @@ -191,7 +183,7 @@ def test_charge_collection_efficiency( ) def test_quantum_efficiency_effective( wavelength: u.Quantity | na.AbstractScalar, - direction: na.AbstractCartesian3dVectorArray, + direction: float | na.AbstractScalar, thickness_oxide: u.Quantity | na.AbstractScalar, thickness_implant: u.Quantity | na.AbstractScalar, thickness_substrate: u.Quantity | na.AbstractScalar, @@ -299,8 +291,8 @@ def test_electrons_measured_approx( @pytest.mark.parametrize( argnames="method", argvalues=[ - "exact", - "approx", + "monte-carlo", + "expected", ], ) def test_signal( @@ -336,19 +328,22 @@ class AbstractTestAbstractSensorMaterial( AbstractTestAbstractMaterial, ): @pytest.mark.parametrize( - argnames="rays", + argnames="photons", argvalues=[ - optika.rays.RayVectorArray( - intensity=1e-6 * u.erg, - wavelength=100 * u.AA, - direction=na.Cartesian3dVectorArray(0, 0, 1), - ), + 1e-6 * u.erg, ], ) @pytest.mark.parametrize( - argnames="normal", + argnames="wavelength", argvalues=[ - na.Cartesian3dVectorArray(0, 0, -1), + 100 * u.AA, + ], + ) + @pytest.mark.parametrize( + argnames="direction", + argvalues=[ + 1, + 0.5, ], ) @pytest.mark.parametrize( @@ -358,13 +353,19 @@ class AbstractTestAbstractSensorMaterial( def test_signal( self, a: optika.sensors.materials.AbstractSensorMaterial, - rays: optika.rays.RayVectorArray, - normal: na.AbstractCartesian3dVectorArray, + photons: u.Quantity | na.AbstractScalar, + wavelength: u.Quantity | na.AbstractScalar, + direction: float | na.AbstractScalar, noise: bool, ): - result = a.signal(rays, normal, noise=noise) - assert isinstance(result, optika.rays.RayVectorArray) - assert np.all(result.intensity >= 0 * u.electron) + result = a.signal( + photons=photons, + wavelength=wavelength, + direction=direction, + noise=noise, + ) + assert isinstance(na.as_named_array(result), na.AbstractScalar) + assert np.all(result >= 0 * u.electron) @pytest.mark.parametrize( argnames="electrons", @@ -407,32 +408,6 @@ def test_photons_incident( assert isinstance(na.as_named_array(result), na.AbstractScalar) assert result.unit.is_equivalent(u.photon) - @pytest.mark.parametrize( - argnames="rays", - argvalues=[ - optika.rays.RayVectorArray( - intensity=10 * u.electron, - wavelength=100 * u.AA, - position=na.Cartesian3dVectorArray() * u.mm, - direction=na.Cartesian3dVectorArray(0, 0, 1), - ), - ], - ) - @pytest.mark.parametrize( - argnames="normal", - argvalues=[ - na.Cartesian3dVectorArray(0, 0, -1), - ], - ) - def test_charge_diffusion( - self, - a: optika.sensors.materials.AbstractSensorMaterial, - rays: optika.rays.RayVectorArray, - normal: na.AbstractCartesian3dVectorArray, - ): - result = a.charge_diffusion(rays, normal) - assert isinstance(result, optika.rays.RayVectorArray) - @pytest.mark.parametrize( argnames="a", @@ -522,154 +497,190 @@ def test_depletion( ) @pytest.mark.parametrize( - argnames="rays", - argvalues=[ - optika.rays.RayVectorArray( - wavelength=100 * u.AA, - position=na.Cartesian3dVectorArray() * u.mm, - direction=na.Cartesian3dVectorArray(0, 0, 1), - ), - ], - ) - @pytest.mark.parametrize( - argnames="normal", + argnames="wavelength", argvalues=[ - na.Cartesian3dVectorArray(0, 0, -1), + 100 * u.AA, ], ) def test_width_charge_diffusion( self, a: optika.sensors.materials.AbstractBackIlluminatedSiliconSensorMaterial, - rays: optika.rays.AbstractRayVectorArray, - normal: na.AbstractCartesian3dVectorArray, + wavelength: u.Quantity | na.AbstractScalar, ): - result = a.width_charge_diffusion(rays, normal) + result = a.width_charge_diffusion(wavelength=wavelength) assert np.all(result >= 0 * u.um) @pytest.mark.parametrize( - argnames="rays", + argnames="wavelength", + argvalues=[ + 100 * u.AA, + ], + ) + @pytest.mark.parametrize( + argnames="direction", argvalues=[ - optika.rays.RayVectorArray( - wavelength=100 * u.AA, - direction=na.Cartesian3dVectorArray(0, 0, 1), - ), + None, ], ) @pytest.mark.parametrize( argnames="normal", argvalues=[ - na.Cartesian3dVectorArray(0, 0, -1), + None, ], ) def test_absorbance( self, a: optika.sensors.materials.AbstractBackIlluminatedSiliconSensorMaterial, - rays: optika.rays.AbstractRayVectorArray, - normal: na.AbstractCartesian3dVectorArray, + wavelength: u.Quantity | na.AbstractScalar, + direction: None | na.AbstractCartesian3dVectorArray, + normal: None | na.AbstractCartesian3dVectorArray, ): - result = a.absorbance(rays, normal) + result = a.absorbance( + wavelength=wavelength, + direction=direction, + normal=normal, + ) assert np.all(result >= 0) assert np.all(result <= 1) @pytest.mark.parametrize( - argnames="rays", + argnames="wavelength", + argvalues=[ + 100 * u.AA, + ], + ) + @pytest.mark.parametrize( + argnames="direction", argvalues=[ - optika.rays.RayVectorArray( - wavelength=100 * u.AA, - direction=na.Cartesian3dVectorArray(0, 0, 1), - ), + None, ], ) @pytest.mark.parametrize( argnames="normal", argvalues=[ - na.Cartesian3dVectorArray(0, 0, -1), + None, ], ) def test_charge_collection_efficiency( self, a: optika.sensors.materials.AbstractBackIlluminatedSiliconSensorMaterial, - rays: optika.rays.AbstractRayVectorArray, - normal: na.AbstractCartesian3dVectorArray, + wavelength: u.Quantity | na.AbstractScalar, + direction: None | na.AbstractCartesian3dVectorArray, + normal: None | na.AbstractCartesian3dVectorArray, ): - result = a.charge_collection_efficiency(rays, normal) + result = a.charge_collection_efficiency( + wavelength=wavelength, + direction=direction, + normal=normal, + ) assert np.all(result >= 0) assert np.all(result <= 1) @pytest.mark.parametrize( - argnames="rays", + argnames="wavelength", argvalues=[ - optika.rays.RayVectorArray( - wavelength=100 * u.AA, - direction=na.Cartesian3dVectorArray(0, 0, 1), - ), + 100 * u.AA, + ], + ) + @pytest.mark.parametrize( + argnames="direction", + argvalues=[ + None, ], ) @pytest.mark.parametrize( argnames="normal", argvalues=[ - na.Cartesian3dVectorArray(0, 0, -1), + None, ], ) def test_quantum_efficiency( self, a: optika.sensors.materials.AbstractBackIlluminatedSiliconSensorMaterial, - rays: optika.rays.AbstractRayVectorArray, - normal: na.AbstractCartesian3dVectorArray, + wavelength: u.Quantity | na.AbstractScalar, + direction: None | na.AbstractCartesian3dVectorArray, + normal: None | na.AbstractCartesian3dVectorArray, ): - result = a.quantum_efficiency(rays, normal) + result = a.quantum_efficiency( + wavelength=wavelength, + direction=direction, + normal=normal, + ) assert result > 0 * u.electron / u.photon @pytest.mark.parametrize( - argnames="rays", + argnames="wavelength", argvalues=[ - optika.rays.RayVectorArray( - wavelength=100 * u.AA, - direction=na.Cartesian3dVectorArray(0, 0, 1), - ), + 100 * u.AA, + ], + ) + @pytest.mark.parametrize( + argnames="direction", + argvalues=[ + None, ], ) @pytest.mark.parametrize( argnames="normal", argvalues=[ - na.Cartesian3dVectorArray(0, 0, -1), + None, ], ) def test_probability_measurement( self, a: optika.sensors.materials.AbstractBackIlluminatedSiliconSensorMaterial, - rays: optika.rays.AbstractRayVectorArray, - normal: na.AbstractCartesian3dVectorArray, + wavelength: u.Quantity | na.AbstractScalar, + direction: None | na.AbstractCartesian3dVectorArray, + normal: None | na.AbstractCartesian3dVectorArray, ): - result = a.probability_measurement(rays, normal) + result = a.probability_measurement( + wavelength=wavelength, + direction=direction, + normal=normal, + ) assert np.all(result >= 0) assert np.all(result <= 1) @pytest.mark.parametrize( - argnames="rays", + argnames="photons_transmitted", + argvalues=[ + (100 * u.photon).astype(int), + ], + ) + @pytest.mark.parametrize( + argnames="wavelength", + argvalues=[ + 100 * u.AA, + ], + ) + @pytest.mark.parametrize( + argnames="direction", argvalues=[ - optika.rays.RayVectorArray( - intensity=100 * u.photon, - wavelength=100 * u.AA, - direction=na.Cartesian3dVectorArray(0, 0, 1), - ), + None, ], ) @pytest.mark.parametrize( argnames="normal", argvalues=[ - na.Cartesian3dVectorArray(0, 0, -1), + None, ], ) def test_electrons_measured( self, a: optika.sensors.materials.AbstractBackIlluminatedSiliconSensorMaterial, - rays: optika.rays.RayVectorArray, - normal: na.AbstractCartesian3dVectorArray, + photons_transmitted: u.Quantity | na.AbstractScalar, + wavelength: u.Quantity | na.AbstractScalar, + direction: None | na.AbstractCartesian3dVectorArray, + normal: None | na.AbstractCartesian3dVectorArray, ): - result = a.electrons_measured(rays, normal) - assert isinstance(result, optika.rays.RayVectorArray) - assert np.all(result.intensity >= 0 * u.electron) + result = a.electrons_measured( + photons_transmitted=photons_transmitted, + wavelength=wavelength, + direction=direction, + normal=normal, + ) + assert isinstance(na.as_named_array(result), na.AbstractScalar) + assert np.all(result >= 0 * u.electron) @pytest.mark.parametrize( diff --git a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py index 795cda98..c8114574 100644 --- a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py +++ b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py @@ -472,7 +472,9 @@ def electrons_measured( thickness_implant: u.Quantity | na.AbstractScalar = _thickness_implant, thickness_depletion: None | u.Quantity | na.AbstractScalar = None, thickness_substrate: u.Quantity | na.AbstractScalar = _thickness_substrate, - width_pixel: u.Quantity | na.ScalarArray = _width_pixel, + width_pixel: ( + u.Quantity | na.AbstractScalar | na.AbstractCartesian2dVectorArray + ) = _width_pixel, cce_backsurface: u.Quantity | na.AbstractScalar = _cce_backsurface, temperature: u.Quantity | na.ScalarArray = 300 * u.K, axis_xy: None | tuple[str, str] = None, @@ -494,7 +496,10 @@ def electrons_measured( wavelength The vacuum wavelength of the absorbed photons. absorption - The absorption coefficient of silicon at the given wavelength. + The absorption coefficient of silicon per unit perpendicular depth. + For oblique incidence, supply the effective coefficient from + :func:`optika.sensors.absorption_effective`, which folds in the + refracted angle, so no separate angle argument is needed. thickness_implant The thickness of the implant layer, where partial-charge collection occurs. thickness_depletion @@ -506,6 +511,10 @@ def electrons_measured( The thickness of the entire light-sensitive region of the device. width_pixel The size of a single pixel on the sensor. + A scalar gives square pixels; a + :class:`named_arrays.AbstractCartesian2dVectorArray` whose ``x``/``y`` + components are the pixel widths along ``axis_xy[0]``/``axis_xy[1]`` + gives rectangular pixels. cce_backsurface The differential charge collection efficiency on the back surface of the sensor. @@ -587,6 +596,12 @@ def electrons_measured( if shape_random is None: shape_random = dict() + if not isinstance(width_pixel, na.AbstractCartesian2dVectorArray): + width_pixel = na.Cartesian2dVectorArray(width_pixel, width_pixel) + + width_pixel_x = width_pixel.x + width_pixel_y = width_pixel.y + shape = na.broadcast_shapes( na.shape(photons_transmitted), na.shape(wavelength), @@ -594,7 +609,8 @@ def electrons_measured( na.shape(thickness_implant), na.shape(thickness_depletion), na.shape(thickness_substrate), - na.shape(width_pixel), + na.shape(width_pixel_x), + na.shape(width_pixel_y), na.shape(cce_backsurface), na.shape(temperature), shape_random, @@ -615,7 +631,8 @@ def electrons_measured( thickness_implant = na.broadcast_to(thickness_implant, shape) thickness_depletion = na.broadcast_to(thickness_depletion, shape) thickness_substrate = na.broadcast_to(thickness_substrate, shape) - width_pixel = na.broadcast_to(width_pixel, shape) + width_pixel_x = na.broadcast_to(width_pixel_x, shape) + width_pixel_y = na.broadcast_to(width_pixel_y, shape) cce_backsurface = na.broadcast_to(cce_backsurface, shape) if not isinstance(cce_backsurface, u.Quantity): @@ -645,7 +662,8 @@ def electrons_measured( thickness_implant=thickness_implant.ndarray, thickness_depletion=thickness_depletion.ndarray, thickness_substrate=thickness_substrate.ndarray, - width_pixel=width_pixel, + width_pixel_x=width_pixel_x.ndarray, + width_pixel_y=width_pixel_y.ndarray, cce_backsurface=cce_backsurface.ndarray, p_n=p_n.ndarray, n=n.ndarray, @@ -653,11 +671,16 @@ def electrons_measured( fano_inf=fano_inf.ndarray, ) - return na.ScalarArray( + result = na.ScalarArray( ndarray=result, axes=tuple(shape), ) + if axis_xy is None: + result = result[{axis_x: 0, axis_y: 0}] + + return result + def _electrons_measured_quantity( photons_transmitted: u.Quantity, @@ -666,7 +689,8 @@ def _electrons_measured_quantity( thickness_implant: u.Quantity, thickness_depletion: u.Quantity, thickness_substrate: u.Quantity, - width_pixel: u.Quantity, + width_pixel_x: u.Quantity, + width_pixel_y: u.Quantity, cce_backsurface: u.Quantity, p_n: np.ndarray, n: np.ndarray, @@ -681,7 +705,8 @@ def _electrons_measured_quantity( thickness_depletion.shape, thickness_substrate.shape, cce_backsurface.shape, - width_pixel.shape, + width_pixel_x.shape, + width_pixel_y.shape, ) num_x = shape[~1] @@ -695,7 +720,8 @@ def _electrons_measured_quantity( thickness_implant = thickness_implant.to_value(unit_length) thickness_depletion = thickness_depletion.to_value(unit_length) thickness_substrate = thickness_substrate.to_value(unit_length) - width_pixel = width_pixel.to_value(unit_length) + width_pixel_x = width_pixel_x.to_value(unit_length) + width_pixel_y = width_pixel_y.to_value(unit_length) cce_backsurface = cce_backsurface.to_value(u.dimensionless_unscaled) energy_pair_inf = energy_pair_inf.to_value(u.eV) fano_inf = fano_inf.to_value(u.electron / u.photon) @@ -707,7 +733,8 @@ def _electrons_measured_quantity( thickness_implant=thickness_implant.reshape(-1, num_x, num_y), thickness_depletion=thickness_depletion.reshape(-1, num_x, num_y), thickness_substrate=thickness_substrate.reshape(-1, num_x, num_y), - width_pixel=width_pixel.reshape(-1, num_x, num_y), + width_pixel_x=width_pixel_x.reshape(-1, num_x, num_y), + width_pixel_y=width_pixel_y.reshape(-1, num_x, num_y), cce_backsurface=cce_backsurface.reshape(-1, num_x, num_y), p_n=p_n.reshape(-1, num_x, num_y, p_n.shape[~0]), n=n.reshape(-1, num_x, num_y, n.shape[~0]), @@ -722,7 +749,11 @@ def _electrons_measured_quantity( return result -@numba.njit(fastmath=True, parallel=True) +@numba.njit( + cache=True, + fastmath=True, + parallel=True, +) def _electrons_measured_numba( photons_transmitted: np.ndarray, energy: np.ndarray, @@ -730,7 +761,8 @@ def _electrons_measured_numba( thickness_implant: np.ndarray, thickness_depletion: np.ndarray, thickness_substrate: np.ndarray, - width_pixel: np.ndarray, + width_pixel_x: np.ndarray, + width_pixel_y: np.ndarray, cce_backsurface: np.ndarray, p_n: np.ndarray, n: np.ndarray, @@ -756,6 +788,8 @@ def _electrons_measured_numba( fano_inf_i = fano_inf[i, x, y] z_substrate = thickness_substrate[i, x, y] z_ff = z_substrate - thickness_depletion[i, x, y] + wp_x = width_pixel_x[i, x, y] + wp_y = width_pixel_y[i, x, y] d = 1 / a @@ -795,11 +829,11 @@ def _electrons_measured_numba( v = random.uniform(-0.5, 0.5) for e in range(m_ij): - if z_ij < z_ff: - w = z_ff * math.sqrt(1 - z_ij / z_ff) / width_pixel + if z_ij < z_ff and wp_x > 0 and wp_y > 0: + w = z_ff * math.sqrt(1 - z_ij / z_ff) - p = random.gauss(u, w) - q = random.gauss(v, w) + p = random.gauss(u, w / wp_x) + q = random.gauss(v, w / wp_y) p = round(p) q = round(q) diff --git a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020_test.py b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020_test.py index 6ef10eb3..59c304f3 100644 --- a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020_test.py +++ b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020_test.py @@ -98,7 +98,7 @@ def test_fano_factor_inf( @pytest.mark.parametrize( - argnames="photons_absorbed", + argnames="photons_transmitted", argvalues=[ 100 * u.photon, ], @@ -125,20 +125,28 @@ def test_fano_factor_inf( 0.5, ], ) +@pytest.mark.parametrize( + argnames="width_pixel", + argvalues=[ + 15 * u.um, + ], +) @pytest.mark.parametrize("temperature", _temperture) def test_electrons_measured( - photons_absorbed: u.Quantity | na.AbstractScalar, + photons_transmitted: u.Quantity | na.AbstractScalar, wavelength: u.Quantity | na.ScalarArray, absorption: u.Quantity | na.AbstractScalar, thickness_implant: u.Quantity | na.AbstractScalar, cce_backsurface: u.Quantity | na.AbstractScalar, temperature: u.Quantity | na.ScalarArray, + width_pixel: u.Quantity | na.AbstractCartesian2dVectorArray, ): result = _ramanathan_2020.electrons_measured( - photons_absorbed=photons_absorbed, + photons_transmitted=photons_transmitted, wavelength=wavelength, absorption=absorption, thickness_implant=thickness_implant, + width_pixel=width_pixel, cce_backsurface=cce_backsurface, temperature=temperature, ) @@ -146,7 +154,7 @@ def test_electrons_measured( assert np.all(result >= 0 * u.electron) shape = na.shape_broadcasted( - photons_absorbed, + photons_transmitted, wavelength, absorption, thickness_implant, From 7057014479a9706c63fd99c635d8abe2b0001e0c Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Sun, 7 Jun 2026 10:39:45 -0600 Subject: [PATCH 04/19] Wire binned detector physics and fix internal consistency of the sensor refactor. Sensor pipeline: - Rename the sensor methods to `collect` -> `expose` -> `measure` (was `bin_rays` -> `expose` -> `readout`), update the `_sequential.py` caller. - `expose` centers the wavelength bin edges and is shared by any per-pixel photon image; tidy the empty-pixel mean-cosine division. Review fixes (internal consistency of the breaking changes): - `fit_eqe` now passes the cosine `direction=1` instead of a 3d vector to the (now cosine-based) standalone `quantum_efficiency_effective`. This was masked by the joblib cache, so a cold cache would have broken all three device materials; verified by running the fit uncached. - Update the stale docstring examples that still used the old API (`quantum_efficiency_effective(rays=...)`, `width_charge_diffusion(rays=...)`, and a vector `direction` to the standalone QE); confirmed by a full doc build. - Revert the `RayVectorArray.intensity` default back to `1` and fix a cosmetic typo in the new `RayVectorArray.n` property. Add a `test_n` exercising the `RayVectorArray.n` complex-index property. Co-Authored-By: Claude Opus 4.8 --- optika/rays/_ray_vectors.py | 4 +++ optika/rays/_tests/test_ray_vectors.py | 15 ++++++++ optika/sensors/_sensors.py | 36 +++++++++++++------ optika/sensors/_sensors_test.py | 4 +-- .../materials/_e2v_ccd203/_e2v_ccd203.py | 12 ++----- .../materials/_e2v_ccd97/_e2v_ccd97.py | 18 +++------- optika/sensors/materials/_materials.py | 10 ++---- .../_tektronix_tk512cb/_tektronix_tk512cb.py | 12 ++----- optika/surfaces.py | 3 +- optika/systems/_sequential.py | 2 +- 10 files changed, 63 insertions(+), 53 deletions(-) diff --git a/optika/rays/_ray_vectors.py b/optika/rays/_ray_vectors.py index cc3176b2..a8f4a9d4 100644 --- a/optika/rays/_ray_vectors.py +++ b/optika/rays/_ray_vectors.py @@ -70,6 +70,10 @@ def unvignetted(self) -> na.AbstractScalar: and :obj:`False` indicates the ray has been blocked by an aperture. """ + @property + def n(self) -> complex | na.AbstractScalar: + return self.index_refraction + self.attenuation * self.wavelength / (4 * np.pi) * 1j + @property def type_abstract(self) -> type[na.AbstractArray]: return AbstractRayVectorArray diff --git a/optika/rays/_tests/test_ray_vectors.py b/optika/rays/_tests/test_ray_vectors.py index c5aa626a..432ef85a 100644 --- a/optika/rays/_tests/test_ray_vectors.py +++ b/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) diff --git a/optika/sensors/_sensors.py b/optika/sensors/_sensors.py index 8d266e28..dbedc940 100644 --- a/optika/sensors/_sensors.py +++ b/optika/sensors/_sensors.py @@ -81,7 +81,7 @@ def aperture(self): half_width=self.width_pixel * self.num_pixel / 2, ) - def bin_rays( + def collect( self, rays: optika.rays.RayVectorArray, wavelength: na.AbstractScalar, @@ -136,7 +136,8 @@ def bin_rays( # flux-weighted mean cosine; empty pixels carry no flux, assume normal incidence nonempty = image.outputs > 0 - direction = np.where(nonempty, moment.outputs / image.outputs, 1) + with np.errstate(invalid="ignore", divide="ignore"): + direction = np.where(nonempty, moment.outputs / image.outputs, 1) return image, direction @@ -147,6 +148,7 @@ def expose( na.AbstractScalar, ], direction: float | na.AbstractScalar = 1, + axis_wavelength: None | str = None, timedelta: None | u.Quantity | na.AbstractScalar = None, noise: bool = True, ) -> na.FunctionArray[ @@ -158,7 +160,7 @@ def expose( This is the detector-physics step shared by every optical system: it operates on a pixel grid (a photon image plus an incidence-cosine map), - so it can be driven by :meth:`bin_rays` for a ray-traced system, or by + 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 @@ -171,8 +173,14 @@ def expose( 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 incidence angle in each pixel. + 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` @@ -180,24 +188,32 @@ def expose( 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: + 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 - sgn = np.sign(photons) electrons = self.material.signal( - photons=np.abs(photons), - wavelength=image.inputs.wavelength, + 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=sgn * electrons) + return dataclasses.replace(image, outputs=electrons) - def readout( + def measure( self, rays: optika.rays.RayVectorArray, wavelength: na.AbstractScalar, @@ -213,7 +229,7 @@ def readout( Bin a set of rays onto the pixel grid and convert them to the electrons measured by the sensor. - This composes :meth:`bin_rays` (gather rays into the pixel grid) with + This composes :meth:`collect` (gather rays into the pixel grid) with :meth:`expose` (apply the detector physics). Parameters @@ -233,7 +249,7 @@ def readout( noise Whether to add shot noise and intrinsic sensor noise to the result. """ - image, direction = self.bin_rays( + image, direction = self.collect( rays=rays, wavelength=wavelength, axis=axis, diff --git a/optika/sensors/_sensors_test.py b/optika/sensors/_sensors_test.py index a5c2995c..a48e8f88 100644 --- a/optika/sensors/_sensors_test.py +++ b/optika/sensors/_sensors_test.py @@ -46,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) diff --git a/optika/sensors/materials/_e2v_ccd203/_e2v_ccd203.py b/optika/sensors/materials/_e2v_ccd203/_e2v_ccd203.py index 8d5511bc..d33193ec 100644 --- a/optika/sensors/materials/_e2v_ccd203/_e2v_ccd203.py +++ b/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 diff --git a/optika/sensors/materials/_e2v_ccd97/_e2v_ccd97.py b/optika/sensors/materials/_e2v_ccd97/_e2v_ccd97.py index c445dbd8..b6e0cb21 100644 --- a/optika/sensors/materials/_e2v_ccd97/_e2v_ccd97.py +++ b/optika/sensors/materials/_e2v_ccd97/_e2v_ccd97.py @@ -95,10 +95,8 @@ def e2v_ccd97( # 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), ) @@ -161,10 +159,8 @@ def e2v_ccd97( material_ccd_aia = optika.sensors.materials.e2v_ccd203() eqe_fit_aia = material_ccd_aia.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), ) @@ -205,11 +201,7 @@ def e2v_ccd97( # 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 diff --git a/optika/sensors/materials/_materials.py b/optika/sensors/materials/_materials.py index f937c393..4ccc3794 100644 --- a/optika/sensors/materials/_materials.py +++ b/optika/sensors/materials/_materials.py @@ -602,13 +602,9 @@ def quantum_efficiency_effective( # Define an array of wavelengths with which to sample the EQE wavelength = na.geomspace(10, 10000, axis="wavelength", num=1001) * u.AA - # Define the incidence directions + # Define the cosines of the incidence angles angle = na.linspace(0, 30, axis="angle", num=2) * u.deg - direction = na.Cartesian3dVectorArray( - x=np.sin(angle), - y=0, - z=np.cos(angle), - ) + direction = np.cos(angle) eqe = optika.sensors.quantum_efficiency_effective( wavelength=wavelength, @@ -2103,7 +2099,7 @@ def eqe_rms_difference(x: tuple[float, float, float, float]): ) = x qe_fit = quantum_efficiency_effective( wavelength=eqe_measured.inputs, - direction=na.Cartesian3dVectorArray(0, 0, 1), + direction=1, thickness_oxide=thickness_oxide << unit_thickness_oxide, thickness_implant=thickness_implant << unit_thickness_implant, thickness_substrate=thickness_substrate, diff --git a/optika/sensors/materials/_tektronix_tk512cb/_tektronix_tk512cb.py b/optika/sensors/materials/_tektronix_tk512cb/_tektronix_tk512cb.py index 7cd263cf..6f167983 100644 --- a/optika/sensors/materials/_tektronix_tk512cb/_tektronix_tk512cb.py +++ b/optika/sensors/materials/_tektronix_tk512cb/_tektronix_tk512cb.py @@ -101,10 +101,8 @@ def tektronix_tk512cb( # 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), ) @@ -167,11 +165,7 @@ def tektronix_tk512cb( # 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 diff --git a/optika/surfaces.py b/optika/surfaces.py index 07ec6a86..ae25f1d5 100644 --- a/optika/surfaces.py +++ b/optika/surfaces.py @@ -123,7 +123,6 @@ def is_stop(self) -> bool: def propagate_rays( self, rays: optika.rays.RayVectorArray, - # material: None | optika.materials.AbstractMaterial = None, ) -> optika.rays.RayVectorArray: r""" Refract, reflect, and/or diffract the given rays off of this surface @@ -131,7 +130,7 @@ def propagate_rays( Parameters ---------- rays - a set of input rays that will interact with this surface + A set of input rays that will interact with this surface. """ sag = self.sag material = self.material diff --git a/optika/systems/_sequential.py b/optika/systems/_sequential.py index a304f3b2..971e17a3 100644 --- a/optika/systems/_sequential.py +++ b/optika/systems/_sequential.py @@ -915,7 +915,7 @@ def image( axis=axis_wavelength, ) - return self.sensor.readout( + return self.sensor.measure( rays=rayfunction.outputs, wavelength=wavelength, axis=(axis_wavelength,) + axis_field + axis_pupil, From 6d2ebd8df5eb99dadab01302064c6f1138c115d0 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Sun, 7 Jun 2026 13:00:20 -0600 Subject: [PATCH 05/19] Add a docstring for `RayVectorArray.n` and apply black formatting. - Document the `n` property (complex index of refraction: real part is the index of refraction, imaginary part is the extinction coefficient k = alpha * lambda / 4 pi). - Run black on the sensor material modules to satisfy the formatting check. Co-Authored-By: Claude Opus 4.8 --- optika/rays/_ray_vectors.py | 14 ++++++++++- optika/sensors/materials/_materials.py | 25 +++++++++++-------- .../_ramanathan_2020/_ramanathan_2020.py | 2 +- 3 files changed, 29 insertions(+), 12 deletions(-) diff --git a/optika/rays/_ray_vectors.py b/optika/rays/_ray_vectors.py index a8f4a9d4..49fb0bfc 100644 --- a/optika/rays/_ray_vectors.py +++ b/optika/rays/_ray_vectors.py @@ -72,7 +72,19 @@ def unvignetted(self) -> na.AbstractScalar: @property def n(self) -> complex | na.AbstractScalar: - return self.index_refraction + self.attenuation * self.wavelength / (4 * np.pi) * 1j + """ + 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]: diff --git a/optika/sensors/materials/_materials.py b/optika/sensors/materials/_materials.py index 4ccc3794..28e5540e 100644 --- a/optika/sensors/materials/_materials.py +++ b/optika/sensors/materials/_materials.py @@ -271,7 +271,7 @@ def absorbance( wavelength=wavelength, method="exact", ) - + absorbance_beer = optika.sensors.absorbance( wavelength=wavelength, method="Beer-Lambert", @@ -361,7 +361,9 @@ def absorbance( direction_substrate=direction_substrate, ) - _transmittance_total = _transmittance * np.exp(-absorption * thickness_substrate) + _transmittance_total = _transmittance * np.exp( + -absorption * thickness_substrate + ) result = _transmittance - _transmittance_total @@ -658,7 +660,7 @@ def quantum_efficiency_effective( roughness_oxide=roughness_oxide, roughness_substrate=roughness_substrate, ) - + n_substrate = chemical_substrate.n(wavelength) direction_substrate = optika.materials.snells_law_scalar( @@ -1401,7 +1403,8 @@ def signal( direction: float | na.AbstractScalar = 1, width_pixel: ( u.Quantity | na.AbstractScalar | na.AbstractCartesian2dVectorArray - ) = 0 * u.um, + ) = 0 + * u.um, axis_xy: None | tuple[str, str] = None, noise: bool = True, ) -> na.AbstractScalar: @@ -1470,7 +1473,8 @@ def signal( direction: float | na.AbstractScalar = 1, width_pixel: ( u.Quantity | na.AbstractScalar | na.AbstractCartesian2dVectorArray - ) = 0 * u.um, + ) = 0 + * u.um, axis_xy: None | tuple[str, str] = None, noise: bool = True, ) -> na.AbstractScalar: @@ -1709,10 +1713,10 @@ def absorbance( """ if direction is None: direction = na.Cartesian3dVectorArray(0, 0, 1) - + if normal is None: normal = na.Cartesian3dVectorArray(0, 0, -1) - + return absorbance( wavelength=wavelength, direction=-direction @ normal, @@ -1891,7 +1895,8 @@ def electrons_measured( normal: None | na.AbstractCartesian3dVectorArray = None, width_pixel: ( u.Quantity | na.AbstractScalar | na.AbstractCartesian2dVectorArray - ) = 0 * u.um, + ) = 0 + * u.um, axis_xy: None | tuple[str, str] = None, ) -> na.AbstractScalar: """ @@ -1940,7 +1945,8 @@ def signal( n: complex | na.AbstractScalar = 1, width_pixel: ( u.Quantity | na.AbstractScalar | na.AbstractCartesian2dVectorArray - ) = 0 * u.um, + ) = 0 + * u.um, axis_xy: None | tuple[str, str] = None, noise: bool = True, ) -> na.AbstractScalar: @@ -1972,7 +1978,6 @@ def signal( axis_xy=axis_xy, ) - def photons_incident( self, electrons: u.Quantity | na.AbstractScalar, diff --git a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py index c8114574..7e6bc1cf 100644 --- a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py +++ b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py @@ -769,7 +769,7 @@ def _electrons_measured_numba( energy_pair_inf: np.ndarray, fano_inf: np.ndarray, ) -> np.ndarray: - + num_i, num_x, num_y, num_n = p_n.shape result = np.zeros((num_i, num_x, num_y)) From f0f356f1302b50c936c58b9f86ed140f13b6a01f Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Sun, 7 Jun 2026 18:05:24 -0600 Subject: [PATCH 06/19] Align `electrons_measured_approx()` with `electrons_measured()`. - Rename its first argument `photons_absorbed` -> `photons_transmitted` and model the transmitted->absorbed step (the fraction `1 - e^(-alpha t)` that is absorbed before escaping out the back), so it is now drop-in comparable with the Monte-Carlo `electrons_measured()`: at 5.9 keV the two agree to ~0.003% on the mean and ~0.6% on the std. - Add the matching `thickness_substrate` parameter and align the shared parameter docstrings/summary with `electrons_measured()`. - Fix two latent bugs: the docstring example called `electrons_measured` with the old `photons_absorbed=` keyword, and the broadcast shape accidentally referenced the module-level `absorbance` function. - Update the MC summary line, which still said "photons absorbed" despite the argument now being `photons_transmitted`. Co-Authored-By: Claude Opus 4.8 --- optika/sensors/materials/_materials.py | 57 ++++++++++++------- optika/sensors/materials/_materials_test.py | 8 +-- .../_ramanathan_2020/_ramanathan_2020.py | 7 ++- 3 files changed, 44 insertions(+), 28 deletions(-) diff --git a/optika/sensors/materials/_materials.py b/optika/sensors/materials/_materials.py index 28e5540e..aed01e06 100644 --- a/optika/sensors/materials/_materials.py +++ b/optika/sensors/materials/_materials.py @@ -820,10 +820,11 @@ def _discrete_gamma( def electrons_measured_approx( - photons_absorbed: u.Quantity | na.AbstractScalar, + photons_transmitted: u.Quantity | na.AbstractScalar, wavelength: u.Quantity | na.ScalarArray, absorption: None | u.Quantity | na.AbstractScalar = None, thickness_implant: u.Quantity | na.AbstractScalar = _thickness_implant, + thickness_substrate: u.Quantity | na.AbstractScalar = _thickness_substrate, cce_backsurface: u.Quantity | na.AbstractScalar = _cce_backsurface, temperature: u.Quantity | na.ScalarArray = 300 * u.K, iqy: None | u.Quantity | na.AbstractScalar = None, @@ -832,30 +833,35 @@ def electrons_measured_approx( ) -> na.AbstractScalar: r""" A random sample from an approximate distribution of measured electrons - given the number of photons absorbed by the light-sensitive layer of the - sensor. + given the number of photons transmitted into the light-sensitive region + of the sensor. - This function accounts for both Fano noise and recombination noise due to + This function accounts for the fraction of photons that escape without being + absorbed, as well as both Fano noise and recombination noise due to partial-charge collection. Parameters ---------- - photons_absorbed - The number of photons absorbed by the light-sensitive layer of the sensor. + photons_transmitted + The number of photons transmitted into the light-sensitive region + of the sensor. wavelength The vacuum wavelength of the absorbed photons. absorption - The absorption coefficient of the light-sensitive material for the - wavelength of interest. + The absorption coefficient of silicon per unit perpendicular depth. + For oblique incidence, supply the effective coefficient from + :func:`optika.sensors.absorption_effective`, which folds in the + refracted angle, so no separate angle argument is needed. thickness_implant - The thickness of the implant layer. - Default is the value given in :cite:t:`Stern1994`. + The thickness of the implant layer, where partial-charge collection occurs. + thickness_substrate + The thickness of the entire light-sensitive region of the device. cce_backsurface The differential charge collection efficiency on the back surface of the sensor. - Default is the value given in :cite:t:`Stern1994`. temperature - The temperature of the light-sensitive silicon layer. + The temperature of the silicon detector. + Default is room temperature. iqy The ideal quantum yield of the sensor in electrons per photon. If :obj:`None` (the default), the result of :func:`quantum_yield_ideal` @@ -888,7 +894,7 @@ def electrons_measured_approx( # Define the expected number of photons # for each experiment - photons_absorbed = (100 * u.photon).astype(int) + photons_transmitted = (100 * u.photon).astype(int) # Define the wavelength at which to sample the distribution wavelength = 5.9 * u.keV @@ -896,14 +902,14 @@ def electrons_measured_approx( # Compute the actual number of electrons measured for each experiment electrons_exact = optika.sensors.electrons_measured( - photons_absorbed=photons_absorbed, + photons_transmitted=photons_transmitted, wavelength=wavelength, shape_random=dict(experiment=num_experiments), ) # Compute the approximate number of electrons measured for each experiment electrons_approx = optika.sensors.electrons_measured_approx( - photons_absorbed=photons_absorbed, + photons_transmitted=photons_transmitted, wavelength=wavelength, shape_random=dict(experiment=num_experiments), ) @@ -962,11 +968,11 @@ def electrons_measured_approx( shape_random = dict() shape = na.shape_broadcasted( - photons_absorbed, + photons_transmitted, absorption, - absorbance, iqy, thickness_implant, + thickness_substrate, cce_backsurface, fano_factor, ) @@ -978,15 +984,24 @@ def electrons_measured_approx( W = thickness_implant n0 = cce_backsurface aW = (a * W).to(u.dimensionless_unscaled).value + aDW = (a * (thickness_substrate - W)).to(u.dimensionless_unscaled).value - fraction_complete = np.exp(-aW) - fraction_partial = 1 - fraction_complete + # A transmitted photon is absorbed in the implant region (partial charge + # collection) with probability 1 - exp(-alpha * W); of those penetrating + # deeper, the fraction absorbed within the remaining substrate (complete + # charge collection) is 1 - exp(-alpha * (D - W)), and the rest escape. + fraction_partial = 1 - np.exp(-aW) photons_absorbed_partial = na.random.binomial( - n=photons_absorbed, + n=photons_transmitted, p=fraction_partial, shape_random=shape, ) - photons_absorbed_complete = photons_absorbed - photons_absorbed_partial + fraction_complete = 1 - np.exp(-aDW) + photons_absorbed_complete = na.random.binomial( + n=(photons_transmitted - photons_absorbed_partial).astype(int), + p=fraction_complete, + shape_random=shape, + ) mean_p = n0 + (1 - n0) / (aW) + (1 - n0) / (1 - np.exp(aW)) var_p = np.square(n0 - 1) * (4 / np.square(aW) - 1 / np.square(np.sinh(aW / 2))) / 4 diff --git a/optika/sensors/materials/_materials_test.py b/optika/sensors/materials/_materials_test.py index ba4154ff..fe5766d8 100644 --- a/optika/sensors/materials/_materials_test.py +++ b/optika/sensors/materials/_materials_test.py @@ -222,7 +222,7 @@ def test_probability_measurement( @pytest.mark.parametrize( - argnames="photons_absorbed", + argnames="photons_transmitted", argvalues=[ (100 * u.photon).astype(int), ], @@ -248,14 +248,14 @@ def test_probability_measurement( ) @pytest.mark.parametrize("temperature", [300 * u.K]) def test_electrons_measured_approx( - photons_absorbed: u.Quantity | na.AbstractScalar, + photons_transmitted: u.Quantity | na.AbstractScalar, wavelength: u.Quantity | na.ScalarArray, thickness_implant: u.Quantity | na.AbstractScalar, cce_backsurface: u.Quantity | na.AbstractScalar, temperature: u.Quantity | na.ScalarArray, ): result = optika.sensors.electrons_measured_approx( - photons_absorbed=photons_absorbed, + photons_transmitted=photons_transmitted, wavelength=wavelength, thickness_implant=thickness_implant, cce_backsurface=cce_backsurface, @@ -265,7 +265,7 @@ def test_electrons_measured_approx( assert np.all(result >= 0 * u.electron) shape = na.shape_broadcasted( - photons_absorbed, + photons_transmitted, wavelength, thickness_implant, cce_backsurface, diff --git a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py index 7e6bc1cf..5a436188 100644 --- a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py +++ b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py @@ -482,10 +482,11 @@ def electrons_measured( ) -> na.AbstractScalar: r""" A random sample from the distribution of measured electrons - given the number of photons absorbed by the light-sensitive layer of the - sensor. + given the number of photons transmitted into the light-sensitive region + of the sensor. - This function accounts for both Fano noise and recombination noise due to + This function accounts for the fraction of photons that escape without being + absorbed, as well as both Fano noise and recombination noise due to partial-charge collection. Parameters From a410c1ec10905ab66a988f48be1ca97bd27650e2 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Sun, 7 Jun 2026 22:50:40 -0600 Subject: [PATCH 07/19] Fix patch coverage on the sensor refactor. - Mark the defensive `ValueError` branches with `# pragma: nocover`, matching the existing convention (e.g. `_sequential.py`): the `expose` ambiguous wavelength-axis check and the `absorbance` method dispatch. - Mark the `@numba.njit` `_electrons_measured_numba` kernel `# pragma: nocover`; JIT-compiled code cannot be traced by coverage.py. - Extend `test_electrons_measured` with an `axis_xy` parametrize and a per-pixel photon grid so the charge-diffusion branch of `electrons_measured` is exercised. `_ramanathan_2020.py` is now at 100% line coverage. Co-Authored-By: Claude Opus 4.8 --- optika/sensors/_sensors.py | 2 +- optika/sensors/materials/_materials.py | 2 +- .../materials/_ramanathan_2020/_ramanathan_2020.py | 2 +- .../_ramanathan_2020/_ramanathan_2020_test.py | 11 ++++++++++- 4 files changed, 13 insertions(+), 4 deletions(-) diff --git a/optika/sensors/_sensors.py b/optika/sensors/_sensors.py index dbedc940..38262cbc 100644 --- a/optika/sensors/_sensors.py +++ b/optika/sensors/_sensors.py @@ -190,7 +190,7 @@ def expose( """ if axis_wavelength is None: shape_wavelength = na.shape(image.inputs.wavelength) - if len(shape_wavelength) != 1: + 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}." diff --git a/optika/sensors/materials/_materials.py b/optika/sensors/materials/_materials.py index aed01e06..0ca42b5e 100644 --- a/optika/sensors/materials/_materials.py +++ b/optika/sensors/materials/_materials.py @@ -367,7 +367,7 @@ def absorbance( result = _transmittance - _transmittance_total - else: + else: # pragma: nocover raise ValueError(f"Method {method} not implemented") return np.real(result) diff --git a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py index 5a436188..fc8d40a4 100644 --- a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py +++ b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py @@ -755,7 +755,7 @@ def _electrons_measured_quantity( fastmath=True, parallel=True, ) -def _electrons_measured_numba( +def _electrons_measured_numba( # pragma: nocover photons_transmitted: np.ndarray, energy: np.ndarray, absorption: np.ndarray, diff --git a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020_test.py b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020_test.py index 59c304f3..02151268 100644 --- a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020_test.py +++ b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020_test.py @@ -100,7 +100,14 @@ def test_fano_factor_inf( @pytest.mark.parametrize( argnames="photons_transmitted", argvalues=[ - 100 * u.photon, + na.broadcast_to((100 * u.photon).astype(int), dict(pixel_x=2, pixel_y=2)), + ], +) +@pytest.mark.parametrize( + argnames="axis_xy", + argvalues=[ + None, + ("pixel_x", "pixel_y"), ], ) @pytest.mark.parametrize( @@ -140,6 +147,7 @@ def test_electrons_measured( cce_backsurface: u.Quantity | na.AbstractScalar, temperature: u.Quantity | na.ScalarArray, width_pixel: u.Quantity | na.AbstractCartesian2dVectorArray, + axis_xy: None | tuple[str, str], ): result = _ramanathan_2020.electrons_measured( photons_transmitted=photons_transmitted, @@ -149,6 +157,7 @@ def test_electrons_measured( width_pixel=width_pixel, cce_backsurface=cce_backsurface, temperature=temperature, + axis_xy=axis_xy, ) assert np.all(result >= 0 * u.electron) From dec99d9c8cca6372fdeec1b01da11c4703159ef0 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Mon, 8 Jun 2026 06:48:46 -0600 Subject: [PATCH 08/19] Cover the remaining new branches in the sensor material tests. - Parametrize `test_absorbance` over `method` ("Beer-Lambert" and "exact") so the exact (`layer_absorbance`) branch is exercised. - Add a `transmittance` method test mirroring the other method tests, covering its `direction`/`normal` None-default handling. Co-Authored-By: Claude Opus 4.8 --- optika/sensors/materials/_materials_test.py | 42 +++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/optika/sensors/materials/_materials_test.py b/optika/sensors/materials/_materials_test.py index fe5766d8..93567629 100644 --- a/optika/sensors/materials/_materials_test.py +++ b/optika/sensors/materials/_materials_test.py @@ -89,12 +89,20 @@ def test_transmittance( 10 * u.um, ], ) +@pytest.mark.parametrize( + argnames="method", + argvalues=[ + "Beer-Lambert", + "exact", + ], +) def test_absorbance( wavelength: u.Quantity | na.AbstractScalar, direction: float | na.AbstractScalar, n: float | na.AbstractScalar, thickness_oxide: u.Quantity | na.AbstractScalar, thickness_substrate: u.Quantity | na.AbstractScalar, + method: str, ): result = optika.sensors.absorbance( wavelength=wavelength, @@ -102,6 +110,7 @@ def test_absorbance( n=n, thickness_oxide=thickness_oxide, thickness_substrate=thickness_substrate, + method=method, ) assert np.all(result >= 0) @@ -510,6 +519,39 @@ def test_width_charge_diffusion( result = a.width_charge_diffusion(wavelength=wavelength) assert np.all(result >= 0 * u.um) + @pytest.mark.parametrize( + argnames="wavelength", + argvalues=[ + 100 * u.AA, + ], + ) + @pytest.mark.parametrize( + argnames="direction", + argvalues=[ + None, + ], + ) + @pytest.mark.parametrize( + argnames="normal", + argvalues=[ + None, + ], + ) + def test_transmittance( + self, + a: optika.sensors.materials.AbstractBackIlluminatedSiliconSensorMaterial, + wavelength: u.Quantity | na.AbstractScalar, + direction: None | na.AbstractCartesian3dVectorArray, + normal: None | na.AbstractCartesian3dVectorArray, + ): + result = a.transmittance( + wavelength=wavelength, + direction=direction, + normal=normal, + ) + assert np.all(result >= 0) + assert np.all(result <= 1) + @pytest.mark.parametrize( argnames="wavelength", argvalues=[ From e4bd0b53720695b4525725f7f6d1c2d9a849315e Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Mon, 8 Jun 2026 09:41:31 -0600 Subject: [PATCH 09/19] Rework sensor-material model to use `photons_absorbed` instead of `photons_transmitted`. The Monte-Carlo kernel now samples the absorption depth from an exponential truncated to the light-sensitive region, so every input photon is treated as absorbed. The escape fraction is applied as Poisson thinning in `signal()` (exactly equivalent). `electrons_measured_approx()` reverts to splitting the absorbed photons into partial/complete charge collection, and `signal()` computes absorbance from the start rather than transmittance. Co-Authored-By: Claude Opus 4.8 --- optika/sensors/materials/_materials.py | 71 ++++++++----------- optika/sensors/materials/_materials_test.py | 14 ++-- .../_ramanathan_2020/_ramanathan_2020.py | 46 ++++++------ .../_ramanathan_2020/_ramanathan_2020_test.py | 8 +-- 4 files changed, 62 insertions(+), 77 deletions(-) diff --git a/optika/sensors/materials/_materials.py b/optika/sensors/materials/_materials.py index 0ca42b5e..6a216530 100644 --- a/optika/sensors/materials/_materials.py +++ b/optika/sensors/materials/_materials.py @@ -820,11 +820,10 @@ def _discrete_gamma( def electrons_measured_approx( - photons_transmitted: u.Quantity | na.AbstractScalar, + photons_absorbed: u.Quantity | na.AbstractScalar, wavelength: u.Quantity | na.ScalarArray, absorption: None | u.Quantity | na.AbstractScalar = None, thickness_implant: u.Quantity | na.AbstractScalar = _thickness_implant, - thickness_substrate: u.Quantity | na.AbstractScalar = _thickness_substrate, cce_backsurface: u.Quantity | na.AbstractScalar = _cce_backsurface, temperature: u.Quantity | na.ScalarArray = 300 * u.K, iqy: None | u.Quantity | na.AbstractScalar = None, @@ -833,17 +832,16 @@ def electrons_measured_approx( ) -> na.AbstractScalar: r""" A random sample from an approximate distribution of measured electrons - given the number of photons transmitted into the light-sensitive region + given the number of photons absorbed by the light-sensitive layer of the sensor. - This function accounts for the fraction of photons that escape without being - absorbed, as well as both Fano noise and recombination noise due to + This function accounts for both Fano noise and recombination noise due to partial-charge collection. Parameters ---------- - photons_transmitted - The number of photons transmitted into the light-sensitive region + photons_absorbed + The number of photons absorbed by the light-sensitive layer of the sensor. wavelength The vacuum wavelength of the absorbed photons. @@ -854,8 +852,6 @@ def electrons_measured_approx( refracted angle, so no separate angle argument is needed. thickness_implant The thickness of the implant layer, where partial-charge collection occurs. - thickness_substrate - The thickness of the entire light-sensitive region of the device. cce_backsurface The differential charge collection efficiency on the back surface of the sensor. @@ -894,7 +890,7 @@ def electrons_measured_approx( # Define the expected number of photons # for each experiment - photons_transmitted = (100 * u.photon).astype(int) + photons_absorbed = (100 * u.photon).astype(int) # Define the wavelength at which to sample the distribution wavelength = 5.9 * u.keV @@ -902,14 +898,14 @@ def electrons_measured_approx( # Compute the actual number of electrons measured for each experiment electrons_exact = optika.sensors.electrons_measured( - photons_transmitted=photons_transmitted, + photons_absorbed=photons_absorbed, wavelength=wavelength, shape_random=dict(experiment=num_experiments), ) # Compute the approximate number of electrons measured for each experiment electrons_approx = optika.sensors.electrons_measured_approx( - photons_transmitted=photons_transmitted, + photons_absorbed=photons_absorbed, wavelength=wavelength, shape_random=dict(experiment=num_experiments), ) @@ -968,11 +964,10 @@ def electrons_measured_approx( shape_random = dict() shape = na.shape_broadcasted( - photons_transmitted, + photons_absorbed, absorption, iqy, thickness_implant, - thickness_substrate, cce_backsurface, fano_factor, ) @@ -984,24 +979,17 @@ def electrons_measured_approx( W = thickness_implant n0 = cce_backsurface aW = (a * W).to(u.dimensionless_unscaled).value - aDW = (a * (thickness_substrate - W)).to(u.dimensionless_unscaled).value - # A transmitted photon is absorbed in the implant region (partial charge - # collection) with probability 1 - exp(-alpha * W); of those penetrating - # deeper, the fraction absorbed within the remaining substrate (complete - # charge collection) is 1 - exp(-alpha * (D - W)), and the rest escape. + # An absorbed photon is absorbed within the implant region (partial charge + # collection) with probability 1 - exp(-alpha * W); otherwise it is absorbed + # deeper in the substrate (complete charge collection). fraction_partial = 1 - np.exp(-aW) photons_absorbed_partial = na.random.binomial( - n=photons_transmitted, + n=photons_absorbed, p=fraction_partial, shape_random=shape, ) - fraction_complete = 1 - np.exp(-aDW) - photons_absorbed_complete = na.random.binomial( - n=(photons_transmitted - photons_absorbed_partial).astype(int), - p=fraction_complete, - shape_random=shape, - ) + photons_absorbed_complete = photons_absorbed - photons_absorbed_partial mean_p = n0 + (1 - n0) / (aW) + (1 - n0) / (1 - np.exp(aW)) var_p = np.square(n0 - 1) * (4 / np.square(aW) - 1 / np.square(np.sinh(aW / 2))) / 4 @@ -1035,6 +1023,7 @@ def electrons_measured_approx( _transmittance = transmittance +_absorbance = absorbance def signal( @@ -1043,7 +1032,7 @@ def signal( direction: float | na.AbstractScalar = 1, n: complex | na.AbstractScalar = 1, n_substrate: None | complex | na.AbstractScalar = None, - transmittance: None | float | na.AbstractScalar = None, + absorbance: None | float | na.AbstractScalar = None, thickness_implant: u.Quantity | na.AbstractScalar = _thickness_implant, thickness_depletion: u.Quantity | na.AbstractScalar = _thickness_substrate, thickness_substrate: None | na.AbstractScalar = _thickness_substrate, @@ -1079,10 +1068,10 @@ def signal( The complex index of refraction of the light-sensitive material If :obj:`None` (the default), the result of :meth:`optika.chemicals.Chemical.n` for silicon will be used. - transmittance - The fraction of incident energy transmitted to the light-sensitive layer + absorbance + The fraction of incident energy absorbed by the light-sensitive region of the detector. - If :obj:`None` (the default), the result of :func:`transmittance` + If :obj:`None` (the default), the result of :func:`absorbance` called with default values will be used. thickness_implant The thickness of the implant layer. @@ -1162,11 +1151,12 @@ def signal( ax.set_ylabel(f"variance-to-mean ratio ({electrons.unit:latex_inline})"); """ - if transmittance is None: - transmittance = _transmittance( + if absorbance is None: + absorbance = _absorbance( wavelength=wavelength, direction=direction, n=n, + thickness_substrate=thickness_substrate, ).average if n_substrate is None: @@ -1190,10 +1180,6 @@ def signal( temperature=temperature, ) - transmittance_total = transmittance * np.exp(-absorption * thickness_substrate) - - absorbance = transmittance - transmittance_total - cce = charge_collection_efficiency( absorption=absorption, thickness_implant=thickness_implant, @@ -1203,14 +1189,13 @@ def signal( elif method == "monte-carlo": - photons_expected = transmittance * photons_expected.to(u.ph) photons = na.random.poisson( - lam=photons_expected, + lam=absorbance * photons_expected.to(u.ph), shape_random=shape_random, ).astype(int) return electrons_measured( - photons_transmitted=photons, + photons_absorbed=photons, wavelength=wavelength, absorption=absorption, thickness_implant=thickness_implant, @@ -1903,7 +1888,7 @@ def probability_measurement( def electrons_measured( self, - photons_transmitted: na.AbstractScalar, + photons_absorbed: na.AbstractScalar, wavelength: u.Quantity | na.AbstractScalar, direction: None | na.AbstractCartesian3dVectorArray = None, n: complex | na.AbstractScalar = 1, @@ -1936,7 +1921,7 @@ def electrons_measured( ) return electrons_measured( - photons_transmitted=photons_transmitted, + photons_absorbed=photons_absorbed, wavelength=wavelength, absorption=absorption_effective( wavelength=wavelength, @@ -1982,7 +1967,7 @@ def signal( direction=direction, n=n, n_substrate=self._chemical.n(wavelength), - transmittance=1, + absorbance=1, thickness_implant=self.thickness_implant, thickness_depletion=self.depletion.thickness, thickness_substrate=self.thickness_substrate, @@ -2033,7 +2018,7 @@ def efficiency( rays: optika.rays.AbstractRayVectorArray, normal: na.AbstractCartesian3dVectorArray, ) -> na.ScalarLike: - result = self.transmittance( + result = self.absorbance( wavelength=rays.wavelength, direction=rays.direction, n=rays.n, diff --git a/optika/sensors/materials/_materials_test.py b/optika/sensors/materials/_materials_test.py index 93567629..856ba537 100644 --- a/optika/sensors/materials/_materials_test.py +++ b/optika/sensors/materials/_materials_test.py @@ -231,7 +231,7 @@ def test_probability_measurement( @pytest.mark.parametrize( - argnames="photons_transmitted", + argnames="photons_absorbed", argvalues=[ (100 * u.photon).astype(int), ], @@ -257,14 +257,14 @@ def test_probability_measurement( ) @pytest.mark.parametrize("temperature", [300 * u.K]) def test_electrons_measured_approx( - photons_transmitted: u.Quantity | na.AbstractScalar, + photons_absorbed: u.Quantity | na.AbstractScalar, wavelength: u.Quantity | na.ScalarArray, thickness_implant: u.Quantity | na.AbstractScalar, cce_backsurface: u.Quantity | na.AbstractScalar, temperature: u.Quantity | na.ScalarArray, ): result = optika.sensors.electrons_measured_approx( - photons_transmitted=photons_transmitted, + photons_absorbed=photons_absorbed, wavelength=wavelength, thickness_implant=thickness_implant, cce_backsurface=cce_backsurface, @@ -274,7 +274,7 @@ def test_electrons_measured_approx( assert np.all(result >= 0 * u.electron) shape = na.shape_broadcasted( - photons_transmitted, + photons_absorbed, wavelength, thickness_implant, cce_backsurface, @@ -684,7 +684,7 @@ def test_probability_measurement( assert np.all(result <= 1) @pytest.mark.parametrize( - argnames="photons_transmitted", + argnames="photons_absorbed", argvalues=[ (100 * u.photon).astype(int), ], @@ -710,13 +710,13 @@ def test_probability_measurement( def test_electrons_measured( self, a: optika.sensors.materials.AbstractBackIlluminatedSiliconSensorMaterial, - photons_transmitted: u.Quantity | na.AbstractScalar, + photons_absorbed: u.Quantity | na.AbstractScalar, wavelength: u.Quantity | na.AbstractScalar, direction: None | na.AbstractCartesian3dVectorArray, normal: None | na.AbstractCartesian3dVectorArray, ): result = a.electrons_measured( - photons_transmitted=photons_transmitted, + photons_absorbed=photons_absorbed, wavelength=wavelength, direction=direction, normal=normal, diff --git a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py index fc8d40a4..8ac3a732 100644 --- a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py +++ b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py @@ -466,7 +466,7 @@ def probability_of_n_pairs( def electrons_measured( - photons_transmitted: u.Quantity | na.AbstractScalar, + photons_absorbed: u.Quantity | na.AbstractScalar, wavelength: u.Quantity | na.ScalarArray, absorption: None | u.Quantity | na.AbstractScalar = None, thickness_implant: u.Quantity | na.AbstractScalar = _thickness_implant, @@ -482,18 +482,17 @@ def electrons_measured( ) -> na.AbstractScalar: r""" A random sample from the distribution of measured electrons - given the number of photons transmitted into the light-sensitive region - of the sensor. + given the number of photons absorbed by the light-sensitive layer of the + sensor. - This function accounts for the fraction of photons that escape without being - absorbed, as well as both Fano noise and recombination noise due to + This function accounts for both Fano noise and recombination noise due to partial-charge collection. Parameters ---------- - photons_transmitted - The number of photons transmitted into the light-sensitive region - of the sensor. + photons_absorbed + The number of photons absorbed by the light-sensitive layer of the + sensor. wavelength The vacuum wavelength of the absorbed photons. absorption @@ -548,7 +547,7 @@ def electrons_measured( # Define the expected number of photons # for each experiment - photons_transmitted = (100 * u.photon).astype(int) + photons_absorbed = (100 * u.photon).astype(int) # Define the wavelength at which to sample the distribution wavelength = 5.9 * u.keV @@ -556,7 +555,7 @@ def electrons_measured( # Compute the actual number of electrons measured for each experiment electrons = optika.sensors.electrons_measured( - photons_transmitted=photons_transmitted, + photons_absorbed=photons_absorbed, wavelength=wavelength, shape_random=dict(experiment=num_experiments), ) @@ -604,7 +603,7 @@ def electrons_measured( width_pixel_y = width_pixel.y shape = na.broadcast_shapes( - na.shape(photons_transmitted), + na.shape(photons_absorbed), na.shape(wavelength), na.shape(absorption), na.shape(thickness_implant), @@ -627,7 +626,7 @@ def electrons_measured( shape[axis_x] = 1 shape[axis_y] = 1 - photons_transmitted = na.broadcast_to(photons_transmitted, shape) + photons_absorbed = na.broadcast_to(photons_absorbed, shape) absorption = na.broadcast_to(absorption, shape) thickness_implant = na.broadcast_to(thickness_implant, shape) thickness_depletion = na.broadcast_to(thickness_depletion, shape) @@ -657,7 +656,7 @@ def electrons_measured( fano_inf = fano_inf.broadcast_to(shape) result = _electrons_measured_quantity( - photons_transmitted=photons_transmitted.ndarray, + photons_absorbed=photons_absorbed.ndarray, wavelength=wavelength.ndarray, absorption=absorption.ndarray, thickness_implant=thickness_implant.ndarray, @@ -684,7 +683,7 @@ def electrons_measured( def _electrons_measured_quantity( - photons_transmitted: u.Quantity, + photons_absorbed: u.Quantity, wavelength: u.Quantity, absorption: u.Quantity, thickness_implant: u.Quantity, @@ -700,7 +699,7 @@ def _electrons_measured_quantity( ) -> u.Quantity: shape = np.broadcast_shapes( - photons_transmitted.shape, + photons_absorbed.shape, absorption.shape, thickness_implant.shape, thickness_depletion.shape, @@ -715,7 +714,7 @@ def _electrons_measured_quantity( unit_length = u.mm - photons_transmitted = photons_transmitted.to_value(u.photon) + photons_absorbed = photons_absorbed.to_value(u.photon) energy = wavelength.to(u.eV, equivalencies=u.spectral()) absorption = absorption.to_value(1 / unit_length) thickness_implant = thickness_implant.to_value(unit_length) @@ -728,7 +727,7 @@ def _electrons_measured_quantity( fano_inf = fano_inf.to_value(u.electron / u.photon) result = _electrons_measured_numba( - photons_transmitted=photons_transmitted.reshape(-1, num_x, num_y), + photons_absorbed=photons_absorbed.reshape(-1, num_x, num_y), energy=energy.reshape(-1, num_x, num_y), absorption=absorption.reshape(-1, num_x, num_y), thickness_implant=thickness_implant.reshape(-1, num_x, num_y), @@ -756,7 +755,7 @@ def _electrons_measured_quantity( parallel=True, ) def _electrons_measured_numba( # pragma: nocover - photons_transmitted: np.ndarray, + photons_absorbed: np.ndarray, energy: np.ndarray, absorption: np.ndarray, thickness_implant: np.ndarray, @@ -778,7 +777,7 @@ def _electrons_measured_numba( # pragma: nocover for i in numba.prange(num_i): for x in range(num_x): for y in range(num_y): - num_photon = int(photons_transmitted[i, x, y]) + num_photon = int(photons_absorbed[i, x, y]) energy_i = energy[i, x, y] a = absorption[i, x, y] W = thickness_implant[i, x, y] @@ -794,6 +793,10 @@ def _electrons_measured_numba( # pragma: nocover d = 1 / a + # every photon is absorbed within the substrate, so sample the + # depth from the exponential truncated to [0, z_substrate) + fraction_absorbed = 1 - math.exp(-a * z_substrate) + mean_inf = energy_i / energy_pair_inf_i std_inf = math.sqrt(fano_inf_i * mean_inf) @@ -817,7 +820,7 @@ def _electrons_measured_numba( # pragma: nocover n_ij = round(n_ij) y_ij = random.uniform(0, 1) - z_ij = -d * math.log(1 - y_ij) + z_ij = -d * math.log(1 - y_ij * fraction_absorbed) if z_ij < W: h_ij = h_0 + (1 - h_0) * z_ij / W @@ -839,9 +842,6 @@ def _electrons_measured_numba( # pragma: nocover p = round(p) q = round(q) - elif z_ij > z_substrate: - continue - else: p = q = 0 diff --git a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020_test.py b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020_test.py index 02151268..dd18a87b 100644 --- a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020_test.py +++ b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020_test.py @@ -98,7 +98,7 @@ def test_fano_factor_inf( @pytest.mark.parametrize( - argnames="photons_transmitted", + argnames="photons_absorbed", argvalues=[ na.broadcast_to((100 * u.photon).astype(int), dict(pixel_x=2, pixel_y=2)), ], @@ -140,7 +140,7 @@ def test_fano_factor_inf( ) @pytest.mark.parametrize("temperature", _temperture) def test_electrons_measured( - photons_transmitted: u.Quantity | na.AbstractScalar, + photons_absorbed: u.Quantity | na.AbstractScalar, wavelength: u.Quantity | na.ScalarArray, absorption: u.Quantity | na.AbstractScalar, thickness_implant: u.Quantity | na.AbstractScalar, @@ -150,7 +150,7 @@ def test_electrons_measured( axis_xy: None | tuple[str, str], ): result = _ramanathan_2020.electrons_measured( - photons_transmitted=photons_transmitted, + photons_absorbed=photons_absorbed, wavelength=wavelength, absorption=absorption, thickness_implant=thickness_implant, @@ -163,7 +163,7 @@ def test_electrons_measured( assert np.all(result >= 0 * u.electron) shape = na.shape_broadcasted( - photons_transmitted, + photons_absorbed, wavelength, absorption, thickness_implant, From f809d829cf2f148c25c95ac340d853fcffb7e9d9 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Mon, 8 Jun 2026 13:34:41 -0600 Subject: [PATCH 10/19] Condition `electrons_measured_approx` partial fraction on absorption within the substrate. The exact `electrons_measured` kernel samples the absorption depth from an exponential truncated to the substrate `[0, D)`, so the approximate model must do the same. The partial-collection fraction is now the implant fraction conditioned on absorption within the substrate, `(1 - exp(-alpha*W)) / (1 - exp(-alpha*D))`, which restores agreement between `electrons_measured_approx` and `electrons_measured`. Co-Authored-By: Claude Opus 4.8 --- optika/sensors/materials/_materials.py | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/optika/sensors/materials/_materials.py b/optika/sensors/materials/_materials.py index 6a216530..9d85c646 100644 --- a/optika/sensors/materials/_materials.py +++ b/optika/sensors/materials/_materials.py @@ -824,6 +824,7 @@ def electrons_measured_approx( wavelength: u.Quantity | na.ScalarArray, absorption: None | u.Quantity | na.AbstractScalar = None, thickness_implant: u.Quantity | na.AbstractScalar = _thickness_implant, + thickness_substrate: u.Quantity | na.AbstractScalar = _thickness_substrate, cce_backsurface: u.Quantity | na.AbstractScalar = _cce_backsurface, temperature: u.Quantity | na.ScalarArray = 300 * u.K, iqy: None | u.Quantity | na.AbstractScalar = None, @@ -852,6 +853,10 @@ def electrons_measured_approx( refracted angle, so no separate angle argument is needed. thickness_implant The thickness of the implant layer, where partial-charge collection occurs. + thickness_substrate + The thickness of the entire light-sensitive region of the device. + The absorbed photons are distributed within this region, which sets + the fraction that land in the implant layer. cce_backsurface The differential charge collection efficiency on the back surface of the sensor. @@ -968,6 +973,7 @@ def electrons_measured_approx( absorption, iqy, thickness_implant, + thickness_substrate, cce_backsurface, fano_factor, ) @@ -977,13 +983,18 @@ def electrons_measured_approx( a = absorption W = thickness_implant + D = thickness_substrate n0 = cce_backsurface aW = (a * W).to(u.dimensionless_unscaled).value - - # An absorbed photon is absorbed within the implant region (partial charge - # collection) with probability 1 - exp(-alpha * W); otherwise it is absorbed - # deeper in the substrate (complete charge collection). - fraction_partial = 1 - np.exp(-aW) + aD = (a * D).to(u.dimensionless_unscaled).value + + # The absorbed photons are distributed within the substrate `[0, D)`, so the + # depth follows an exponential truncated to that region (matching the exact + # `electrons_measured` kernel). Conditioned on absorption within the + # substrate, a photon lands in the implant region (partial charge collection) + # with probability (1 - exp(-alpha * W)) / (1 - exp(-alpha * D)); otherwise + # it is absorbed deeper in the substrate (complete charge collection). + fraction_partial = (1 - np.exp(-aW)) / (1 - np.exp(-aD)) photons_absorbed_partial = na.random.binomial( n=photons_absorbed, p=fraction_partial, From 47c9fc25e99228b0834070dafb0b197b8d26d257 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Mon, 8 Jun 2026 13:59:53 -0600 Subject: [PATCH 11/19] changed examples to use 20 photons --- optika/sensors/materials/_materials.py | 4 ++-- optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/optika/sensors/materials/_materials.py b/optika/sensors/materials/_materials.py index 9d85c646..54706fb2 100644 --- a/optika/sensors/materials/_materials.py +++ b/optika/sensors/materials/_materials.py @@ -879,7 +879,7 @@ def electrons_measured_approx( Examples -------- - Plot the energy spectrum of 100 6 keV photons emitted from an Fe-55 + Plot the energy spectrum of twenty 6 keV photons emitted from an Fe-55 radioactive source and compare it to the exact spectrum .. jupyter-execute:: @@ -895,7 +895,7 @@ def electrons_measured_approx( # Define the expected number of photons # for each experiment - photons_absorbed = (100 * u.photon).astype(int) + photons_absorbed = (20 * u.photon).astype(int) # Define the wavelength at which to sample the distribution wavelength = 5.9 * u.keV diff --git a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py index 8ac3a732..37b75fda 100644 --- a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py +++ b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py @@ -531,7 +531,7 @@ def electrons_measured( Examples -------- - Plot the energy spectrum of 100 6 keV photons emitted from an Fe-55 + Plot the energy spectrum of twenty 6 keV photons emitted from an Fe-55 radioactive source. .. jupyter-execute:: @@ -547,7 +547,7 @@ def electrons_measured( # Define the expected number of photons # for each experiment - photons_absorbed = (100 * u.photon).astype(int) + photons_absorbed = (20 * u.photon).astype(int) # Define the wavelength at which to sample the distribution wavelength = 5.9 * u.keV From 8d28f6c24ca1b24367a2019814a981986b516ff7 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Mon, 8 Jun 2026 17:01:33 -0600 Subject: [PATCH 12/19] Fix `signal` method docstring and add a charge-diffusion test. The `signal()` `method` parameter docstring described nonexistent `exact`/ `approx` methods; it now documents the actual `monte-carlo` and `expected` options, and notes that the per-pixel `expected` method applies no charge diffusion. Adds `test_electrons_measured_diffusion`, which injects photons into a single pixel and checks that the spatial spread of the diffused charge matches the analytic width from `optika.sensors.charge_diffusion`. Co-Authored-By: Claude Opus 4.8 --- optika/sensors/materials/_materials.py | 14 ++--- .../_ramanathan_2020/_ramanathan_2020_test.py | 57 +++++++++++++++++++ 2 files changed, 64 insertions(+), 7 deletions(-) diff --git a/optika/sensors/materials/_materials.py b/optika/sensors/materials/_materials.py index 54706fb2..bd3af815 100644 --- a/optika/sensors/materials/_materials.py +++ b/optika/sensors/materials/_materials.py @@ -1104,13 +1104,13 @@ def signal( temperature The temperature of the light-sensitive silicon layer. method - The method used to generate random samples of measured electrons. - The `exact` method simulates every photon so it is accurate for low - photon counts but slow for high photon counts. - The `approx` method is much faster, but is only accurate if the number - of photons is high. - The `expected` method does not add any noise to the signal and just - returns the expected number of electrons. + The method used to generate samples of measured electrons. + The `monte-carlo` method draws a random sample by simulating every + photon using :func:`electrons_measured`, including shot, Fano, and + recombination noise as well as charge diffusion. + The `expected` method adds no noise and just returns the expected + number of electrons in each pixel; since it is a per-pixel expectation, + it does not apply charge diffusion. axis_xy The two logical axes corresponding to the pixel grid of the sensor along which electrons will diffuse. diff --git a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020_test.py b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020_test.py index dd18a87b..8b099f08 100644 --- a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020_test.py +++ b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020_test.py @@ -2,6 +2,7 @@ import numpy as np import astropy.units as u import named_arrays as na +import optika from . import _ramanathan_2020 _wavelength = [ @@ -172,3 +173,59 @@ def test_electrons_measured( ) assert result.shape == shape + + +def test_electrons_measured_diffusion(): + """ + Inject many photons into a single pixel and check that the spatial spread + of the diffused charge matches the analytic charge-diffusion width given by + :func:`optika.sensors.charge_diffusion`. + """ + num = 41 + axis_xy = ("pixel_x", "pixel_y") + + absorption = 1 / u.um + thickness_substrate = 14 * u.um + thickness_depletion = 0 * u.um + width_pixel = 3 * u.um + + # Place all of the photons in the central pixel of an otherwise empty grid. + photons = np.zeros((num, num)) + photons[num // 2, num // 2] = 20000 + photons = na.ScalarArray(photons << u.photon, axes=axis_xy).astype(int) + + electrons = _ramanathan_2020.electrons_measured( + photons_absorbed=photons, + wavelength=500 * u.nm, + absorption=absorption, + thickness_implant=0 * u.um, + thickness_depletion=thickness_depletion, + thickness_substrate=thickness_substrate, + width_pixel=width_pixel, + cce_backsurface=1, + axis_xy=axis_xy, + ) + + # The diffused charge should spread beyond the central pixel. + assert electrons[{axis_xy[0]: num // 2, axis_xy[1]: num // 2}] < electrons.sum( + axis_xy + ) + + # Physical offset of each pixel from the center of the grid. + offset_x = (na.arange(0, num, axis=axis_xy[0]) - num // 2) * width_pixel + offset_y = (na.arange(0, num, axis=axis_xy[1]) - num // 2) * width_pixel + + total = electrons.sum(axis_xy) + mean_x = (electrons * offset_x).sum(axis_xy) / total + mean_y = (electrons * offset_y).sum(axis_xy) / total + var_x = (electrons * np.square(offset_x - mean_x)).sum(axis_xy) / total + var_y = (electrons * np.square(offset_y - mean_y)).sum(axis_xy) / total + std_measured = np.sqrt((var_x + var_y) / 2) + + std_expected = optika.sensors.charge_diffusion( + absorption=absorption, + thickness_substrate=thickness_substrate, + thickness_depletion=thickness_depletion, + ) + + assert np.allclose(std_measured, std_expected, rtol=0.05) From 151a2163214e8beb3902f70c2b01dfb76804e38b Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 9 Jun 2026 10:43:13 -0600 Subject: [PATCH 13/19] Refract incident rays into the substrate inside the sensor's `collect` method. The per-pixel charge-collection model needs the cosine of the refracted angle inside the light-sensitive region, which depends on the ambient index of refraction. Rather than threading an ambient-index argument through `expose`/`signal`, `collect` now refracts each ray with its own `rays.n` via the new `AbstractSensorMaterial.direction_refracted` method and flux-weights the resulting (complex) substrate-side cosine per pixel. `material.signal` consumes this already-refracted cosine and forwards it to `signal` with equal ambient/substrate indices so the internal Snell step is a no-op while the effective absorption still uses the correct substrate index; its `n` parameter is removed. `AbstractBackIlluminatedSiliconSensorMaterial.electrons_measured` is refactored to reuse `direction_refracted`. Co-Authored-By: Claude Opus 4.8 --- optika/sensors/_sensors.py | 34 +++++++-- optika/sensors/materials/_materials.py | 96 +++++++++++++++++++++----- 2 files changed, 105 insertions(+), 25 deletions(-) diff --git a/optika/sensors/_sensors.py b/optika/sensors/_sensors.py index 38262cbc..bd4ed762 100644 --- a/optika/sensors/_sensors.py +++ b/optika/sensors/_sensors.py @@ -95,9 +95,12 @@ def collect( 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 incidence angle in each pixel: the two - quantities :meth:`expose` needs. Performing the ray-to-cosine projection - here is what lets :meth:`expose` be shared with systems that have no rays. + 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 ---------- @@ -113,7 +116,16 @@ def collect( where = where & rays.unvignetted normal = self.sag.normal(rays.position) - direction = -(rays.direction @ normal) + + # 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, + ) flux = rays.intensity * where @@ -132,12 +144,20 @@ def collect( ) image = na.histogram(a, bins=bins, axis=axis, weights=flux) - moment = na.histogram(a, bins=bins, axis=axis, weights=flux * direction) + 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 cosine; empty pixels carry no flux, assume normal incidence + # 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 = np.where(nonempty, moment.outputs / image.outputs, 1) + 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 diff --git a/optika/sensors/materials/_materials.py b/optika/sensors/materials/_materials.py index bd3af815..1992074e 100644 --- a/optika/sensors/materials/_materials.py +++ b/optika/sensors/materials/_materials.py @@ -1441,6 +1441,36 @@ def signal( Whether to add noise to the result. """ + @abc.abstractmethod + def direction_refracted( + self, + wavelength: u.Quantity | na.AbstractScalar, + direction: None | na.AbstractCartesian3dVectorArray = None, + n: complex | na.AbstractScalar = 1, + normal: None | na.AbstractCartesian3dVectorArray = None, + ) -> complex | na.AbstractScalar: + """ + The cosine of the refracted propagation angle inside the light-sensitive + region of the sensor. + + This is the quantity :meth:`signal` expects as its ``direction`` + argument. Performing the refraction here lets the sensor's ``collect`` + method fold the ambient index of refraction into the per-pixel cosine, + so :meth:`signal` does not need a separate ambient-index argument. + + Parameters + ---------- + wavelength + The wavelength of the incident light in vacuum. + direction + The propagation direction of the incident light in the ambient medium. + If :obj:`None` (default), normal incidence (:math:`\\hat{z}`) is assumed. + n + The complex index of refraction of the ambient medium. + normal + The vector perpendicular to the surface of the sensor. + """ + @abc.abstractmethod def photons_incident( self, @@ -1503,6 +1533,21 @@ def signal( return electrons + def direction_refracted( + self, + wavelength: u.Quantity | na.AbstractScalar, + direction: None | na.AbstractCartesian3dVectorArray = None, + n: complex | na.AbstractScalar = 1, + normal: None | na.AbstractCartesian3dVectorArray = None, + ) -> complex | na.AbstractScalar: + # An ideal material does not refract, so the refracted cosine is just + # the cosine of the incidence angle. + if direction is None: + direction = na.Cartesian3dVectorArray(0, 0, 1) + if normal is None: + normal = na.Cartesian3dVectorArray(0, 0, -1) + return -direction @ normal + def photons_incident( self, electrons: u.Quantity | na.AbstractScalar, @@ -1897,6 +1942,25 @@ def probability_measurement( ), ) + def direction_refracted( + self, + wavelength: u.Quantity | na.AbstractScalar, + direction: None | na.AbstractCartesian3dVectorArray = None, + n: complex | na.AbstractScalar = 1, + normal: None | na.AbstractCartesian3dVectorArray = None, + ) -> complex | na.AbstractScalar: + if direction is None: + direction = na.Cartesian3dVectorArray(0, 0, 1) + + if normal is None: + normal = na.Cartesian3dVectorArray(0, 0, -1) + + return optika.materials.snells_law_scalar( + cos_incidence=-direction @ normal, + index_refraction=n, + index_refraction_new=self._chemical.n(wavelength), + ) + def electrons_measured( self, photons_absorbed: na.AbstractScalar, @@ -1915,20 +1979,11 @@ def electrons_measured( absorbed photons using :func:`electrons_measured`. """ - if direction is None: - direction = na.Cartesian3dVectorArray(0, 0, 1) - - if normal is None: - normal = na.Cartesian3dVectorArray(0, 0, -1) - - direction = -direction @ normal - - n_substrate = self._chemical.n(wavelength) - - direction_substrate = optika.materials.snells_law_scalar( - cos_incidence=direction, - index_refraction=n, - index_refraction_new=n_substrate, + direction_substrate = self.direction_refracted( + wavelength=wavelength, + direction=direction, + n=n, + normal=normal, ) return electrons_measured( @@ -1936,7 +1991,7 @@ def electrons_measured( wavelength=wavelength, absorption=absorption_effective( wavelength=wavelength, - n_substrate=n_substrate, + n_substrate=self._chemical.n(wavelength), direction_substrate=direction_substrate, ), thickness_implant=self.thickness_implant, @@ -1953,7 +2008,6 @@ def signal( photons: u.Quantity | na.AbstractScalar, wavelength: u.Quantity | na.AbstractScalar, direction: float | na.AbstractScalar = 1, - n: complex | na.AbstractScalar = 1, width_pixel: ( u.Quantity | na.AbstractScalar | na.AbstractCartesian2dVectorArray ) = 0 @@ -1972,12 +2026,18 @@ def signal( else: method = "expected" + # `direction` is the cosine of the refracted angle *inside* the substrate + # (computed by the sensor's `collect` method), so use equal ambient and + # substrate indices to make the Snell refraction inside `signal` a no-op + # while still computing the effective absorption with the correct index. + n_substrate = self._chemical.n(wavelength) + return signal( photons_expected=photons, wavelength=wavelength, direction=direction, - n=n, - n_substrate=self._chemical.n(wavelength), + n=n_substrate, + n_substrate=n_substrate, absorbance=1, thickness_implant=self.thickness_implant, thickness_depletion=self.depletion.thickness, From 4a7e084441a82c4f0e36941699a140287e3379c3 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 9 Jun 2026 18:37:32 -0600 Subject: [PATCH 14/19] Make off-grid charge diffusion lossy by default, with an opt-in toroidal `wrap`. The Monte-Carlo charge-diffusion kernel previously wrapped charge that diffused past the edge of the pixel grid back onto the opposite edge (a toroidal boundary), which is unphysical for a finite sensor. By default, charge that diffuses off the grid is now lost, as it would be at the physical edge of a sensor. The old toroidal behavior is still available via a new `wrap` keyword argument, threaded from the kernel up through `electrons_measured`, `signal`, and the sensor-material `signal`/`electrons_measured` methods (but not the sensor-level `expose`/`measure`, which use the default). Adds a test covering both boundary behaviors. Co-Authored-By: Claude Opus 4.8 --- optika/sensors/materials/_materials.py | 20 ++++++++++++ .../_ramanathan_2020/_ramanathan_2020.py | 21 ++++++++++++- .../_ramanathan_2020/_ramanathan_2020_test.py | 31 +++++++++++++++++++ 3 files changed, 71 insertions(+), 1 deletion(-) diff --git a/optika/sensors/materials/_materials.py b/optika/sensors/materials/_materials.py index 1992074e..03ea3a81 100644 --- a/optika/sensors/materials/_materials.py +++ b/optika/sensors/materials/_materials.py @@ -1054,6 +1054,7 @@ def signal( temperature: u.Quantity | na.ScalarArray = 300 * u.K, method: Literal["monte-carlo", "expected"] = "monte-carlo", axis_xy: None | tuple[str, str] = None, + wrap: bool = False, shape_random: None | dict[str, int] = None, ) -> na.AbstractScalar: r""" @@ -1115,6 +1116,12 @@ def signal( The two logical axes corresponding to the pixel grid of the sensor along which electrons will diffuse. If :obj:`None` (the default), there is no charge diffusion. + wrap + Controls how diffused charge is treated at the edges of the pixel grid. + If :obj:`False` (the default), charge that diffuses past the edge of the + grid is lost, as it would be at the physical edge of a sensor. + If :obj:`True`, the grid is treated as periodic and the charge re-enters + the opposite edge (a toroidal boundary). shape_random Additional shape used to specify the number of samples to draw. @@ -1216,6 +1223,7 @@ def signal( cce_backsurface=cce_backsurface, temperature=temperature, axis_xy=axis_xy, + wrap=wrap, shape_random=shape_random, ) @@ -1418,6 +1426,7 @@ def signal( * u.um, axis_xy: None | tuple[str, str] = None, noise: bool = True, + wrap: bool = False, ) -> na.AbstractScalar: """ Given the photons incident on each pixel, compute the number of @@ -1439,6 +1448,12 @@ def signal( If :obj:`None` (the default), no diffusion is performed. noise Whether to add noise to the result. + wrap + Controls how diffused charge is treated at the edges of the pixel grid. + If :obj:`False` (the default), charge that diffuses past the edge of + the grid is lost, as it would be at the physical edge of a sensor. + If :obj:`True`, the grid is treated as periodic and the charge + re-enters the opposite edge (a toroidal boundary). """ @abc.abstractmethod @@ -1518,6 +1533,7 @@ def signal( * u.um, axis_xy: None | tuple[str, str] = None, noise: bool = True, + wrap: bool = False, ) -> na.AbstractScalar: if not photons.unit.is_equivalent(u.photon): @@ -1973,6 +1989,7 @@ def electrons_measured( ) = 0 * u.um, axis_xy: None | tuple[str, str] = None, + wrap: bool = False, ) -> na.AbstractScalar: """ Randomly sample the number of measured electrons given the number of @@ -2000,6 +2017,7 @@ def electrons_measured( width_pixel=width_pixel, cce_backsurface=self.cce_backsurface, axis_xy=axis_xy, + wrap=wrap, temperature=self.temperature, ) @@ -2014,6 +2032,7 @@ def signal( * u.um, axis_xy: None | tuple[str, str] = None, noise: bool = True, + wrap: bool = False, ) -> na.AbstractScalar: if not photons.unit.is_equivalent(u.photon): @@ -2047,6 +2066,7 @@ def signal( temperature=self.temperature, method=method, axis_xy=axis_xy, + wrap=wrap, ) def photons_incident( diff --git a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py index 37b75fda..e99a9720 100644 --- a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py +++ b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py @@ -478,6 +478,7 @@ def electrons_measured( cce_backsurface: u.Quantity | na.AbstractScalar = _cce_backsurface, temperature: u.Quantity | na.ScalarArray = 300 * u.K, axis_xy: None | tuple[str, str] = None, + wrap: bool = False, shape_random: None | dict[str, int] = None, ) -> na.AbstractScalar: r""" @@ -525,6 +526,12 @@ def electrons_measured( The two logical axes corresponding to the pixel grid of the sensor along which electrons will diffuse. If :obj:`None` (the default), there is no charge diffusion. + wrap + Controls how diffused charge is treated at the edges of the pixel grid. + If :obj:`False` (the default), charge that diffuses past the edge of the + grid is lost, as it would be at the physical edge of a sensor. + If :obj:`True`, the grid is treated as periodic and the charge re-enters + the opposite edge (a toroidal boundary). shape_random Additional shape used to specify the number of samples to draw. @@ -669,6 +676,7 @@ def electrons_measured( n=n.ndarray, energy_pair_inf=energy_pair_inf.ndarray, fano_inf=fano_inf.ndarray, + wrap=wrap, ) result = na.ScalarArray( @@ -696,6 +704,7 @@ def _electrons_measured_quantity( n: np.ndarray, energy_pair_inf: u.Quantity, fano_inf: u.Quantity, + wrap: bool, ) -> u.Quantity: shape = np.broadcast_shapes( @@ -740,6 +749,7 @@ def _electrons_measured_quantity( n=n.reshape(-1, num_x, num_y, n.shape[~0]), energy_pair_inf=energy_pair_inf.reshape(-1, num_x, num_y), fano_inf=fano_inf.reshape(-1, num_x, num_y), + wrap=wrap, ) result = result.reshape(shape) @@ -768,6 +778,7 @@ def _electrons_measured_numba( # pragma: nocover n: np.ndarray, energy_pair_inf: np.ndarray, fano_inf: np.ndarray, + wrap: bool, ) -> np.ndarray: num_i, num_x, num_y, num_n = p_n.shape @@ -845,6 +856,14 @@ def _electrons_measured_numba( # pragma: nocover else: p = q = 0 - result[i, (x + p) % num_x, (y + q) % num_y] += 1 + x_e = x + p + y_e = y + q + + if wrap: + # toroidal boundary: charge re-enters the opposite edge + result[i, x_e % num_x, y_e % num_y] += 1 + elif (0 <= x_e < num_x) and (0 <= y_e < num_y): + result[i, x_e, y_e] += 1 + # otherwise the electron diffused off the sensor and is lost return result diff --git a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020_test.py b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020_test.py index 8b099f08..0dda2fc0 100644 --- a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020_test.py +++ b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020_test.py @@ -229,3 +229,34 @@ def test_electrons_measured_diffusion(): ) assert np.allclose(std_measured, std_expected, rtol=0.05) + + +def test_electrons_measured_wrap(): + """ + On a grid small compared to the diffusion width, charge that diffuses off + the grid is lost when ``wrap=False`` but retained (toroidally) when + ``wrap=True``, so the wrapped grid holds strictly more charge. + """ + num = 3 + axis_xy = ("pixel_x", "pixel_y") + + photons = np.zeros((num, num)) + photons[num // 2, num // 2] = 5000 + photons = na.ScalarArray(photons << u.photon, axes=axis_xy).astype(int) + + kwargs = dict( + photons_absorbed=photons, + wavelength=500 * u.nm, + absorption=1 / u.um, + thickness_implant=0 * u.um, + thickness_depletion=0 * u.um, + thickness_substrate=14 * u.um, + width_pixel=2 * u.um, + cce_backsurface=1, + axis_xy=axis_xy, + ) + + total_drop = _ramanathan_2020.electrons_measured(**kwargs, wrap=False).sum(axis_xy) + total_wrap = _ramanathan_2020.electrons_measured(**kwargs, wrap=True).sum(axis_xy) + + assert total_wrap > total_drop From 423461cda7c89582989ffad988d5f56294fb9f35 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 9 Jun 2026 18:43:54 -0600 Subject: [PATCH 15/19] Guard the charge-collection kernel against a non-absorbing medium. When the absorption coefficient is zero, the reciprocal absorption length `1/a` and the truncated-exponential depth sampling are undefined and previously produced NaN depths. The kernel now samples the absorption depth from the uniform distribution over the substrate, which is the `a -> 0` limit of the truncated exponential, so a direct call at a non-absorbing wavelength returns finite results. Co-Authored-By: Claude Opus 4.8 --- .../materials/_ramanathan_2020/_ramanathan_2020.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py index e99a9720..b7432f33 100644 --- a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py +++ b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py @@ -802,7 +802,9 @@ def _electrons_measured_numba( # pragma: nocover wp_x = width_pixel_x[i, x, y] wp_y = width_pixel_y[i, x, y] - d = 1 / a + # guard against a non-absorbing medium (a == 0), where the + # reciprocal absorption length is undefined + d = 1 / a if a > 0 else 0.0 # every photon is absorbed within the substrate, so sample the # depth from the exponential truncated to [0, z_substrate) @@ -831,7 +833,12 @@ def _electrons_measured_numba( # pragma: nocover n_ij = round(n_ij) y_ij = random.uniform(0, 1) - z_ij = -d * math.log(1 - y_ij * fraction_absorbed) + if a > 0: + z_ij = -d * math.log(1 - y_ij * fraction_absorbed) + else: + # with no absorption the truncated-exponential depth + # distribution reduces to a uniform one over [0, z_substrate) + z_ij = y_ij * z_substrate if z_ij < W: h_ij = h_0 + (1 - h_0) * z_ij / W From d760bddabec513566f51be1ec250c64ad1ebf0a9 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 9 Jun 2026 19:02:27 -0600 Subject: [PATCH 16/19] Make sensor docstrings consistent with the refraction and signal changes. After moving the refraction into `collect`, `direction` carries the cosine of the refracted angle inside the light-sensitive region, not the incidence angle. Update the `expose` and `AbstractSensorMaterial.signal` docstrings accordingly, and drop the stale "approximate" from the standalone `signal` summary (its monte-carlo method simulates every photon, matching `electrons_measured`). Co-Authored-By: Claude Opus 4.8 --- optika/sensors/_sensors.py | 5 +++-- optika/sensors/materials/_materials.py | 5 +++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/optika/sensors/_sensors.py b/optika/sensors/_sensors.py index bd4ed762..67cd83e7 100644 --- a/optika/sensors/_sensors.py +++ b/optika/sensors/_sensors.py @@ -179,7 +179,7 @@ def expose( 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 an incidence-cosine map), + 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. @@ -196,7 +196,8 @@ def expose( The wavelength inputs (``image.inputs.wavelength``) must be the bin *edges*, not the centers. direction - The cosine of the incidence angle in each pixel. + 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 diff --git a/optika/sensors/materials/_materials.py b/optika/sensors/materials/_materials.py index 03ea3a81..3753318e 100644 --- a/optika/sensors/materials/_materials.py +++ b/optika/sensors/materials/_materials.py @@ -1058,7 +1058,7 @@ def signal( shape_random: None | dict[str, int] = None, ) -> na.AbstractScalar: r""" - A random sample from the approximate distribution of measured electrons + A random sample from the distribution of measured electrons given the expected number of photons incident on the front surface of the sensor. @@ -1439,7 +1439,8 @@ def signal( wavelength An assumed grid of wavelengths for the incident photons. direction - The cosine of the incidence angle. + The cosine of the refracted angle inside the light-sensitive region, + as produced by :meth:`direction_refracted`. width_pixel The physical size of each pixel, used by the charge-diffusion model. axis_xy From 3bcd3417165f8201750012bef78243548c1bea6b Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 9 Jun 2026 19:40:44 -0600 Subject: [PATCH 17/19] Cover the `direction_refracted` default-argument branches. `collect` always passes an explicit direction and normal, so the normal-incidence defaults in `direction_refracted` were never exercised, failing patch coverage. Add `test_direction_refracted` to the shared sensor material test base so the defaults run for every material type. Co-Authored-By: Claude Opus 4.8 --- optika/sensors/materials/_materials_test.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/optika/sensors/materials/_materials_test.py b/optika/sensors/materials/_materials_test.py index 856ba537..4367359f 100644 --- a/optika/sensors/materials/_materials_test.py +++ b/optika/sensors/materials/_materials_test.py @@ -417,6 +417,23 @@ def test_photons_incident( assert isinstance(na.as_named_array(result), na.AbstractScalar) assert result.unit.is_equivalent(u.photon) + @pytest.mark.parametrize( + argnames="wavelength", + argvalues=[ + 100 * u.AA, + ], + ) + def test_direction_refracted( + self, + a: optika.sensors.materials.AbstractSensorMaterial, + wavelength: u.Quantity | na.AbstractScalar, + ): + # with the default direction and normal (normal incidence) the + # refracted cosine is unity + result = a.direction_refracted(wavelength=wavelength) + assert isinstance(na.as_named_array(result), na.AbstractScalar) + assert np.all(np.real(result) > 0) + @pytest.mark.parametrize( argnames="a", From e269beccbd05c183fd6a0ec966431d11b054be90 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 9 Jun 2026 20:30:27 -0600 Subject: [PATCH 18/19] fix ordering --- optika/sensors/__init__.py | 4 ++-- optika/sensors/materials/_materials.py | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/optika/sensors/__init__.py b/optika/sensors/__init__.py index b4139dd7..d5d89bba 100644 --- a/optika/sensors/__init__.py +++ b/optika/sensors/__init__.py @@ -16,8 +16,8 @@ fano_factor, fano_factor_inf, transmittance, - absorption_effective, absorbance, + absorption_effective, charge_collection_efficiency, quantum_efficiency_effective, probability_measurement, @@ -43,8 +43,8 @@ "fano_factor", "fano_factor_inf", "transmittance", - "absorption_effective", "absorbance", + "absorption_effective", "charge_collection_efficiency", "quantum_efficiency_effective", "probability_measurement", diff --git a/optika/sensors/materials/_materials.py b/optika/sensors/materials/_materials.py index 3753318e..770905b8 100644 --- a/optika/sensors/materials/_materials.py +++ b/optika/sensors/materials/_materials.py @@ -35,6 +35,7 @@ "fano_factor", "fano_factor_inf", "absorption_effective", + "transmittance", "absorbance", "charge_collection_efficiency", "quantum_efficiency_effective", From 092eaa3af7ae07b8bcda918222eedb3322f4e41e Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 9 Jun 2026 20:32:38 -0600 Subject: [PATCH 19/19] doc tweaks --- optika/sensors/_sensors.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/optika/sensors/_sensors.py b/optika/sensors/_sensors.py index 67cd83e7..43586b3a 100644 --- a/optika/sensors/_sensors.py +++ b/optika/sensors/_sensors.py @@ -99,7 +99,8 @@ def collect( 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, + material's :meth:`~optika.sensors.materials.AbstractSensorMaterial.signal` + model be shared with systems that have no rays, without threading a separate ambient-index argument through them. Parameters