From 2f7719a3013d26032d0d1b19523ca98645107166 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Fri, 28 Aug 2026 22:52:42 -0600 Subject: [PATCH 1/2] Integrate the charge-diffusion correction by Gauss-Legendre quadrature The average in the charge-diffusion term of `vmr_signal` has no closed form. The obstruction is the coincidence probability, which contains erf(1 / 2 sigma) with sigma proportional to sqrt(z_f - z); squaring it to cover the two axes of the sensor gives a product of two error functions integrated against the exponential absorption profile. Sympy returns NonElementaryIntegral for the underlying antiderivative, while happily evaluating the same integrand without the reciprocal-argument exponential, so the obstruction is that combination specifically. Routing the calculation through the separation instead gives an incomplete Bessel integral, which is no better and is not in scipy. It was previously evaluated with a 1001-node midpoint rule in the cumulative absorption probability. That converges only algebraically, because sigma(z) puts a square-root branch point at the edge of the field-free region, and it reached only 2.2e-4 across 1 to 10000 angstroms for the e2v CCD97. Substituting s = 1 - r^2 removes that branch point while keeping the nodes distributed by where photons are actually absorbed, which matters because the optical depth of the field-free region spans 0.02 to 1271 over that wavelength range for the same sensor. The interval is split where the implant ends and the differential CCE saturates, and Gauss-Legendre is applied to each part. The optical depth is evaluated as -log(exp(-a z_f) + (1 - exp(-a z_f)) r^2) rather than through s, so it stays finite once exp(-a z_f) underflows, which it does beyond an optical depth of about 745. With 32 nodes per part the worst error over the same wavelength range is 1.6e-7, against 2.2e-4 from 1001 midpoints: three orders of magnitude more accurate with fifteen times fewer evaluations. Convergence was checked separately against optical depths from 0.02 to 3000 to confirm the node count is not tuned to one sensor. A plain sqrt substitution in depth was tried first and rejected: it is spectrally accurate at low optical depth but needs nodes scaling as sqrt(alpha z_f), and fails at 9e-3 with 16 nodes where a real sensor reaches 1271. Results are unchanged to within the quadrature error, and the no-diffusion path is untouched. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_0192abyvBf2Zw5rq5CQ62JFb --- optika/sensors/materials/_materials.py | 118 ++++++++++++++++++-- optika/sensors/materials/_materials_test.py | 30 +++++ 2 files changed, 136 insertions(+), 12 deletions(-) diff --git a/optika/sensors/materials/_materials.py b/optika/sensors/materials/_materials.py index 94f8099..881425a 100644 --- a/optika/sensors/materials/_materials.py +++ b/optika/sensors/materials/_materials.py @@ -1,4 +1,4 @@ -from typing import Literal +from typing import Callable, Literal from typing_extensions import Self import abc import functools @@ -1268,6 +1268,53 @@ def _probability_same_pixel( return np.where(where, result, 1) +_num_gauss_legendre = 32 +""" +The number of Gauss-Legendre nodes used on each subinterval by +:func:`_integrate_gauss_legendre`. +Chosen so that the charge-diffusion integral of :func:`vmr_signal` is accurate +to better than one part in :math:`10^6` over the full range of optical depths +encountered by a silicon sensor between 1 and 10000 angstroms. +""" + + +def _integrate_gauss_legendre( + integrand: Callable[[na.AbstractScalar], na.AbstractScalar], + lower: float | na.AbstractScalar, + upper: float | na.AbstractScalar, + axis: str, +) -> na.AbstractScalar: + """ + Integrate `integrand` between `lower` and `upper` using Gauss-Legendre + quadrature with :obj:`_num_gauss_legendre` nodes. + + The limits may be arrays, in which case a separate quadrature rule is + applied to every element, and `axis` is the logical axis along which the + nodes are placed. + + Parameters + ---------- + integrand + The function to integrate. + lower + The lower limit of integration. + upper + The upper limit of integration. + axis + The logical axis along which to place the quadrature nodes. + Consumed by the sum, so it does not appear in the result. + """ + nodes, weights = scipy.special.roots_legendre(_num_gauss_legendre) + + nodes = na.ScalarArray(nodes, axes=(axis,)) + weights = na.ScalarArray(weights, axes=(axis,)) + + half = (upper - lower) / 2 + center = (upper + lower) / 2 + + return half * (weights * integrand(half * nodes + center)).sum(axis) + + def vmr_signal( wavelength: u.Quantity | na.ScalarArray, direction: float | na.AbstractScalar = 1, @@ -1535,9 +1582,44 @@ def vmr_signal( F_\text{diffusion} = - \left( \overline{n} + \mathcal{F} - 1 \right) \frac{\left\langle \left[ 1 - D(z) \right] \eta^2(z) \right\rangle}{\left\langle \eta(z) \right\rangle} - from Equation :eq:`vmr-compact`, - where the average is computed numerically over the field-free region - using the absorption-depth distribution truncated to the substrate. + from Equation :eq:`vmr-compact`. + The integrand vanishes outside the field-free region, since :math:`D = 1` + where the charge does not spread, so only that region is integrated. + + The average has no closed form. + The obstruction is :math:`D`, which contains + :math:`\text{erf}(1 / 2 \sigma)` with :math:`\sigma \propto \sqrt{z_f - z}`, + and squaring it to cover the two axes of the sensor produces a product of + two error functions integrated against the exponential absorption profile. + It is therefore evaluated by quadrature, in a variable chosen to make that + quadrature converge quickly. + Writing the average as an integral over the cumulative absorption + probability :math:`s` within the field-free region distributes the nodes + according to where photons are actually absorbed, which matters because the + optical depth of that region spans four orders of magnitude across the + wavelengths of interest. + Substituting :math:`s = 1 - r^2` then removes the square-root branch point + which :math:`\sigma(z)` places at :math:`z = z_f`, since the diffusion width + becomes proportional to :math:`r` near that end of the range. + In terms of :math:`r` the optical depth is + + .. math:: + + \alpha z = -\log \left( e^{-\alpha z_f} + \left( 1 - e^{-\alpha z_f} \right) r^2 \right), + + which is evaluated in this form rather than through :math:`s` so that it + stays finite when :math:`e^{-\alpha z_f}` underflows. + The remaining integrand is smooth apart from the point where the implant + ends and :math:`\eta` saturates, so the interval is split there and + Gauss-Legendre quadrature is applied to each part. + + With :obj:`_num_gauss_legendre` nodes per part this is accurate to better + than one part in :math:`10^6` for optical depths + :math:`\alpha z_f` between :math:`0.02` and :math:`3000`, + which covers silicon between 1 and 10000 angstroms; + a midpoint rule in :math:`s` needs some thirty times as many nodes to reach + a hundred times worse accuracy. + This result assumes uniform illumination and a periodic pixel grid, so it corresponds to ``wrap=True`` in :func:`signal`, and is a good approximation away from the edges of a sensor @@ -1621,19 +1703,31 @@ def ratio_ff(width: u.Quantity | na.AbstractScalar) -> na.AbstractScalar: r_y = ratio_ff(width_pixel.y) fraction_absorbed = -np.expm1(-az_substrate) - t_ff = -np.expm1(-az_ff) / fraction_absorbed + fraction_ff = -np.expm1(-az_ff) + t_ff = fraction_ff / fraction_absorbed az_ff_safe = np.where(az_ff > 0, az_ff, 1) axis_z = "_vmr_signal_depth" - num_z = 1001 - s = (na.arange(0, num_z, axis=axis_z) + 0.5) / num_z - az = -np.log1p(-t_ff * s * fraction_absorbed) - eta = np.minimum(n0 + (1 - n0) * az / aW, 1) - w = np.sqrt(np.maximum(1 - az / az_ff_safe, 0)) - D = _probability_same_pixel(r_x * w) * _probability_same_pixel(r_y * w) + def integrand(r: na.AbstractScalar) -> na.AbstractScalar: + az = -np.log(np.exp(-az_ff) + fraction_ff * np.square(r)) + eta = np.minimum(n0 + (1 - n0) * az / aW, 1) + w = np.sqrt(np.maximum(1 - az / az_ff_safe, 0)) + D = _probability_same_pixel(r_x * w) * _probability_same_pixel(r_y * w) + return (1 - D) * np.square(eta) * 2 * r + + # The integrand is smooth except where the implant ends and the + # differential CCE saturates, so integrate up to that point and beyond + # it separately. The breakpoint collapses to zero if the implant is + # thicker than the field-free region, leaving a single interval. + r_implant = np.sqrt( + np.maximum((np.exp(-aW) - np.exp(-az_ff)) / fraction_ff, 0), + ) - integral = t_ff * ((1 - D) * np.square(eta)).mean(axis_z) + integral = t_ff * ( + _integrate_gauss_legendre(integrand, 0, r_implant, axis_z) + + _integrate_gauss_legendre(integrand, r_implant, 1, axis_z) + ) unit = u.electron / u.photon diff --git a/optika/sensors/materials/_materials_test.py b/optika/sensors/materials/_materials_test.py index b432945..355a2f3 100644 --- a/optika/sensors/materials/_materials_test.py +++ b/optika/sensors/materials/_materials_test.py @@ -376,6 +376,36 @@ def test_vmr_signal_diffusion(): assert np.all(result_default == result_default_no_diffusion) +def test_vmr_signal_quadrature(): + """ + The charge-diffusion integral should be converged at the default number of + quadrature nodes, across the full range of optical depths that silicon + spans between 1 and 10000 angstroms. + """ + ccd = optika.sensors.materials.e2v_ccd97() + + kwargs = dict( + wavelength=na.geomspace(1, 10000, axis="wavelength", num=101) * u.AA, + thickness_implant=ccd.thickness_implant, + thickness_depletion=ccd.depletion.thickness, + thickness_substrate=ccd.thickness_substrate, + width_pixel=16 * u.um, + cce_backsurface=ccd.cce_backsurface, + temperature=ccd.temperature, + ) + + result = optika.sensors.vmr_signal(**kwargs) + + num = optika.sensors.materials._materials._num_gauss_legendre + try: + optika.sensors.materials._materials._num_gauss_legendre = 8 * num + expected = optika.sensors.vmr_signal(**kwargs) + finally: + optika.sensors.materials._materials._num_gauss_legendre = num + + assert np.all(np.abs(result / expected - 1) < 1e-5) + + class AbstractTestAbstractSensorMaterial( AbstractTestAbstractMaterial, ): From d3bace984f2323eb3d0c5e59d0f439cd37bed359 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Sat, 29 Aug 2026 10:03:02 -0600 Subject: [PATCH 2/2] Record why the quadrature has a breakpoint to split at The split exists only because the differential CCE of Stern 1994 is piecewise linear and therefore kinks where the implant ends. A smooth profile would need no split and half as many nodes, so the split reads like an accident worth optimizing away. It is not: the piecewise-linear form is kept deliberately, so the fitted implant thickness and back-surface CCE stay comparable with the values published by Stern and by Boerner et al. Note that no choice of CCE would produce a closed form anyway, since setting eta = 1 still leaves the erf-squared obstruction. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_0192abyvBf2Zw5rq5CQ62JFb --- optika/sensors/materials/_materials.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/optika/sensors/materials/_materials.py b/optika/sensors/materials/_materials.py index 881425a..7d54008 100644 --- a/optika/sensors/materials/_materials.py +++ b/optika/sensors/materials/_materials.py @@ -1720,6 +1720,14 @@ def integrand(r: na.AbstractScalar) -> na.AbstractScalar: # differential CCE saturates, so integrate up to that point and beyond # it separately. The breakpoint collapses to zero if the implant is # thicker than the field-free region, leaving a single interval. + # + # The kink is inherited from the piecewise-linear differential CCE of + # `charge_collection_efficiency`, which follows :cite:t:`Stern1994`. + # A smooth profile would need no split and half as many nodes, but the + # piecewise-linear form is kept deliberately so that the fitted implant + # thickness and back-surface CCE remain comparable with the values + # published by :cite:t:`Stern1994` and :cite:t:`Boerner2012`. Do not + # trade that away for the quadrature. r_implant = np.sqrt( np.maximum((np.exp(-aW) - np.exp(-az_ff)) / fraction_ff, 0), )