Integrate the charge-diffusion correction by Gauss-Legendre quadrature

[roytsmart]

The charge-diffusion average in vmr_signal was evaluated with a 1001-node midpoint rule. This replaces it with 64 Gauss-Legendre nodes in a better variable: three orders of magnitude more accurate with fifteen times fewer evaluations.

Why there is no closed form

The obstruction is the coincidence probability $\mathcal{P}$, which contains $\operatorname{erf}(1/2\sigma)$ with $\sigma \propto \sqrt{z_f - z}$. Squaring it to cover the two axes gives a product of two error functions integrated against the exponential absorption profile. Substituting $x = 1/2\sigma$ reduces everything to three families, and sympy's verdict is:

integrand	result
$x^{-3}e^{c/x^2}$, $x^{-5}e^{c/x^2}$	closed form
$x^{-3}e^{-x^2}e^{c/x^2}$	NonElementaryIntegral
$\operatorname{erf}(x),x^{-3}e^{c/x^2}$	unevaluated
$\operatorname{erf}(x)^2x^{-3}e^{c/x^2}$	unevaluated
control: $\operatorname{erf}(x)^2e^{-x^2}$	closed form, $\sqrt{\pi}\operatorname{erf}(x)^3/6$

The control row matters: sympy handles $\operatorname{erf}^2$ integrals fine in general, so the obstruction is specifically its combination with the reciprocal-argument exponential. Integrating over the separation instead yields an incomplete Bessel integral, which is no better and is not in scipy.

Why the old scheme was slow to converge

$\sigma(z) \propto \sqrt{z_f - z}$ puts a square-root branch point at the edge of the field-free region, so the midpoint rule converges only algebraically. Over 1–10000 Å for the e2v CCD97 it reached only 2.2e-4 with 1001 nodes.

The substitution

Writing the average over the cumulative absorption probability $s$ keeps nodes where photons are actually absorbed — necessary, since $\alpha z_f$ ranges over 0.02 to 1271 across that wavelength range for one sensor. Then $s = 1-r^2$ removes the branch point, since the diffusion width becomes linear in $r$ there. In $r$,

$$\alpha z = -\log!\big(e^{-\alpha z_f} + (1 - e^{-\alpha z_f}),r^2\big)$$

evaluated in this form rather than through $s$ so it stays finite once $e^{-\alpha z_f}$ underflows, which happens beyond optical depth ~745 — well inside the range a real sensor spans. The interval is split where the implant ends and $\eta$ saturates, and Gauss-Legendre is applied to each part.

Accuracy

Against a 400-node reference, over 1–10000 Å:

nodes/segment	total	max rel. error
8	16	4.1e-05
16	32	2.8e-06
32	64	1.6e-07
48	96	2.7e-08

versus 2.2e-04 from the previous 1001 midpoints. Convergence was also checked directly against optical depths from 0.02 to 3000, so the node count is not tuned to a single sensor.

A rejected alternative

Substituting $t = \sqrt{1-z/z_f}$ in depth is spectrally accurate at low optical depth (1e-13 with 48 nodes on a toy case) but needs nodes scaling as $\sqrt{\alpha z_f}$, and errs by 9e-3 with 16 nodes at $\alpha z_f = 500$ — while real sensors reach 1271. The CDF-based version is uniformly robust because it adapts to the exponential.

Checks

199 tests pass, including the existing Monte Carlo comparison. New test_vmr_signal_quadrature asserts the default node count is converged to 1e-5 against an 8× denser rule across the full wavelength range. The no-diffusion path is bit-identical. black, ruff, and the docs build are clean, and the new Notes text was verified to parse into the doctree with no stray markup.

🤖 Generated with Claude Code

https://claude.ai/code/session_0192abyvBf2Zw5rq5CQ62JFb

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 100.00%. Comparing base (472fcd1) to head (d3bace9).

Additional details and impacted files

@@            Coverage Diff            @@
##              main      #217   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files          121       121           
  Lines         7424      7445   +21     
=========================================
+ Hits          7424      7445   +21     

Flag	Coverage Δ
unittests	100.00% <100.00%> (ø)

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.