warehack.ing/pg_orrery
2
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
pg_orrery/test/sql/orbital_elements.sql
Ryan Malloy b33d63034b Add v0.9.0 apparent position features: equatorial type, refraction, proper motion, light-time 
New equatorial type (24 bytes: RA/Dec/distance) captures apparent coordinates
of date — what the observation pipeline computes at precession step 3 but was
discarding before hour angle conversion. Matches telescope GoTo mount conventions.

24 new SQL functions (82 → 106 total):
- equatorial type I/O + 3 accessors (eq_ra, eq_dec, eq_distance)
- Satellite RA/Dec: eci_to_equatorial (topocentric), eci_to_equatorial_geo (geocentric)
- Solar system equatorial: planet/sun/moon/small_body_equatorial
- Atmospheric refraction: Bennett (1982) with domain clamp at -1 deg
- Refracted pass prediction: predict_passes_refracted (horizon at -0.569 deg)
- Stellar proper motion: star_observe_pm, star_equatorial_pm (Hipparcos/Gaia convention)
- Light-time correction: planet/sun/small_body_observe_apparent, *_equatorial_apparent
- DE equatorial variants: planet_equatorial_de, moon_equatorial_de

Also includes v0.8.0 orbital_elements type (MPC parser, small_body_observe),
GiST 0-based indexing fix, llms.txt updates, and doc improvements.

All 18 regression suites pass. Zero build warnings (GCC + Clang).
2026-02-21 15:31:46 -07:00

209 lines
11 KiB
SQL
Raw Blame History

