Ryan Malloy
024c0c1e0c
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.
406 lines
13 KiB
Plaintext
406 lines
13 KiB
Plaintext
-- 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)
|
|
\set t '''2000-01-01 12:00:00+00'''
|
|
-- ============================================================
|
|
-- Sun-Earth L1/L2: should be ~0.01 AU from Earth (~1.5 million km)
|
|
-- SOHO is at L1, JWST at L2.
|
|
-- ============================================================
|
|
-- L1 heliocentric: should be close to Earth's heliocentric (~1 AU from Sun)
|
|
SELECT
|
|
round(helio_distance(lagrange_heliocentric(3, 1, :t ::timestamptz))::numeric, 2) AS sun_dist_au;
|
|
sun_dist_au
|
|
-------------
|
|
0.97
|
|
(1 row)
|
|
|
|
-- L2 heliocentric: also ~1 AU from Sun, slightly further than L1
|
|
SELECT
|
|
round(helio_distance(lagrange_heliocentric(3, 2, :t ::timestamptz))::numeric, 2) AS sun_dist_au;
|
|
sun_dist_au
|
|
-------------
|
|
0.99
|
|
(1 row)
|
|
|
|
-- L1 between Sun and Earth (closer to Sun than L2)
|
|
SELECT
|
|
helio_distance(lagrange_heliocentric(3, 1, :t ::timestamptz))
|
|
<
|
|
helio_distance(lagrange_heliocentric(3, 2, :t ::timestamptz))
|
|
AS l1_closer_than_l2;
|
|
l1_closer_than_l2
|
|
-------------------
|
|
t
|
|
(1 row)
|
|
|
|
-- ============================================================
|
|
-- Sun-Jupiter L4/L5: ~60 degrees from Jupiter, ~5.2 AU from Sun
|
|
-- These are the Trojan asteroid zones.
|
|
-- ============================================================
|
|
-- L4/L5 should be ~5.2 AU from Sun
|
|
SELECT
|
|
round(helio_distance(lagrange_heliocentric(5, 4, :t ::timestamptz))::numeric, 1) AS l4_sun_dist;
|
|
l4_sun_dist
|
|
-------------
|
|
5.0
|
|
(1 row)
|
|
|
|
SELECT
|
|
round(helio_distance(lagrange_heliocentric(5, 5, :t ::timestamptz))::numeric, 1) AS l5_sun_dist;
|
|
l5_sun_dist
|
|
-------------
|
|
5.0
|
|
(1 row)
|
|
|
|
-- L4 and L5 equidistant from Sun (within 0.001 AU)
|
|
SELECT
|
|
abs(
|
|
helio_distance(lagrange_heliocentric(5, 4, :t ::timestamptz))
|
|
-
|
|
helio_distance(lagrange_heliocentric(5, 5, :t ::timestamptz))
|
|
) < 0.001 AS l4_l5_equidistant;
|
|
l4_l5_equidistant
|
|
-------------------
|
|
t
|
|
(1 row)
|
|
|
|
-- ============================================================
|
|
-- Earth-Moon L1: ~326,000 km from Earth
|
|
-- ============================================================
|
|
-- lunar_lagrange_equatorial returns distance in km
|
|
SELECT
|
|
round(eq_distance(lunar_lagrange_equatorial(1, :t ::timestamptz))::numeric, -3)
|
|
BETWEEN 300000 AND 360000 AS em_l1_in_range;
|
|
em_l1_in_range
|
|
----------------
|
|
t
|
|
(1 row)
|
|
|
|
-- ============================================================
|
|
-- lagrange_observe returns valid az/el
|
|
-- ============================================================
|
|
SELECT
|
|
topo_elevation(lagrange_observe(3, 2, :obs ::observer, :t ::timestamptz))
|
|
BETWEEN -90 AND 90 AS valid_elevation;
|
|
valid_elevation
|
|
-----------------
|
|
t
|
|
(1 row)
|
|
|
|
-- lagrange_equatorial returns valid RA/Dec
|
|
SELECT
|
|
eq_ra(lagrange_equatorial(3, 1, :t ::timestamptz)) BETWEEN 0 AND 24 AS valid_ra,
|
|
eq_dec(lagrange_equatorial(3, 1, :t ::timestamptz)) BETWEEN -90 AND 90 AS valid_dec;
|
|
valid_ra | valid_dec
|
|
----------+-----------
|
|
t | t
|
|
(1 row)
|
|
|
|
-- ============================================================
|
|
-- lagrange_distance self-test: L-point distance to itself ≈ 0
|
|
-- ============================================================
|
|
SELECT
|
|
round(lagrange_distance(
|
|
5, 4,
|
|
lagrange_heliocentric(5, 4, :t ::timestamptz),
|
|
:t ::timestamptz
|
|
)::numeric, 10) AS self_distance;
|
|
self_distance
|
|
---------------
|
|
0.0000000000
|
|
(1 row)
|
|
|
|
-- ============================================================
|
|
-- Hill radius
|
|
-- ============================================================
|
|
-- Jupiter Hill radius ~0.35 AU
|
|
SELECT
|
|
round(hill_radius(5, :t ::timestamptz)::numeric, 2)
|
|
BETWEEN 0.30 AND 0.40 AS jupiter_hill_ok;
|
|
jupiter_hill_ok
|
|
-----------------
|
|
t
|
|
(1 row)
|
|
|
|
-- Earth Hill radius ~0.01 AU
|
|
SELECT
|
|
round(hill_radius(3, :t ::timestamptz)::numeric, 3)
|
|
BETWEEN 0.008 AND 0.012 AS earth_hill_ok;
|
|
earth_hill_ok
|
|
---------------
|
|
t
|
|
(1 row)
|
|
|
|
-- Lunar Hill radius (much smaller, AU)
|
|
SELECT
|
|
hill_radius_lunar(:t ::timestamptz) > 0 AS lunar_hill_positive;
|
|
lunar_hill_positive
|
|
---------------------
|
|
t
|
|
(1 row)
|
|
|
|
-- ============================================================
|
|
-- Zone radius
|
|
-- ============================================================
|
|
SELECT
|
|
lagrange_zone_radius(5, 4, :t ::timestamptz) > 0 AS jup_l4_zone_positive;
|
|
jup_l4_zone_positive
|
|
----------------------
|
|
t
|
|
(1 row)
|
|
|
|
SELECT
|
|
lagrange_zone_radius(5, 1, :t ::timestamptz) > 0 AS jup_l1_zone_positive;
|
|
jup_l1_zone_positive
|
|
----------------------
|
|
t
|
|
(1 row)
|
|
|
|
-- ============================================================
|
|
-- Convenience functions
|
|
-- ============================================================
|
|
-- lagrange_mass_ratio returns small positive number
|
|
SELECT
|
|
lagrange_mass_ratio(5) > 0 AND lagrange_mass_ratio(5) < 0.01 AS jupiter_mu_ok;
|
|
jupiter_mu_ok
|
|
---------------
|
|
t
|
|
(1 row)
|
|
|
|
SELECT
|
|
lagrange_mass_ratio(3) > 0 AND lagrange_mass_ratio(3) < 0.001 AS earth_mu_ok;
|
|
earth_mu_ok
|
|
-------------
|
|
t
|
|
(1 row)
|
|
|
|
-- lagrange_point_name
|
|
SELECT lagrange_point_name(1) AS l1_name;
|
|
l1_name
|
|
---------
|
|
L1
|
|
(1 row)
|
|
|
|
SELECT lagrange_point_name(5) AS l5_name;
|
|
l5_name
|
|
---------
|
|
L5
|
|
(1 row)
|
|
|
|
-- ============================================================
|
|
-- All planets produce valid results
|
|
-- ============================================================
|
|
SELECT body_id,
|
|
round(helio_distance(lagrange_heliocentric(body_id, 1, :t ::timestamptz))::numeric, 2) AS sun_dist_au
|
|
FROM generate_series(1, 8) AS body_id
|
|
ORDER BY body_id;
|
|
body_id | sun_dist_au
|
|
---------+-------------
|
|
1 | 0.46
|
|
2 | 0.71
|
|
3 | 0.97
|
|
4 | 1.38
|
|
5 | 4.63
|
|
6 | 8.77
|
|
7 | 19.44
|
|
8 | 29.35
|
|
(8 rows)
|
|
|
|
-- ============================================================
|
|
-- Planetary moon Lagrange points
|
|
-- ============================================================
|
|
-- Galilean: Io L4 (body=0, point=4)
|
|
SELECT
|
|
eq_ra(galilean_lagrange_equatorial(0, 4, :t ::timestamptz)) BETWEEN 0 AND 24
|
|
AS io_l4_valid_ra;
|
|
io_l4_valid_ra
|
|
----------------
|
|
t
|
|
(1 row)
|
|
|
|
-- Saturn: Titan L1 (body=5, point=1)
|
|
SELECT
|
|
eq_ra(saturn_moon_lagrange_equatorial(5, 1, :t ::timestamptz)) BETWEEN 0 AND 24
|
|
AS titan_l1_valid_ra;
|
|
titan_l1_valid_ra
|
|
-------------------
|
|
t
|
|
(1 row)
|
|
|
|
-- Uranus: Titania L2 (body=3, point=2)
|
|
SELECT
|
|
eq_ra(uranus_moon_lagrange_equatorial(3, 2, :t ::timestamptz)) BETWEEN 0 AND 24
|
|
AS titania_l2_valid_ra;
|
|
titania_l2_valid_ra
|
|
---------------------
|
|
t
|
|
(1 row)
|
|
|
|
-- Mars: Phobos L5 (body=0, point=5)
|
|
SELECT
|
|
eq_ra(mars_moon_lagrange_equatorial(0, 5, :t ::timestamptz)) BETWEEN 0 AND 24
|
|
AS phobos_l5_valid_ra;
|
|
phobos_l5_valid_ra
|
|
--------------------
|
|
t
|
|
(1 row)
|
|
|
|
-- Galilean observe returns valid topocentric
|
|
SELECT
|
|
topo_elevation(galilean_lagrange_observe(2, 3, :obs ::observer, :t ::timestamptz))
|
|
BETWEEN -90 AND 90 AS ganymede_l3_valid_el;
|
|
ganymede_l3_valid_el
|
|
----------------------
|
|
t
|
|
(1 row)
|
|
|
|
-- ============================================================
|
|
-- DE fallback (no DE loaded, should produce same results as VSOP87)
|
|
-- ============================================================
|
|
SELECT
|
|
round(helio_distance(lagrange_heliocentric_de(3, 1, :t ::timestamptz))::numeric, 4) AS de_l1_dist;
|
|
de_l1_dist
|
|
------------
|
|
0.9735
|
|
(1 row)
|
|
|
|
SELECT
|
|
round(helio_distance(lagrange_heliocentric(3, 1, :t ::timestamptz))::numeric, 4) AS vsop_l1_dist;
|
|
vsop_l1_dist
|
|
--------------
|
|
0.9735
|
|
(1 row)
|
|
|
|
-- DE hill_radius fallback
|
|
SELECT
|
|
round(hill_radius_de(5, :t ::timestamptz)::numeric, 4) =
|
|
round(hill_radius(5, :t ::timestamptz)::numeric, 4)
|
|
AS hill_de_matches_vsop;
|
|
hill_de_matches_vsop
|
|
----------------------
|
|
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
|
|
-- ============================================================
|
|
-- Bad body_id
|
|
SELECT lagrange_heliocentric(0, 1, :t ::timestamptz); -- Sun not valid
|
|
ERROR: lagrange_heliocentric: body_id 0 must be 1-8 (Mercury-Neptune)
|
|
SELECT lagrange_heliocentric(9, 1, :t ::timestamptz); -- body 9 invalid
|
|
ERROR: lagrange_heliocentric: body_id 9 must be 1-8 (Mercury-Neptune)
|
|
-- Bad point_id
|
|
SELECT lagrange_heliocentric(3, 0, :t ::timestamptz); -- point 0 invalid
|
|
ERROR: lagrange_heliocentric: point_id 0 must be 1-5 (L1-L5)
|
|
SELECT lagrange_heliocentric(3, 6, :t ::timestamptz); -- point 6 invalid
|
|
ERROR: lagrange_heliocentric: point_id 6 must be 1-5 (L1-L5)
|
|
-- Bad lunar point_id
|
|
SELECT lunar_lagrange_equatorial(0, :t ::timestamptz); -- point 0 invalid
|
|
ERROR: lunar_lagrange_equatorial: point_id 0 must be 1-5
|
|
SELECT lunar_lagrange_equatorial(6, :t ::timestamptz); -- point 6 invalid
|
|
ERROR: lunar_lagrange_equatorial: point_id 6 must be 1-5
|
|
-- Bad planetary moon body_id
|
|
SELECT galilean_lagrange_equatorial(4, 1, :t ::timestamptz); -- Galilean 4 invalid
|
|
ERROR: galilean_lagrange_equatorial: body_id 4 must be 0-3
|
|
SELECT saturn_moon_lagrange_equatorial(8, 1, :t ::timestamptz); -- Saturn 8 invalid
|
|
ERROR: saturn_moon_lagrange_equatorial: body_id 8 must be 0-7
|
|
SELECT uranus_moon_lagrange_equatorial(5, 1, :t ::timestamptz); -- Uranus 5 invalid
|
|
ERROR: uranus_moon_lagrange_equatorial: body_id 5 must be 0-4
|
|
SELECT mars_moon_lagrange_equatorial(2, 1, :t ::timestamptz); -- Mars 2 invalid
|
|
ERROR: mars_moon_lagrange_equatorial: body_id 2 must be 0-1
|
|
-- lagrange_mass_ratio bad body
|
|
SELECT lagrange_mass_ratio(0);
|
|
ERROR: lagrange_mass_ratio: body_id 0 must be 1-8
|
|
SELECT lagrange_mass_ratio(9);
|
|
ERROR: lagrange_mass_ratio: body_id 9 must be 1-8
|
|
-- lagrange_point_name bad id
|
|
SELECT lagrange_point_name(0);
|
|
ERROR: lagrange_point_name: point_id 0 must be 1-5
|
|
SELECT lagrange_point_name(6);
|
|
ERROR: lagrange_point_name: point_id 6 must be 1-5
|