Cross-check duplicate numerics, fix scalar high-e Kepler crash, canonicalize geodesics (item 6)

[wiyth00]

Summary

Work-plan item 6: the repo carries duplicate implementations of three numeric pieces. This PR pins them against each other with tests, fixes a real latent bug the cross-check flushed out, and reduces the geodesics to one canonical implementation.

• Bug fix — orbital_dynamics.solve_kepler crashed with TypeError for scalar input when eccentricity > 0.8: NumPy ufuncs collapse 0-d inputs to scalars, so np.mod(np.asarray(scalar), 2π) returns np.float64 and the high-eccentricity starting guess eccentric[...] = np.pi failed on item assignment. This means elements_to_state with e > 0.8 crashed on main. Fixed with np.full_like / explicit array copy; behavior for array inputs and e ≤ 0.8 is bit-identical.
• One canonical geodesic — validation.wgs84_geodesic_distance_km now delegates to pyproj.Geod(ellps="WGS84").inv (Karney), the same implementation search.py and sensitivity.py already use (pyproj became a declared dependency in fix: declare pyproj as a direct dependency #5). Signature and name unchanged. The retired hand-rolled Vincenty inverse moves verbatim into the test module as an independent reference.
• New tests/test_numeric_crosschecks.py:
  • both Kepler solvers agree to < 1e-12 rad (wrapped) over an 11×720 (e, M) grid — the two solvers normalize M into different ranges, so the comparison wraps the difference;
  • both element→state converters (orbital_dynamics.elements_to_state vs stability.orbital_elements_to_icrf_state, which is deliberately standalone) agree to < 1 mm and < 1 µm/s over a 7,200-case element grid;
  • Geod vs the Vincenty reference: < 1 mm at validation scales, < 1 cm on continental-scale pairs between the four NASA reference sites;
  • regression tests that fail on main for the scalar high-eccentricity crash (per ground rule: every bug fix ships a test that fails on main).

Measured agreement in a NumPy-only replica before encoding tolerances: 8.9e-15 rad (Kepler, wrapped), 2.9e-10 km and 2.0e-14 km/s (element→state grid) — the asserted tolerances have 3+ orders of margin.

stability._solve_kepler and its converter are intentionally kept (the module documents standalone capability); the cross-check tests now guard against divergence instead of deleting the redundancy.

Scientific impact

• No change to computed results
• Changes computed results (explain below)

Effectively none, but stated honestly:

• The Karney geodesic and the old Vincenty inverse agree to sub-millimetre over the validation pairs (proven by the new tests), so validation_report.json coordinate-error values may shift in far decimal places, orders of magnitude below the 25 km target. No pass/fail status can change.
• The Kepler fix only affects a path that previously crashed (scalar input, e > 0.8); no currently-computable result changes. The baseline design (e = 0.04) is untouched.
• The 2026-08-12 anchor event is unaffected.

Checklist

• Tests added/updated (new cross-check + regression module)
• Regression test fails on main (scalar high-e Kepler tests)
• CI green (first integrated run happens in CI — sandbox here has no scipy/pyproj/skyfield; numerics verified in a NumPy-only replica)
• No results-changing edits to the immovable 2026-08-12 event