Fast analytic alpha-beta engine for MetSimErosion

[AstroEloy]

Summary

• Adds wmpl/MetSim/MetSimErosionAlphaBeta.py: a tolerance-matched, drop-in-compatible replacement for MetSimErosion.py's per-timestep RK4 simulation, using the closed-form alpha-beta ablation solution as a fast analytic propagator inside the same erosion/disruption/complex-fragmentation framework — not a different model, just a fast engine for the same physics.
• Full feature parity with the reference: continuous two-phase erosion with grain spawning, compressive-strength disruption, and all five complex-fragmentation entry types (M/F/EF/D/A), including full retroactive re-splitting of already-alive daughters/grains for "A" events.
• Real, measured speedups up to  21x across representative scenarios (see table below), from closed-form main/daughter dynamics (no per-tick RK4 stepping at all for those), a new Cython-accelerated grain-population stepper, and several profiling-driven fixes (shared spline construction, O(N) aggregation, elimination of redundant spline queries, skipping full Fragment object construction for grain spawning).
• Adds a new companion Cython module, wmpl/MetSim/MetSimErosionAlphaBetaCyTools.pyx, for the grain-population stepping hot path.
• Adds wmpl/MetSim/Tests/test_MetSimErosionAlphaBeta.py (49 tests), validating both the underlying math and end-to-end fidelity against the real RK4 engine across every mechanism.
• wmpl/Utils/AlphaBeta.py and wmpl/MetSim/MetSimErosion.py are unchanged — used throughout only as read-only reference/ground truth. Wiring wmpl/Dynesty/DynestyMetSim.py to actually use this new engine is intentionally not part of this PR.

Motivation

Every single log-likelihood evaluation in DynestyMetSim.py's dynesty-based Bayesian fits runs one complete RK4 simulation of MetSimErosion.py. Per wmpl/Dynesty/README.md, this makes a full fit cost 1–2 days on 96 cores per meteor (up to 6 days for fireballs) — one logged run took 2.88M likelihood calls over ~60 hours. The bottleneck isn't the RK4 arithmetic itself (already Cython); it's the outer Python-level loop cost of stepping every fragment every 5ms, and of re-binning grain populations from scratch every step during erosion. This engine replaces that per-step stepping with closed-form per-segment evaluation for the main body and every daughter, collapsing the dominant cost from O(timesteps × fragments) down to O(segments), while keeping the RK4-numerics-matching behavior for the parts (small grains) where the reference's own coarse-dt behavior is deliberately reproduced rather than "corrected."

Architecture

