Ryan Malloy
0544a78276
Phase 1 — Stars, comets, Keplerian propagation: - star_observe() / star_observe_safe(): fixed star alt/az via IAU 1976 precession, equatorial-to-horizontal transform - kepler_propagate(): two-body Keplerian orbit propagation for elliptic, parabolic, and hyperbolic orbits - comet_observe(): observe comets/asteroids from orbital elements - heliocentric type: ecliptic J2000 position (x, y, z in AU) Phase 2 — VSOP87 planets, ELP82B Moon, Sun: - planet_heliocentric(): VSOP87 heliocentric ecliptic J2000 positions for Mercury through Neptune (Bretagnon & Francou, MIT) - planet_observe(): full observation pipeline for any planet - sun_observe(): Sun position from negated Earth VSOP87 - moon_observe(): ELP2000-82B lunar position (Chapront-Touzé, MIT) - Clean-room precession (IAU 2006) and sidereal time (IERS 2010) - elliptic_to_rectangular utility (Stellarium, MIT) All Stellarium extractions are MIT-licensed, thread-safe (static caching removed for PARALLEL SAFE), zero external data files. All 9 regression tests pass (90ms total).
105 lines
5.2 KiB
SQL
105 lines
5.2 KiB
SQL
-- kepler_comet regression tests
|
|
--
|
|
-- Tests Keplerian propagation and comet observation.
|
|
-- Verifies Kepler equation solver for elliptic, parabolic,
|
|
-- and hyperbolic orbits.
|
|
|
|
\set boulder '''40.015N 105.270W 1655m'''::observer
|
|
|
|
-- ============================================================
|
|
-- Test 1: Circular orbit (e=0) at 1 AU
|
|
-- A body in a circular orbit at 1 AU with i=0, omega=0, Omega=0,
|
|
-- perihelion at J2000.0 should complete one orbit in ~365.25 days.
|
|
-- At T=0 (perihelion), position should be (1, 0, 0).
|
|
-- ============================================================
|
|
SELECT 'circular_t0' AS test,
|
|
round(helio_x(kepler_propagate(1.0, 0.0, 0.0, 0.0, 0.0,
|
|
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS x,
|
|
round(helio_y(kepler_propagate(1.0, 0.0, 0.0, 0.0, 0.0,
|
|
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS y,
|
|
round(helio_z(kepler_propagate(1.0, 0.0, 0.0, 0.0, 0.0,
|
|
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS z;
|
|
|
|
-- ============================================================
|
|
-- Test 2: Circular orbit at 1 AU, quarter orbit (~91.3 days)
|
|
-- Should be at approximately (0, 1, 0) for i=0 prograde orbit.
|
|
-- ============================================================
|
|
SELECT 'circular_quarter' AS test,
|
|
round(helio_x(kepler_propagate(1.0, 0.0, 0.0, 0.0, 0.0,
|
|
2451545.0, '2000-04-01 18:00:00+00'))::numeric, 1) AS x_approx,
|
|
round(helio_y(kepler_propagate(1.0, 0.0, 0.0, 0.0, 0.0,
|
|
2451545.0, '2000-04-01 18:00:00+00'))::numeric, 1) AS y_approx;
|
|
|
|
-- ============================================================
|
|
-- Test 3: Eccentric elliptic orbit (e=0.5)
|
|
-- q=0.5 AU, e=0.5 -> a=1.0 AU, same period as Earth.
|
|
-- At perihelion (t=0), position should be (0.5, 0, 0).
|
|
-- ============================================================
|
|
SELECT 'eccentric_t0' AS test,
|
|
round(helio_x(kepler_propagate(0.5, 0.5, 0.0, 0.0, 0.0,
|
|
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS x,
|
|
round(helio_y(kepler_propagate(0.5, 0.5, 0.0, 0.0, 0.0,
|
|
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS y;
|
|
|
|
-- ============================================================
|
|
-- Test 4: Inclined orbit (i=90 deg)
|
|
-- With i=90, the orbit is polar. At perihelion the position
|
|
-- is still in the orbital plane but rotated by inclination.
|
|
-- q=1.0 AU, e=0, i=90, omega=0, Omega=0 at t=0:
|
|
-- position should be (1, 0, 0) regardless of inclination
|
|
-- (perihelion direction is in the node line).
|
|
-- ============================================================
|
|
SELECT 'inclined_t0' AS test,
|
|
round(helio_x(kepler_propagate(1.0, 0.0, 90.0, 0.0, 0.0,
|
|
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS x,
|
|
round(helio_z(kepler_propagate(1.0, 0.0, 90.0, 0.0, 0.0,
|
|
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS z;
|
|
|
|
-- ============================================================
|
|
-- Test 5: Heliocentric distance conservation for circular orbit
|
|
-- Distance should remain ~1.0 AU at all times.
|
|
-- ============================================================
|
|
SELECT 'distance_conservation' AS test,
|
|
round(helio_distance(kepler_propagate(1.0, 0.0, 0.0, 0.0, 0.0,
|
|
2451545.0, '2000-07-01 12:00:00+00'))::numeric, 6) AS dist_au;
|
|
|
|
-- ============================================================
|
|
-- Test 6: Hyperbolic orbit (e=1.5)
|
|
-- q=1.0 AU, e=1.5, i=0, omega=0, Omega=0
|
|
-- At perihelion, position = (1, 0, 0). Should propagate without error.
|
|
-- ============================================================
|
|
SELECT 'hyperbolic_t0' AS test,
|
|
round(helio_x(kepler_propagate(1.0, 1.5, 0.0, 0.0, 0.0,
|
|
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS x,
|
|
round(helio_distance(kepler_propagate(1.0, 1.5, 0.0, 0.0, 0.0,
|
|
2451545.0, '2000-07-01 12:00:00+00'))::numeric, 2) AS dist_later;
|
|
|
|
-- ============================================================
|
|
-- Test 7: Near-parabolic orbit (e=1.0)
|
|
-- Comet on parabolic trajectory, q=1.0 AU.
|
|
-- At perihelion, position = (1, 0, 0).
|
|
-- ============================================================
|
|
SELECT 'parabolic_t0' AS test,
|
|
round(helio_x(kepler_propagate(1.0, 1.0, 0.0, 0.0, 0.0,
|
|
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS x;
|
|
|
|
-- ============================================================
|
|
-- Test 8: Comet observation with known Earth position
|
|
-- Place a "comet" at (2, 0, 0) AU heliocentric ecliptic.
|
|
-- Earth at (1, 0, 0) AU. Geocentric comet is at (1, 0, 0) AU,
|
|
-- which in equatorial J2000 is RA=0h, Dec=0 (on ecliptic plane
|
|
-- after rotation, dec ~ -obliquity... actually RA=0, Dec=0
|
|
-- because ecliptic-to-equatorial of (1,0,0) = (1,0,0)).
|
|
-- Geocentric distance should be 1 AU = 149597870.7 km.
|
|
-- ============================================================
|
|
SELECT 'comet_obs_range' AS test,
|
|
round(topo_range(comet_observe(
|
|
2.0, 0.0, 0.0, 0.0, 0.0, 2451545.0,
|
|
1.0, 0.0, 0.0,
|
|
:boulder, '2000-01-01 12:00:00+00'))::numeric, 0) AS range_km;
|
|
|
|
-- ============================================================
|
|
-- Test 9: Input validation
|
|
-- ============================================================
|
|
SELECT 'negative_q' AS test, kepler_propagate(-1.0, 0.5, 0.0, 0.0, 0.0, 2451545.0, now());
|