From 024c0c1e0c942a0b946893b9a75b9ae15bcf8f67 Mon Sep 17 00:00:00 2001 From: Ryan Malloy Date: Sat, 28 Feb 2026 19:09:08 -0700 Subject: [PATCH] Harden Newton-Raphson gamma bounds, improve v0.20.0 test coverage Add positive-gamma clamp in L1/L2/L3 Newton-Raphson iterations to prevent divergence on extreme mass ratios. Add missing CREATE EXTENSION, tighter L1/L2 precision checks (4 decimal places), lagrange_distance_oe test with Ceres, L1-Earth-L2 ordering verification, and DE fallback tests for planetary moon Lagrange functions. --- src/lagrange.h | 6 +++ test/expected/v020_features.out | 82 +++++++++++++++++++++++++++++++++ test/sql/v020_features.sql | 57 +++++++++++++++++++++++ 3 files changed, 145 insertions(+) diff --git a/src/lagrange.h b/src/lagrange.h index efc1d3c..ea908c1 100644 --- a/src/lagrange.h +++ b/src/lagrange.h @@ -149,6 +149,8 @@ lagrange_corotating(double mu, int point_id, double *x, double *y) return -1; gamma_new = gamma - f / fp; + if (gamma_new <= 0.0) + gamma_new = gamma * 0.5; /* keep gamma positive */ if (fabs(gamma_new - gamma) < 1e-15) break; gamma = gamma_new; @@ -184,6 +186,8 @@ lagrange_corotating(double mu, int point_id, double *x, double *y) return -1; gamma_new = gamma - f / fp; + if (gamma_new <= 0.0) + gamma_new = gamma * 0.5; /* keep gamma positive */ if (fabs(gamma_new - gamma) < 1e-15) break; gamma = gamma_new; @@ -223,6 +227,8 @@ lagrange_corotating(double mu, int point_id, double *x, double *y) return -1; gamma_new = gamma - f / fp; + if (gamma_new <= 0.0) + gamma_new = gamma * 0.5; /* keep gamma positive */ if (fabs(gamma_new - gamma) < 1e-15) break; gamma = gamma_new; diff --git a/test/expected/v020_features.out b/test/expected/v020_features.out index a2a934a..bb2d386 100644 --- a/test/expected/v020_features.out +++ b/test/expected/v020_features.out @@ -1,6 +1,8 @@ -- v020_features: Lagrange point support -- Tests Sun-planet, Earth-Moon, planetary moon Lagrange points, -- Hill radius, zone radius, DE fallback, and input validation. +CREATE EXTENSION IF NOT EXISTS pg_orrery; +NOTICE: extension "pg_orrery" already exists, skipping -- Reference observer: Greenwich, UK \set obs '''(51.4769,-0.0005,0)''' -- Reference time: J2000 epoch (2000-01-01 12:00:00 UTC) @@ -284,6 +286,86 @@ SELECT t (1 row) +-- ============================================================ +-- Tighter L1/L2 precision (4 decimal places) +-- ============================================================ +SELECT + round(helio_distance(lagrange_heliocentric(3, 1, :t ::timestamptz))::numeric, 4) AS earth_l1_dist; + earth_l1_dist +--------------- + 0.9735 +(1 row) + +SELECT + round(helio_distance(lagrange_heliocentric(3, 2, :t ::timestamptz))::numeric, 4) AS earth_l2_dist; + earth_l2_dist +--------------- + 0.9932 +(1 row) + +-- L1 and L2 bracket Earth's heliocentric distance +SELECT + helio_distance(lagrange_heliocentric(3, 1, :t ::timestamptz)) + < + helio_distance(planet_heliocentric(3, :t ::timestamptz)) + AND + helio_distance(planet_heliocentric(3, :t ::timestamptz)) + < + helio_distance(lagrange_heliocentric(3, 2, :t ::timestamptz)) + AS l1_earth_l2_ordering; + l1_earth_l2_ordering +---------------------- + t +(1 row) + +-- ============================================================ +-- lagrange_distance_oe — Ceres distance from Jupiter L4 +-- Ceres orbits at ~2.77 AU, Jupiter L4 at ~5.2 AU, so distance > 2 AU +-- ============================================================ +SELECT + round(lagrange_distance_oe( + 5, 4, + oe_from_mpc('00001 3.33 0.12 K24AM 60.07966 73.42937 80.26860 10.58664 0.0789126 0.21406048 2.7660961 0 MPO838504 8738 115 1801-2024 0.65 M-v 30k MPCLINUX 0000 (1) Ceres 20240825'), + :t ::timestamptz + )::numeric, 2) AS ceres_jup_l4_dist; + ceres_jup_l4_dist +------------------- + 3.03 +(1 row) + +-- Distance should be positive and > 2 AU (main belt vs Trojan zone) +SELECT + lagrange_distance_oe( + 5, 4, + oe_from_mpc('00001 3.33 0.12 K24AM 60.07966 73.42937 80.26860 10.58664 0.0789126 0.21406048 2.7660961 0 MPO838504 8738 115 1801-2024 0.65 M-v 30k MPCLINUX 0000 (1) Ceres 20240825'), + :t ::timestamptz + ) > 2.0 AS ceres_far_from_trojan; + ceres_far_from_trojan +----------------------- + t +(1 row) + +-- ============================================================ +-- DE fallback for planetary moon Lagrange functions +-- ============================================================ +SELECT + round(eq_ra(galilean_lagrange_equatorial_de(0, 4, :t ::timestamptz))::numeric, 4) = + round(eq_ra(galilean_lagrange_equatorial(0, 4, :t ::timestamptz))::numeric, 4) + AS galilean_de_fallback_matches; + galilean_de_fallback_matches +------------------------------ + t +(1 row) + +SELECT + round(eq_ra(saturn_moon_lagrange_equatorial_de(5, 1, :t ::timestamptz))::numeric, 4) = + round(eq_ra(saturn_moon_lagrange_equatorial(5, 1, :t ::timestamptz))::numeric, 4) + AS saturn_de_fallback_matches; + saturn_de_fallback_matches +---------------------------- + t +(1 row) + -- ============================================================ -- Input validation -- ============================================================ diff --git a/test/sql/v020_features.sql b/test/sql/v020_features.sql index a578ba3..e2bde8b 100644 --- a/test/sql/v020_features.sql +++ b/test/sql/v020_features.sql @@ -2,6 +2,8 @@ -- Tests Sun-planet, Earth-Moon, planetary moon Lagrange points, -- Hill radius, zone radius, DE fallback, and input validation. +CREATE EXTENSION IF NOT EXISTS pg_orrery; + -- Reference observer: Greenwich, UK \set obs '''(51.4769,-0.0005,0)''' @@ -178,6 +180,61 @@ SELECT round(hill_radius(5, :t ::timestamptz)::numeric, 4) AS hill_de_matches_vsop; +-- ============================================================ +-- Tighter L1/L2 precision (4 decimal places) +-- ============================================================ + +SELECT + round(helio_distance(lagrange_heliocentric(3, 1, :t ::timestamptz))::numeric, 4) AS earth_l1_dist; + +SELECT + round(helio_distance(lagrange_heliocentric(3, 2, :t ::timestamptz))::numeric, 4) AS earth_l2_dist; + +-- L1 and L2 bracket Earth's heliocentric distance +SELECT + helio_distance(lagrange_heliocentric(3, 1, :t ::timestamptz)) + < + helio_distance(planet_heliocentric(3, :t ::timestamptz)) + AND + helio_distance(planet_heliocentric(3, :t ::timestamptz)) + < + helio_distance(lagrange_heliocentric(3, 2, :t ::timestamptz)) + AS l1_earth_l2_ordering; + +-- ============================================================ +-- lagrange_distance_oe — Ceres distance from Jupiter L4 +-- Ceres orbits at ~2.77 AU, Jupiter L4 at ~5.2 AU, so distance > 2 AU +-- ============================================================ + +SELECT + round(lagrange_distance_oe( + 5, 4, + oe_from_mpc('00001 3.33 0.12 K24AM 60.07966 73.42937 80.26860 10.58664 0.0789126 0.21406048 2.7660961 0 MPO838504 8738 115 1801-2024 0.65 M-v 30k MPCLINUX 0000 (1) Ceres 20240825'), + :t ::timestamptz + )::numeric, 2) AS ceres_jup_l4_dist; + +-- Distance should be positive and > 2 AU (main belt vs Trojan zone) +SELECT + lagrange_distance_oe( + 5, 4, + oe_from_mpc('00001 3.33 0.12 K24AM 60.07966 73.42937 80.26860 10.58664 0.0789126 0.21406048 2.7660961 0 MPO838504 8738 115 1801-2024 0.65 M-v 30k MPCLINUX 0000 (1) Ceres 20240825'), + :t ::timestamptz + ) > 2.0 AS ceres_far_from_trojan; + +-- ============================================================ +-- DE fallback for planetary moon Lagrange functions +-- ============================================================ + +SELECT + round(eq_ra(galilean_lagrange_equatorial_de(0, 4, :t ::timestamptz))::numeric, 4) = + round(eq_ra(galilean_lagrange_equatorial(0, 4, :t ::timestamptz))::numeric, 4) + AS galilean_de_fallback_matches; + +SELECT + round(eq_ra(saturn_moon_lagrange_equatorial_de(5, 1, :t ::timestamptz))::numeric, 4) = + round(eq_ra(saturn_moon_lagrange_equatorial(5, 1, :t ::timestamptz))::numeric, 4) + AS saturn_de_fallback_matches; + -- ============================================================ -- Input validation -- ============================================================