From f9835fe31446921fbcbf64d41d7241cb51845c8b Mon Sep 17 00:00:00 2001 From: Josh <162485031+wiyth00@users.noreply.github.com> Date: Sat, 25 Jul 2026 14:15:46 -0400 Subject: [PATCH] Cross-check duplicate numerics, fix scalar high-e Kepler, canonicalize geodesics - Fix orbital_dynamics.solve_kepler: scalar input with eccentricity > 0.8 raised TypeError because NumPy ufuncs collapse 0-d inputs to scalars and the high-eccentricity starting guess used in-place item assignment. - Make pyproj.Geod (Karney) the canonical WGS84 geodesic in validation.py, matching search.py and sensitivity.py; the retired Vincenty inverse is preserved verbatim in the new test module as an independent reference. - Add tests/test_numeric_crosschecks.py pinning the two Kepler solvers (< 1e-12 rad wrapped), the two element->state converters (< 1 mm, < 1 um/s), and Geod vs the Vincenty reference (sub-mm), plus a regression test for the scalar high-eccentricity crash. --- chained_eclipse/orbital_dynamics.py | 9 +- chained_eclipse/validation.py | 100 +++--------- tests/test_numeric_crosschecks.py | 234 ++++++++++++++++++++++++++++ 3 files changed, 260 insertions(+), 83 deletions(-) create mode 100644 tests/test_numeric_crosschecks.py diff --git a/chained_eclipse/orbital_dynamics.py b/chained_eclipse/orbital_dynamics.py index 2b116de..d9ba917 100644 --- a/chained_eclipse/orbital_dynamics.py +++ b/chained_eclipse/orbital_dynamics.py @@ -35,9 +35,14 @@ def solve_kepler(mean_anomaly_rad: float | np.ndarray, eccentricity: float) -> n """Solve Kepler's equation with a vectorized Newton iteration.""" mean = np.mod(np.asarray(mean_anomaly_rad, dtype=float), 2.0 * np.pi) - eccentric = mean.copy() + # NumPy ufuncs collapse 0-d array inputs to scalars, so a scalar mean + # anomaly reaches this point as np.float64; the previous in-place item + # assignment of the high-eccentricity starting guess therefore raised + # TypeError for scalar inputs with eccentricity > 0.8. if eccentricity > 0.8: - eccentric[...] = np.pi + eccentric = np.full_like(np.asarray(mean), np.pi) + else: + eccentric = np.asarray(mean).copy() for _ in range(20): residual = eccentric - eccentricity * np.sin(eccentric) - mean delta = residual / (1.0 - eccentricity * np.cos(eccentric)) diff --git a/chained_eclipse/validation.py b/chained_eclipse/validation.py index 7f12219..46231f2 100644 --- a/chained_eclipse/validation.py +++ b/chained_eclipse/validation.py @@ -11,10 +11,11 @@ from dataclasses import asdict, dataclass from datetime import datetime, timedelta, timezone -import math from typing import Sequence -from .constants import SECONDS_PER_DAY, WGS84_A_KM, WGS84_F +from pyproj import Geod + +from .constants import SECONDS_PER_DAY from .eclipse_geometry import solve_local_circumstances from .ephemeris import EphemerisContext, enumerate_real_solar_eclipses from .models import RealEclipse @@ -22,6 +23,8 @@ TIMING_TARGET_SECONDS = 60.0 POSITION_TARGET_KM = 25.0 +_WGS84_GEOD = Geod(ellps="WGS84") + @dataclass(frozen=True, slots=True) class PublishedEclipseReference: @@ -162,90 +165,25 @@ def wgs84_geodesic_distance_km( latitude_2_deg: float, longitude_2_deg: float, ) -> float: - """Vincenty inverse distance on the WGS84 reference ellipsoid. - - Validation points are close together, well away from Vincenty's - near-antipodal non-convergence case. + """Geodesic inverse distance on the WGS84 reference ellipsoid. + + ``pyproj.Geod`` (Karney's algorithm) is the project's canonical geodesic + implementation, shared with ``search.py`` and ``sensitivity.py``. The + hand-rolled Vincenty inverse that previously lived here is preserved in + ``tests/test_numeric_crosschecks.py`` as an independent reference, where + the two implementations are asserted to agree to sub-millimetre level + over validation-scale point pairs. """ if latitude_1_deg == latitude_2_deg and longitude_1_deg == longitude_2_deg: return 0.0 - - semi_major = WGS84_A_KM - flattening = WGS84_F - semi_minor = semi_major * (1.0 - flattening) - phi_1 = math.radians(latitude_1_deg) - phi_2 = math.radians(latitude_2_deg) - longitude_delta = math.radians( - (longitude_2_deg - longitude_1_deg + 180.0) % 360.0 - 180.0 - ) - reduced_1 = math.atan((1.0 - flattening) * math.tan(phi_1)) - reduced_2 = math.atan((1.0 - flattening) * math.tan(phi_2)) - sin_u1, cos_u1 = math.sin(reduced_1), math.cos(reduced_1) - sin_u2, cos_u2 = math.sin(reduced_2), math.cos(reduced_2) - longitude = longitude_delta - - for _ in range(100): - sin_longitude = math.sin(longitude) - cos_longitude = math.cos(longitude) - sin_sigma = math.hypot( - cos_u2 * sin_longitude, - cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos_longitude, - ) - if sin_sigma == 0.0: - return 0.0 - cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos_longitude - sigma = math.atan2(sin_sigma, cos_sigma) - sin_alpha = cos_u1 * cos_u2 * sin_longitude / sin_sigma - cos_sq_alpha = 1.0 - sin_alpha * sin_alpha - cos_two_sigma_m = ( - cos_sigma - 2.0 * sin_u1 * sin_u2 / cos_sq_alpha - if cos_sq_alpha > 1e-16 - else 0.0 - ) - coefficient = flattening / 16.0 * cos_sq_alpha * ( - 4.0 + flattening * (4.0 - 3.0 * cos_sq_alpha) - ) - previous = longitude - longitude = longitude_delta + (1.0 - coefficient) * flattening * sin_alpha * ( - sigma - + coefficient - * sin_sigma - * ( - cos_two_sigma_m - + coefficient - * cos_sigma - * (-1.0 + 2.0 * cos_two_sigma_m * cos_two_sigma_m) - ) - ) - if abs(longitude - previous) < 1e-12: - break - else: - raise RuntimeError("WGS84 Vincenty inverse failed to converge") - - reduced_sq = cos_sq_alpha * ( - (semi_major * semi_major - semi_minor * semi_minor) / (semi_minor * semi_minor) - ) - series_a = 1.0 + reduced_sq / 16_384.0 * ( - 4_096.0 + reduced_sq * (-768.0 + reduced_sq * (320.0 - 175.0 * reduced_sq)) - ) - series_b = reduced_sq / 1_024.0 * ( - 256.0 + reduced_sq * (-128.0 + reduced_sq * (74.0 - 47.0 * reduced_sq)) - ) - delta_sigma = series_b * sin_sigma * ( - cos_two_sigma_m - + series_b - / 4.0 - * ( - cos_sigma * (-1.0 + 2.0 * cos_two_sigma_m * cos_two_sigma_m) - - series_b - / 6.0 - * cos_two_sigma_m - * (-3.0 + 4.0 * sin_sigma * sin_sigma) - * (-3.0 + 4.0 * cos_two_sigma_m * cos_two_sigma_m) - ) + _, _, distance_m = _WGS84_GEOD.inv( + longitude_1_deg, + latitude_1_deg, + longitude_2_deg, + latitude_2_deg, ) - return semi_minor * series_a * (sigma - delta_sigma) + return float(distance_m) / 1_000.0 def _validate_one( diff --git a/tests/test_numeric_crosschecks.py b/tests/test_numeric_crosschecks.py new file mode 100644 index 0000000..a593961 --- /dev/null +++ b/tests/test_numeric_crosschecks.py @@ -0,0 +1,234 @@ +"""Cross-checks between the project's duplicate numeric implementations. + +The repository intentionally carries two Kepler solvers (vectorized in +``orbital_dynamics``, scalar in ``stability``, which documents standalone +capability) and two element->state converters. These tests pin the pairs to +each other so any future divergence fails CI, and they pin the canonical +pyproj/Karney geodesic to the hand-rolled Vincenty inverse that previously +lived in ``validation.py``, preserved below as an independent reference. +""" + +from __future__ import annotations + +import itertools +import math + +import numpy as np +import pytest + +from chained_eclipse.constants import MU_EARTH_KM3_S2, WGS84_A_KM, WGS84_F +from chained_eclipse.models import OrbitalElements +from chained_eclipse.orbital_dynamics import elements_to_state, solve_kepler +from chained_eclipse.stability import _solve_kepler, orbital_elements_to_icrf_state +from chained_eclipse.validation import NASA_GSFC_REFERENCES, wgs84_geodesic_distance_km + +_REFERENCE_SITES = tuple( + (reference.latitude_deg, reference.longitude_deg) for reference in NASA_GSFC_REFERENCES +) + + +def _vincenty_reference_km( + latitude_1_deg: float, + longitude_1_deg: float, + latitude_2_deg: float, + longitude_2_deg: float, +) -> float: + """Vincenty inverse distance on the WGS84 reference ellipsoid. + + This is the implementation removed from ``validation.py`` when + ``pyproj.Geod`` became the project's single canonical geodesic. It is + kept verbatim as an independent reference for the cross-check below. + Validation points are well away from Vincenty's near-antipodal + non-convergence case. + """ + + if latitude_1_deg == latitude_2_deg and longitude_1_deg == longitude_2_deg: + return 0.0 + + semi_major = WGS84_A_KM + flattening = WGS84_F + semi_minor = semi_major * (1.0 - flattening) + phi_1 = math.radians(latitude_1_deg) + phi_2 = math.radians(latitude_2_deg) + longitude_delta = math.radians( + (longitude_2_deg - longitude_1_deg + 180.0) % 360.0 - 180.0 + ) + reduced_1 = math.atan((1.0 - flattening) * math.tan(phi_1)) + reduced_2 = math.atan((1.0 - flattening) * math.tan(phi_2)) + sin_u1, cos_u1 = math.sin(reduced_1), math.cos(reduced_1) + sin_u2, cos_u2 = math.sin(reduced_2), math.cos(reduced_2) + longitude = longitude_delta + + for _ in range(100): + sin_longitude = math.sin(longitude) + cos_longitude = math.cos(longitude) + sin_sigma = math.hypot( + cos_u2 * sin_longitude, + cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos_longitude, + ) + if sin_sigma == 0.0: + return 0.0 + cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos_longitude + sigma = math.atan2(sin_sigma, cos_sigma) + sin_alpha = cos_u1 * cos_u2 * sin_longitude / sin_sigma + cos_sq_alpha = 1.0 - sin_alpha * sin_alpha + cos_two_sigma_m = ( + cos_sigma - 2.0 * sin_u1 * sin_u2 / cos_sq_alpha + if cos_sq_alpha > 1e-16 + else 0.0 + ) + coefficient = flattening / 16.0 * cos_sq_alpha * ( + 4.0 + flattening * (4.0 - 3.0 * cos_sq_alpha) + ) + previous = longitude + longitude = longitude_delta + (1.0 - coefficient) * flattening * sin_alpha * ( + sigma + + coefficient + * sin_sigma + * ( + cos_two_sigma_m + + coefficient + * cos_sigma + * (-1.0 + 2.0 * cos_two_sigma_m * cos_two_sigma_m) + ) + ) + if abs(longitude - previous) < 1e-12: + break + else: + raise RuntimeError("WGS84 Vincenty inverse failed to converge") + + reduced_sq = cos_sq_alpha * ( + (semi_major * semi_major - semi_minor * semi_minor) / (semi_minor * semi_minor) + ) + series_a = 1.0 + reduced_sq / 16_384.0 * ( + 4_096.0 + reduced_sq * (-768.0 + reduced_sq * (320.0 - 175.0 * reduced_sq)) + ) + series_b = reduced_sq / 1_024.0 * ( + 256.0 + reduced_sq * (-128.0 + reduced_sq * (74.0 - 47.0 * reduced_sq)) + ) + delta_sigma = series_b * sin_sigma * ( + cos_two_sigma_m + + series_b + / 4.0 + * ( + cos_sigma * (-1.0 + 2.0 * cos_two_sigma_m * cos_two_sigma_m) + - series_b + / 6.0 + * cos_two_sigma_m + * (-3.0 + 4.0 * sin_sigma * sin_sigma) + * (-3.0 + 4.0 * cos_two_sigma_m * cos_two_sigma_m) + ) + ) + return semi_minor * series_a * (sigma - delta_sigma) + + +def _wrapped_difference_rad(first: float, second: float) -> float: + return abs((first - second + math.pi) % (2.0 * math.pi) - math.pi) + + +def test_kepler_solvers_agree_on_grid() -> None: + """Both Kepler solvers must agree to < 1e-12 rad across (e, M) space.""" + + eccentricities = (0.0, 0.01, 0.04, 0.2, 0.5, 0.79, 0.8, 0.85, 0.9, 0.95, 0.99) + worst = 0.0 + for eccentricity in eccentricities: + for mean_deg in np.linspace(0.0, 359.5, 720): + mean = math.radians(float(mean_deg)) + vectorized = float(solve_kepler(mean, eccentricity)) + scalar = _solve_kepler(mean, eccentricity) + worst = max(worst, _wrapped_difference_rad(vectorized, scalar)) + assert worst < 1e-12 + + +def test_solve_kepler_accepts_scalar_high_eccentricity() -> None: + """Regression: scalar input with e > 0.8 raised TypeError before the fix. + + NumPy ufuncs collapse 0-d inputs to scalars, so the high-eccentricity + starting-guess item assignment failed for scalar mean anomalies. + """ + + eccentric = float(solve_kepler(1.0, 0.9)) + assert abs(eccentric - 0.9 * math.sin(eccentric) - 1.0) < 1e-12 + + +def test_elements_to_state_accepts_high_eccentricity() -> None: + """Regression: elements_to_state crashed for e > 0.8 before the fix.""" + + elements = OrbitalElements(eccentricity=0.9, mean_anomaly_deg=57.0) + state = elements_to_state(elements) + assert state.shape == (6,) + assert bool(np.all(np.isfinite(state))) + + +@pytest.mark.parametrize("eccentricity", [0.04, 0.9]) +def test_solve_kepler_vectorized_matches_scalar_calls(eccentricity: float) -> None: + means = np.linspace(0.0, 2.0 * np.pi, 97) + vectorized = np.asarray(solve_kepler(means, eccentricity)) + one_by_one = np.asarray([float(solve_kepler(mean, eccentricity)) for mean in means]) + assert float(np.max(np.abs(vectorized - one_by_one))) < 1e-13 + + +def test_element_to_state_paths_agree() -> None: + """The two element->state converters must agree to < 1 mm and < 1 um/s.""" + + worst_position_km = 0.0 + worst_velocity_km_s = 0.0 + for a, e, inc, node, argp, mean in itertools.product( + (8_000.0, 42_164.0, 180_000.0, 384_400.0), + (0.0, 0.04, 0.3, 0.7, 0.9), + (0.0, 5.0, 28.5, 63.4, 90.0, 150.0), + (0.0, 45.0, 123.4, 270.0), + (0.0, 90.0, 251.0), + (0.0, 10.0, 179.0, 181.0, 350.0), + ): + elements = OrbitalElements( + semimajor_axis_km=a, + eccentricity=e, + inclination_deg=inc, + longitude_ascending_node_deg=node, + argument_periapsis_deg=argp, + mean_anomaly_deg=mean, + ) + primary = elements_to_state(elements) + secondary = orbital_elements_to_icrf_state(elements, MU_EARTH_KM3_S2) + worst_position_km = max( + worst_position_km, float(np.max(np.abs(primary[:3] - secondary[:3]))) + ) + worst_velocity_km_s = max( + worst_velocity_km_s, float(np.max(np.abs(primary[3:] - secondary[3:]))) + ) + assert worst_position_km < 1e-6 + assert worst_velocity_km_s < 1e-9 + + +def test_geod_matches_vincenty_reference_for_validation_scale_pairs() -> None: + """Canonical Geod vs the Vincenty reference: < 1 mm at validation scales.""" + + offsets_deg = (0.0004, 0.01, 0.05, 0.2) + for latitude, longitude in _REFERENCE_SITES: + for offset in offsets_deg: + for delta_lat, delta_lon in ((offset, 0.0), (0.0, offset), (offset, -offset)): + expected = _vincenty_reference_km( + latitude, longitude, latitude + delta_lat, longitude + delta_lon + ) + actual = wgs84_geodesic_distance_km( + latitude, longitude, latitude + delta_lat, longitude + delta_lon + ) + assert abs(actual - expected) < 1e-6 + + +def test_geod_matches_vincenty_reference_for_long_pairs() -> None: + """Continental-scale pairs between reference sites agree to < 1 cm.""" + + for (lat_1, lon_1), (lat_2, lon_2) in itertools.combinations(_REFERENCE_SITES, 2): + expected = _vincenty_reference_km(lat_1, lon_1, lat_2, lon_2) + actual = wgs84_geodesic_distance_km(lat_1, lon_1, lat_2, lon_2) + assert abs(actual - expected) < 1e-5 + + +def test_geodesic_zero_and_symmetry() -> None: + latitude, longitude = _REFERENCE_SITES[0] + assert wgs84_geodesic_distance_km(latitude, longitude, latitude, longitude) == 0.0 + forward = wgs84_geodesic_distance_km(latitude, longitude, -12.5, 101.25) + backward = wgs84_geodesic_distance_km(-12.5, 101.25, latitude, longitude) + assert abs(forward - backward) < 1e-9