pg_orrery/test/expected/equatorial.out
Ryan Malloy b0741c553b Add v0.10.0: aberration, DE apparent, angular separation, stellar parallax
Annual stellar aberration (~20 arcsec) added to all 6 existing _apparent()
functions via classical first-order v/c projection (Ron & Vondrak). Earth
velocity sourced from VSOP87 xyz[3..5] (analytic) or DE numerical
differentiation.

New functions (106 -> 114):
- eq_angular_distance(): Vincenty formula, stable at 0 and 180 deg
- eq_within_cone(): cosine shortcut for fast cone-search predicate
- <-> operator on equatorial type
- 6 DE apparent variants with VSOP87 fallback:
  planet/sun/moon_observe_apparent_de(),
  planet/moon_equatorial_apparent_de(),
  small_body_observe_apparent_de()

Stellar parallax now functional in star_observe_pm() and
star_equatorial_pm() — Green (1985) Eq. 11.3 displacement using
Earth heliocentric position from VSOP87.

All 19 regression suites pass (18 existing + new aberration suite).
2026-02-21 21:47:42 -07:00

281 lines
13 KiB
Plaintext

-- equatorial type regression tests
--
-- Tests equatorial type I/O, accessor functions, satellite RA/Dec,
-- planet/Sun/Moon/small body equatorial coordinates, stellar
-- precession, proper motion, light-time correction, and DE variants.
\set boulder '''40.015N 105.270W 1655m'''::observer
-- ============================================================
-- Test 1: equatorial type I/O round-trip
-- ============================================================
SELECT 'eq_io' AS test,
'(12.00000000,45.00000000,0.000)'::equatorial::text AS val;
test | val
-------+---------------------------------
eq_io | (12.00000000,45.00000000,0.000)
(1 row)
-- ============================================================
-- Test 2: equatorial accessor functions
-- ============================================================
SELECT 'eq_accessors' AS test,
round(eq_ra(e)::numeric, 4) AS ra_hours,
round(eq_dec(e)::numeric, 4) AS dec_deg,
round(eq_distance(e)::numeric, 1) AS dist_km
FROM (SELECT '(6.50000000,-23.00000000,149597870.700)'::equatorial AS e) sub;
test | ra_hours | dec_deg | dist_km
--------------+----------+----------+-------------
eq_accessors | 6.5000 | -23.0000 | 149597870.7
(1 row)
-- ============================================================
-- Test 3: equatorial input validation - RA out of range
-- ============================================================
SELECT 'eq_bad_ra' AS test, '(25.0,0.0,0.0)'::equatorial;
ERROR: right ascension out of range: 25.000000
LINE 1: SELECT 'eq_bad_ra' AS test, '(25.0,0.0,0.0)'::equatorial;
^
HINT: RA must be in [0, 24) hours.
-- ============================================================
-- Test 4: equatorial input validation - Dec out of range
-- ============================================================
SELECT 'eq_bad_dec' AS test, '(12.0,91.0,0.0)'::equatorial;
ERROR: declination out of range: 91.000000
LINE 1: SELECT 'eq_bad_dec' AS test, '(12.0,91.0,0.0)'::equatorial;
^
HINT: Declination must be between -90 and +90 degrees.
-- ============================================================
-- Test 5: Sun equatorial RA/Dec at J2000.0
-- At 2000-01-01 12:00:00 UTC the Sun should be near RA ~18.75h, Dec ~-23 deg
-- (winter solstice was ~10 days earlier)
-- ============================================================
SELECT 'sun_eq' AS test,
round(eq_ra(sun_equatorial('2000-01-01 12:00:00+00'))::numeric, 1) AS ra_h,
round(eq_dec(sun_equatorial('2000-01-01 12:00:00+00'))::numeric, 0) AS dec_deg;
test | ra_h | dec_deg
--------+------+---------
sun_eq | 18.8 | -23
(1 row)
-- ============================================================
-- Test 6: Sun equatorial at summer solstice ~2024
-- RA should be ~6h, Dec ~+23.4 deg
-- ============================================================
SELECT 'sun_eq_solstice' AS test,
round(eq_ra(sun_equatorial('2024-06-21 12:00:00+00'))::numeric, 0) AS ra_h,
round(eq_dec(sun_equatorial('2024-06-21 12:00:00+00'))::numeric, 0) AS dec_deg;
test | ra_h | dec_deg
-----------------+------+---------
sun_eq_solstice | 6 | 23
(1 row)
-- ============================================================
-- Test 7: Moon equatorial returns valid coordinates
-- Range should be 356k-407k km
-- ============================================================
SELECT 'moon_eq' AS test,
eq_ra(moon_equatorial('2024-06-21 12:00:00+00')) BETWEEN 0 AND 24 AS ra_valid,
eq_dec(moon_equatorial('2024-06-21 12:00:00+00')) BETWEEN -90 AND 90 AS dec_valid,
eq_distance(moon_equatorial('2024-06-21 12:00:00+00')) BETWEEN 350000 AND 410000 AS range_valid;
test | ra_valid | dec_valid | range_valid
---------+----------+-----------+-------------
moon_eq | t | t | t
(1 row)
-- ============================================================
-- Test 8: Jupiter equatorial returns valid coordinates
-- Distance should be ~588M-968M km (3.93-6.47 AU)
-- ============================================================
SELECT 'jupiter_eq' AS test,
eq_ra(planet_equatorial(5, '2024-06-21 12:00:00+00')) BETWEEN 0 AND 24 AS ra_valid,
eq_dec(planet_equatorial(5, '2024-06-21 12:00:00+00')) BETWEEN -90 AND 90 AS dec_valid,
eq_distance(planet_equatorial(5, '2024-06-21 12:00:00+00')) BETWEEN 500000000 AND 1000000000 AS range_valid;
test | ra_valid | dec_valid | range_valid
------------+----------+-----------+-------------
jupiter_eq | t | t | t
(1 row)
-- ============================================================
-- Test 9: planet_equatorial error - cannot observe Earth
-- ============================================================
SELECT 'earth_eq_error' AS test, planet_equatorial(3, now());
ERROR: cannot observe Earth from Earth
-- ============================================================
-- Test 10: star_equatorial at J2000.0 (precession = identity)
-- Polaris RA 2.530303h Dec 89.2641 deg should be unchanged
-- ============================================================
SELECT 'star_eq_j2000' AS test,
round(eq_ra(star_equatorial(2.530303, 89.2641, '2000-01-01 12:00:00+00'))::numeric, 4) AS ra_h,
round(eq_dec(star_equatorial(2.530303, 89.2641, '2000-01-01 12:00:00+00'))::numeric, 4) AS dec_deg,
round(eq_distance(star_equatorial(2.530303, 89.2641, '2000-01-01 12:00:00+00'))::numeric, 1) AS dist;
test | ra_h | dec_deg | dist
---------------+--------+---------+------
star_eq_j2000 | 2.5303 | 89.2641 | 0.0
(1 row)
-- ============================================================
-- Test 11: star_equatorial with precession (25 years from J2000)
-- RA should shift slightly due to precession
-- ============================================================
SELECT 'star_eq_precessed' AS test,
round(eq_ra(star_equatorial(2.530303, 89.2641, '2025-06-15 12:00:00+00'))::numeric, 4) AS ra_h,
round(eq_dec(star_equatorial(2.530303, 89.2641, '2025-06-15 12:00:00+00'))::numeric, 4) AS dec_deg;
test | ra_h | dec_deg
-------------------+--------+---------
star_eq_precessed | 3.0836 | 89.3695
(1 row)
-- ============================================================
-- Test 12: star_equatorial_pm for Barnard's Star
-- Barnard's Star: RA 17.963472h, Dec 4.6933 deg
-- pm_ra = -798.58 mas/yr, pm_dec = 10328.12 mas/yr
-- parallax = 545.4 mas, rv = -110.51 km/s
-- After 25 years, Dec should shift by ~0.26 deg
-- ============================================================
SELECT 'barnard_pm' AS test,
round(eq_dec(star_equatorial_pm(
17.963472, 4.6933, -798.58, 10328.12, 545.4, -110.51,
'2025-01-01 12:00:00+00'))::numeric, 2) AS dec_deg,
round(eq_distance(star_equatorial_pm(
17.963472, 4.6933, -798.58, 10328.12, 545.4, -110.51,
'2025-01-01 12:00:00+00'))::numeric, -6) AS dist_km;
test | dec_deg | dist_km
------------+---------+----------------
barnard_pm | 4.76 | 56576466000000
(1 row)
-- ============================================================
-- Test 13: star_observe_pm returns valid topocentric
-- ============================================================
SELECT 'barnard_topo' AS test,
round(topo_azimuth(star_observe_pm(
17.963472, 4.6933, -798.58, 10328.12, 545.4, -110.51,
:boulder, '2024-07-15 04:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(star_observe_pm(
17.963472, 4.6933, -798.58, 10328.12, 545.4, -110.51,
:boulder, '2024-07-15 04:00:00+00'))::numeric, 0) AS el_deg;
test | az_deg | el_deg
--------------+--------+--------
barnard_topo | 146 | 50
(1 row)
-- ============================================================
-- Test 14: Satellite equatorial (geocentric) from ECI
-- ISS at known ECI: should get valid RA/Dec
-- ============================================================
SELECT 'sat_eq_geo' AS test,
eq_ra(eci_to_equatorial_geo(
'(-4400.594,1586.943,4772.175,2.604,7.004,2.126)'::eci_position,
'2024-06-21 12:00:00+00')) BETWEEN 0 AND 24 AS ra_valid,
eq_dec(eci_to_equatorial_geo(
'(-4400.594,1586.943,4772.175,2.604,7.004,2.126)'::eci_position,
'2024-06-21 12:00:00+00')) BETWEEN -90 AND 90 AS dec_valid;
test | ra_valid | dec_valid
------------+----------+-----------
sat_eq_geo | t | t
(1 row)
-- ============================================================
-- Test 15: Satellite equatorial (topocentric vs geocentric)
-- Topocentric should differ from geocentric by ~0.5-1 deg for LEO
-- ============================================================
SELECT 'sat_parallax' AS test,
round(abs(
eq_ra(eci_to_equatorial(
'(-4400.594,1586.943,4772.175,2.604,7.004,2.126)'::eci_position,
:boulder, '2024-06-21 12:00:00+00'))
- eq_ra(eci_to_equatorial_geo(
'(-4400.594,1586.943,4772.175,2.604,7.004,2.126)'::eci_position,
'2024-06-21 12:00:00+00'))
)::numeric, 2) AS parallax_hours;
test | parallax_hours
--------------+----------------
sat_parallax | 0.16
(1 row)
-- ============================================================
-- Test 16: All planets return valid equatorial coordinates
-- ============================================================
SELECT 'all_planets_eq' AS test,
body_id,
round(eq_ra(planet_equatorial(body_id, '2024-06-21 12:00:00+00'))::numeric, 2) AS ra_h,
round(eq_dec(planet_equatorial(body_id, '2024-06-21 12:00:00+00'))::numeric, 1) AS dec_deg
FROM generate_series(1, 8) AS body_id
WHERE body_id != 3;
test | body_id | ra_h | dec_deg
----------------+---------+-------+---------
all_planets_eq | 1 | 6.65 | 24.9
all_planets_eq | 2 | 6.38 | 23.9
all_planets_eq | 4 | 2.47 | 13.5
all_planets_eq | 5 | 4.29 | 20.6
all_planets_eq | 6 | 23.40 | -6.0
all_planets_eq | 7 | 3.53 | 18.8
all_planets_eq | 8 | 0.03 | -1.2
(7 rows)
-- ============================================================
-- Test 17: planet_equatorial_apparent (light-time corrected)
-- Jupiter's light-time is ~35-52 min; apparent RA should differ
-- from geometric by a small amount
-- ============================================================
SELECT 'jupiter_apparent' AS test,
round(abs(
eq_ra(planet_equatorial_apparent(5, '2024-06-21 12:00:00+00'))
- eq_ra(planet_equatorial(5, '2024-06-21 12:00:00+00'))
)::numeric, 4) AS ra_diff_hours;
test | ra_diff_hours
------------------+---------------
jupiter_apparent | 0.0005
(1 row)
-- ============================================================
-- Test 18: moon_equatorial_apparent (light-time ~1.3 sec)
-- Difference from geometric should be tiny
-- ============================================================
SELECT 'moon_apparent' AS test,
round(abs(
eq_ra(moon_equatorial_apparent('2024-06-21 12:00:00+00'))
- eq_ra(moon_equatorial('2024-06-21 12:00:00+00'))
)::numeric, 6) AS ra_diff_hours;
test | ra_diff_hours
---------------+---------------
moon_apparent | 0.000405
(1 row)
-- ============================================================
-- Test 19: Small body equatorial with Ceres
-- Ceres orbital elements (approximate)
-- ============================================================
SELECT 'ceres_eq' AS test,
eq_ra(small_body_equatorial(
'(2460400.5,2.5577,0.0785,0.1849,1.2836,1.4013,2460500.0,3.53,0.12)'::orbital_elements,
'2024-06-21 12:00:00+00')) BETWEEN 0 AND 24 AS ra_valid;
test | ra_valid
----------+----------
ceres_eq | t
(1 row)
-- ============================================================
-- Test 20: DE equatorial variants (fall back to VSOP87)
-- Without DE configured, these should produce same results as VSOP87
-- ============================================================
SELECT 'de_planet_eq' AS test,
round(eq_ra(planet_equatorial_de(5, '2024-06-21 12:00:00+00'))::numeric, 4) AS ra_h,
round(eq_ra(planet_equatorial(5, '2024-06-21 12:00:00+00'))::numeric, 4) AS ra_h_vsop,
round(eq_ra(planet_equatorial_de(5, '2024-06-21 12:00:00+00'))::numeric, 4) =
round(eq_ra(planet_equatorial(5, '2024-06-21 12:00:00+00'))::numeric, 4) AS match;
test | ra_h | ra_h_vsop | match
--------------+--------+-----------+-------
de_planet_eq | 4.2922 | 4.2922 | t
(1 row)
SELECT 'de_moon_eq' AS test,
round(eq_ra(moon_equatorial_de('2024-06-21 12:00:00+00'))::numeric, 4) AS ra_h,
round(eq_ra(moon_equatorial('2024-06-21 12:00:00+00'))::numeric, 4) AS ra_h_elp,
round(eq_ra(moon_equatorial_de('2024-06-21 12:00:00+00'))::numeric, 4) =
round(eq_ra(moon_equatorial('2024-06-21 12:00:00+00'))::numeric, 4) AS match;
test | ra_h | ra_h_elp | match
------------+---------+----------+-------
de_moon_eq | 17.5281 | 17.5281 | t
(1 row)