-- orbital_elements regression tests
--
-- Tests orbital_elements type I/O, accessors, MPC parser,
-- small_body_heliocentric(), and small_body_observe().
CREATE EXTENSION IF NOT EXISTS pg_orrery;
\set boulder '''40.015N 105.270W 1655m'''::observer
-- ============================================================
-- Test 1: Type I/O round-trip
-- Construct orbital_elements from text literal, verify output matches.
-- Angles in degrees, stored as radians internally.
-- ============================================================
SELECT 'io_roundtrip' AS test,
'(2460600.000000,2.5478000000,0.0789126000,10.586640,73.429370,80.268600,2460319.000000,3.33,0.12)'::orbital_elements AS oe;
-- ============================================================
-- Test 2: Accessor functions (direct field access)
-- ============================================================
SELECT 'accessor_epoch' AS test,
round(oe_epoch('(2460600.0,2.5478,0.0789126,10.58664,73.42937,80.2686,2460319.0,3.33,0.12)'::orbital_elements)::numeric, 1) AS epoch_jd;
SELECT 'accessor_q' AS test,
round(oe_perihelion('(2460600.0,2.5478,0.0789126,10.58664,73.42937,80.2686,2460319.0,3.33,0.12)'::orbital_elements)::numeric, 4) AS q_au;
SELECT 'accessor_e' AS test,
round(oe_eccentricity('(2460600.0,2.5478,0.0789126,10.58664,73.42937,80.2686,2460319.0,3.33,0.12)'::orbital_elements)::numeric, 7) AS ecc;
SELECT 'accessor_inc' AS test,
round(oe_inclination('(2460600.0,2.5478,0.0789126,10.58664,73.42937,80.2686,2460319.0,3.33,0.12)'::orbital_elements)::numeric, 5) AS inc_deg;
SELECT 'accessor_omega' AS test,
round(oe_arg_perihelion('(2460600.0,2.5478,0.0789126,10.58664,73.42937,80.2686,2460319.0,3.33,0.12)'::orbital_elements)::numeric, 5) AS omega_deg;
SELECT 'accessor_Omega' AS test,
round(oe_raan('(2460600.0,2.5478,0.0789126,10.58664,73.42937,80.2686,2460319.0,3.33,0.12)'::orbital_elements)::numeric, 4) AS Omega_deg;
SELECT 'accessor_tp' AS test,
round(oe_tp('(2460600.0,2.5478,0.0789126,10.58664,73.42937,80.2686,2460319.0,3.33,0.12)'::orbital_elements)::numeric, 1) AS tp_jd;
SELECT 'accessor_h' AS test,
round(oe_h_mag('(2460600.0,2.5478,0.0789126,10.58664,73.42937,80.2686,2460319.0,3.33,0.12)'::orbital_elements)::numeric, 2) AS h_mag;
SELECT 'accessor_g' AS test,
round(oe_g_slope('(2460600.0,2.5478,0.0789126,10.58664,73.42937,80.2686,2460319.0,3.33,0.12)'::orbital_elements)::numeric, 2) AS g_slope;
-- ============================================================
-- Test 3: Computed accessors
-- ============================================================
-- Semi-major axis: a = q / (1 - e) = 2.5478 / (1 - 0.0789126) ~ 2.766 AU
SELECT 'sma_elliptic' AS test,
round(oe_semi_major_axis('(2460600.0,2.5478,0.0789126,10.58664,73.42937,80.2686,2460319.0,3.33,0.12)'::orbital_elements)::numeric, 3) AS a_au;
-- Period: a^1.5 years. a ~ 2.766, period ~ 4.599 years
SELECT 'period_elliptic' AS test,
round(oe_period_years('(2460600.0,2.5478,0.0789126,10.58664,73.42937,80.2686,2460319.0,3.33,0.12)'::orbital_elements)::numeric, 2) AS period_yr;
-- Hyperbolic orbit (e=1.5): semi-major axis and period should be NULL
SELECT 'sma_hyperbolic' AS test,
oe_semi_major_axis('(2460600.0,1.0,1.5,10.0,73.0,80.0,2460319.0,NaN,NaN)'::orbital_elements) IS NULL AS is_null;
SELECT 'period_hyperbolic' AS test,
oe_period_years('(2460600.0,1.0,1.5,10.0,73.0,80.0,2460319.0,NaN,NaN)'::orbital_elements) IS NULL AS is_null;
-- ============================================================
-- Test 4: MPC parser -- (1) Ceres
--
-- MPCORB.DAT line for Ceres (epoch 2024 Oct 17.0 = K24AM):
-- Packed epoch K24AM -> K=2000, 24, A=Oct, M=22 -> 2024-10-22
-- (Actually: K=2000, year=24, month=A=10, day=M=22)
-- JD for 2024-10-22 = 2460605.5
--
-- a = 2.7660961 AU, e = 0.0789126
-- q = a*(1-e) = 2.7660961*(1-0.0789126) = 2.5478...
-- M = 60.07966 deg
-- ============================================================
SELECT 'mpc_ceres_sma' AS test,
round(oe_semi_major_axis(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'
))::numeric, 4) AS a_au;
SELECT 'mpc_ceres_q' AS test,
round(oe_perihelion(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'
))::numeric, 4) AS q_au;
SELECT 'mpc_ceres_ecc' AS test,
round(oe_eccentricity(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'
))::numeric, 7) AS ecc;
SELECT 'mpc_ceres_inc' AS test,
round(oe_inclination(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'
))::numeric, 5) AS inc_deg;
SELECT 'mpc_ceres_h' AS test,
round(oe_h_mag(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'
))::numeric, 2) AS h_mag;
-- ============================================================
-- Test 5: small_body_heliocentric -- compare to kepler_propagate
-- Both should produce the same heliocentric position for Ceres.
-- ============================================================
SELECT 'helio_vs_kepler' AS test,
round(helio_x(small_body_heliocentric(
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'),
'2025-01-01 00:00:00+00'
))::numeric, 6) AS x_au,
round(helio_y(small_body_heliocentric(
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'),
'2025-01-01 00:00:00+00'
))::numeric, 6) AS y_au,
round(helio_z(small_body_heliocentric(
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'),
'2025-01-01 00:00:00+00'
))::numeric, 6) AS z_au;
-- Verify heliocentric distance is reasonable (Ceres orbits ~2.5-3.0 AU)
SELECT 'helio_distance' AS test,
round(helio_distance(small_body_heliocentric(
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'),
'2025-01-01 00:00:00+00'
))::numeric, 2) AS dist_au,
helio_distance(small_body_heliocentric(
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'),
'2025-01-01 00:00:00+00'
)) BETWEEN 2.5 AND 3.0 AS in_range;
-- ============================================================
-- Test 6: small_body_observe -- topocentric observation
-- Verify we get reasonable az/el/range for Ceres from Boulder.
-- Range should be on the order of 1.5-4.5 AU in km.
-- ============================================================
SELECT 'observe_range' AS test,
topo_range(small_body_observe(
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'),
:boulder, '2025-01-01 00:00:00+00'
)) BETWEEN 1.5 * 149597870.7 AND 4.5 * 149597870.7 AS range_reasonable;
-- Elevation should be a finite number
SELECT 'observe_elevation' AS test,
topo_elevation(small_body_observe(
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'),
:boulder, '2025-01-01 00:00:00+00'
)) BETWEEN -90 AND 90 AS el_reasonable;
-- ============================================================
-- Test 7: Parabolic/hyperbolic elements via text I/O
-- Verify type handles e >= 1 without error.
-- ============================================================
SELECT 'parabolic_io' AS test,
oe_eccentricity('(2460600.0,1.0,1.0,45.0,90.0,180.0,2460500.0,12.0,0.15)'::orbital_elements) AS ecc;
SELECT 'hyperbolic_io' AS test,
oe_eccentricity('(2460600.0,0.5,2.5,30.0,60.0,120.0,2460400.0,8.0,0.04)'::orbital_elements) AS ecc;
-- ============================================================
-- Test 8: Error paths
-- ============================================================
-- Invalid text input (wrong number of fields)
SELECT 'bad_text' AS test, '(1.0,2.0,3.0)'::orbital_elements;
-- Negative perihelion distance
SELECT 'negative_q' AS test, '(2460600.0,-1.0,0.5,10.0,73.0,80.0,2460319.0,3.33,0.12)'::orbital_elements;
-- Negative eccentricity
SELECT 'negative_e' AS test, '(2460600.0,1.0,-0.1,10.0,73.0,80.0,2460319.0,3.33,0.12)'::orbital_elements;
-- MPC line too short
SELECT 'mpc_short' AS test, oe_from_mpc('too short');
-- ============================================================
-- Test 9: MPC packed date decoding
-- Verify K24AM -> 2024-10-22 -> JD 2460605.5
-- ============================================================
SELECT 'mpc_epoch' AS test,
round(oe_epoch(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'
))::numeric, 1) AS epoch_jd;
-- ============================================================
-- Test 10: Consistency -- small_body_observe vs comet_observe
-- Both pipelines should produce the same topocentric result
-- when fed the same elements and Earth position.
-- ============================================================
WITH ceres AS (
SELECT 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'
) AS oe
), earth AS (
SELECT planet_heliocentric(3, '2025-06-15 12:00:00+00') AS pos
)
SELECT 'pipeline_match' AS test,
round(abs(
topo_elevation(small_body_observe(c.oe, :boulder, '2025-06-15 12:00:00+00'))
- topo_elevation(comet_observe(
oe_perihelion(c.oe), oe_eccentricity(c.oe),
oe_inclination(c.oe), oe_arg_perihelion(c.oe), oe_raan(c.oe),
oe_tp(c.oe),
helio_x(e.pos), helio_y(e.pos), helio_z(e.pos),
:boulder, '2025-06-15 12:00:00+00'
))
)::numeric, 6) AS el_diff_deg
FROM ceres c, earth e;
Reference in New Issue View Git Blame Copy Permalink