Coefficients to Tables

How the C++ engine in orbitals.h/pointcloud.h turns three quantum numbers into the sampling tables every point draw consumes — the associated Legendre and Laguerre polynomials, and the inverse-CDF construction that turns each into an O(1) lookup.

Companion

This is the build side of Tables to Point Clouds, whose pipelines start from the OrbitalSampler / OrbitalAngularTables this one produces. Same three-axis split (radial, colatitude, azimuthal); this diagram covers how each axis's density is sampled and inverted.

Building one OrbitalSampler(n, ℓ, m)

src/physics/orbitals.h · pointcloud.h — one run per library entry (36 hydrogen presets)

Given — three quantum numbers

(n, , m),   0 ≤ n−1,   −m
table resolution N = 1001 samples per axis
kOrbitalTableSize

1 · Associated Legendre coefficients

Pm(θ) = (Σk ck uk) · sin|m|θ,   u = cosθ
seed c0 = (−1); downward recurrence i = … |m|+1

Di = √((+1) − i(i−1)); each step folds one degree of the ladder relation into the coefficient array: ck−1 += k·ck/Di, then ck+1 = −(k+2ick/Di — a closed-form polynomial, not a per-angle recursion, so it's built once and evaluated cheaply afterward. The ck−1 update is skipped at k = 0, and wrong-parity entries stay zeroed (c0 is re-zeroed whenever the parity offset kStart = 1).

legendreCoeffs()
c_k (unnormalized)

2 · Normalize to unit peak

max = maxθ∈[0,π/2), step π/100 |Pm(θ)|
ckck / max,   for every k

Sampled at the exact same 50 samples (step π/100) in every port (C++/MicroPython/JS), so all three land on the identical normalization constant instead of three slightly different ones.

legendreCoeffs() / computePLM()
c_k (normalized)

3 · Associated Laguerre coefficients

Rnℓ(r) = (Σk bk rk) · r · er/n,   degree = n−1
b0 = 1;   bk+1 = −2(1 − (+k+1)/n)(k+1)(2+k+2) · bk

Forward recurrence this time, k = 0 up to degree−1 — the polynomial part of the radial wavefunction, expressed directly in n and ℓ rather than the generic Laguerre index.

laguerreCoeffs()
b_k, c_k

4 · Sample the radial density

wr(ri) = [ri · Rnℓ(ri)]²,   ri = i·Δr,   Δr = maxR/(N−1)
maxR = 6n²   (heuristic — covers essentially all of the wavefunction)

Squaring r·R rather than R alone folds in the r² Jacobian of the spherical volume element, so wr is already the correctly-weighted radial density, not just |R|².

buildOrbitalSamplerConstexpr() radial weight loop

5 · Sample the colatitude density

wθ(θi) = Pm(θi)² · sinθi,   θi = i·Δθ,   Δθ = π/(N−1)

sinθ is the colatitude factor of the solid-angle element dΩ = sinθ dθ dφ — without it, points would bunch up near the poles.

buildOrbitalSamplerConstexpr() θ-loop

6 · Sample the azimuthal density

wφ(φi) = Φm(φi)²,   φi = i·Δφ,   Δφ = 2π/(N−1)
Φm(φ) =
cos()   m ≥ 0
sin(|m|φ)   m < 0

m = 0 makes Φ ≡ 1 everywhere, so wφ is flat — an s orbital's azimuthal symmetry falls out for free, no special case needed.

buildOrbitalSamplerConstexpr() φ-loop
w_r, w_θ, w_φ

7 · Invert each into a quantile table

cumulativei = cumulativei−1 + wi  (running sum, cumulative−1=0)
Ci = cumulativei / cumulativeN−1  (normalized CDF)
F−1(uk) = (j₀ + t·(j₁−j₀))·Δ,   t = uk − CjCj − Cj
for uk = k/(N−1), where Cj < uk ≤ Cj bracket it (j₀ = 0 on the first bin; a flat CDF region collapses t to 0)

One monotonic forward sweep per axis, O(N) — the CDF only ever increases, and the target quantile only ever increases, so no index is ever revisited. This is the exact mirror of the runtime lookup in Tables to Point Clouds's step 2: there, an index interpolates a value; here, a probability interpolates an index. A degenerate all-zero density (total ≤ 0) is guarded to a uniform table instead of dividing by zero.

buildInverseCdf()
invRTable, invThetaTable, invPhiTable

Deliverable — OrbitalSampler(n, ℓ, m)

{ n, , m, maxR, invRTable[N], invThetaTable[N], invPhiTable[N] }

This is exactly the dashed "Given" box at the top of the hydrogen pipeline in Tables to Point Clouds.

buildOrbitalSamplerConstexpr()
Reduced variant

Angular-only tables, reused by the atom model

src/physics/pointcloud.h · angular_library.h — one run per (ℓ, m), independent of n and Z

Steps 1, 2, 5, 6, 7 only

Legendre coefficients → normalize → sample wθ, wφ → invert both

Steps 3–4 (the Laguerre/radial half) are skipped entirely: an orbital's angular shape never depends on n or the effective nuclear charge Zeff, only its radial extent does (see the atom pipeline's own r → Zeff·r substitution) — so the θ and φ tables alone are reusable across every element and every n that shares an (ℓ, m).

buildAngularTablesConstexpr()

Deliverable — OrbitalAngularTables(ℓ, m)

{ invThetaTable[N], invPhiTable[N] }

16 entries, ℓ ≤ 3 — one of the "Given" inputs in the atom pipeline in Tables to Point Clouds.

kAngularLibrary
Where this actually executes

16 angular tables

Run inside the C++ compiler itself, building the firmware image

Genuinely compile-time: kAngularLibrary is a constexpr array, baked straight into .rodata — zero runtime cost to build, ~128KB of flash.

36 full hydrogen samplers

Run host-side, in Python, when the firmware image is packaged

Same math, ported to micropython/orbitals.py + pointcloud.py — kept off the compiled image (would cost ~423KB) and written to data/orbital_samplers.bin instead, flashed to SPIFFS, read on demand at runtime. Regenerate with tools/orbital_table_gen.py whenever kOrbitalLibrary changes.

Bit-identical either way — the C++, MicroPython, and JS implementations are cross-checked against each other by tools/orbitals_host/run_crosscheck.sh.

tools/orbital_table_gen.py