nulltracer

Kerr-Newman geodesic ray tracing rebuilt from a code audit, on three pillars:

1. Mino-time separated formulation. With d/dlam = Sigma d/dsigma the geodesic equations decouple into polynomial, velocity-independent second-order ODEs r'' = R'(r)/2, u'' = U'(u)/2 (u = cos theta). Turning points and poles are ordinary points; no square roots, no sign bookkeeping. Tracing uses the inverted radius w = 1/r, under which the radial equation is form-invariant (Rw(w) = w^4 R(1/w) is the reversed-coefficient quartic) and every dynamical variable is O(1) — this removes a conditioning failure (relative tolerances on dr/dlam ~ E r^2 admit absolute errors that scramble near-critical photons).

2. Adaptive Dormand-Prince 8(7) — vectorised over rays, PI step controller, NaN-hardened rejection, cubic-Hermite event localisation. The tableau is the audit-verified coefficient set, re-verified against order conditions at machine precision in the test suite. No FSAL reuse.

3. Two genuinely symplectic integrators.

   • MinoSymplectic: the Mino Hamiltonian h = |v|^2/2 - (Rw + U)/2 is exactly separable, so explicit kick-drift splitting applies legitimately; composed to order 4 (Yoshida) or 8 (Kahan-Li s15odr8). phi and t ride inside the splitting via conjugate-pair augmentation with zero momenta (exact, order-preserving — not bolted-on quadrature).
   • TaoSymplectic: Tao (2016) explicit extended-phase-space method for the nonseparable Boyer-Lindquist Hamiltonian — phase space is doubled to (q, p, x, y) with an omega-binding; flows A, B, C are exact; symplecticity of the composed map is verified numerically via J^T Omega J = Omega in the tests.

Disk physics (disk.py) replaces the audited-wrong chain: circular-orbit kinematics from the metric-derivative quadratic (exact for Kerr-Newman), Page-Thorne flux from the defining integral and the independent 1974 closed form (required to agree by the tests), Cunningham plunge inside the ISCO, full kinematic redshift.

Quick start

import numpy as np
from nulltracer import KerrNewman, render

bh = KerrNewman(M=1.0, a=0.9)
out = render(bh, n=300, fov=12.0, theta0=np.deg2rad(75), method="adaptive")
# out["image"], out["g"], out["r_hit"], out["status"]

See notebooks/ for shadow/disk imaging, the integrator comparison (constraint boundedness vs RK drift, Tao convergence), and disk physics.

Testing

python -m pytest tests/        # 18 tests: symbolic-verified formulas,
                               # cross-formulation agreement vs scipy DOP853,
                               # analytic benchmarks (3*sqrt(3) shadow,
                               # Kerr photon rings, 4M/b deflection),
                               # empirical integrator orders, symplecticity,
                               # Page-Thorne numeric vs closed form

Every closed-form expression was derived/verified symbolically first (derivations/derive.py).