• Segments, not timesteps. Each fragment (main body, disruption daughter, complex-fragmentation daughter) is represented as one or more segments of constant (K, sigma_eff), each with a closed-form v(t)/h(t)/m(t) solution, built once and queried many times.
• Two corrections beyond AlphaBeta.py's existing formulas were needed to make the closed form exactly reproduce this model's own ODEs: beta = sigma_eff*v_start**2/6 (not /2 — a mu=2/3 self-similar-shrinking-body correction, not the naive mu=0 you'd get by pattern-matching the mass formula alone), and a finite-starting-height generalization of the standard (vacuum-start) closed form, needed for any fragment spawned mid-flight.
• Atmosphere reconciliation: the closed form assumes an exponential atmosphere; the reference uses a real polynomial density fit. Reconciled via a cumulative-column-density-matched equivalent-height map (AtmEquivHeightMap), built once per simulation and reused for every fragment.
• Erosion: replaces per-5ms grain re-binning with epochs spaced uniformly in cumulative eroded mass along a segment (closed-form, no per-step tracking), with epoch budget allocated proportionally to each segment's own share of the total eroded mass across the whole simulation.
• Grains: two evolution modes. Default (grain_evolution_analytic=False) steps grains via RK4, deliberately reproducing the reference's own coarse-dt numerical behavior (confirmed the reference itself is under-resolved for small grains — matching it, not "fixing" it, is the point of a tolerance-matched substitute). An opt-in grain_evolution_analytic=True mode uses the true continuous-physics closed form instead, trading reference-numerics-matching for physical exactness. The main body and every daughter fragment are always evaluated via the closed-form alpha-beta segment chain, in both configurations — that's where the speedup comes from, and it is not a toggle. Only the grains differ: the default ("hybrid", grain_evolution_analytic=False) steps grains via RK4, deliberately reproducing the reference's own coarse-dt numerical behavior (confirmed the reference itself is under-resolved for small grains — matching it, not "fixing" it, is the point of a tolerance-matched substitute). The opt-in "100% alpha-beta" mode (grain_evolution_analytic=True) evolves grains via the true continuous-physics closed form instead — no RK4 stepping anywhere in the engine in this mode — trading reference-numerics-matching for physical exactness. These are the two configurations named in the performance table below.
• Disruption & complex fragmentation: tick-exact dynamic-pressure trigger for disruption; all five fragmentation entry types implemented, including full retroactive re-splitting for type "A" (which overwrites sigma/gamma for every already-alive fragment, not just future ones) — genuinely the hardest piece, requiring a two-phase (segment re-split, then grain-batch tick-level re-splice) implementation.
• A near-singular mass-blowup-near-death artifact in the closed form (a real, characterized numerical regime, not a bug) is handled via automatic detection and a short exact-RK4 "tail" splice for just the final few ticks of an affected segment — a hybrid, not a wholesale RK4 fallback, so it doesn't reintroduce the cost this engine exists to eliminate.

Accuracy

Tolerance-matched, not bit-identical, against MetSimErosion.runSimulation() on the same Constants, across single-body, erosion, disruption, and complex-fragmentation scenarios:

• Main/daughter dynamics: sub-meter height error, <0.01% velocity error over the first 90% of a fragment's own lifetime, in the well-validated cases.
• Frame-averaged (30fps) light curve: worst-case |Δmag| ranging from ~0.01–0.15 for single-mechanism scenarios up to ~0.75 for the most demanding combined scenario (erosion + all three complex-fragmentation types + disruption in one flight) — the residual there is root-caused to an ordinary peak-timing-alignment sensitivity, not a structural gap.
• Daughter death times match the reference exactly (0.0ms) after fixing a real upstream-sync bug (a disruption daughter's own erosion coefficient wasn't staying fixed for its whole life, matching a fix landed in the reference tool itself).
• A handful of deliberate, explicitly documented deviations from the reference exist (e.g. a single-tick simplification in when a disrupting fragment's daughters start contributing to mass_total_active) — each is a real, small, characterized effect, not an unknown one.

Performance

Best-of-5 wall-clock, same Constants object fed to both engines, measured via the new benchmark_speed_across_scenarios() (see Testing below):

Scenario	RK4 (reference)	Hybrid alpha-beta	Speedup	100% alpha-beta	Speedup
Erosion only	4593 ms	218 ms	21.1x	512 ms	9.0x
Disruption only	11.6 ms	23.5 ms	0.49x	32.4 ms	0.36x
Disruption + erosion	7129 ms	406 ms	17.5x	844 ms	8.4x
Complex (erosion + EF/A/D + disruption)	842 ms	219 ms	3.8x	526 ms	1.6x

"Disruption only" is the one scenario that remains slower than the reference: it's a deliberately tiny scenario (milliseconds either way — the main body disrupts almost immediately into a handful of non-eroding daughters) where this engine's own fixed per-call overhead (building the atmosphere map, constructing per-daughter trajectory splines) outweighs the reference's own near-trivial workload. It's a real, bounded, already-diagnosed gap, not a regression risk for any realistic (erosion-heavy) fit.

Testing

wmpl/MetSim/Tests/test_MetSimErosionAlphaBeta.py, 49 tests, runnable via pytest or standalone:

python -m pytest wmpl/MetSim/Tests/test_MetSimErosionAlphaBeta.py -v
python -m wmpl.MetSim.Tests.test_MetSimErosionAlphaBeta              # asserts only
python -m wmpl.MetSim.Tests.test_MetSimErosionAlphaBeta --plot       # + diagnostic comparison figures
python -m wmpl.MetSim.Tests.test_MetSimErosionAlphaBeta --benchmark  # + the speed table above

Coverage spans three levels:

• Math-level unit tests: the two closed-form corrections (mu=2/3, finite-start), the atmosphere-equivalent-height map, the velocity-inversion spline, gravity-drop/curvature geometry — each checked against an independent from-scratch RK4 mirror, not just against the engine's own other pieces.
• End-to-end accuracy tests: single-body, erosion, disruption, and every complex-fragmentation entry type, each checked against real MetSimErosion.runSimulation() output on matched Constants, plus dedicated tests for results_list aggregation correctness (the vectorized "brightest fragment"/"leading fragment" selection logic) and structural invariants that must hold regardless of exact grain-population alignment between the two engines.
• Regression/bit-exactness tests for internal rewrites: several pieces were rewritten for speed during development (the grain-population stepper fused into Cython, cross-daughter trajectory-spline batching, the mass-binning fast path) — each has a dedicated test comparing the rewrite's output against the original, unoptimized implementation (or an independent scalar reference) to near machine precision, so a future change that accidentally alters behavior while "just" optimizing gets caught immediately.

Two explicit NotImplementedError cases are locked in by their own tests rather than left as silent wrong-output risks: extremely grazing entries (where reported height is provably non-monotonic, a regime this engine's atmosphere-lookup machinery doesn't yet support) and upward_only fragmentation entries (the same underlying geometric limitation).

Known limitations / explicitly out of scope

• Wake computation (compute_wake=True) is not implemented — confirmed no current caller in this codebase actually triggers it (DynestyMetSim.py always passes compute_wake=False), so this raises a clean NotImplementedError rather than silently producing wrong output.
• The "disruption only" scenario's own performance gap (above) is real, small, and not yet fixed.
• Wiring DynestyMetSim.py (or GUI.py/AutoRefineFit.py/other consumers of MetSimErosion.runSimulation()) to actually import and use this new engine is a natural, low-risk follow-up — deliberately not included here, since the public interface (Constants, Fragment, runSimulation()) was kept compatible specifically to make that a later, separate, easy change.

🤖 Co-Authored-By: Claude Fable 5 & Claude Opus 4.8