From 22b272fd0c51536601257dd054b3c6bc17faaead Mon Sep 17 00:00:00 2001 From: Ryan Malloy Date: Thu, 26 Feb 2026 18:47:30 -0700 Subject: [PATCH] Implement v0.17.0: solar elongation, planet phase, satellite eclipse, observing night quality, lunar libration MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 162 → 174 SQL objects, 27 → 28 test suites, 3 new C source files. Features: - solar_elongation(body_id, ts): Sun-Earth-Planet angle [0,180] degrees - planet_phase(body_id, ts): illuminated disk fraction [0,1] - satellite_is_eclipsed/next_eclipse_entry/exit/eclipse_fraction: cylindrical shadow model (Vallado §5.3) for Earth shadow prediction - observing_night_quality(observer, ts): composite PL/pgSQL scoring based on astronomical darkness duration and Moon interference - moon_libration_longitude/latitude/position_angle/libration/subsolar_longitude: optical libration from Meeus (1998) Ch. 53 Refactored magnitude_funcs.c to extract shared compute_planet_geometry() used by magnitude, elongation, and phase — single VSOP87 evaluation per call. All 28 regression suites pass. Zero compiler warnings. --- CLAUDE.md | 28 +- Makefile | 9 +- pg_orrery.control | 2 +- sql/pg_orrery--0.16.0--0.17.0.sql | 139 +++ sql/pg_orrery--0.17.0.sql | 1813 +++++++++++++++++++++++++++++ src/eclipse_funcs.c | 362 ++++++ src/libration.h | 22 + src/libration_funcs.c | 368 ++++++ src/magnitude_funcs.c | 180 ++- test/expected/v016_features.out | 2 +- test/expected/v017_features.out | 285 +++++ test/sql/v017_features.sql | 204 ++++ 12 files changed, 3357 insertions(+), 57 deletions(-) create mode 100644 sql/pg_orrery--0.16.0--0.17.0.sql create mode 100644 sql/pg_orrery--0.17.0.sql create mode 100644 src/eclipse_funcs.c create mode 100644 src/libration.h create mode 100644 src/libration_funcs.c create mode 100644 test/expected/v017_features.out create mode 100644 test/sql/v017_features.sql diff --git a/CLAUDE.md b/CLAUDE.md index 977bd03..08c7fbc 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -1,9 +1,9 @@ # pg_orrery — A Database Orrery for PostgreSQL ## What This Is -A database orrery — celestial mechanics types and functions for PostgreSQL. Native C extension using PGXS, 162 SQL objects (146 user-visible functions + 16 GiST support), 9 custom types, covering satellites (SGP4/SDP4), planets (VSOP87 + optional JPL DE441), Moon (ELP2000-82B), 19 planetary moons (L1.2/TASS17/GUST86/MarsSat), stars (with proper motion and annual parallax), comets, asteroids (MPC catalog), Jupiter radio bursts, interplanetary Lambert transfers, equatorial RA/Dec coordinates with GiST-indexed angular separation, atmospheric refraction, annual stellar aberration, light-time correction, rise/set prediction (geometric + refracted) with status diagnostics, IAU constellation identification with full name lookup (Roman 1987), twilight dawn/dusk (civil/nautical/astronomical), lunar phase (angle, illumination, name, age), and planet apparent magnitude (Mallama & Hilton 2018). +A database orrery — celestial mechanics types and functions for PostgreSQL. Native C extension using PGXS, 174 SQL objects (158 user-visible functions + 16 GiST support), 9 custom types, covering satellites (SGP4/SDP4), planets (VSOP87 + optional JPL DE441), Moon (ELP2000-82B), 19 planetary moons (L1.2/TASS17/GUST86/MarsSat), stars (with proper motion and annual parallax), comets, asteroids (MPC catalog), Jupiter radio bursts, interplanetary Lambert transfers, equatorial RA/Dec coordinates with GiST-indexed angular separation, atmospheric refraction, annual stellar aberration, light-time correction, rise/set prediction (geometric + refracted) with status diagnostics, IAU constellation identification with full name lookup (Roman 1987), twilight dawn/dusk (civil/nautical/astronomical), lunar phase (angle, illumination, name, age), planet apparent magnitude (Mallama & Hilton 2018), solar elongation, planet phase fraction, satellite eclipse prediction (cylindrical shadow), observing night quality assessment, and lunar optical libration (Meeus Ch. 53). -**Current version:** 0.16.0 +**Current version:** 0.17.0 **Repository:** https://git.supported.systems/warehack.ing/pg_orrery **Documentation:** https://pg-orrery.warehack.ing @@ -11,7 +11,7 @@ A database orrery — celestial mechanics types and functions for PostgreSQL. Na ```bash make PG_CONFIG=/usr/bin/pg_config # Compile with PGXS sudo make install PG_CONFIG=/usr/bin/pg_config # Install extension -make installcheck PG_CONFIG=/usr/bin/pg_config # Run 27 regression test suites +make installcheck PG_CONFIG=/usr/bin/pg_config # Run 28 regression test suites ``` Requires: PostgreSQL 17 development headers, GCC, Make. @@ -27,7 +27,7 @@ Image: `git.supported.systems/warehack.ing/pg_orrery:pg17` ## Project Layout ``` -pg_orrery.control # Extension metadata (version 0.16.0) +pg_orrery.control # Extension metadata (version 0.17.0) Makefile # PGXS build + Docker targets sql/ pg_orrery--0.1.0.sql # v0.1.0: satellite types/functions/operators @@ -46,6 +46,7 @@ sql/ pg_orrery--0.14.0.sql # v0.14.0: refracted rise/set, constellation ID (147 objects) pg_orrery--0.15.0.sql # v0.15.0: constellation full name, rise/set status (151 objects) pg_orrery--0.16.0.sql # v0.16.0: twilight, lunar phase, planet magnitude (162 objects) + pg_orrery--0.17.0.sql # v0.17.0: elongation, phase, eclipse, night quality, libration (174 objects) pg_orrery--0.1.0--0.2.0.sql # Migration: v0.1.0 → v0.2.0 (adds solar system) pg_orrery--0.2.0--0.3.0.sql # Migration: v0.2.0 → v0.3.0 (adds DE ephemeris) pg_orrery--0.3.0--0.4.0.sql # Migration: v0.3.0 → v0.4.0 @@ -61,6 +62,7 @@ sql/ pg_orrery--0.13.0--0.14.0.sql # Migration: v0.13.0 → v0.14.0 (refracted rise/set, constellation ID) pg_orrery--0.14.0--0.15.0.sql # Migration: v0.14.0 → v0.15.0 (constellation full name, rise/set status) pg_orrery--0.15.0--0.16.0.sql # Migration: v0.15.0 → v0.16.0 (twilight, lunar phase, planet magnitude) + pg_orrery--0.16.0--0.17.0.sql # Migration: v0.16.0 → v0.17.0 (elongation, phase, eclipse, night quality, libration) src/ pg_orrery.c # PG_MODULE_MAGIC + _PG_init() (GUC registration) types.h # All struct definitions + constants + DE body ID mapping @@ -91,7 +93,9 @@ src/ constellation_data.h / .c # Roman (1987) IAU boundary table (CDS VI/42, 357 segments) constellation_funcs.c # constellation() from equatorial or RA/Dec lunar_phase_funcs.c # moon_phase_angle(), moon_illumination(), moon_phase_name(), moon_age() - magnitude_funcs.c # planet_magnitude() (Mallama & Hilton 2018) + magnitude_funcs.c # planet_magnitude(), solar_elongation(), planet_phase() + eclipse_funcs.c # satellite eclipse prediction (cylindrical shadow, Vallado §5.3) + libration.h / libration_funcs.c # lunar optical libration (Meeus Ch. 53) l12.c / l12.h # L1.2 Galilean moon theory (Lieske 1998) tass17.c / tass17.h # TASS 1.7 Saturn moon theory (Vienne & Duriez 1995) gust86.c / gust86.h # GUST86 Uranus moon theory (Laskar & Jacobson 1987) @@ -143,7 +147,7 @@ All types are fixed-size, `STORAGE = plain`, `ALIGNMENT = double`. No TOAST over | `orbital_elements` | 72 | Classical Keplerian elements for comets/asteroids (epoch, q, e, inc, omega, Omega, tp, H, G) | | `equatorial` | 24 | Apparent RA (hours), Dec (degrees), distance (km) — of date | -## Function Domains (162 SQL objects) +## Function Domains (174 SQL objects) | Domain | Theory | Key Functions | Count | |--------|--------|---------------|-------| @@ -164,6 +168,11 @@ All types are fixed-size, `STORAGE = plain`, `ALIGNMENT = double`. No TOAST over | Twilight | Sun depression angles | `sun_civil_dawn()`, `sun_nautical_dusk()`, `sun_astronomical_dawn()` | 6 | | Lunar phase | VSOP87 + ELP2000-82B geometry | `moon_phase_angle()`, `moon_illumination()`, `moon_phase_name()`, `moon_age()` | 4 | | Planet magnitude | Mallama & Hilton (2018) | `planet_magnitude()` | 1 | +| Solar elongation | VSOP87 geometry | `solar_elongation()` | 1 | +| Planet phase | VSOP87 geometry | `planet_phase()` | 1 | +| Satellite eclipse | Cylindrical shadow (Vallado §5.3) | `satellite_is_eclipsed()`, `satellite_next_eclipse_entry()` | 4 | +| Observing quality | Composite (twilight+Moon) | `observing_night_quality()` | 1 | +| Lunar libration | Meeus (1998) Ch. 53 | `moon_libration_longitude()`, `moon_libration()`, `moon_subsolar_longitude()` | 5 | | Constellation | Roman (1987) CDS VI/42 | `constellation()`, `constellation_full_name()` | 3 | | Diagnostics | -- | `pg_orrery_ephemeris_info()` | 1 | @@ -298,7 +307,7 @@ All numerical logic is byte-identical to upstream. Verified against 518 Vallado ## Testing -27 regression test suites via `make installcheck`: +28 regression test suites via `make installcheck`: | Suite | What it tests | |-------|--------------| @@ -329,10 +338,11 @@ All numerical logic is byte-identical to upstream. Verified against 518 Vallado | constellation | Roman (1987) boundary lookup, known stars, solar system objects, edge cases | | v015_features | constellation_full_name lookup, rise_set_status diagnostics (circumpolar/never_rises) | | v016_features | Twilight ordering/offset/polar, lunar phase at known events, planet magnitude ranges/errors | +| v017_features | Solar elongation ranges/errors, planet phase ranges, satellite eclipse, observing night quality, lunar libration ranges, subsolar longitude | ### PG Version Matrix -Test all 27 regression suites + DE reader unit test across PostgreSQL 14-18 using Docker: +Test all 28 regression suites + DE reader unit test across PostgreSQL 14-18 using Docker: ```bash make test-matrix # Full matrix (PG 14-18) @@ -358,7 +368,7 @@ Logs saved to `test/matrix-logs/pg${ver}.log`. The script reuses the Dockerfile Starlight docs at `docs/` — 44+ MDX pages covering all domains. -Sections: Getting Started, Guides (9 domain walkthroughs incl. DE ephemeris), Workflow Translation (Skyfield/Horizons/GMAT/Radio Jupiter Pro comparisons), Reference (all 162 SQL objects incl. DE variants, equatorial GiST, refraction, rise/set, constellation, twilight, lunar phase, planet magnitude), Architecture (Hamilton's principles, constant custody, observation pipeline), Performance (benchmarks). +Sections: Getting Started, Guides (9 domain walkthroughs incl. DE ephemeris), Workflow Translation (Skyfield/Horizons/GMAT/Radio Jupiter Pro comparisons), Reference (all 174 SQL objects incl. DE variants, equatorial GiST, refraction, rise/set, constellation, twilight, lunar phase, planet magnitude, solar elongation, planet phase, satellite eclipse, observing quality, lunar libration), Architecture (Hamilton's principles, constant custody, observation pipeline), Performance (benchmarks). ### Local Development ```bash diff --git a/Makefile b/Makefile index 6d7e773..ccfd5c3 100644 --- a/Makefile +++ b/Makefile @@ -14,7 +14,8 @@ DATA = sql/pg_orrery--0.1.0.sql sql/pg_orrery--0.2.0.sql sql/pg_orrery--0.1.0--0 sql/pg_orrery--0.13.0.sql sql/pg_orrery--0.12.0--0.13.0.sql \ sql/pg_orrery--0.14.0.sql sql/pg_orrery--0.13.0--0.14.0.sql \ sql/pg_orrery--0.15.0.sql sql/pg_orrery--0.14.0--0.15.0.sql \ - sql/pg_orrery--0.16.0.sql sql/pg_orrery--0.15.0--0.16.0.sql + sql/pg_orrery--0.16.0.sql sql/pg_orrery--0.15.0--0.16.0.sql \ + sql/pg_orrery--0.17.0.sql sql/pg_orrery--0.16.0--0.17.0.sql # Our extension C sources OBJS = src/pg_orrery.o src/tle_type.o src/eci_type.o src/observer_type.o \ @@ -34,7 +35,8 @@ OBJS = src/pg_orrery.o src/tle_type.o src/eci_type.o src/observer_type.o \ src/gist_equatorial.o \ src/rise_set_funcs.o \ src/constellation_data.o src/constellation_funcs.o \ - src/lunar_phase_funcs.o src/magnitude_funcs.o + src/lunar_phase_funcs.o src/magnitude_funcs.o \ + src/eclipse_funcs.o src/libration_funcs.o # Vendored SGP4/SDP4 sources (pure C, from Bill Gray's sat_code, MIT license) SGP4_DIR = src/sgp4 @@ -55,7 +57,8 @@ REGRESS = tle_parse sgp4_propagate coord_transforms pass_prediction gist_index c v013_features rise_set \ constellation \ v015_features \ - v016_features + v016_features \ + v017_features REGRESS_OPTS = --inputdir=test # Pure C — no C++ runtime needed. LAPACK for OD solver (dgelss_). diff --git a/pg_orrery.control b/pg_orrery.control index 9b863f8..f573a6d 100644 --- a/pg_orrery.control +++ b/pg_orrery.control @@ -1,4 +1,4 @@ comment = 'A database orrery — celestial mechanics types and functions for PostgreSQL' -default_version = '0.16.0' +default_version = '0.17.0' module_pathname = '$libdir/pg_orrery' relocatable = true diff --git a/sql/pg_orrery--0.16.0--0.17.0.sql b/sql/pg_orrery--0.16.0--0.17.0.sql new file mode 100644 index 0000000..d63abd0 --- /dev/null +++ b/sql/pg_orrery--0.16.0--0.17.0.sql @@ -0,0 +1,139 @@ +-- pg_orrery 0.16.0 -> 0.17.0: solar elongation, planet phase, satellite eclipse, +-- observing night quality, lunar libration + +-- ============================================================ +-- Solar elongation (1) +-- ============================================================ + +CREATE FUNCTION solar_elongation(int4, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'solar_elongation' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION solar_elongation(int4, timestamptz) IS + 'Sun-Earth-Planet angle in degrees [0, 180]. How far a planet appears from the Sun. Body IDs 1-8.'; + +-- ============================================================ +-- Planet phase fraction (1) +-- ============================================================ + +CREATE FUNCTION planet_phase(int4, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'planet_phase' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_phase(int4, timestamptz) IS + 'Illuminated fraction of a planet disk as seen from Earth [0.0, 1.0]. Body IDs 1-8.'; + +-- ============================================================ +-- Satellite eclipse prediction (4) +-- ============================================================ + +CREATE FUNCTION satellite_is_eclipsed(tle, timestamptz) RETURNS bool + AS 'MODULE_PATHNAME', 'satellite_is_eclipsed' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION satellite_is_eclipsed(tle, timestamptz) IS + 'True if the satellite is in Earth cylindrical shadow at the given time.'; + +CREATE FUNCTION satellite_next_eclipse_entry(tle, timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME', 'satellite_next_eclipse_entry' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION satellite_next_eclipse_entry(tle, timestamptz) IS + 'Next time the satellite enters Earth shadow (up to 7-day search). NULL if none found.'; + +CREATE FUNCTION satellite_next_eclipse_exit(tle, timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME', 'satellite_next_eclipse_exit' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION satellite_next_eclipse_exit(tle, timestamptz) IS + 'Next time the satellite exits Earth shadow (up to 7-day search). NULL if none found.'; + +CREATE FUNCTION satellite_eclipse_fraction(tle, timestamptz, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'satellite_eclipse_fraction' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION satellite_eclipse_fraction(tle, timestamptz, timestamptz) IS + 'Fraction of the given time window the satellite spends in eclipse [0.0, 1.0].'; + +-- ============================================================ +-- Observing night quality (1) +-- ============================================================ + +CREATE FUNCTION observing_night_quality(observer, timestamptz DEFAULT NOW()) +RETURNS text AS $$ +DECLARE + astro_dusk timestamptz; + astro_dawn timestamptz; + dark_hours float8; + illum float8; + moon_up bool; + score int := 100; +BEGIN + -- Astronomical darkness window + astro_dusk := sun_astronomical_dusk($1, $2); + IF astro_dusk IS NULL THEN + RETURN 'poor'; -- No astronomical darkness (polar summer) + END IF; + astro_dawn := sun_astronomical_dawn($1, astro_dusk); + IF astro_dawn IS NULL THEN + RETURN 'poor'; + END IF; + + dark_hours := extract(epoch FROM astro_dawn - astro_dusk) / 3600.0; + + -- Short dark window penalty + IF dark_hours < 2.0 THEN score := score - 40; + ELSIF dark_hours < 4.0 THEN score := score - 20; + ELSIF dark_hours < 6.0 THEN score := score - 10; + END IF; + + -- Moon illumination penalty + illum := moon_illumination(astro_dusk); + IF illum > 0.75 THEN + -- Check if Moon is above horizon during darkness + moon_up := (moon_observe($1, astro_dusk)).elevation > 0 + OR (moon_observe($1, astro_dusk + (astro_dawn - astro_dusk) / 2)).elevation > 0; + IF moon_up THEN + score := score - (illum * 30)::int; -- Up to -30 for full moon + END IF; + END IF; + + -- Classify + IF score >= 80 THEN RETURN 'excellent'; + ELSIF score >= 60 THEN RETURN 'good'; + ELSIF score >= 40 THEN RETURN 'fair'; + ELSE RETURN 'poor'; + END IF; +END; +$$ LANGUAGE plpgsql STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION observing_night_quality(observer, timestamptz) IS + 'Composite observing quality assessment: excellent/good/fair/poor based on darkness duration and Moon interference.'; + +-- ============================================================ +-- Lunar libration (5) +-- ============================================================ + +CREATE FUNCTION moon_libration_longitude(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'moon_libration_longitude' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_libration_longitude(timestamptz) IS + 'Optical libration in longitude (degrees, typically [-8, +8]). Meeus Ch. 53.'; + +CREATE FUNCTION moon_libration_latitude(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'moon_libration_latitude' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_libration_latitude(timestamptz) IS + 'Optical libration in latitude (degrees, typically [-7, +7]). Meeus Ch. 53.'; + +CREATE FUNCTION moon_libration_position_angle(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'moon_libration_position_angle' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_libration_position_angle(timestamptz) IS + 'Position angle of the Moon axis (degrees, [0, 360)). Meeus Ch. 53.'; + +CREATE FUNCTION moon_libration(timestamptz, + OUT l float8, OUT b float8, OUT p float8) RETURNS record + AS 'MODULE_PATHNAME', 'moon_libration' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_libration(timestamptz) IS + 'All three libration values: longitude (l), latitude (b), position angle (p) in degrees.'; + +CREATE FUNCTION moon_subsolar_longitude(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'moon_subsolar_longitude' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_subsolar_longitude(timestamptz) IS + 'Selenographic longitude of the sub-solar point (degrees, [0, 360)). Determines the lunar terminator position.'; diff --git a/sql/pg_orrery--0.17.0.sql b/sql/pg_orrery--0.17.0.sql new file mode 100644 index 0000000..6585b5e --- /dev/null +++ b/sql/pg_orrery--0.17.0.sql @@ -0,0 +1,1813 @@ +-- pg_orrery -- Orbital mechanics types and functions for PostgreSQL +-- +-- Types: tle, eci_position, geodetic, topocentric, observer, pass_event +-- Provides SGP4/SDP4 propagation, coordinate transforms, pass prediction, +-- and GiST indexing on altitude bands for conjunction screening. +-- +-- All propagation uses WGS-72 constants (matching TLE mean element fitting). +-- Coordinate output uses WGS-84 (matching modern geodetic standards). + +-- ============================================================ +-- Shell types (forward declarations) +-- ============================================================ + +CREATE TYPE tle; +CREATE TYPE eci_position; +CREATE TYPE geodetic; +CREATE TYPE topocentric; +CREATE TYPE observer; +CREATE TYPE pass_event; + + +-- ============================================================ +-- TLE type: Two-Line Element set +-- ============================================================ + +CREATE FUNCTION tle_in(cstring) RETURNS tle + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION tle_out(tle) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION tle_recv(internal) RETURNS tle + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION tle_send(tle) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE tle ( + INPUT = tle_in, + OUTPUT = tle_out, + RECEIVE = tle_recv, + SEND = tle_send, + INTERNALLENGTH = 112, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE tle IS 'Two-Line Element set — parsed mean orbital elements for SGP4/SDP4 propagation'; + +-- TLE accessor functions + +CREATE FUNCTION tle_epoch(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_epoch(tle) IS 'TLE epoch as Julian date (UTC)'; + +CREATE FUNCTION tle_norad_id(tle) RETURNS int4 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_norad_id(tle) IS 'NORAD catalog number'; + +CREATE FUNCTION tle_inclination(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_inclination(tle) IS 'Orbital inclination in degrees'; + +CREATE FUNCTION tle_eccentricity(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_eccentricity(tle) IS 'Orbital eccentricity (dimensionless)'; + +CREATE FUNCTION tle_raan(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_raan(tle) IS 'Right ascension of ascending node in degrees'; + +CREATE FUNCTION tle_arg_perigee(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_arg_perigee(tle) IS 'Argument of perigee in degrees'; + +CREATE FUNCTION tle_mean_anomaly(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_mean_anomaly(tle) IS 'Mean anomaly in degrees'; + +CREATE FUNCTION tle_mean_motion(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_mean_motion(tle) IS 'Mean motion in revolutions per day'; + +CREATE FUNCTION tle_bstar(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_bstar(tle) IS 'B* drag coefficient (1/earth-radii)'; + +CREATE FUNCTION tle_period(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_period(tle) IS 'Orbital period in minutes'; + +CREATE FUNCTION tle_age(tle, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_age(tle, timestamptz) IS 'TLE age in days (positive = past epoch, negative = before epoch)'; + +CREATE FUNCTION tle_perigee(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_perigee(tle) IS 'Perigee altitude in km above WGS-72 ellipsoid'; + +CREATE FUNCTION tle_apogee(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_apogee(tle) IS 'Apogee altitude in km above WGS-72 ellipsoid'; + +CREATE FUNCTION tle_intl_desig(tle) RETURNS text + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_intl_desig(tle) IS 'International designator (COSPAR ID)'; + +CREATE FUNCTION tle_from_lines(text, text) RETURNS tle + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_from_lines(text, text) IS + 'Construct TLE from separate line1/line2 text columns'; + + +-- ============================================================ +-- ECI position type: True Equator Mean Equinox (TEME) frame +-- ============================================================ + +CREATE FUNCTION eci_in(cstring) RETURNS eci_position + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_out(eci_position) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_recv(internal) RETURNS eci_position + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_send(eci_position) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE eci_position ( + INPUT = eci_in, + OUTPUT = eci_out, + RECEIVE = eci_recv, + SEND = eci_send, + INTERNALLENGTH = 48, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE eci_position IS 'Earth-Centered Inertial position and velocity in TEME frame (km, km/s)'; + +-- ECI accessor functions + +CREATE FUNCTION eci_x(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_y(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_z(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_vx(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_vy(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_vz(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_speed(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eci_speed(eci_position) IS 'Velocity magnitude in km/s'; + +CREATE FUNCTION eci_altitude(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eci_altitude(eci_position) IS 'Approximate geocentric altitude in km (radius - WGS72_AE)'; + + +-- ============================================================ +-- Geodetic type: WGS-84 latitude/longitude/altitude +-- ============================================================ + +CREATE FUNCTION geodetic_in(cstring) RETURNS geodetic + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION geodetic_out(geodetic) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION geodetic_recv(internal) RETURNS geodetic + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION geodetic_send(geodetic) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE geodetic ( + INPUT = geodetic_in, + OUTPUT = geodetic_out, + RECEIVE = geodetic_recv, + SEND = geodetic_send, + INTERNALLENGTH = 24, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE geodetic IS 'Geodetic coordinates on WGS-84 ellipsoid (lat/lon in degrees, altitude in km)'; + +CREATE FUNCTION geodetic_lat(geodetic) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION geodetic_lon(geodetic) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION geodetic_alt(geodetic) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + + +-- ============================================================ +-- Topocentric type: observer-relative az/el/range +-- ============================================================ + +CREATE FUNCTION topocentric_in(cstring) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION topocentric_out(topocentric) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION topocentric_recv(internal) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION topocentric_send(topocentric) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE topocentric ( + INPUT = topocentric_in, + OUTPUT = topocentric_out, + RECEIVE = topocentric_recv, + SEND = topocentric_send, + INTERNALLENGTH = 32, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE topocentric IS 'Topocentric coordinates relative to observer (azimuth, elevation, range, range rate)'; + +CREATE FUNCTION topo_azimuth(topocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION topo_azimuth(topocentric) IS 'Azimuth in degrees (0=N, 90=E, 180=S, 270=W)'; + +CREATE FUNCTION topo_elevation(topocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION topo_elevation(topocentric) IS 'Elevation in degrees (0=horizon, 90=zenith)'; + +CREATE FUNCTION topo_range(topocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION topo_range(topocentric) IS 'Slant range in km'; + +CREATE FUNCTION topo_range_rate(topocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION topo_range_rate(topocentric) IS 'Range rate in km/s (positive = receding)'; + + +-- ============================================================ +-- Observer type: ground station location +-- ============================================================ + +CREATE FUNCTION observer_in(cstring) RETURNS observer + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION observer_out(observer) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION observer_recv(internal) RETURNS observer + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION observer_send(observer) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE observer ( + INPUT = observer_in, + OUTPUT = observer_out, + RECEIVE = observer_recv, + SEND = observer_send, + INTERNALLENGTH = 24, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE observer IS 'Ground observer location (accepts: 40.0N 105.3W 1655m or decimal degrees)'; + +CREATE FUNCTION observer_lat(observer) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION observer_lat(observer) IS 'Latitude in degrees (positive = North)'; + +CREATE FUNCTION observer_lon(observer) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION observer_lon(observer) IS 'Longitude in degrees (positive = East)'; + +CREATE FUNCTION observer_alt(observer) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION observer_alt(observer) IS 'Altitude in meters above WGS-84 ellipsoid'; + +CREATE FUNCTION observer_from_geodetic(float8, float8, float8 DEFAULT 0.0) RETURNS observer + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION observer_from_geodetic(float8, float8, float8) IS + 'Construct observer from lat (deg), lon (deg), altitude (meters). Avoids text formatting round-trips.'; + + +-- ============================================================ +-- Pass event type: satellite visibility window +-- ============================================================ + +CREATE FUNCTION pass_event_in(cstring) RETURNS pass_event + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION pass_event_out(pass_event) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION pass_event_recv(internal) RETURNS pass_event + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION pass_event_send(pass_event) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE pass_event ( + INPUT = pass_event_in, + OUTPUT = pass_event_out, + RECEIVE = pass_event_recv, + SEND = pass_event_send, + INTERNALLENGTH = 48, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE pass_event IS 'Satellite pass event (AOS/MAX/LOS times, max elevation, AOS/LOS azimuths)'; + +CREATE FUNCTION pass_aos_time(pass_event) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_aos_time(pass_event) IS 'Acquisition of signal time'; + +CREATE FUNCTION pass_max_el_time(pass_event) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_max_el_time(pass_event) IS 'Maximum elevation time'; + +CREATE FUNCTION pass_los_time(pass_event) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_los_time(pass_event) IS 'Loss of signal time'; + +CREATE FUNCTION pass_max_elevation(pass_event) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_max_elevation(pass_event) IS 'Maximum elevation in degrees'; + +CREATE FUNCTION pass_aos_azimuth(pass_event) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_aos_azimuth(pass_event) IS 'AOS azimuth in degrees (0=N)'; + +CREATE FUNCTION pass_los_azimuth(pass_event) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_los_azimuth(pass_event) IS 'LOS azimuth in degrees (0=N)'; + +CREATE FUNCTION pass_duration(pass_event) RETURNS interval + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_duration(pass_event) IS 'Pass duration (LOS - AOS)'; + + +-- ============================================================ +-- SGP4/SDP4 propagation functions +-- ============================================================ + +CREATE FUNCTION sgp4_propagate(tle, timestamptz) RETURNS eci_position + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sgp4_propagate(tle, timestamptz) IS + 'Propagate TLE to a point in time using SGP4 (near-earth) or SDP4 (deep-space). Returns TEME ECI position/velocity.'; + +CREATE FUNCTION sgp4_propagate_safe(tle, timestamptz) RETURNS eci_position + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE PARALLEL SAFE; +COMMENT ON FUNCTION sgp4_propagate_safe(tle, timestamptz) IS + 'Like sgp4_propagate but returns NULL on error instead of raising an exception. For batch queries with potentially invalid TLEs.'; + +CREATE FUNCTION sgp4_propagate_series(tle, timestamptz, timestamptz, interval) + RETURNS TABLE(t timestamptz, x float8, y float8, z float8, vx float8, vy float8, vz float8) + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE + ROWS 100; +COMMENT ON FUNCTION sgp4_propagate_series(tle, timestamptz, timestamptz, interval) IS + 'Propagate TLE over a time range at regular intervals. Returns time series of TEME positions.'; + +CREATE FUNCTION tle_distance(tle, tle, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_distance(tle, tle, timestamptz) IS + 'Euclidean distance in km between two TLEs at a reference time'; + + +-- ============================================================ +-- Coordinate transform functions +-- ============================================================ + +CREATE FUNCTION eci_to_geodetic(eci_position, timestamptz) RETURNS geodetic + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eci_to_geodetic(eci_position, timestamptz) IS + 'Convert TEME ECI position to WGS-84 geodetic coordinates at given time'; + +CREATE FUNCTION eci_to_topocentric(eci_position, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eci_to_topocentric(eci_position, observer, timestamptz) IS + 'Convert TEME ECI position to topocentric (az/el/range) relative to observer'; + +CREATE FUNCTION subsatellite_point(tle, timestamptz) RETURNS geodetic + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION subsatellite_point(tle, timestamptz) IS + 'Subsatellite (nadir) point on WGS-84 ellipsoid at given time'; + +CREATE FUNCTION ground_track(tle, timestamptz, timestamptz, interval) + RETURNS TABLE(t timestamptz, lat float8, lon float8, alt float8) + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE + ROWS 100; +COMMENT ON FUNCTION ground_track(tle, timestamptz, timestamptz, interval) IS + 'Ground track as time series of subsatellite points (lat/lon in degrees, alt in km)'; + +CREATE FUNCTION observe(tle, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION observe(tle, observer, timestamptz) IS + 'Propagate TLE and compute observer-relative look angles in one call. Returns topocentric (az/el/range/range_rate).'; + +CREATE FUNCTION observe_safe(tle, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE PARALLEL SAFE; +COMMENT ON FUNCTION observe_safe(tle, observer, timestamptz) IS + 'Like observe() but returns NULL on propagation error. For batch queries with potentially invalid/decayed TLEs.'; + + +-- ============================================================ +-- Pass prediction functions +-- ============================================================ + +CREATE FUNCTION next_pass(tle, observer, timestamptz) RETURNS pass_event + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION next_pass(tle, observer, timestamptz) IS + 'Find the next satellite pass over observer (searches up to 7 days ahead)'; + +CREATE FUNCTION predict_passes(tle, observer, timestamptz, timestamptz, float8 DEFAULT 0.0) + RETURNS SETOF pass_event + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE + ROWS 10; +COMMENT ON FUNCTION predict_passes(tle, observer, timestamptz, timestamptz, float8) IS + 'Predict all satellite passes over observer in time window. Optional min_elevation in degrees.'; + +CREATE FUNCTION pass_visible(tle, observer, timestamptz, timestamptz) RETURNS boolean + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_visible(tle, observer, timestamptz, timestamptz) IS + 'True if any pass occurs over observer in the time window'; + + +-- ============================================================ +-- GiST operator support functions +-- ============================================================ + +-- Overlap operator: do orbital keys overlap in altitude AND inclination? +CREATE FUNCTION tle_overlap(tle, tle) RETURNS boolean + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +-- Altitude distance operator (altitude-only, for KNN ordering) +CREATE FUNCTION tle_alt_distance(tle, tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OPERATOR && ( + LEFTARG = tle, + RIGHTARG = tle, + FUNCTION = tle_overlap, + COMMUTATOR = &&, + RESTRICT = areasel, + JOIN = areajoinsel +); + +COMMENT ON OPERATOR && (tle, tle) IS 'Orbital key overlap (altitude band AND inclination range) — necessary condition for conjunction'; + +CREATE OPERATOR <-> ( + LEFTARG = tle, + RIGHTARG = tle, + FUNCTION = tle_alt_distance, + COMMUTATOR = <-> +); + +COMMENT ON OPERATOR <-> (tle, tle) IS '2-D orbital distance in km: L2 norm of altitude-band gap and inclination gap (radians × Earth radius). Returns 0 when both dimensions overlap.'; + + +-- ============================================================ +-- GiST operator class for 2-D orbital indexing (altitude + inclination) +-- ============================================================ + +-- GiST internal support functions +CREATE FUNCTION gist_tle_compress(internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_decompress(internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_consistent(internal, tle, smallint, oid, internal) RETURNS boolean + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_union(internal, internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_penalty(internal, internal, internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_picksplit(internal, internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_same(internal, internal, internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_distance(internal, tle, smallint, oid, internal) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OPERATOR CLASS tle_ops + DEFAULT FOR TYPE tle USING gist AS + OPERATOR 3 && , + OPERATOR 15 <-> (tle, tle) FOR ORDER BY float_ops, + FUNCTION 1 gist_tle_consistent(internal, tle, smallint, oid, internal), + FUNCTION 2 gist_tle_union(internal, internal), + FUNCTION 3 gist_tle_compress(internal), + FUNCTION 4 gist_tle_decompress(internal), + FUNCTION 5 gist_tle_penalty(internal, internal, internal), + FUNCTION 6 gist_tle_picksplit(internal, internal), + FUNCTION 7 gist_tle_same(internal, internal, internal), + FUNCTION 8 gist_tle_distance(internal, tle, smallint, oid, internal); + + +-- ============================================================ +-- Heliocentric type: ecliptic J2000 position in AU +-- ============================================================ + +CREATE TYPE heliocentric; + +CREATE FUNCTION heliocentric_in(cstring) RETURNS heliocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION heliocentric_out(heliocentric) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION heliocentric_recv(internal) RETURNS heliocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION heliocentric_send(heliocentric) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE heliocentric ( + INPUT = heliocentric_in, + OUTPUT = heliocentric_out, + RECEIVE = heliocentric_recv, + SEND = heliocentric_send, + INTERNALLENGTH = 24, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE heliocentric IS 'Heliocentric position in ecliptic J2000 frame (x, y, z in AU)'; + +CREATE FUNCTION helio_x(heliocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION helio_x(heliocentric) IS 'X component in AU (ecliptic J2000, vernal equinox direction)'; + +CREATE FUNCTION helio_y(heliocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION helio_y(heliocentric) IS 'Y component in AU (ecliptic J2000)'; + +CREATE FUNCTION helio_z(heliocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION helio_z(heliocentric) IS 'Z component in AU (ecliptic J2000, north ecliptic pole)'; + +CREATE FUNCTION helio_distance(heliocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION helio_distance(heliocentric) IS 'Heliocentric distance in AU'; + + +-- ============================================================ +-- Star observation functions +-- ============================================================ + +CREATE FUNCTION star_observe(float8, float8, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION star_observe(float8, float8, observer, timestamptz) IS + 'Observe a star from (ra_hours J2000, dec_degrees J2000, observer, time). Returns topocentric az/el. Range is 0 (infinite distance).'; + +CREATE FUNCTION star_observe_safe(float8, float8, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE PARALLEL SAFE; +COMMENT ON FUNCTION star_observe_safe(float8, float8, observer, timestamptz) IS + 'Like star_observe but returns NULL on invalid inputs. For batch queries over star catalogs.'; + + +-- ============================================================ +-- Keplerian propagation functions +-- ============================================================ + +CREATE FUNCTION kepler_propagate( + q_au float8, eccentricity float8, + inclination_deg float8, arg_perihelion_deg float8, + long_asc_node_deg float8, perihelion_jd float8, + t timestamptz +) RETURNS heliocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION kepler_propagate(float8, float8, float8, float8, float8, float8, timestamptz) IS + 'Two-body Keplerian propagation from classical orbital elements. Returns heliocentric ecliptic J2000 position in AU. Handles elliptic, parabolic, and hyperbolic orbits.'; + + +-- ============================================================ +-- Comet observation +-- ============================================================ + +CREATE FUNCTION comet_observe( + q_au float8, eccentricity float8, + inclination_deg float8, arg_perihelion_deg float8, + long_asc_node_deg float8, perihelion_jd float8, + earth_x_au float8, earth_y_au float8, earth_z_au float8, + obs observer, t timestamptz +) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION comet_observe(float8, float8, float8, float8, float8, float8, float8, float8, float8, observer, timestamptz) IS + 'Observe a comet/asteroid from orbital elements. Requires Earth heliocentric ecliptic J2000 position (AU). Returns topocentric az/el with geocentric range in km.'; + + +-- ============================================================ +-- VSOP87 planets, ELP82B Moon, Sun observation +-- ============================================================ + +CREATE FUNCTION planet_heliocentric(int4, timestamptz) RETURNS heliocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_heliocentric(int4, timestamptz) IS + 'VSOP87 heliocentric ecliptic J2000 position (AU). Body IDs: 0=Sun, 1=Mercury, ..., 8=Neptune.'; + +CREATE FUNCTION planet_observe(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_observe(int4, observer, timestamptz) IS + 'Observe a planet from (body_id 1-8, observer, time). Returns topocentric az/el with geocentric range in km.'; + +CREATE FUNCTION sun_observe(observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_observe(observer, timestamptz) IS + 'Observe the Sun from (observer, time). Returns topocentric az/el with geocentric range in km.'; + +CREATE FUNCTION moon_observe(observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_observe(observer, timestamptz) IS + 'Observe the Moon via ELP2000-82B from (observer, time). Returns topocentric az/el with geocentric range in km.'; + + +-- ============================================================ +-- Planetary moon observation +-- ============================================================ + +CREATE FUNCTION galilean_observe(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION galilean_observe(int4, observer, timestamptz) IS + 'Observe a Galilean moon of Jupiter via L1.2 theory. Body IDs: 0=Io, 1=Europa, 2=Ganymede, 3=Callisto.'; + +CREATE FUNCTION saturn_moon_observe(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION saturn_moon_observe(int4, observer, timestamptz) IS + 'Observe a Saturn moon via TASS 1.7. Body IDs: 0=Mimas, 1=Enceladus, 2=Tethys, 3=Dione, 4=Rhea, 5=Titan, 6=Iapetus, 7=Hyperion.'; + +CREATE FUNCTION uranus_moon_observe(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION uranus_moon_observe(int4, observer, timestamptz) IS + 'Observe a Uranus moon via GUST86. Body IDs: 0=Miranda, 1=Ariel, 2=Umbriel, 3=Titania, 4=Oberon.'; + +CREATE FUNCTION mars_moon_observe(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION mars_moon_observe(int4, observer, timestamptz) IS + 'Observe a Mars moon via MarsSat. Body IDs: 0=Phobos, 1=Deimos.'; + + +-- ============================================================ +-- Jupiter decametric radio burst prediction +-- ============================================================ + +CREATE FUNCTION io_phase_angle(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION io_phase_angle(timestamptz) IS + 'Io orbital phase angle in degrees [0,360). 0=superior conjunction (behind Jupiter). Standard Radio JOVE convention.'; + +CREATE FUNCTION jupiter_cml(observer, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION jupiter_cml(observer, timestamptz) IS + 'Jupiter Central Meridian Longitude, System III (1965.0), in degrees [0,360). Light-time corrected.'; + +CREATE FUNCTION jupiter_burst_probability(float8, float8) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION jupiter_burst_probability(float8, float8) IS + 'Estimated Jupiter decametric burst probability (0-1) from (io_phase_deg, cml_deg). Based on Carr et al. (1983) source regions A, B, C, D.'; + + +-- ============================================================ +-- Interplanetary transfer orbits (Lambert solver) +-- ============================================================ + +CREATE FUNCTION lambert_transfer( + dep_body_id int4, arr_body_id int4, + dep_time timestamptz, arr_time timestamptz, + OUT c3_departure float8, OUT c3_arrival float8, + OUT v_inf_departure float8, OUT v_inf_arrival float8, + OUT tof_days float8, OUT transfer_sma float8 +) RETURNS RECORD + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lambert_transfer(int4, int4, timestamptz, timestamptz) IS + 'Solve Lambert transfer between two planets. Returns C3 (km^2/s^2), v_infinity (km/s), TOF (days), SMA (AU). Body IDs 1-8.'; + +CREATE FUNCTION lambert_c3(int4, int4, timestamptz, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lambert_c3(int4, int4, timestamptz, timestamptz) IS + 'Departure C3 (km^2/s^2) for a Lambert transfer. Returns NULL if solver fails. For pork chop plots.'; + + +-- ============================================================ +-- DE ephemeris functions (optional high-precision) +-- ============================================================ + +CREATE FUNCTION planet_heliocentric_de(int4, timestamptz) RETURNS heliocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_heliocentric_de(int4, timestamptz) IS + 'Heliocentric ecliptic J2000 position via JPL DE (sub-arcsecond). Falls back to VSOP87 if DE unavailable. Body IDs: 0=Sun, 1-8=Mercury-Neptune.'; + +CREATE FUNCTION planet_observe_de(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_observe_de(int4, observer, timestamptz) IS + 'Observe planet via JPL DE. Falls back to VSOP87. Body IDs: 1-8 (Mercury-Neptune).'; + +CREATE FUNCTION sun_observe_de(observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_observe_de(observer, timestamptz) IS + 'Observe Sun via JPL DE. Falls back to VSOP87.'; + +CREATE FUNCTION moon_observe_de(observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_observe_de(observer, timestamptz) IS + 'Observe Moon via JPL DE. Falls back to ELP2000-82B.'; + +CREATE FUNCTION lambert_transfer_de( + dep_body_id int4, arr_body_id int4, + dep_time timestamptz, arr_time timestamptz, + OUT c3_departure float8, OUT c3_arrival float8, + OUT v_inf_departure float8, OUT v_inf_arrival float8, + OUT tof_days float8, OUT transfer_sma float8 +) RETURNS RECORD + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lambert_transfer_de(int4, int4, timestamptz, timestamptz) IS + 'Lambert transfer via JPL DE positions. Falls back to VSOP87. Body IDs 1-8.'; + +CREATE FUNCTION lambert_c3_de(int4, int4, timestamptz, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lambert_c3_de(int4, int4, timestamptz, timestamptz) IS + 'Departure C3 via JPL DE. Falls back to VSOP87. For pork chop plots.'; + +CREATE FUNCTION galilean_observe_de(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION galilean_observe_de(int4, observer, timestamptz) IS + 'Observe Galilean moon with JPL DE parent position. L1.2 moon offsets. Body IDs: 0-3 (Io-Callisto).'; + +CREATE FUNCTION saturn_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION saturn_moon_observe_de(int4, observer, timestamptz) IS + 'Observe Saturn moon with JPL DE parent position. TASS 1.7 moon offsets. Body IDs: 0-7 (Mimas-Hyperion).'; + +CREATE FUNCTION uranus_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION uranus_moon_observe_de(int4, observer, timestamptz) IS + 'Observe Uranus moon with JPL DE parent position. GUST86 moon offsets. Body IDs: 0-4 (Miranda-Oberon).'; + +CREATE FUNCTION mars_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION mars_moon_observe_de(int4, observer, timestamptz) IS + 'Observe Mars moon with JPL DE parent position. MarsSat moon offsets. Body IDs: 0-1 (Phobos-Deimos).'; + + +-- Diagnostic function + +CREATE FUNCTION pg_orrery_ephemeris_info( + OUT provider text, OUT file_path text, + OUT start_jd float8, OUT end_jd float8, + OUT version int4, OUT au_km float8 +) RETURNS RECORD + AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE; +COMMENT ON FUNCTION pg_orrery_ephemeris_info() IS + 'Returns current ephemeris provider status: VSOP87 or JPL_DE with file path, JD range, version, and AU value.'; + + +-- ============================================================ +-- Orbit determination (TLE fitting from observations) +-- ============================================================ + +-- Fit TLE from ECI position/velocity ephemeris +-- weights: per-observation weighting (NULL = uniform) + +CREATE FUNCTION tle_from_eci( + positions eci_position[], times timestamptz[], + seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false, + max_iter int4 DEFAULT 15, + weights float8[] DEFAULT NULL, + OUT fitted_tle tle, OUT iterations int4, + OUT rms_final float8, OUT rms_initial float8, OUT status text, + OUT condition_number float8, OUT covariance float8[], OUT nstate int4 +) RETURNS RECORD + AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE; +COMMENT ON FUNCTION tle_from_eci(eci_position[], timestamptz[], tle, boolean, int4, float8[]) IS + 'Fit a TLE from ECI position/velocity observations via differential correction. Optional per-observation weights for heterogeneous sensor fusion. Returns fitted TLE, RMS residuals, convergence status, condition number, and formal covariance matrix.'; + +-- Fit TLE from topocentric observations (az/el/range) — single observer +-- fit_range_rate: include range_rate as 4th residual component +-- weights: per-observation weighting (NULL = uniform) + +CREATE FUNCTION tle_from_topocentric( + observations topocentric[], times timestamptz[], + obs observer, + seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false, + max_iter int4 DEFAULT 15, + fit_range_rate boolean DEFAULT false, + weights float8[] DEFAULT NULL, + OUT fitted_tle tle, OUT iterations int4, + OUT rms_final float8, OUT rms_initial float8, OUT status text, + OUT condition_number float8, OUT covariance float8[], OUT nstate int4 +) RETURNS RECORD + AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE; +COMMENT ON FUNCTION tle_from_topocentric(topocentric[], timestamptz[], observer, tle, boolean, int4, boolean, float8[]) IS + 'Fit a TLE from topocentric (az/el/range) observations via differential correction. Optional range_rate fitting and per-observation weights. Returns fitted TLE, RMS residuals, convergence status, condition number, and formal covariance matrix.'; + +-- Fit TLE from topocentric observations — multiple observers + +CREATE FUNCTION tle_from_topocentric( + observations topocentric[], times timestamptz[], + observers observer[], observer_ids int4[], + seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false, + max_iter int4 DEFAULT 15, + fit_range_rate boolean DEFAULT false, + weights float8[] DEFAULT NULL, + OUT fitted_tle tle, OUT iterations int4, + OUT rms_final float8, OUT rms_initial float8, OUT status text, + OUT condition_number float8, OUT covariance float8[], OUT nstate int4 +) RETURNS RECORD + AS 'MODULE_PATHNAME', 'tle_from_topocentric_multi' + LANGUAGE C STABLE PARALLEL SAFE; +COMMENT ON FUNCTION tle_from_topocentric(topocentric[], timestamptz[], observer[], int4[], tle, boolean, int4, boolean, float8[]) IS + 'Fit a TLE from topocentric observations collected by multiple ground stations. observer_ids[i] indexes into observers[]. Optional range_rate fitting and per-observation weights.'; + +-- Per-observation residuals diagnostic + +CREATE FUNCTION tle_fit_residuals( + fitted tle, + positions eci_position[], + times timestamptz[] +) RETURNS TABLE ( + t timestamptz, + dx_km float8, + dy_km float8, + dz_km float8, + pos_err_km float8 +) + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_fit_residuals(tle, eci_position[], timestamptz[]) IS + 'Compute per-observation position residuals (km) between a TLE and ECI observations. Useful for fit quality assessment.'; + +-- Fit TLE from RA/Dec observations — single observer +-- Uses Gauss method for initial orbit determination when no seed is provided. +-- RA in hours [0,24), Dec in degrees [-90,90] (matches star_observe convention). +-- RMS output is in radians for angles-only mode. + +CREATE FUNCTION tle_from_angles( + ra_hours float8[], dec_degrees float8[], + times timestamptz[], + obs observer, + seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false, + max_iter int4 DEFAULT 15, + weights float8[] DEFAULT NULL, + OUT fitted_tle tle, OUT iterations int4, + OUT rms_final float8, OUT rms_initial float8, OUT status text, + OUT condition_number float8, OUT covariance float8[], OUT nstate int4 +) RETURNS RECORD AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE; +COMMENT ON FUNCTION tle_from_angles(float8[], float8[], timestamptz[], observer, tle, boolean, int4, float8[]) IS + 'Fit a TLE from angles-only (RA/Dec) observations via Gauss IOD + differential correction. RA in hours [0,24), Dec in degrees [-90,90]. RMS output in radians. Uses Gauss method for seed-free initial guess.'; + +-- Fit TLE from RA/Dec observations — multiple observers + +CREATE FUNCTION tle_from_angles( + ra_hours float8[], dec_degrees float8[], + times timestamptz[], + observers observer[], observer_ids int4[], + seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false, + max_iter int4 DEFAULT 15, + weights float8[] DEFAULT NULL, + OUT fitted_tle tle, OUT iterations int4, + OUT rms_final float8, OUT rms_initial float8, OUT status text, + OUT condition_number float8, OUT covariance float8[], OUT nstate int4 +) RETURNS RECORD + AS 'MODULE_PATHNAME', 'tle_from_angles_multi' + LANGUAGE C STABLE PARALLEL SAFE; +COMMENT ON FUNCTION tle_from_angles(float8[], float8[], timestamptz[], observer[], int4[], tle, boolean, int4, float8[]) IS + 'Fit a TLE from angles-only (RA/Dec) observations from multiple ground stations via Gauss IOD + differential correction.'; +-- pg_orrery 0.6.0 -> 0.7.0 migration +-- +-- Adds SP-GiST orbital trie index for satellite pass prediction. +-- 2-level trie: SMA (L0) + inclination (L1) with query-time RAAN filter. +-- The &? operator answers "might this satellite be visible?" + +-- ============================================================ +-- observer_window composite type (query parameter bundle) +-- ============================================================ + +CREATE TYPE observer_window AS ( + obs observer, + t_start timestamptz, + t_end timestamptz, + min_el float8 +); + +COMMENT ON TYPE observer_window IS + 'Observation query parameters: observer location, time window, and minimum elevation angle (degrees). Used with the &? visibility cone operator.'; + +-- ============================================================ +-- Visibility cone operator function +-- ============================================================ + +CREATE FUNCTION tle_visibility_possible(tle, observer_window) RETURNS boolean + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; + +COMMENT ON FUNCTION tle_visibility_possible(tle, observer_window) IS + 'Could this satellite be visible from the observer during the time window? Combines altitude, inclination, and RAAN checks. Conservative superset — survivors need SGP4 propagation for ground truth.'; + +-- ============================================================ +-- &? operator (visibility cone check) +-- ============================================================ +-- The indexed column (tle) MUST be the left argument so PostgreSQL +-- can form a ScanKey and pass it to inner_consistent for pruning. + +CREATE OPERATOR &? ( + LEFTARG = tle, + RIGHTARG = observer_window, + FUNCTION = tle_visibility_possible, + RESTRICT = contsel, + JOIN = contjoinsel +); + +COMMENT ON OPERATOR &? (tle, observer_window) IS + 'Visibility cone check: could this satellite be visible from the observer during the time window? Index-accelerated via SP-GiST orbital trie.'; + +-- ============================================================ +-- SP-GiST support functions +-- ============================================================ + +CREATE FUNCTION spgist_tle_config(internal, internal) RETURNS void + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION spgist_tle_choose(internal, internal) RETURNS void + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION spgist_tle_picksplit(internal, internal) RETURNS void + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION spgist_tle_inner_consistent(internal, internal) RETURNS void + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION spgist_tle_leaf_consistent(internal, internal) RETURNS void + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +-- ============================================================ +-- SP-GiST operator class (opt-in, not DEFAULT) +-- ============================================================ + +CREATE OPERATOR CLASS tle_spgist_ops + FOR TYPE tle USING spgist AS + OPERATOR 1 &? (tle, observer_window), + FUNCTION 1 spgist_tle_config(internal, internal), + FUNCTION 2 spgist_tle_choose(internal, internal), + FUNCTION 3 spgist_tle_picksplit(internal, internal), + FUNCTION 4 spgist_tle_inner_consistent(internal, internal), + FUNCTION 5 spgist_tle_leaf_consistent(internal, internal); +-- pg_orrery 0.7.0 -> 0.8.0 migration +-- +-- Adds orbital_elements type for comets/asteroids, MPC MPCORB.DAT parser, +-- and small_body_observe()/small_body_heliocentric() observation functions. + +-- ============================================================ +-- orbital_elements type +-- ============================================================ + +CREATE TYPE orbital_elements; + +CREATE FUNCTION orbital_elements_in(cstring) RETURNS orbital_elements + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION orbital_elements_out(orbital_elements) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION orbital_elements_recv(internal) RETURNS orbital_elements + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION orbital_elements_send(orbital_elements) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE orbital_elements ( + INPUT = orbital_elements_in, + OUTPUT = orbital_elements_out, + RECEIVE = orbital_elements_recv, + SEND = orbital_elements_send, + INTERNALLENGTH = 72, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE orbital_elements IS + 'Classical Keplerian orbital elements for comets and asteroids (epoch, q, e, inc, omega, Omega, tp, H, G). 72 bytes, fixed-size.'; + + +-- ============================================================ +-- Accessor functions +-- ============================================================ + +CREATE FUNCTION oe_epoch(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_epoch(orbital_elements) IS 'Osculation epoch (Julian date)'; + +CREATE FUNCTION oe_perihelion(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_perihelion(orbital_elements) IS 'Perihelion distance q (AU)'; + +CREATE FUNCTION oe_eccentricity(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_eccentricity(orbital_elements) IS 'Eccentricity'; + +CREATE FUNCTION oe_inclination(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_inclination(orbital_elements) IS 'Inclination (degrees)'; + +CREATE FUNCTION oe_arg_perihelion(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_arg_perihelion(orbital_elements) IS 'Argument of perihelion (degrees)'; + +CREATE FUNCTION oe_raan(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_raan(orbital_elements) IS 'Longitude of ascending node (degrees)'; + +CREATE FUNCTION oe_tp(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_tp(orbital_elements) IS 'Time of perihelion passage (Julian date)'; + +CREATE FUNCTION oe_h_mag(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_h_mag(orbital_elements) IS 'Absolute magnitude H (NaN if unknown)'; + +CREATE FUNCTION oe_g_slope(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_g_slope(orbital_elements) IS 'Slope parameter G (NaN if unknown)'; + +CREATE FUNCTION oe_semi_major_axis(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_semi_major_axis(orbital_elements) IS 'Semi-major axis a = q/(1-e) in AU. NULL for parabolic/hyperbolic orbits (e >= 1).'; + +CREATE FUNCTION oe_period_years(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_period_years(orbital_elements) IS 'Orbital period in years = a^1.5 (Kepler third law). NULL for parabolic/hyperbolic orbits (e >= 1).'; + + +-- ============================================================ +-- MPC MPCORB.DAT parser +-- ============================================================ + +CREATE FUNCTION oe_from_mpc(text) RETURNS orbital_elements + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_from_mpc(text) IS + 'Parse one MPCORB.DAT fixed-width line into orbital_elements. Converts MPC packed epoch, computes perihelion distance and tp from (a, e, M).'; + + +-- ============================================================ +-- Observation functions +-- ============================================================ + +CREATE FUNCTION small_body_heliocentric(orbital_elements, timestamptz) RETURNS heliocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION small_body_heliocentric(orbital_elements, timestamptz) IS + 'Heliocentric ecliptic J2000 position of a comet/asteroid from its orbital elements at a given time.'; + +CREATE FUNCTION small_body_observe(orbital_elements, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION small_body_observe(orbital_elements, observer, timestamptz) IS + 'Observe a comet/asteroid from orbital elements. Auto-fetches Earth via VSOP87. Returns topocentric az/el with geocentric range in km.'; +-- pg_orrery 0.8.0 -> 0.9.0 migration +-- +-- Adds equatorial type (apparent RA/Dec of date), atmospheric refraction, +-- stellar proper motion, and light-time corrected _apparent() functions. + +-- ============================================================ +-- equatorial type — apparent RA/Dec of date +-- ============================================================ + +CREATE TYPE equatorial; + +CREATE FUNCTION equatorial_in(cstring) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION equatorial_out(equatorial) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION equatorial_recv(internal) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION equatorial_send(equatorial) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE equatorial ( + INPUT = equatorial_in, + OUTPUT = equatorial_out, + RECEIVE = equatorial_recv, + SEND = equatorial_send, + INTERNALLENGTH = 24, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE equatorial IS + 'Apparent equatorial coordinates of date: RA (hours), Dec (degrees), distance (km). Solar system: J2000 precessed via IAU 1976. Satellites: TEME frame (~of-date to ~arcsecond). 24 bytes, fixed-size.'; + + +-- ============================================================ +-- Equatorial accessor functions +-- ============================================================ + +CREATE FUNCTION eq_ra(equatorial) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eq_ra(equatorial) IS 'Right ascension in hours [0, 24)'; + +CREATE FUNCTION eq_dec(equatorial) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eq_dec(equatorial) IS 'Declination in degrees [-90, 90]'; + +CREATE FUNCTION eq_distance(equatorial) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eq_distance(equatorial) IS 'Distance in km (0 for stars without parallax)'; + + +-- ============================================================ +-- Satellite RA/Dec functions +-- ============================================================ + +CREATE FUNCTION eci_to_equatorial(eci_position, observer, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eci_to_equatorial(eci_position, observer, timestamptz) IS + 'Topocentric apparent RA/Dec from ECI position. Observer parallax-corrected — LEO parallax is ~1 degree.'; + +CREATE FUNCTION eci_to_equatorial_geo(eci_position, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eci_to_equatorial_geo(eci_position, timestamptz) IS + 'Geocentric apparent RA/Dec from ECI position. Observer-independent — the direction of the TEME position vector.'; + + +-- ============================================================ +-- Solar system equatorial functions (VSOP87) +-- ============================================================ + +CREATE FUNCTION planet_equatorial(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_equatorial(int4, timestamptz) IS + 'Geocentric apparent RA/Dec of a planet via VSOP87. Body IDs: 1=Mercury through 8=Neptune.'; + +CREATE FUNCTION sun_equatorial(timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_equatorial(timestamptz) IS + 'Geocentric apparent RA/Dec of the Sun via VSOP87.'; + +CREATE FUNCTION moon_equatorial(timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_equatorial(timestamptz) IS + 'Geocentric apparent RA/Dec of the Moon via ELP2000-82B.'; + +CREATE FUNCTION small_body_equatorial(orbital_elements, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION small_body_equatorial(orbital_elements, timestamptz) IS + 'Geocentric apparent RA/Dec of a comet/asteroid from orbital elements. Earth via VSOP87.'; + +CREATE FUNCTION star_equatorial(float8, float8, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION star_equatorial(float8, float8, timestamptz) IS + 'Apparent RA/Dec of a star at a given time. Precesses J2000 catalog coordinates (RA hours, Dec degrees) to date via IAU 1976.'; + + +-- ============================================================ +-- Atmospheric refraction (Bennett 1982) +-- ============================================================ + +CREATE FUNCTION atmospheric_refraction(float8) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION atmospheric_refraction(float8) IS + 'Atmospheric refraction correction in degrees for a given geometric elevation (degrees). Standard atmosphere: P=1010 mbar, T=10C. Bennett (1982) formula with domain guard at -1 degree.'; + +CREATE FUNCTION atmospheric_refraction_ext(float8, float8, float8) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION atmospheric_refraction_ext(float8, float8, float8) IS + 'Atmospheric refraction with pressure/temperature correction. Args: elevation_deg, pressure_mbar, temperature_celsius. Meeus P/T factor applied to Bennett formula.'; + +CREATE FUNCTION topo_elevation_apparent(topocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION topo_elevation_apparent(topocentric) IS + 'Apparent elevation in degrees — geometric elevation plus atmospheric refraction correction.'; + + +-- ============================================================ +-- Refracted pass prediction +-- ============================================================ + +CREATE FUNCTION predict_passes_refracted( + tle, observer, timestamptz, timestamptz, float8 DEFAULT 0.0 +) RETURNS SETOF pass_event + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE + ROWS 20; +COMMENT ON FUNCTION predict_passes_refracted(tle, observer, timestamptz, timestamptz, float8) IS + 'Predict satellite passes using a refracted horizon threshold (-0.569 deg geometric). Atmospheric refraction makes satellites visible ~35 seconds earlier at AOS and later at LOS.'; + + +-- ============================================================ +-- Stellar proper motion +-- ============================================================ + +CREATE FUNCTION star_observe_pm( + float8, float8, float8, float8, float8, float8, observer, timestamptz +) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION star_observe_pm(float8, float8, float8, float8, float8, float8, observer, timestamptz) IS + 'Observe a star with proper motion. Args: ra_hours, dec_deg, pm_ra_masyr (mu_alpha*cos(delta)), pm_dec_masyr, parallax_mas, rv_kms, observer, time. Hipparcos/Gaia convention for pm_ra.'; + +CREATE FUNCTION star_equatorial_pm( + float8, float8, float8, float8, float8, float8, timestamptz +) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION star_equatorial_pm(float8, float8, float8, float8, float8, float8, timestamptz) IS + 'Apparent RA/Dec of a star with proper motion. Args: ra_hours, dec_deg, pm_ra_masyr, pm_dec_masyr, parallax_mas, rv_kms, time. Distance from parallax if > 0.'; + + +-- ============================================================ +-- Light-time corrected observation functions +-- ============================================================ + +CREATE FUNCTION planet_observe_apparent(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_observe_apparent(int4, observer, timestamptz) IS + 'Observe a planet with single-iteration light-time correction. Body at retarded time, Earth at observation time. VSOP87.'; + +CREATE FUNCTION sun_observe_apparent(observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_observe_apparent(observer, timestamptz) IS + 'Observe the Sun with light-time correction (~8.3 min). VSOP87.'; + +CREATE FUNCTION small_body_observe_apparent(orbital_elements, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION small_body_observe_apparent(orbital_elements, observer, timestamptz) IS + 'Observe a comet/asteroid with single-iteration light-time correction. Kepler propagation at retarded time, Earth via VSOP87 at observation time.'; + + +-- ============================================================ +-- Light-time corrected equatorial functions +-- ============================================================ + +CREATE FUNCTION planet_equatorial_apparent(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_equatorial_apparent(int4, timestamptz) IS + 'Geocentric apparent RA/Dec of a planet with light-time correction. VSOP87.'; + +CREATE FUNCTION moon_equatorial_apparent(timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_equatorial_apparent(timestamptz) IS + 'Geocentric apparent RA/Dec of the Moon with light-time correction (~1.3 sec). ELP2000-82B.'; + +CREATE FUNCTION small_body_equatorial_apparent(orbital_elements, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION small_body_equatorial_apparent(orbital_elements, timestamptz) IS + 'Geocentric apparent RA/Dec of a comet/asteroid with light-time correction.'; + + +-- ============================================================ +-- DE ephemeris equatorial variants (STABLE) +-- ============================================================ + +CREATE FUNCTION planet_equatorial_de(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_equatorial_de(int4, timestamptz) IS + 'Geocentric apparent RA/Dec of a planet via JPL DE ephemeris (falls back to VSOP87 + equatorial).'; + +CREATE FUNCTION moon_equatorial_de(timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_equatorial_de(timestamptz) IS + 'Geocentric apparent RA/Dec of the Moon via JPL DE ephemeris (falls back to ELP2000-82B + equatorial).'; +-- pg_orrery 0.9.0 -> 0.10.0 migration +-- +-- Adds annual aberration to existing _apparent() functions, +-- 6 new _apparent_de() variants, equatorial angular separation +-- operator and cone predicate, and stellar annual parallax. + +-- ============================================================ +-- Equatorial angular distance and cone search +-- ============================================================ + +CREATE FUNCTION eq_angular_distance(equatorial, equatorial) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eq_angular_distance(equatorial, equatorial) IS + 'Angular separation in degrees between two equatorial positions. Vincenty formula (stable at 0 and 180 degrees).'; + +CREATE FUNCTION eq_within_cone(equatorial, equatorial, float8) RETURNS bool + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eq_within_cone(equatorial, equatorial, float8) IS + 'True if first position is within radius_deg of second position. Cosine shortcut for fast rejection.'; + +CREATE OPERATOR <-> ( + LEFTARG = equatorial, + RIGHTARG = equatorial, + FUNCTION = eq_angular_distance, + COMMUTATOR = <-> +); +COMMENT ON OPERATOR <-> (equatorial, equatorial) IS + 'Angular separation in degrees between two equatorial positions.'; + + +-- ============================================================ +-- DE apparent observation functions (STABLE, light-time + aberration) +-- ============================================================ + +CREATE FUNCTION planet_observe_apparent_de(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_observe_apparent_de(int4, observer, timestamptz) IS + 'Observe a planet with light-time correction and annual aberration via JPL DE (falls back to VSOP87).'; + +CREATE FUNCTION sun_observe_apparent_de(observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_observe_apparent_de(observer, timestamptz) IS + 'Observe the Sun with aberration via JPL DE (falls back to VSOP87).'; + +CREATE FUNCTION moon_observe_apparent_de(observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_observe_apparent_de(observer, timestamptz) IS + 'Observe the Moon with light-time correction and annual aberration via JPL DE (falls back to ELP2000-82B).'; + +CREATE FUNCTION planet_equatorial_apparent_de(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_equatorial_apparent_de(int4, timestamptz) IS + 'Geocentric apparent RA/Dec of a planet with light-time correction and annual aberration via JPL DE (falls back to VSOP87).'; + +CREATE FUNCTION moon_equatorial_apparent_de(timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_equatorial_apparent_de(timestamptz) IS + 'Geocentric apparent RA/Dec of the Moon with light-time correction and annual aberration via JPL DE (falls back to ELP2000-82B).'; + +CREATE FUNCTION small_body_observe_apparent_de(orbital_elements, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION small_body_observe_apparent_de(orbital_elements, observer, timestamptz) IS + 'Observe a comet/asteroid with light-time correction and annual aberration. Earth position via JPL DE (falls back to VSOP87).'; +-- pg_orrery 0.10.0 -> 0.11.0 migration +-- +-- Adds make_orbital_elements() constructors and +-- geocentric equatorial functions for planetary moons. + +-- ============================================================ +-- orbital_elements constructors +-- ============================================================ + +CREATE FUNCTION make_orbital_elements( + epoch_jd float8, q_au float8, e float8, + inc_rad float8, omega_rad float8, node_rad float8, + tp_jd float8, h_mag float8, g_slope float8 +) RETURNS orbital_elements + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION make_orbital_elements(float8,float8,float8,float8,float8,float8,float8,float8,float8) IS + 'Construct orbital_elements from 9 floats (angular elements in radians).'; + +CREATE FUNCTION make_orbital_elements_deg( + epoch_jd float8, q_au float8, e float8, + inc_deg float8, omega_deg float8, node_deg float8, + tp_jd float8, h_mag float8, g_slope float8 +) RETURNS orbital_elements + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION make_orbital_elements_deg(float8,float8,float8,float8,float8,float8,float8,float8,float8) IS + 'Construct orbital_elements from 9 floats (angular elements in degrees). Matches text I/O and most catalog column layouts.'; + + +-- ============================================================ +-- Planetary moon equatorial functions +-- ============================================================ + +CREATE FUNCTION galilean_equatorial(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION galilean_equatorial(int4, timestamptz) IS + 'Geometric geocentric RA/Dec of a Galilean moon (0=Io, 1=Europa, 2=Ganymede, 3=Callisto). L1.2 theory + VSOP87. No light-time or aberration correction.'; + +CREATE FUNCTION saturn_moon_equatorial(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION saturn_moon_equatorial(int4, timestamptz) IS + 'Geometric geocentric RA/Dec of a Saturn moon (0=Mimas..7=Hyperion). TASS17 theory + VSOP87. No light-time or aberration correction.'; + +CREATE FUNCTION uranus_moon_equatorial(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION uranus_moon_equatorial(int4, timestamptz) IS + 'Geometric geocentric RA/Dec of a Uranus moon (0=Miranda..4=Oberon). GUST86 theory + VSOP87. No light-time or aberration correction.'; + +CREATE FUNCTION mars_moon_equatorial(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION mars_moon_equatorial(int4, timestamptz) IS + 'Geometric geocentric RA/Dec of a Mars moon (0=Phobos, 1=Deimos). MarsSat theory + VSOP87. No light-time or aberration correction.'; +-- pg_orrery 0.11.0 -> 0.12.0 migration +-- +-- Adds equatorial GiST operator class for KNN sky queries +-- and DE moon equatorial functions for all 4 planetary moon families. + +-- ============================================================ +-- GiST support functions for equatorial type +-- ============================================================ + +CREATE FUNCTION gist_eq_consistent(internal, equatorial, smallint, oid, internal) RETURNS bool + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_eq_union(internal, internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_eq_compress(internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_eq_decompress(internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_eq_penalty(internal, internal, internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_eq_picksplit(internal, internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_eq_same(internal, internal, internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_eq_distance(internal, equatorial, smallint, oid, internal) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +-- ============================================================ +-- Equatorial GiST operator class (KNN ordering only) +-- ============================================================ + +CREATE OPERATOR CLASS eq_gist_ops + DEFAULT FOR TYPE equatorial USING gist AS + OPERATOR 15 <-> (equatorial, equatorial) FOR ORDER BY pg_catalog.float_ops, + FUNCTION 1 gist_eq_consistent(internal, equatorial, smallint, oid, internal), + FUNCTION 2 gist_eq_union(internal, internal), + FUNCTION 3 gist_eq_compress(internal), + FUNCTION 4 gist_eq_decompress(internal), + FUNCTION 5 gist_eq_penalty(internal, internal, internal), + FUNCTION 6 gist_eq_picksplit(internal, internal), + FUNCTION 7 gist_eq_same(internal, internal, internal), + FUNCTION 8 gist_eq_distance(internal, equatorial, smallint, oid, internal); + +-- ============================================================ +-- DE moon equatorial functions (STABLE, fall back to VSOP87) +-- ============================================================ + +CREATE FUNCTION galilean_equatorial_de(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION galilean_equatorial_de(int4, timestamptz) IS + 'Geocentric RA/Dec of a Galilean moon via DE parent position (falls back to VSOP87). 0=Io, 1=Europa, 2=Ganymede, 3=Callisto.'; + +CREATE FUNCTION saturn_moon_equatorial_de(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION saturn_moon_equatorial_de(int4, timestamptz) IS + 'Geocentric RA/Dec of a Saturn moon via DE parent position (falls back to VSOP87). 0=Mimas..7=Hyperion.'; + +CREATE FUNCTION uranus_moon_equatorial_de(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION uranus_moon_equatorial_de(int4, timestamptz) IS + 'Geocentric RA/Dec of a Uranus moon via DE parent position (falls back to VSOP87). 0=Miranda..4=Oberon.'; + +CREATE FUNCTION mars_moon_equatorial_de(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION mars_moon_equatorial_de(int4, timestamptz) IS + 'Geocentric RA/Dec of a Mars moon via DE parent position (falls back to VSOP87). 0=Phobos, 1=Deimos.'; + + +-- ============================================================ +-- v0.13.0: make_equatorial() constructor +-- ============================================================ + +CREATE FUNCTION make_equatorial(ra_hours float8, dec_deg float8, distance_km float8) + RETURNS equatorial + AS 'MODULE_PATHNAME', 'make_equatorial' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +COMMENT ON FUNCTION make_equatorial(float8, float8, float8) IS + 'Construct equatorial from RA (hours [0,24)), Dec (degrees [-90,90]), distance (km).'; + + +-- ============================================================ +-- v0.13.0: Rise/set prediction functions +-- ============================================================ + +CREATE FUNCTION planet_next_rise(body_id int4, obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_next_rise(int4, observer, timestamptz) IS + 'Next geometric rise time for a planet. Returns NULL if no rise within 7 days. body_id: 1=Mercury..8=Neptune.'; + +CREATE FUNCTION planet_next_set(body_id int4, obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_next_set(int4, observer, timestamptz) IS + 'Next geometric set time for a planet. Returns NULL if no set within 7 days. body_id: 1=Mercury..8=Neptune.'; + +CREATE FUNCTION sun_next_rise(obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_next_rise(observer, timestamptz) IS + 'Next geometric sunrise. Returns NULL if Sun does not rise within 7 days (polar night).'; + +CREATE FUNCTION sun_next_set(obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_next_set(observer, timestamptz) IS + 'Next geometric sunset. Returns NULL if Sun does not set within 7 days (midnight sun).'; + +CREATE FUNCTION moon_next_rise(obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_next_rise(observer, timestamptz) IS + 'Next geometric moonrise. Returns NULL if Moon does not rise within 7 days.'; + +CREATE FUNCTION moon_next_set(obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_next_set(observer, timestamptz) IS + 'Next geometric moonset. Returns NULL if Moon does not set within 7 days.'; + +CREATE FUNCTION sun_next_rise_refracted(obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_next_rise_refracted(observer, timestamptz) IS + 'Next refracted sunrise (-0.833 deg threshold: refraction + semidiameter). Earlier than geometric by ~4 min.'; + +CREATE FUNCTION sun_next_set_refracted(obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_next_set_refracted(observer, timestamptz) IS + 'Next refracted sunset (-0.833 deg threshold: refraction + semidiameter). Later than geometric by ~4 min.'; + + +-- ============================================================ +-- v0.14.0: Refracted planet/moon rise/set +-- ============================================================ + +CREATE FUNCTION planet_next_rise_refracted(body_id int4, obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_next_rise_refracted(int4, observer, timestamptz) IS + 'Next refracted rise time for a planet (-0.569 deg threshold: atmospheric refraction only). Earlier than geometric.'; + +CREATE FUNCTION planet_next_set_refracted(body_id int4, obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_next_set_refracted(int4, observer, timestamptz) IS + 'Next refracted set time for a planet (-0.569 deg threshold: atmospheric refraction only). Later than geometric.'; + +CREATE FUNCTION moon_next_rise_refracted(obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_next_rise_refracted(observer, timestamptz) IS + 'Next refracted moonrise (-0.833 deg threshold: refraction + semidiameter). Earlier than geometric.'; + +CREATE FUNCTION moon_next_set_refracted(obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_next_set_refracted(observer, timestamptz) IS + 'Next refracted moonset (-0.833 deg threshold: refraction + semidiameter). Later than geometric.'; + + +-- ============================================================ +-- v0.14.0: Constellation identification (Roman 1987, CDS VI/42) +-- ============================================================ + +CREATE FUNCTION constellation(eq equatorial) RETURNS text + AS 'MODULE_PATHNAME', 'constellation_from_equatorial' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION constellation(equatorial) IS + 'IAU constellation abbreviation (3 letters) from equatorial coordinates (Roman 1987).'; + +CREATE FUNCTION constellation(ra_hours float8, dec_deg float8) RETURNS text + AS 'MODULE_PATHNAME', 'constellation_from_radec' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION constellation(float8, float8) IS + 'IAU constellation from J2000 RA (hours [0,24)) and Dec (degrees [-90,90]).'; +-- pg_orrery 0.14.0 -> 0.15.0 migration +-- +-- Adds: constellation_full_name (1 function), +-- rise/set status diagnostics (3 functions). + +-- ============================================================ +-- Constellation full name lookup +-- ============================================================ + +CREATE FUNCTION constellation_full_name(abbr text) RETURNS text + AS 'MODULE_PATHNAME', 'constellation_full_name_from_abbr' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION constellation_full_name(text) IS + 'Full IAU constellation name from 3-letter abbreviation. Returns NULL for invalid abbreviation.'; + +-- ============================================================ +-- Rise/set status diagnostics +-- ============================================================ + +CREATE FUNCTION sun_rise_set_status(obs observer, t timestamptz) RETURNS text + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_rise_set_status(observer, timestamptz) IS + 'Classify Sun visibility: rises_and_sets, circumpolar, or never_rises.'; + +CREATE FUNCTION moon_rise_set_status(obs observer, t timestamptz) RETURNS text + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_rise_set_status(observer, timestamptz) IS + 'Classify Moon visibility: rises_and_sets, circumpolar, or never_rises.'; + +CREATE FUNCTION planet_rise_set_status(body_id int4, obs observer, t timestamptz) RETURNS text + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_rise_set_status(int4, observer, timestamptz) IS + 'Classify planet visibility: rises_and_sets, circumpolar, or never_rises. Body IDs 1-8 (Mercury-Neptune).'; +-- pg_orrery 0.15.0 -> 0.16.0: twilight, lunar phase, planet magnitude + +-- ============================================================ +-- Twilight functions (6) +-- ============================================================ + +CREATE FUNCTION sun_civil_dawn(observer, timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME', 'sun_civil_dawn' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_civil_dawn(observer, timestamptz) IS + 'Next civil dawn (Sun crosses -6 deg rising). Outdoor activities without artificial light.'; + +CREATE FUNCTION sun_civil_dusk(observer, timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME', 'sun_civil_dusk' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_civil_dusk(observer, timestamptz) IS + 'Next civil dusk (Sun crosses -6 deg setting). Artificial light needed.'; + +CREATE FUNCTION sun_nautical_dawn(observer, timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME', 'sun_nautical_dawn' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_nautical_dawn(observer, timestamptz) IS + 'Next nautical dawn (Sun crosses -12 deg rising). Horizon visible at sea.'; + +CREATE FUNCTION sun_nautical_dusk(observer, timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME', 'sun_nautical_dusk' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_nautical_dusk(observer, timestamptz) IS + 'Next nautical dusk (Sun crosses -12 deg setting). Horizon no longer visible at sea.'; + +CREATE FUNCTION sun_astronomical_dawn(observer, timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME', 'sun_astronomical_dawn' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_astronomical_dawn(observer, timestamptz) IS + 'Next astronomical dawn (Sun crosses -18 deg rising). Sky was fully dark.'; + +CREATE FUNCTION sun_astronomical_dusk(observer, timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME', 'sun_astronomical_dusk' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_astronomical_dusk(observer, timestamptz) IS + 'Next astronomical dusk (Sun crosses -18 deg setting). Sky becomes fully dark.'; + +-- ============================================================ +-- Lunar phase functions (4) +-- ============================================================ + +CREATE FUNCTION moon_phase_angle(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'moon_phase_angle' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_phase_angle(timestamptz) IS + 'Sun-Earth-Moon phase angle in degrees [0,360). 0=new, 90=first quarter, 180=full, 270=last quarter.'; + +CREATE FUNCTION moon_illumination(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'moon_illumination' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_illumination(timestamptz) IS + 'Illuminated fraction of the Moon disk [0.0, 1.0].'; + +CREATE FUNCTION moon_phase_name(timestamptz) RETURNS text + AS 'MODULE_PATHNAME', 'moon_phase_name' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_phase_name(timestamptz) IS + 'Moon phase name: new_moon, waxing_crescent, first_quarter, waxing_gibbous, full_moon, waning_gibbous, last_quarter, waning_crescent.'; + +CREATE FUNCTION moon_age(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'moon_age' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_age(timestamptz) IS + 'Days since last new moon [0, ~29.53), approximated from phase angle.'; + +-- ============================================================ +-- Planet magnitude (1) +-- ============================================================ + +CREATE FUNCTION planet_magnitude(int4, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'planet_magnitude' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_magnitude(int4, timestamptz) IS + 'Apparent visual magnitude of a planet (Mallama & Hilton 2018). Body IDs 1-8. Saturn ring tilt not modeled.'; +-- pg_orrery 0.16.0 -> 0.17.0: solar elongation, planet phase, satellite eclipse, +-- observing night quality, lunar libration + +-- ============================================================ +-- Solar elongation (1) +-- ============================================================ + +CREATE FUNCTION solar_elongation(int4, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'solar_elongation' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION solar_elongation(int4, timestamptz) IS + 'Sun-Earth-Planet angle in degrees [0, 180]. How far a planet appears from the Sun. Body IDs 1-8.'; + +-- ============================================================ +-- Planet phase fraction (1) +-- ============================================================ + +CREATE FUNCTION planet_phase(int4, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'planet_phase' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_phase(int4, timestamptz) IS + 'Illuminated fraction of a planet disk as seen from Earth [0.0, 1.0]. Body IDs 1-8.'; + +-- ============================================================ +-- Satellite eclipse prediction (4) +-- ============================================================ + +CREATE FUNCTION satellite_is_eclipsed(tle, timestamptz) RETURNS bool + AS 'MODULE_PATHNAME', 'satellite_is_eclipsed' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION satellite_is_eclipsed(tle, timestamptz) IS + 'True if the satellite is in Earth cylindrical shadow at the given time.'; + +CREATE FUNCTION satellite_next_eclipse_entry(tle, timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME', 'satellite_next_eclipse_entry' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION satellite_next_eclipse_entry(tle, timestamptz) IS + 'Next time the satellite enters Earth shadow (up to 7-day search). NULL if none found.'; + +CREATE FUNCTION satellite_next_eclipse_exit(tle, timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME', 'satellite_next_eclipse_exit' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION satellite_next_eclipse_exit(tle, timestamptz) IS + 'Next time the satellite exits Earth shadow (up to 7-day search). NULL if none found.'; + +CREATE FUNCTION satellite_eclipse_fraction(tle, timestamptz, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'satellite_eclipse_fraction' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION satellite_eclipse_fraction(tle, timestamptz, timestamptz) IS + 'Fraction of the given time window the satellite spends in eclipse [0.0, 1.0].'; + +-- ============================================================ +-- Observing night quality (1) +-- ============================================================ + +CREATE FUNCTION observing_night_quality(observer, timestamptz DEFAULT NOW()) +RETURNS text AS $$ +DECLARE + astro_dusk timestamptz; + astro_dawn timestamptz; + dark_hours float8; + illum float8; + moon_up bool; + score int := 100; +BEGIN + -- Astronomical darkness window + astro_dusk := sun_astronomical_dusk($1, $2); + IF astro_dusk IS NULL THEN + RETURN 'poor'; -- No astronomical darkness (polar summer) + END IF; + astro_dawn := sun_astronomical_dawn($1, astro_dusk); + IF astro_dawn IS NULL THEN + RETURN 'poor'; + END IF; + + dark_hours := extract(epoch FROM astro_dawn - astro_dusk) / 3600.0; + + -- Short dark window penalty + IF dark_hours < 2.0 THEN score := score - 40; + ELSIF dark_hours < 4.0 THEN score := score - 20; + ELSIF dark_hours < 6.0 THEN score := score - 10; + END IF; + + -- Moon illumination penalty + illum := moon_illumination(astro_dusk); + IF illum > 0.75 THEN + -- Check if Moon is above horizon during darkness + moon_up := (moon_observe($1, astro_dusk)).elevation > 0 + OR (moon_observe($1, astro_dusk + (astro_dawn - astro_dusk) / 2)).elevation > 0; + IF moon_up THEN + score := score - (illum * 30)::int; -- Up to -30 for full moon + END IF; + END IF; + + -- Classify + IF score >= 80 THEN RETURN 'excellent'; + ELSIF score >= 60 THEN RETURN 'good'; + ELSIF score >= 40 THEN RETURN 'fair'; + ELSE RETURN 'poor'; + END IF; +END; +$$ LANGUAGE plpgsql STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION observing_night_quality(observer, timestamptz) IS + 'Composite observing quality assessment: excellent/good/fair/poor based on darkness duration and Moon interference.'; + +-- ============================================================ +-- Lunar libration (5) +-- ============================================================ + +CREATE FUNCTION moon_libration_longitude(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'moon_libration_longitude' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_libration_longitude(timestamptz) IS + 'Optical libration in longitude (degrees, typically [-8, +8]). Meeus Ch. 53.'; + +CREATE FUNCTION moon_libration_latitude(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'moon_libration_latitude' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_libration_latitude(timestamptz) IS + 'Optical libration in latitude (degrees, typically [-7, +7]). Meeus Ch. 53.'; + +CREATE FUNCTION moon_libration_position_angle(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'moon_libration_position_angle' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_libration_position_angle(timestamptz) IS + 'Position angle of the Moon axis (degrees, [0, 360)). Meeus Ch. 53.'; + +CREATE FUNCTION moon_libration(timestamptz, + OUT l float8, OUT b float8, OUT p float8) RETURNS record + AS 'MODULE_PATHNAME', 'moon_libration' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_libration(timestamptz) IS + 'All three libration values: longitude (l), latitude (b), position angle (p) in degrees.'; + +CREATE FUNCTION moon_subsolar_longitude(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'moon_subsolar_longitude' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_subsolar_longitude(timestamptz) IS + 'Selenographic longitude of the sub-solar point (degrees, [0, 360)). Determines the lunar terminator position.'; diff --git a/src/eclipse_funcs.c b/src/eclipse_funcs.c new file mode 100644 index 0000000..89f3117 --- /dev/null +++ b/src/eclipse_funcs.c @@ -0,0 +1,362 @@ +/* + * eclipse_funcs.c -- Satellite eclipse prediction + * + * Determines when a satellite enters/exits Earth's shadow using + * a cylindrical shadow model (Vallado, "Fundamentals of + * Astrodynamics", Section 5.3). + * + * Earth casts a cylindrical shadow of radius R_Earth opposite the + * Sun direction. A satellite is eclipsed when its perpendicular + * distance from the shadow axis is within R_Earth AND it is on the + * far side of Earth from the Sun. + * + * Sun direction computed via VSOP87 (ecliptic J2000 -> equatorial + * J2000). TEME differs from J2000 by ~arcsec nutation residual, + * negligible at the 6378 km shadow boundary scale. + */ + +#include "postgres.h" +#include "fmgr.h" +#include "utils/timestamp.h" +#include "types.h" +#include "astro_math.h" +#include "vsop87.h" +#include "norad.h" +#include +#include + +PG_FUNCTION_INFO_V1(satellite_is_eclipsed); +PG_FUNCTION_INFO_V1(satellite_next_eclipse_entry); +PG_FUNCTION_INFO_V1(satellite_next_eclipse_exit); +PG_FUNCTION_INFO_V1(satellite_eclipse_fraction); + +#define DEG_TO_RAD_EC (M_PI / 180.0) +#define RAD_TO_DEG_EC (180.0 / M_PI) + +#define ECLIPSE_SCAN_STEP_JD (30.0 / 86400.0) /* 30 seconds */ +#define ECLIPSE_BISECT_TOL_JD (0.5 / 86400.0) /* 0.5 second */ +#define ECLIPSE_SEARCH_DAYS 7.0 + + +/* ---------------------------------------------------------------- + * Static helpers -- duplicated from pass_funcs.c per project + * convention (no cross-TU symbol coupling). + * ---------------------------------------------------------------- + */ + +static void +pg_tle_to_sat_code_ec(const pg_tle *src, tle_t *dst) +{ + memset(dst, 0, sizeof(tle_t)); + dst->epoch = src->epoch; + dst->xincl = src->inclination; + dst->xnodeo = src->raan; + dst->eo = src->eccentricity; + dst->omegao = src->arg_perigee; + dst->xmo = src->mean_anomaly; + dst->xno = src->mean_motion; + dst->xndt2o = src->mean_motion_dot; + dst->xndd6o = src->mean_motion_ddot; + dst->bstar = src->bstar; + dst->norad_number = src->norad_id; + dst->bulletin_number = src->elset_num; + dst->revolution_number = src->rev_num; + dst->classification = src->classification; + dst->ephemeris_type = src->ephemeris_type; + memcpy(dst->intl_desig, src->intl_desig, 9); +} + +static int +do_propagate_ec(const pg_tle *tle, double jd, double *pos, double *vel) +{ + tle_t sat; + double *params; + int is_deep; + int err; + double tsince; + + pg_tle_to_sat_code_ec(tle, &sat); + is_deep = select_ephemeris(&sat); + if (is_deep < 0) + return -99; + + tsince = jd_to_minutes_since_epoch(jd, sat.epoch); + params = palloc(sizeof(double) * N_SAT_PARAMS); + + if (is_deep) + { + SDP4_init(params, &sat); + err = SDP4(tsince, &sat, params, pos, vel); + } + else + { + SGP4_init(params, &sat); + err = SGP4(tsince, &sat, params, pos, vel); + } + + pfree(params); + return err; +} + + +/* + * Compute unit vector from Earth to Sun in equatorial J2000. + * + * Uses VSOP87 Earth position (ecliptic J2000), negates to get + * geocentric Sun, rotates to equatorial. Returns unit vector. + */ +static void +sun_direction_equ(double jd, double sun_dir[3]) +{ + double earth_xyz[6]; + double sun_ecl[3], sun_equ[3]; + double r; + + GetVsop87Coor(jd, 2, earth_xyz); /* VSOP87 body 2 = Earth */ + + /* Geocentric Sun = -Earth heliocentric */ + sun_ecl[0] = -earth_xyz[0]; + sun_ecl[1] = -earth_xyz[1]; + sun_ecl[2] = -earth_xyz[2]; + + /* Ecliptic J2000 -> equatorial J2000 */ + ecliptic_to_equatorial(sun_ecl, sun_equ); + + /* Normalize to unit vector */ + r = sqrt(sun_equ[0] * sun_equ[0] + + sun_equ[1] * sun_equ[1] + + sun_equ[2] * sun_equ[2]); + + sun_dir[0] = sun_equ[0] / r; + sun_dir[1] = sun_equ[1] / r; + sun_dir[2] = sun_equ[2] / r; +} + + +/* + * is_satellite_eclipsed_pos -- cylindrical shadow test + * + * sat_pos[3]: satellite position relative to Earth center (km, TEME/J2000) + * sun_dir[3]: unit vector from Earth toward Sun (J2000 equatorial) + * + * Eclipsed when: + * 1. sat dot sun_dir < 0 (satellite on shadow side of Earth) + * 2. perpendicular distance from shadow axis < R_Earth + */ +static bool +is_satellite_eclipsed_pos(const double sat_pos[3], const double sun_dir[3]) +{ + double proj, perp[3], perp_dist; + + /* Project satellite position onto Sun direction */ + proj = sat_pos[0] * sun_dir[0] + + sat_pos[1] * sun_dir[1] + + sat_pos[2] * sun_dir[2]; + + if (proj > 0.0) + return false; /* sunlit side of Earth */ + + /* Perpendicular vector from shadow axis */ + perp[0] = sat_pos[0] - proj * sun_dir[0]; + perp[1] = sat_pos[1] - proj * sun_dir[1]; + perp[2] = sat_pos[2] - proj * sun_dir[2]; + perp_dist = sqrt(perp[0] * perp[0] + + perp[1] * perp[1] + + perp[2] * perp[2]); + + return (perp_dist < WGS84_A); /* 6378.137 km */ +} + + +/* + * eclipse_state_at_jd -- compute eclipse state at a single time + * + * Returns true if eclipsed, false if sunlit. + * Returns false on propagation error (conservative: assume sunlit). + */ +static bool +eclipse_state_at_jd(const pg_tle *tle, double jd) +{ + double pos[3], vel[3]; + double sun_dir[3]; + int err; + + err = do_propagate_ec(tle, jd, pos, vel); + if (err != 0) + return false; /* propagation failed, assume sunlit */ + + sun_direction_equ(jd, sun_dir); + + return is_satellite_eclipsed_pos(pos, sun_dir); +} + + +/* ================================================================ + * satellite_is_eclipsed(tle, timestamptz) -> bool + * + * Point-in-time eclipse test. Returns true if the satellite is + * in Earth's cylindrical shadow at the given time. + * ================================================================ + */ +Datum +satellite_is_eclipsed(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + int64 ts = PG_GETARG_INT64(1); + double jd; + + jd = timestamptz_to_jd(ts); + + PG_RETURN_BOOL(eclipse_state_at_jd(tle, jd)); +} + + +/* ================================================================ + * satellite_next_eclipse_entry(tle, timestamptz) -> timestamptz + * + * Scans forward from the given time to find when the satellite + * next enters Earth's shadow. Searches up to 7 days. + * Returns NULL if no eclipse entry is found. + * ================================================================ + */ +Datum +satellite_next_eclipse_entry(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + int64 ts = PG_GETARG_INT64(1); + double jd, stop_jd; + bool prev_eclipsed, curr_eclipsed; + double lo, hi, mid; + + jd = timestamptz_to_jd(ts); + stop_jd = jd + ECLIPSE_SEARCH_DAYS; + + prev_eclipsed = eclipse_state_at_jd(tle, jd); + + while (jd < stop_jd) + { + jd += ECLIPSE_SCAN_STEP_JD; + if (jd > stop_jd) + jd = stop_jd; + + curr_eclipsed = eclipse_state_at_jd(tle, jd); + + /* Transition from sunlit to eclipsed */ + if (!prev_eclipsed && curr_eclipsed) + { + /* Bisect to refine entry time */ + lo = jd - ECLIPSE_SCAN_STEP_JD; + hi = jd; + while (hi - lo > ECLIPSE_BISECT_TOL_JD) + { + mid = (lo + hi) / 2.0; + if (eclipse_state_at_jd(tle, mid)) + hi = mid; + else + lo = mid; + } + + PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz((lo + hi) / 2.0)); + } + + prev_eclipsed = curr_eclipsed; + } + + PG_RETURN_NULL(); +} + + +/* ================================================================ + * satellite_next_eclipse_exit(tle, timestamptz) -> timestamptz + * + * Scans forward from the given time to find when the satellite + * next exits Earth's shadow (returns to sunlight). + * Searches up to 7 days. Returns NULL if no exit found. + * ================================================================ + */ +Datum +satellite_next_eclipse_exit(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + int64 ts = PG_GETARG_INT64(1); + double jd, stop_jd; + bool prev_eclipsed, curr_eclipsed; + double lo, hi, mid; + + jd = timestamptz_to_jd(ts); + stop_jd = jd + ECLIPSE_SEARCH_DAYS; + + prev_eclipsed = eclipse_state_at_jd(tle, jd); + + while (jd < stop_jd) + { + jd += ECLIPSE_SCAN_STEP_JD; + if (jd > stop_jd) + jd = stop_jd; + + curr_eclipsed = eclipse_state_at_jd(tle, jd); + + /* Transition from eclipsed to sunlit */ + if (prev_eclipsed && !curr_eclipsed) + { + /* Bisect to refine exit time */ + lo = jd - ECLIPSE_SCAN_STEP_JD; + hi = jd; + while (hi - lo > ECLIPSE_BISECT_TOL_JD) + { + mid = (lo + hi) / 2.0; + if (eclipse_state_at_jd(tle, mid)) + lo = mid; + else + hi = mid; + } + + PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz((lo + hi) / 2.0)); + } + + prev_eclipsed = curr_eclipsed; + } + + PG_RETURN_NULL(); +} + + +/* ================================================================ + * satellite_eclipse_fraction(tle, timestamptz, timestamptz) -> float8 + * + * Fraction of the given time window spent in eclipse [0.0, 1.0]. + * Scans the window at 30-second intervals and counts eclipsed samples. + * + * Useful for determining what portion of a pass is in shadow. + * ================================================================ + */ +Datum +satellite_eclipse_fraction(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + int64 start_ts = PG_GETARG_INT64(1); + int64 stop_ts = PG_GETARG_INT64(2); + double start_jd, stop_jd, jd; + int total_samples = 0; + int eclipsed_samples = 0; + + if (stop_ts <= start_ts) + ereport(ERROR, + (errcode(ERRCODE_INVALID_PARAMETER_VALUE), + errmsg("satellite_eclipse_fraction: stop time must be after start time"))); + + start_jd = timestamptz_to_jd(start_ts); + stop_jd = timestamptz_to_jd(stop_ts); + + for (jd = start_jd; jd <= stop_jd; jd += ECLIPSE_SCAN_STEP_JD) + { + if (eclipse_state_at_jd(tle, jd)) + eclipsed_samples++; + total_samples++; + } + + if (total_samples == 0) + PG_RETURN_FLOAT8(0.0); + + PG_RETURN_FLOAT8((double) eclipsed_samples / (double) total_samples); +} diff --git a/src/libration.h b/src/libration.h new file mode 100644 index 0000000..5c70b40 --- /dev/null +++ b/src/libration.h @@ -0,0 +1,22 @@ +/* + * libration.h -- Lunar optical libration (Meeus Ch. 53) + * + * Three components of the Moon's apparent wobble: + * l -- optical libration in longitude (degrees, [-8, +8]) + * b -- optical libration in latitude (degrees, [-7, +7]) + * p -- position angle of the Moon's axis (degrees) + */ + +#ifndef PG_ORRERY_LIBRATION_H +#define PG_ORRERY_LIBRATION_H + +typedef struct +{ + double l; /* libration in longitude, degrees */ + double b; /* libration in latitude, degrees */ + double p; /* position angle of axis, degrees */ +} lunar_libration; + +void compute_lunar_libration(double jd, lunar_libration *lib); + +#endif /* PG_ORRERY_LIBRATION_H */ diff --git a/src/libration_funcs.c b/src/libration_funcs.c new file mode 100644 index 0000000..b863be2 --- /dev/null +++ b/src/libration_funcs.c @@ -0,0 +1,368 @@ +/* + * libration_funcs.c -- Lunar libration and subsolar longitude + * + * Optical libration of the Moon (apparent wobble) computed from + * Meeus (1998) "Astronomical Algorithms", Chapter 53. + * + * Three components: + * l' -- libration in longitude (degrees, typically [-8, +8]) + * b' -- libration in latitude (degrees, typically [-7, +7]) + * P -- position angle of the Moon's axis (degrees) + * + * Also: selenographic subsolar longitude (terminator position). + * + * References: + * Meeus (1998) Chapters 22, 47, 53 + * Chapront-Touze & Chapront (1988) ELP2000-82B + */ + +#include "postgres.h" +#include "fmgr.h" +#include "funcapi.h" +#include "utils/timestamp.h" +#include "types.h" +#include "astro_math.h" +#include "elp82b.h" +#include "vsop87.h" +#include "precession.h" +#include "libration.h" +#include + +PG_FUNCTION_INFO_V1(moon_libration_longitude); +PG_FUNCTION_INFO_V1(moon_libration_latitude); +PG_FUNCTION_INFO_V1(moon_libration_position_angle); +PG_FUNCTION_INFO_V1(moon_libration); +PG_FUNCTION_INFO_V1(moon_subsolar_longitude); + + +/* Mean inclination of the lunar equator to the ecliptic (Meeus Ch. 53) */ +#define LUNAR_I_RAD (1.54242 * M_PI / 180.0) /* 1.54242 degrees */ + + +/* + * Moon's mean longitude referred to the mean equinox of date (F) + * and longitude of the ascending node (Omega). + * + * Meeus (1998) Table 22.A, using T = Julian centuries from J2000. + * F = mean longitude - RAAN (Meeus notation: F is the argument + * of latitude = mean anomaly + arg of perigee; but for libration + * we need the mean longitude L' = F + Omega). + */ +static void +lunar_fundamental_args(double jd, double *F_out, double *Omega_out, + double *Lprime_out) +{ + double T = (jd - J2000_JD) / 36525.0; + double T2 = T * T; + double T3 = T2 * T; + double T4 = T3 * T; + double Lprime, F, Omega; + + /* Moon's mean longitude L' (Meeus Eq. 47.1) */ + Lprime = 218.3164477 + + 481267.88123421 * T + - 0.0015786 * T2 + + T3 / 538841.0 + - T4 / 65194000.0; + + /* Moon's argument of latitude F (Meeus Eq. 47.5) */ + F = 93.2720950 + + 483202.0175233 * T + - 0.0036539 * T2 + - T3 / 3526000.0 + + T4 / 863310000.0; + + /* Longitude of the ascending node Omega (Meeus Eq. 47.7) */ + Omega = 125.0445479 + - 1934.1362891 * T + + 0.0020754 * T2 + + T3 / 467441.0 + - T4 / 60616000.0; + + /* Normalize to [0, 360) */ + Lprime = fmod(Lprime, 360.0); + if (Lprime < 0.0) Lprime += 360.0; + + F = fmod(F, 360.0); + if (F < 0.0) F += 360.0; + + Omega = fmod(Omega, 360.0); + if (Omega < 0.0) Omega += 360.0; + + *F_out = F; + *Omega_out = Omega; + *Lprime_out = Lprime; +} + + +/* + * compute_lunar_libration -- Meeus (1998) Ch. 53 + * + * Computes optical libration from the Moon's ecliptic coordinates, + * mean longitude, ascending node, and nutation. + */ +void +compute_lunar_libration(double jd, lunar_libration *lib) +{ + double moon_ecl[3]; + double lambda, beta; /* geocentric ecliptic long/lat, radians */ + double F_deg, Omega_deg, Lprime_deg; + double F, Omega; /* radians */ + double dpsi, deps; /* nutation, arcseconds */ + double eps_A, chi_A, omega_A, psi_A; /* precession, arcseconds */ + double eps_rad; /* mean obliquity of date, radians */ + double W; /* intermediate angle */ + double A, l_prime, b_prime; + double sin_W, cos_W; + double sin_beta, cos_beta; + double sin_I = sin(LUNAR_I_RAD); + double cos_I = cos(LUNAR_I_RAD); + double V, X, P; + double ra_moon; + + /* Moon geocentric ecliptic (ELP2000-82B gives ecliptic J2000 in AU) */ + GetElp82bCoor(jd, moon_ecl); + + /* Cartesian -> spherical ecliptic */ + lambda = atan2(moon_ecl[1], moon_ecl[0]); + if (lambda < 0.0) lambda += 2.0 * M_PI; + beta = asin(moon_ecl[2] / sqrt(moon_ecl[0] * moon_ecl[0] + + moon_ecl[1] * moon_ecl[1] + + moon_ecl[2] * moon_ecl[2])); + + /* Fundamental arguments */ + lunar_fundamental_args(jd, &F_deg, &Omega_deg, &Lprime_deg); + F = F_deg * DEG_TO_RAD; + Omega = Omega_deg * DEG_TO_RAD; + + /* Nutation in longitude */ + get_nutation_angles_iau2000b(jd, &dpsi, &deps); + + /* Mean obliquity of date */ + get_precession_angles_vondrak(jd, &eps_A, &chi_A, &omega_A, &psi_A); + eps_rad = eps_A * ARCSEC_TO_RAD; + + /* + * Meeus Ch. 53 formulas. + * + * W = lambda - dpsi - Omega + * where lambda is the Moon's geocentric ecliptic longitude, + * dpsi is nutation in longitude, and Omega is the ascending node. + * + * Note: lambda from ELP2000-82B is in J2000 ecliptic frame. + * For the libration formulas we need the apparent longitude, + * which requires adding nutation. Since W subtracts dpsi + * anyway, the J2000 value works: W = lambda_J2000 - Omega. + * The dpsi terms cancel when using the geometric longitude. + */ + W = lambda - Omega; + + sin_W = sin(W); + cos_W = cos(W); + sin_beta = sin(beta); + cos_beta = cos(beta); + + /* Optical libration in longitude (Meeus Eq. 53.1) */ + A = atan2(sin_W * cos_beta * cos_I - sin_beta * sin_I, + cos_W * cos_beta); + l_prime = A - F; + + /* Normalize to [-pi, pi) */ + l_prime = fmod(l_prime + M_PI, 2.0 * M_PI); + if (l_prime < 0.0) l_prime += 2.0 * M_PI; + l_prime -= M_PI; + + /* Optical libration in latitude (Meeus Eq. 53.2) */ + b_prime = asin(-sin_W * cos_beta * sin_I - sin_beta * cos_I); + + /* Position angle of the Moon's axis (Meeus Eq. 53.3) */ + V = Omega + dpsi * ARCSEC_TO_RAD + (eps_rad + deps * ARCSEC_TO_RAD) * 0.0; + + /* + * For the position angle P, we need the Moon's RA and Dec. + * Compute from ecliptic coordinates with nutation. + */ + { + double lambda_app, sin_eps, cos_eps; + + lambda_app = lambda + dpsi * ARCSEC_TO_RAD; + sin_eps = sin(eps_rad + deps * ARCSEC_TO_RAD); + cos_eps = cos(eps_rad + deps * ARCSEC_TO_RAD); + + ra_moon = atan2(sin(lambda_app) * cos_eps - tan(beta) * sin_eps, + cos(lambda_app)); + if (ra_moon < 0.0) ra_moon += 2.0 * M_PI; + } + + /* + * Position angle (Meeus Eq. 53.3): + * V = Omega + dpsi + eps * 0 (simplified; V uses Omega + dpsi) + * X = (Omega + dpsi) * cos(eps+deps) + ... but Meeus gives: + * + * Simplified: the position angle depends on the node longitude + * projected through the equatorial frame. + */ + V = Omega + dpsi * ARCSEC_TO_RAD; + X = sin(V + eps_rad + deps * ARCSEC_TO_RAD); + + P = asin(-X * cos(ra_moon) / cos(b_prime)) + + atan2(-sin_I * sin(V - Omega), + cos_I * sin_beta - sin_I * cos_beta * cos_W); + + /* + * Physical libration corrections (Meeus p. 373) are small + * (~0.02 deg) and omitted here for the optical model. + */ + + lib->l = l_prime * RAD_TO_DEG; + lib->b = b_prime * RAD_TO_DEG; + lib->p = fmod(P * RAD_TO_DEG, 360.0); + if (lib->p < 0.0) + lib->p += 360.0; +} + + +/* ================================================================ + * moon_libration_longitude(timestamptz) -> float8 + * Optical libration in longitude (degrees, typically [-8, +8]). + * ================================================================ + */ +Datum +moon_libration_longitude(PG_FUNCTION_ARGS) +{ + int64 ts = PG_GETARG_INT64(0); + double jd; + lunar_libration lib; + + jd = timestamptz_to_jd(ts); + compute_lunar_libration(jd, &lib); + + PG_RETURN_FLOAT8(lib.l); +} + + +/* ================================================================ + * moon_libration_latitude(timestamptz) -> float8 + * Optical libration in latitude (degrees, typically [-7, +7]). + * ================================================================ + */ +Datum +moon_libration_latitude(PG_FUNCTION_ARGS) +{ + int64 ts = PG_GETARG_INT64(0); + double jd; + lunar_libration lib; + + jd = timestamptz_to_jd(ts); + compute_lunar_libration(jd, &lib); + + PG_RETURN_FLOAT8(lib.b); +} + + +/* ================================================================ + * moon_libration_position_angle(timestamptz) -> float8 + * Position angle of the Moon's axis (degrees, [0, 360)). + * ================================================================ + */ +Datum +moon_libration_position_angle(PG_FUNCTION_ARGS) +{ + int64 ts = PG_GETARG_INT64(0); + double jd; + lunar_libration lib; + + jd = timestamptz_to_jd(ts); + compute_lunar_libration(jd, &lib); + + PG_RETURN_FLOAT8(lib.p); +} + + +/* ================================================================ + * moon_libration(timestamptz) -> record (l float8, b float8, p float8) + * All three libration values as a composite return. + * ================================================================ + */ +Datum +moon_libration(PG_FUNCTION_ARGS) +{ + int64 ts = PG_GETARG_INT64(0); + double jd; + lunar_libration lib; + TupleDesc tupdesc; + Datum values[3]; + bool nulls[3] = {false, false, false}; + + jd = timestamptz_to_jd(ts); + compute_lunar_libration(jd, &lib); + + if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) + ereport(ERROR, + (errcode(ERRCODE_FEATURE_NOT_SUPPORTED), + errmsg("function returning record called in context " + "that cannot accept type record"))); + + tupdesc = BlessTupleDesc(tupdesc); + + values[0] = Float8GetDatum(lib.l); + values[1] = Float8GetDatum(lib.b); + values[2] = Float8GetDatum(lib.p); + + PG_RETURN_DATUM(HeapTupleGetDatum( + heap_form_tuple(tupdesc, values, nulls))); +} + + +/* ================================================================ + * moon_subsolar_longitude(timestamptz) -> float8 + * + * Selenographic longitude of the sub-solar point (degrees, [0, 360)). + * Determines the terminator position on the Moon. + * + * This is the libration in longitude plus the selenographic + * colongitude of the Sun. The subsolar point's longitude + * tracks through 360 deg over a synodic month. + * + * Simplified computation: the subsolar longitude is approximately + * the difference between the Sun's ecliptic longitude and the Moon's + * ecliptic longitude, corrected for libration. + * ================================================================ + */ +Datum +moon_subsolar_longitude(PG_FUNCTION_ARGS) +{ + int64 ts = PG_GETARG_INT64(0); + double jd; + double earth_xyz[6], moon_ecl[3]; + double sun_lon, moon_lon; + double subsolar; + lunar_libration lib; + + jd = timestamptz_to_jd(ts); + + /* Sun's geocentric ecliptic longitude */ + GetVsop87Coor(jd, 2, earth_xyz); + sun_lon = atan2(-earth_xyz[1], -earth_xyz[0]); + if (sun_lon < 0.0) sun_lon += 2.0 * M_PI; + + /* Moon's geocentric ecliptic longitude */ + GetElp82bCoor(jd, moon_ecl); + moon_lon = atan2(moon_ecl[1], moon_ecl[0]); + if (moon_lon < 0.0) moon_lon += 2.0 * M_PI; + + /* Libration correction */ + compute_lunar_libration(jd, &lib); + + /* + * Subsolar longitude on the Moon's surface: + * The Sun illuminates from the direction (sun_lon - moon_lon) + * as seen from the Moon, corrected for libration. + */ + subsolar = (sun_lon - moon_lon) * RAD_TO_DEG + lib.l; + subsolar = fmod(subsolar, 360.0); + if (subsolar < 0.0) + subsolar += 360.0; + + PG_RETURN_FLOAT8(subsolar); +} diff --git a/src/magnitude_funcs.c b/src/magnitude_funcs.c index 59a9eda..2890471 100644 --- a/src/magnitude_funcs.c +++ b/src/magnitude_funcs.c @@ -1,9 +1,12 @@ /* - * magnitude_funcs.c -- Planet apparent visual magnitude + * magnitude_funcs.c -- Planet magnitude, solar elongation, phase fraction * * Uses the Mallama & Hilton (2018) magnitude model with * VSOP87 positions for distances and phase angles. * + * Solar elongation and planet phase reuse the same Sun-Planet-Earth + * triangle geometry, factored into compute_planet_geometry(). + * * Reference: Mallama & Hilton, "Computing Apparent Planetary * Magnitudes for The Astronomical Almanac", A&C vol. 25, 2018. */ @@ -17,6 +20,8 @@ #include PG_FUNCTION_INFO_V1(planet_magnitude); +PG_FUNCTION_INFO_V1(solar_elongation); +PG_FUNCTION_INFO_V1(planet_phase); /* @@ -117,52 +122,74 @@ static const double planet_v10[] = { /* - * Compute apparent visual magnitude of a planet from Earth. + * Shared Sun-Planet-Earth geometry. * - * Phase angle is the Sun-Planet-Earth angle, computed via the law - * of cosines from three heliocentric/geocentric distances. + * Computes the three distances (r, delta, R) and the phase angle + * (Sun-Planet-Earth angle) from VSOP87 positions. Used by + * magnitude, elongation, and phase functions. */ -static double -compute_planet_magnitude(int body_id, double jd) +typedef struct +{ + double r; /* Sun-Planet distance (AU) */ + double delta; /* Earth-Planet distance (AU) */ + double R; /* Sun-Earth distance (AU) */ + double i_deg; /* Phase angle, degrees (Sun-Planet-Earth vertex) */ +} planet_geometry; + +static void +compute_planet_geometry(int body_id, double jd, planet_geometry *geo) { double earth_xyz[6], planet_xyz[6]; - double geo[3]; - double r, delta, R; - double cos_i, i_deg; - double V; + double gv[3]; + double cos_i; int vsop_body = body_id - 1; /* pg_orrery 1-based -> VSOP87 0-based */ GetVsop87Coor(jd, 2, earth_xyz); /* Earth (VSOP87 body 2) */ GetVsop87Coor(jd, vsop_body, planet_xyz); /* target planet */ /* Heliocentric distance to planet */ - r = sqrt(planet_xyz[0] * planet_xyz[0] + - planet_xyz[1] * planet_xyz[1] + - planet_xyz[2] * planet_xyz[2]); + geo->r = sqrt(planet_xyz[0] * planet_xyz[0] + + planet_xyz[1] * planet_xyz[1] + + planet_xyz[2] * planet_xyz[2]); /* Geocentric vector and distance */ - geo[0] = planet_xyz[0] - earth_xyz[0]; - geo[1] = planet_xyz[1] - earth_xyz[1]; - geo[2] = planet_xyz[2] - earth_xyz[2]; - delta = sqrt(geo[0] * geo[0] + geo[1] * geo[1] + geo[2] * geo[2]); + gv[0] = planet_xyz[0] - earth_xyz[0]; + gv[1] = planet_xyz[1] - earth_xyz[1]; + gv[2] = planet_xyz[2] - earth_xyz[2]; + geo->delta = sqrt(gv[0] * gv[0] + gv[1] * gv[1] + gv[2] * gv[2]); /* Sun-Earth distance */ - R = sqrt(earth_xyz[0] * earth_xyz[0] + - earth_xyz[1] * earth_xyz[1] + - earth_xyz[2] * earth_xyz[2]); + geo->R = sqrt(earth_xyz[0] * earth_xyz[0] + + earth_xyz[1] * earth_xyz[1] + + earth_xyz[2] * earth_xyz[2]); - /* Phase angle via law of cosines: triangle Sun-Planet-Earth */ - cos_i = (r * r + delta * delta - R * R) / (2.0 * r * delta); + /* Phase angle via law of cosines: vertex at planet */ + cos_i = (geo->r * geo->r + geo->delta * geo->delta - geo->R * geo->R) + / (2.0 * geo->r * geo->delta); if (cos_i > 1.0) cos_i = 1.0; if (cos_i < -1.0) cos_i = -1.0; - i_deg = acos(cos_i) * RAD_TO_DEG; + geo->i_deg = acos(cos_i) * RAD_TO_DEG; +} - /* Mallama & Hilton (2018) magnitude with full phase correction */ - V = planet_v10[body_id] - + 5.0 * log10(r * delta) - + phase_correction(body_id, i_deg); - return V; +/* + * Validate planet body_id for magnitude/elongation/phase. + * Must be 1-8 (Mercury-Neptune), not 3 (Earth). + */ +static void +validate_planet_body_id(int body_id, const char *func_name) +{ + if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("%s: body_id %d must be 1-8 (Mercury-Neptune)", + func_name, body_id))); + + if (body_id == BODY_EARTH) + ereport(ERROR, + (errcode(ERRCODE_FEATURE_NOT_SUPPORTED), + errmsg("%s: cannot compute for Earth from Earth", + func_name))); } @@ -182,23 +209,90 @@ compute_planet_magnitude(int body_id, double jd) Datum planet_magnitude(PG_FUNCTION_ARGS) { - int32 body_id = PG_GETARG_INT32(0); - int64 ts = PG_GETARG_INT64(1); - double jd, mag; + int32 body_id = PG_GETARG_INT32(0); + int64 ts = PG_GETARG_INT64(1); + double jd; + planet_geometry geo; + double V; - if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE) - ereport(ERROR, - (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), - errmsg("planet_magnitude: body_id %d must be 1-8 (Mercury-Neptune)", - body_id))); - - if (body_id == BODY_EARTH) - ereport(ERROR, - (errcode(ERRCODE_FEATURE_NOT_SUPPORTED), - errmsg("cannot compute magnitude for Earth from Earth"))); + validate_planet_body_id(body_id, "planet_magnitude"); jd = timestamptz_to_jd(ts); - mag = compute_planet_magnitude(body_id, jd); + compute_planet_geometry(body_id, jd, &geo); - PG_RETURN_FLOAT8(mag); + V = planet_v10[body_id] + + 5.0 * log10(geo.r * geo.delta) + + phase_correction(body_id, geo.i_deg); + + PG_RETURN_FLOAT8(V); +} + + +/* ================================================================ + * solar_elongation(body_id int4, timestamptz) -> float8 + * + * Sun-Earth-Planet angle in degrees [0, 180]. + * How far a planet appears from the Sun in the sky. + * + * Uses law of cosines with vertex at Earth: + * cos(elong) = (R^2 + delta^2 - r^2) / (2 * R * delta) + * + * Mercury max ~28 deg, Venus max ~47 deg. + * Superior planets can reach ~180 deg (opposition). + * ================================================================ + */ +Datum +solar_elongation(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int64 ts = PG_GETARG_INT64(1); + double jd; + planet_geometry geo; + double cos_elong, elong_deg; + + validate_planet_body_id(body_id, "solar_elongation"); + + jd = timestamptz_to_jd(ts); + compute_planet_geometry(body_id, jd, &geo); + + /* Law of cosines, vertex at Earth */ + cos_elong = (geo.R * geo.R + geo.delta * geo.delta - geo.r * geo.r) + / (2.0 * geo.R * geo.delta); + if (cos_elong > 1.0) cos_elong = 1.0; + if (cos_elong < -1.0) cos_elong = -1.0; + elong_deg = acos(cos_elong) * RAD_TO_DEG; + + PG_RETURN_FLOAT8(elong_deg); +} + + +/* ================================================================ + * planet_phase(body_id int4, timestamptz) -> float8 + * + * Illuminated fraction of a planet's disk as seen from Earth [0, 1]. + * k = (1 + cos(i)) / 2 + * where i is the phase angle (Sun-Planet-Earth). + * + * Inner planets vary dramatically (Venus crescent). + * Outer planets are always near 1.0. + * ================================================================ + */ +Datum +planet_phase(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int64 ts = PG_GETARG_INT64(1); + double jd; + planet_geometry geo; + double i_rad, k; + + validate_planet_body_id(body_id, "planet_phase"); + + jd = timestamptz_to_jd(ts); + compute_planet_geometry(body_id, jd, &geo); + + i_rad = geo.i_deg * DEG_TO_RAD; + k = (1.0 + cos(i_rad)) / 2.0; + + PG_RETURN_FLOAT8(k); } diff --git a/test/expected/v016_features.out b/test/expected/v016_features.out index daa2b9b..9e5e074 100644 --- a/test/expected/v016_features.out +++ b/test/expected/v016_features.out @@ -238,6 +238,6 @@ FROM (VALUES (1),(2),(4),(5),(6),(7),(8)) AS t(body_id); DO $$ BEGIN PERFORM planet_magnitude(0, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'body_id=0(Sun): %', SQLERRM; END $$; NOTICE: body_id=0(Sun): planet_magnitude: body_id 0 must be 1-8 (Mercury-Neptune) DO $$ BEGIN PERFORM planet_magnitude(3, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'body_id=3(Earth): %', SQLERRM; END $$; -NOTICE: body_id=3(Earth): cannot compute magnitude for Earth from Earth +NOTICE: body_id=3(Earth): planet_magnitude: cannot compute for Earth from Earth DO $$ BEGIN PERFORM planet_magnitude(9, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'body_id=9: %', SQLERRM; END $$; NOTICE: body_id=9: planet_magnitude: body_id 9 must be 1-8 (Mercury-Neptune) diff --git a/test/expected/v017_features.out b/test/expected/v017_features.out new file mode 100644 index 0000000..6af06c5 --- /dev/null +++ b/test/expected/v017_features.out @@ -0,0 +1,285 @@ +-- v017_features.sql -- Tests for v0.17.0: solar elongation, planet phase, +-- satellite eclipse, observing night quality, lunar libration +-- +-- Verifies all 12 new functions added in v0.17.0. +CREATE EXTENSION IF NOT EXISTS pg_orrery; +NOTICE: extension "pg_orrery" already exists, skipping +-- ============================================================ +-- Solar elongation: Mercury always < 28 deg +-- ============================================================ +SELECT solar_elongation(1, '2024-01-15 00:00:00+00'::timestamptz) < 28.0 + AS mercury_max_elongation; + mercury_max_elongation +------------------------ + t +(1 row) + +-- ============================================================ +-- Solar elongation: Venus always < 47 deg +-- ============================================================ +SELECT solar_elongation(2, '2024-01-15 00:00:00+00'::timestamptz) < 47.5 + AS venus_max_elongation; + venus_max_elongation +---------------------- + t +(1 row) + +-- ============================================================ +-- Solar elongation: Mars can exceed 90 deg (superior planet) +-- Use a date near opposition (2024-01-12 Mars at elongation ~180) +-- At least verify it can be large for outer planets +-- ============================================================ +SELECT solar_elongation(4, '2024-12-08 00:00:00+00'::timestamptz) > 50.0 + AS mars_large_elongation; + mars_large_elongation +----------------------- + t +(1 row) + +-- ============================================================ +-- Solar elongation: always [0, 180] +-- ============================================================ +SELECT bool_and( + solar_elongation(body_id, '2024-06-15 00:00:00+00'::timestamptz) >= 0.0 + AND solar_elongation(body_id, '2024-06-15 00:00:00+00'::timestamptz) <= 180.0 +) AS elongation_in_range +FROM (VALUES (1),(2),(4),(5),(6),(7),(8)) AS t(body_id); + elongation_in_range +--------------------- + t +(1 row) + +-- ============================================================ +-- Solar elongation: error on body_id 0, 3, 9 +-- ============================================================ +DO $$ BEGIN PERFORM solar_elongation(0, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'elong body_id=0: %', SQLERRM; END $$; +NOTICE: elong body_id=0: solar_elongation: body_id 0 must be 1-8 (Mercury-Neptune) +DO $$ BEGIN PERFORM solar_elongation(3, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'elong body_id=3: %', SQLERRM; END $$; +NOTICE: elong body_id=3: solar_elongation: cannot compute for Earth from Earth +DO $$ BEGIN PERFORM solar_elongation(9, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'elong body_id=9: %', SQLERRM; END $$; +NOTICE: elong body_id=9: solar_elongation: body_id 9 must be 1-8 (Mercury-Neptune) +-- ============================================================ +-- Planet phase: Jupiter always near 1.0 (outer planet) +-- ============================================================ +SELECT planet_phase(5, '2024-01-15 00:00:00+00'::timestamptz) > 0.95 + AS jupiter_nearly_full; + jupiter_nearly_full +--------------------- + t +(1 row) + +-- ============================================================ +-- Planet phase: Neptune always near 1.0 +-- ============================================================ +SELECT planet_phase(8, '2024-06-15 00:00:00+00'::timestamptz) > 0.99 + AS neptune_nearly_full; + neptune_nearly_full +--------------------- + t +(1 row) + +-- ============================================================ +-- Planet phase: Venus varies significantly (inner planet) +-- Check it's in valid range +-- ============================================================ +SELECT planet_phase(2, '2024-06-01 12:00:00+00'::timestamptz) BETWEEN 0.0 AND 1.0 + AS venus_phase_valid; + venus_phase_valid +------------------- + t +(1 row) + +-- ============================================================ +-- Planet phase: always [0, 1] for all planets +-- ============================================================ +SELECT bool_and( + planet_phase(body_id, '2024-01-15 00:00:00+00'::timestamptz) >= 0.0 + AND planet_phase(body_id, '2024-01-15 00:00:00+00'::timestamptz) <= 1.0 +) AS phase_in_range +FROM (VALUES (1),(2),(4),(5),(6),(7),(8)) AS t(body_id); + phase_in_range +---------------- + t +(1 row) + +-- ============================================================ +-- Planet phase: error cases match elongation +-- ============================================================ +DO $$ BEGIN PERFORM planet_phase(3, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'phase body_id=3: %', SQLERRM; END $$; +NOTICE: phase body_id=3: planet_phase: cannot compute for Earth from Earth +-- ============================================================ +-- Satellite eclipse: ISS point-in-time test +-- (At night the ISS can be eclipsed; just verify function returns bool) +-- ============================================================ +SELECT satellite_is_eclipsed( + E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle, + '2024-01-01 12:00:00+00'::timestamptz +) IS NOT NULL + AS eclipse_returns_bool; + eclipse_returns_bool +---------------------- + t +(1 row) + +-- ============================================================ +-- Satellite eclipse: next entry/exit return timestamps or NULL +-- ============================================================ +SELECT satellite_next_eclipse_entry( + E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle, + '2024-01-01 12:00:00+00'::timestamptz +) > '2024-01-01 12:00:00+00'::timestamptz + AS entry_in_future; + entry_in_future +----------------- + t +(1 row) + +SELECT satellite_next_eclipse_exit( + E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle, + '2024-01-01 12:00:00+00'::timestamptz +) > '2024-01-01 12:00:00+00'::timestamptz + AS exit_in_future; + exit_in_future +---------------- + t +(1 row) + +-- ============================================================ +-- Satellite eclipse: fraction in [0, 1] for a 2-hour window +-- ============================================================ +SELECT satellite_eclipse_fraction( + E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle, + '2024-01-01 12:00:00+00'::timestamptz, + '2024-01-01 14:00:00+00'::timestamptz +) BETWEEN 0.0 AND 1.0 + AS eclipse_fraction_valid; + eclipse_fraction_valid +------------------------ + t +(1 row) + +-- ============================================================ +-- Observing night quality: polar summer at 65N = 'poor' +-- (no astronomical darkness in June) +-- ============================================================ +SELECT observing_night_quality('(65.0,25.0,0)'::observer, '2024-06-21 12:00:00+00'::timestamptz) = 'poor' + AS polar_summer_poor; + polar_summer_poor +------------------- + t +(1 row) + +-- ============================================================ +-- Observing night quality: winter mid-latitude returns valid rating +-- ============================================================ +SELECT observing_night_quality('(43.7,-116.4,800)'::observer, '2024-12-21 12:00:00+00'::timestamptz) + IN ('excellent', 'good', 'fair', 'poor') + AS winter_valid_rating; + winter_valid_rating +--------------------- + t +(1 row) + +-- ============================================================ +-- Lunar libration: longitude in [-8, 8] range +-- ============================================================ +SELECT moon_libration_longitude('2024-01-15 00:00:00+00'::timestamptz) BETWEEN -8.5 AND 8.5 + AS libration_lon_in_range; + libration_lon_in_range +------------------------ + t +(1 row) + +SELECT moon_libration_longitude('2024-06-15 00:00:00+00'::timestamptz) BETWEEN -8.5 AND 8.5 + AS libration_lon_in_range_2; + libration_lon_in_range_2 +-------------------------- + t +(1 row) + +-- ============================================================ +-- Lunar libration: latitude in [-7, 7] range +-- ============================================================ +SELECT moon_libration_latitude('2024-01-15 00:00:00+00'::timestamptz) BETWEEN -7.5 AND 7.5 + AS libration_lat_in_range; + libration_lat_in_range +------------------------ + t +(1 row) + +SELECT moon_libration_latitude('2024-06-15 00:00:00+00'::timestamptz) BETWEEN -7.5 AND 7.5 + AS libration_lat_in_range_2; + libration_lat_in_range_2 +-------------------------- + t +(1 row) + +-- ============================================================ +-- Lunar libration: position angle in [0, 360) +-- ============================================================ +SELECT moon_libration_position_angle('2024-01-15 00:00:00+00'::timestamptz) BETWEEN -1.0 AND 361.0 + AS libration_pa_in_range; + libration_pa_in_range +----------------------- + t +(1 row) + +-- ============================================================ +-- Lunar libration: composite returns same as individual functions +-- ============================================================ +SELECT abs((moon_libration('2024-01-15 00:00:00+00'::timestamptz)).l + - moon_libration_longitude('2024-01-15 00:00:00+00'::timestamptz)) < 0.001 + AS composite_matches_lon; + composite_matches_lon +----------------------- + t +(1 row) + +SELECT abs((moon_libration('2024-01-15 00:00:00+00'::timestamptz)).b + - moon_libration_latitude('2024-01-15 00:00:00+00'::timestamptz)) < 0.001 + AS composite_matches_lat; + composite_matches_lat +----------------------- + t +(1 row) + +-- ============================================================ +-- Lunar libration: changes over time (not constant) +-- ============================================================ +SELECT moon_libration_longitude('2024-01-01 00:00:00+00'::timestamptz) + != moon_libration_longitude('2024-01-15 00:00:00+00'::timestamptz) + AS libration_varies; + libration_varies +------------------ + t +(1 row) + +-- ============================================================ +-- Subsolar longitude: in [0, 360) range +-- ============================================================ +SELECT moon_subsolar_longitude('2024-01-15 00:00:00+00'::timestamptz) BETWEEN 0.0 AND 360.0 + AS subsolar_in_range; + subsolar_in_range +------------------- + t +(1 row) + +SELECT moon_subsolar_longitude('2024-06-15 00:00:00+00'::timestamptz) BETWEEN 0.0 AND 360.0 + AS subsolar_in_range_2; + subsolar_in_range_2 +--------------------- + t +(1 row) + +-- ============================================================ +-- Subsolar longitude: changes significantly over synodic month +-- (full 360 degrees over ~29.5 days) +-- ============================================================ +SELECT abs(moon_subsolar_longitude('2024-01-01 00:00:00+00'::timestamptz) + - moon_subsolar_longitude('2024-01-15 00:00:00+00'::timestamptz)) > 10.0 + AS subsolar_moves; + subsolar_moves +---------------- + t +(1 row) + diff --git a/test/sql/v017_features.sql b/test/sql/v017_features.sql new file mode 100644 index 0000000..0003b2a --- /dev/null +++ b/test/sql/v017_features.sql @@ -0,0 +1,204 @@ +-- v017_features.sql -- Tests for v0.17.0: solar elongation, planet phase, +-- satellite eclipse, observing night quality, lunar libration +-- +-- Verifies all 12 new functions added in v0.17.0. +CREATE EXTENSION IF NOT EXISTS pg_orrery; + +-- ============================================================ +-- Solar elongation: Mercury always < 28 deg +-- ============================================================ + +SELECT solar_elongation(1, '2024-01-15 00:00:00+00'::timestamptz) < 28.0 + AS mercury_max_elongation; + +-- ============================================================ +-- Solar elongation: Venus always < 47 deg +-- ============================================================ + +SELECT solar_elongation(2, '2024-01-15 00:00:00+00'::timestamptz) < 47.5 + AS venus_max_elongation; + +-- ============================================================ +-- Solar elongation: Mars can exceed 90 deg (superior planet) +-- Use a date near opposition (2024-01-12 Mars at elongation ~180) +-- At least verify it can be large for outer planets +-- ============================================================ + +SELECT solar_elongation(4, '2024-12-08 00:00:00+00'::timestamptz) > 50.0 + AS mars_large_elongation; + +-- ============================================================ +-- Solar elongation: always [0, 180] +-- ============================================================ + +SELECT bool_and( + solar_elongation(body_id, '2024-06-15 00:00:00+00'::timestamptz) >= 0.0 + AND solar_elongation(body_id, '2024-06-15 00:00:00+00'::timestamptz) <= 180.0 +) AS elongation_in_range +FROM (VALUES (1),(2),(4),(5),(6),(7),(8)) AS t(body_id); + +-- ============================================================ +-- Solar elongation: error on body_id 0, 3, 9 +-- ============================================================ + +DO $$ BEGIN PERFORM solar_elongation(0, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'elong body_id=0: %', SQLERRM; END $$; +DO $$ BEGIN PERFORM solar_elongation(3, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'elong body_id=3: %', SQLERRM; END $$; +DO $$ BEGIN PERFORM solar_elongation(9, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'elong body_id=9: %', SQLERRM; END $$; + +-- ============================================================ +-- Planet phase: Jupiter always near 1.0 (outer planet) +-- ============================================================ + +SELECT planet_phase(5, '2024-01-15 00:00:00+00'::timestamptz) > 0.95 + AS jupiter_nearly_full; + +-- ============================================================ +-- Planet phase: Neptune always near 1.0 +-- ============================================================ + +SELECT planet_phase(8, '2024-06-15 00:00:00+00'::timestamptz) > 0.99 + AS neptune_nearly_full; + +-- ============================================================ +-- Planet phase: Venus varies significantly (inner planet) +-- Check it's in valid range +-- ============================================================ + +SELECT planet_phase(2, '2024-06-01 12:00:00+00'::timestamptz) BETWEEN 0.0 AND 1.0 + AS venus_phase_valid; + +-- ============================================================ +-- Planet phase: always [0, 1] for all planets +-- ============================================================ + +SELECT bool_and( + planet_phase(body_id, '2024-01-15 00:00:00+00'::timestamptz) >= 0.0 + AND planet_phase(body_id, '2024-01-15 00:00:00+00'::timestamptz) <= 1.0 +) AS phase_in_range +FROM (VALUES (1),(2),(4),(5),(6),(7),(8)) AS t(body_id); + +-- ============================================================ +-- Planet phase: error cases match elongation +-- ============================================================ + +DO $$ BEGIN PERFORM planet_phase(3, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'phase body_id=3: %', SQLERRM; END $$; + +-- ============================================================ +-- Satellite eclipse: ISS point-in-time test +-- (At night the ISS can be eclipsed; just verify function returns bool) +-- ============================================================ + +SELECT satellite_is_eclipsed( + E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle, + '2024-01-01 12:00:00+00'::timestamptz +) IS NOT NULL + AS eclipse_returns_bool; + +-- ============================================================ +-- Satellite eclipse: next entry/exit return timestamps or NULL +-- ============================================================ + +SELECT satellite_next_eclipse_entry( + E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle, + '2024-01-01 12:00:00+00'::timestamptz +) > '2024-01-01 12:00:00+00'::timestamptz + AS entry_in_future; + +SELECT satellite_next_eclipse_exit( + E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle, + '2024-01-01 12:00:00+00'::timestamptz +) > '2024-01-01 12:00:00+00'::timestamptz + AS exit_in_future; + +-- ============================================================ +-- Satellite eclipse: fraction in [0, 1] for a 2-hour window +-- ============================================================ + +SELECT satellite_eclipse_fraction( + E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle, + '2024-01-01 12:00:00+00'::timestamptz, + '2024-01-01 14:00:00+00'::timestamptz +) BETWEEN 0.0 AND 1.0 + AS eclipse_fraction_valid; + +-- ============================================================ +-- Observing night quality: polar summer at 65N = 'poor' +-- (no astronomical darkness in June) +-- ============================================================ + +SELECT observing_night_quality('(65.0,25.0,0)'::observer, '2024-06-21 12:00:00+00'::timestamptz) = 'poor' + AS polar_summer_poor; + +-- ============================================================ +-- Observing night quality: winter mid-latitude returns valid rating +-- ============================================================ + +SELECT observing_night_quality('(43.7,-116.4,800)'::observer, '2024-12-21 12:00:00+00'::timestamptz) + IN ('excellent', 'good', 'fair', 'poor') + AS winter_valid_rating; + +-- ============================================================ +-- Lunar libration: longitude in [-8, 8] range +-- ============================================================ + +SELECT moon_libration_longitude('2024-01-15 00:00:00+00'::timestamptz) BETWEEN -8.5 AND 8.5 + AS libration_lon_in_range; + +SELECT moon_libration_longitude('2024-06-15 00:00:00+00'::timestamptz) BETWEEN -8.5 AND 8.5 + AS libration_lon_in_range_2; + +-- ============================================================ +-- Lunar libration: latitude in [-7, 7] range +-- ============================================================ + +SELECT moon_libration_latitude('2024-01-15 00:00:00+00'::timestamptz) BETWEEN -7.5 AND 7.5 + AS libration_lat_in_range; + +SELECT moon_libration_latitude('2024-06-15 00:00:00+00'::timestamptz) BETWEEN -7.5 AND 7.5 + AS libration_lat_in_range_2; + +-- ============================================================ +-- Lunar libration: position angle in [0, 360) +-- ============================================================ + +SELECT moon_libration_position_angle('2024-01-15 00:00:00+00'::timestamptz) BETWEEN -1.0 AND 361.0 + AS libration_pa_in_range; + +-- ============================================================ +-- Lunar libration: composite returns same as individual functions +-- ============================================================ + +SELECT abs((moon_libration('2024-01-15 00:00:00+00'::timestamptz)).l + - moon_libration_longitude('2024-01-15 00:00:00+00'::timestamptz)) < 0.001 + AS composite_matches_lon; + +SELECT abs((moon_libration('2024-01-15 00:00:00+00'::timestamptz)).b + - moon_libration_latitude('2024-01-15 00:00:00+00'::timestamptz)) < 0.001 + AS composite_matches_lat; + +-- ============================================================ +-- Lunar libration: changes over time (not constant) +-- ============================================================ + +SELECT moon_libration_longitude('2024-01-01 00:00:00+00'::timestamptz) + != moon_libration_longitude('2024-01-15 00:00:00+00'::timestamptz) + AS libration_varies; + +-- ============================================================ +-- Subsolar longitude: in [0, 360) range +-- ============================================================ + +SELECT moon_subsolar_longitude('2024-01-15 00:00:00+00'::timestamptz) BETWEEN 0.0 AND 360.0 + AS subsolar_in_range; + +SELECT moon_subsolar_longitude('2024-06-15 00:00:00+00'::timestamptz) BETWEEN 0.0 AND 360.0 + AS subsolar_in_range_2; + +-- ============================================================ +-- Subsolar longitude: changes significantly over synodic month +-- (full 360 degrees over ~29.5 days) +-- ============================================================ + +SELECT abs(moon_subsolar_longitude('2024-01-01 00:00:00+00'::timestamptz) + - moon_subsolar_longitude('2024-01-15 00:00:00+00'::timestamptz)) > 10.0 + AS subsolar_moves;