diff --git a/Makefile b/Makefile index 8b6b5d9..71630f4 100644 --- a/Makefile +++ b/Makefile @@ -4,7 +4,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.3.0.sql sql/pg_orrery--0.2.0--0.3.0.sql \ sql/pg_orrery--0.4.0.sql sql/pg_orrery--0.3.0--0.4.0.sql \ sql/pg_orrery--0.5.0.sql sql/pg_orrery--0.4.0--0.5.0.sql \ - sql/pg_orrery--0.6.0.sql sql/pg_orrery--0.5.0--0.6.0.sql + sql/pg_orrery--0.6.0.sql sql/pg_orrery--0.5.0--0.6.0.sql \ + sql/pg_orrery--0.7.0.sql sql/pg_orrery--0.6.0--0.7.0.sql # Our extension C sources OBJS = src/pg_orrery.o src/tle_type.o src/eci_type.o src/observer_type.o \ @@ -16,7 +17,8 @@ OBJS = src/pg_orrery.o src/tle_type.o src/eci_type.o src/observer_type.o \ src/moon_funcs.o src/radio_funcs.o \ src/lambert.o src/transfer_funcs.o \ src/de_reader.o src/eph_provider.o src/de_funcs.o \ - src/od_math.o src/od_iod.o src/od_solver.o src/od_funcs.o + src/od_math.o src/od_iod.o src/od_solver.o src/od_funcs.o \ + src/spgist_tle.o # Vendored SGP4/SDP4 sources (pure C, from Bill Gray's sat_code, MIT license) SGP4_DIR = src/sgp4 @@ -31,7 +33,7 @@ OBJS += $(SGP4_OBJS) # Regression tests REGRESS = tle_parse sgp4_propagate coord_transforms pass_prediction gist_index convenience \ star_observe kepler_comet planet_observe moon_observe lambert_transfer \ - de_ephemeris od_fit vallado_518 + de_ephemeris od_fit spgist_tle vallado_518 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 bd38ca9..af75a62 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.6.0' +default_version = '0.7.0' module_pathname = '$libdir/pg_orrery' relocatable = true diff --git a/sql/pg_orrery--0.6.0--0.7.0.sql b/sql/pg_orrery--0.6.0--0.7.0.sql new file mode 100644 index 0000000..0f344a5 --- /dev/null +++ b/sql/pg_orrery--0.6.0--0.7.0.sql @@ -0,0 +1,76 @@ +-- 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(observer_window, tle) RETURNS boolean + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; + +COMMENT ON FUNCTION tle_visibility_possible(observer_window, tle) 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) +-- ============================================================ + +CREATE OPERATOR &? ( + LEFTARG = observer_window, + RIGHTARG = tle, + FUNCTION = tle_visibility_possible, + RESTRICT = contsel, + JOIN = contjoinsel +); + +COMMENT ON OPERATOR &? (observer_window, tle) 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 &? (observer_window, tle), + 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); diff --git a/sql/pg_orrery--0.7.0.sql b/sql/pg_orrery--0.7.0.sql new file mode 100644 index 0000000..77176c5 --- /dev/null +++ b/sql/pg_orrery--0.7.0.sql @@ -0,0 +1,961 @@ +-- 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 'Minimum altitude-band separation in km (0 if overlapping). Altitude-only — does not account for inclination. Use && for 2-D filtering.'; + + +-- ============================================================ +-- 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(observer_window, tle) RETURNS boolean + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; + +COMMENT ON FUNCTION tle_visibility_possible(observer_window, tle) 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) +-- ============================================================ + +CREATE OPERATOR &? ( + LEFTARG = observer_window, + RIGHTARG = tle, + FUNCTION = tle_visibility_possible, + RESTRICT = contsel, + JOIN = contjoinsel +); + +COMMENT ON OPERATOR &? (observer_window, tle) 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 &? (observer_window, tle), + 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); diff --git a/src/spgist_tle.c b/src/spgist_tle.c new file mode 100644 index 0000000..18ef774 --- /dev/null +++ b/src/spgist_tle.c @@ -0,0 +1,841 @@ +/* + * spgist_tle.c -- SP-GiST operator class for orbital trie on TLE + * + * 2-level space-partitioning trie: + * L0: Semi-major axis (altitude regime, from Kepler's 3rd law) + * L1: Inclination (ground-track latitude coverage) + * + * Query-time RAAN filter at leaf level projects ascending node to + * the query midpoint via J2 secular precession and checks alignment + * with the observer's local sidereal position. + * + * The &? operator answers "could this satellite be visible from this + * observer during this time window?" -- a conservative superset of + * the actual answer. Survivors go through SGP4 propagation. + * + * Equal-population splits: picksplit sorts by the level's element + * and divides into floor(sqrt(n)) bins, clamped [2,16]. Dense LEO + * gets finer SMA bins than sparse MEO/GEO. + */ + +#include "postgres.h" +#include "fmgr.h" +#include "access/spgist.h" +#include "access/htup_details.h" +#include "catalog/pg_type_d.h" +#include "utils/builtins.h" +#include "utils/timestamp.h" +#include "executor/executor.h" +#include "types.h" +#include "astro_math.h" +#include +#include + +PG_FUNCTION_INFO_V1(spgist_tle_config); +PG_FUNCTION_INFO_V1(spgist_tle_choose); +PG_FUNCTION_INFO_V1(spgist_tle_picksplit); +PG_FUNCTION_INFO_V1(spgist_tle_inner_consistent); +PG_FUNCTION_INFO_V1(spgist_tle_leaf_consistent); +PG_FUNCTION_INFO_V1(tle_visibility_possible); + +/* Max trie depth: L0 (SMA) + L1 (inclination) */ +#define SPGIST_TLE_MAX_LEVEL 2 + +/* Earth angular rotation rate in radians/day */ +#define EARTH_ROT_RAD_PER_DAY (2.0 * M_PI) + +/* Seconds per day */ +#define SECONDS_PER_DAY 86400.0 + +/* Minutes per day */ +#define MINUTES_PER_DAY 1440.0 + + +/* ---------------------------------------------------------------- + * Helper: semi-major axis in km from mean motion + * + * Kepler's 3rd law with WGS-72: a = (KE / n)^(2/3) * AE + * where n is in radians/minute (TLE internal units). + * ---------------------------------------------------------------- + */ +static inline double +tle_sma_km(const pg_tle *tle) +{ + double n = tle->mean_motion; + + if (n <= 0.0) + return 0.0; + return pow(WGS72_KE / n, 2.0 / 3.0) * WGS72_AE; +} + + +/* ---------------------------------------------------------------- + * Helper: perigee altitude in km above Earth's surface + * ---------------------------------------------------------------- + */ +static inline double +tle_perigee_alt_km(const pg_tle *tle) +{ + double a = tle_sma_km(tle); + + return a * (1.0 - tle->eccentricity) - WGS72_AE; +} + + +/* ---------------------------------------------------------------- + * Helper: apogee altitude in km above Earth's surface + * ---------------------------------------------------------------- + */ +static inline double +tle_apogee_alt_km(const pg_tle *tle) +{ + double a = tle_sma_km(tle); + + return a * (1.0 + tle->eccentricity) - WGS72_AE; +} + + +/* ---------------------------------------------------------------- + * Helper: maximum satellite altitude visible at a given min elevation + * + * At min_el degrees elevation, the observer can see a satellite at + * most this far above the surface. Conservative upper bound using + * the Earth-center angle geometry: + * h_max = Re * (1/cos(el) - 1) roughly, but for a safe upper + * bound we use the slant range limit. + * + * For min_el = 10 deg, practical limit is ~2500 km for LEO passes. + * For min_el = 0 deg, theoretical limit extends to GEO+. + * We use a generous bound: if min_el < 5, return 50000 (no filter). + * Otherwise compute from geometry. + * ---------------------------------------------------------------- + */ +static inline double +max_visible_altitude_km(double min_el_deg) +{ + double el_rad, sin_el; + double rho, h_max; + + if (min_el_deg < 5.0) + return 50000.0; /* effectively no altitude filter */ + + el_rad = min_el_deg * DEG_TO_RAD; + sin_el = sin(el_rad); + + /* + * Maximum slant range rho where a satellite at altitude h is visible + * at elevation el. From the geometry: + * rho = Re * (sqrt((h/Re + 1)^2 - cos^2(el)) - sin(el)) + * Invert for h given a max practical rho. We take rho_max = 5000 km + * (well beyond any LEO pass) and solve for h. + * + * Simpler conservative bound: h_max = rho_max / sin(el) for el > 0. + */ + rho = 5000.0; + h_max = sqrt(rho * rho + 2.0 * WGS72_AE * rho * sin_el + + WGS72_AE * WGS72_AE) - WGS72_AE; + + return h_max; +} + + +/* ---------------------------------------------------------------- + * Helper: angular radius of ground visibility footprint + * + * For a satellite at altitude h km, the half-angle of the visibility + * cone (Earth-center angle) at min_el elevation is: + * lambda = arccos(Re / (Re + h) * cos(el)) - el + * Returns degrees. + * ---------------------------------------------------------------- + */ +static inline double +ground_footprint_deg(double sma_km, double min_el_deg) +{ + double h_km = sma_km - WGS72_AE; + double el_rad, cos_ratio, lambda; + + if (h_km <= 0.0) + return 0.0; + + el_rad = min_el_deg * DEG_TO_RAD; + cos_ratio = WGS72_AE / (WGS72_AE + h_km) * cos(el_rad); + + if (cos_ratio >= 1.0) + return 0.0; + + lambda = acos(cos_ratio) - el_rad; + + return lambda * RAD_TO_DEG; +} + + +/* ---------------------------------------------------------------- + * Helper: J2 secular RAAN precession rate in radians/day + * + * dOmega/dt = -1.5 * n * J2 * (Re/a)^2 * cos(i) + * + * where n is mean motion in rad/s, J2 and Re are WGS-72. + * Result in rad/day (multiply rad/s by 86400). + * ---------------------------------------------------------------- + */ +static inline double +j2_raan_rate(double sma_km, double inc_rad) +{ + double a = sma_km; + double ratio, n_rad_s; + + if (a <= 0.0) + return 0.0; + + ratio = WGS72_AE / a; + n_rad_s = sqrt(WGS72_MU / (a * a * a)); + + return -1.5 * n_rad_s * WGS72_J2 * ratio * ratio * cos(inc_rad) + * SECONDS_PER_DAY; +} + + +/* ---------------------------------------------------------------- + * Traversal state carried down the tree during index scans. + * Accumulates SMA and inclination ranges from L0 and L1. + * ---------------------------------------------------------------- + */ +typedef struct OrbitalTraversal +{ + double sma_low; + double sma_high; + double inc_low; + double inc_high; +} OrbitalTraversal; + + +/* ---------------------------------------------------------------- + * Query parameter extraction from observer_window composite type + * + * observer_window is: (obs observer, t_start timestamptz, + * t_end timestamptz, min_el float8) + * ---------------------------------------------------------------- + */ +typedef struct ObserverWindow +{ + pg_observer obs; + double jd_start; + double jd_end; + double jd_mid; + double min_el_deg; +} ObserverWindow; + + +static void +extract_observer_window(HeapTupleHeader composite, ObserverWindow *win) +{ + bool isnull; + Datum val; + int64 ts_start, ts_end; + + /* + * Field 1: obs (observer -- fixed-size 24B, pass-by-reference). + * GetAttributeByNum returns a Datum that is a direct pointer to + * the pg_observer bytes in the composite tuple. No varlena header. + */ + val = GetAttributeByNum(composite, 1, &isnull); + if (isnull) + ereport(ERROR, + (errcode(ERRCODE_NULL_VALUE_NOT_ALLOWED), + errmsg("observer_window.obs must not be NULL"))); + + memcpy(&win->obs, DatumGetPointer(val), sizeof(pg_observer)); + + /* Field 2: t_start (timestamptz) */ + val = GetAttributeByNum(composite, 2, &isnull); + if (isnull) + ereport(ERROR, + (errcode(ERRCODE_NULL_VALUE_NOT_ALLOWED), + errmsg("observer_window.t_start must not be NULL"))); + ts_start = DatumGetTimestampTz(val); + win->jd_start = timestamptz_to_jd(ts_start); + + /* Field 3: t_end (timestamptz) */ + val = GetAttributeByNum(composite, 3, &isnull); + if (isnull) + ereport(ERROR, + (errcode(ERRCODE_NULL_VALUE_NOT_ALLOWED), + errmsg("observer_window.t_end must not be NULL"))); + ts_end = DatumGetTimestampTz(val); + win->jd_end = timestamptz_to_jd(ts_end); + + /* Validate time window ordering */ + if (win->jd_end < win->jd_start) + ereport(ERROR, + (errcode(ERRCODE_INVALID_PARAMETER_VALUE), + errmsg("observer_window.t_end must not precede t_start"))); + + win->jd_mid = (win->jd_start + win->jd_end) / 2.0; + + /* Field 4: min_el (float8), clamped to [0, 90] */ + val = GetAttributeByNum(composite, 4, &isnull); + if (isnull) + win->min_el_deg = 10.0; /* default */ + else + { + win->min_el_deg = DatumGetFloat8(val); + if (win->min_el_deg < 0.0) + win->min_el_deg = 0.0; + if (win->min_el_deg > 90.0) + win->min_el_deg = 90.0; + } +} + + +/* ---------------------------------------------------------------- + * Comparison function for picksplit qsort. + * Sorts by extracted key value (SMA or inclination). + * ---------------------------------------------------------------- + */ +typedef struct +{ + int orig_index; + double key; +} PicksplitEntry; + +static int +picksplit_entry_cmp(const void *a, const void *b) +{ + double ka = ((const PicksplitEntry *) a)->key; + double kb = ((const PicksplitEntry *) b)->key; + + if (ka < kb) + return -1; + if (ka > kb) + return 1; + return 0; +} + + +/* ================================================================ + * SP-GiST support functions (5 required) + * ================================================================ + */ + + +/* + * spgist_tle_config -- declare trie type system + * + * No prefix data (VOIDOID): bin ranges encoded entirely in node labels. + * Labels are float8 (bin boundary values, sorted ascending). + * Leaves store full TLE unchanged. + */ +Datum +spgist_tle_config(PG_FUNCTION_ARGS) +{ + spgConfigIn *in = (spgConfigIn *) PG_GETARG_POINTER(0); + spgConfigOut *out = (spgConfigOut *) PG_GETARG_POINTER(1); + + out->prefixType = VOIDOID; + out->labelType = FLOAT8OID; + out->leafType = in->attType; /* tle type */ + out->canReturnData = true; + out->longValuesOK = false; + + PG_RETURN_VOID(); +} + + +/* + * spgist_tle_choose -- route a new TLE to the correct child node + * + * L0: route by SMA. L1: route by inclination. + * restDatum = leafDatum (full TLE unchanged), matching the quad-tree + * precedent where the tree terminates by depth, not value exhaustion. + * + * At level >= 2, all nodes are allTheSame (no further partitioning). + */ +Datum +spgist_tle_choose(PG_FUNCTION_ARGS) +{ + spgChooseIn *in = (spgChooseIn *) PG_GETARG_POINTER(0); + spgChooseOut *out = (spgChooseOut *) PG_GETARG_POINTER(1); + + pg_tle *tle = (pg_tle *) DatumGetPointer(in->leafDatum); + int level = in->level; + double val; + + /* Extract the routing key for this level */ + if (level == 0) + val = tle_sma_km(tle); + else + val = tle->inclination; + + if (in->allTheSame || level >= SPGIST_TLE_MAX_LEVEL) + { + out->resultType = spgMatchNode; + out->result.matchNode.nodeN = 0; + out->result.matchNode.levelAdd = 1; + out->result.matchNode.restDatum = in->leafDatum; + PG_RETURN_VOID(); + } + + /* + * Find the bin: labels are sorted ascending, each label marks the + * lower bound of a bin. We want the last bin whose label <= val. + */ + { + int best = 0; + int i; + + for (i = 1; i < in->nNodes; i++) + { + if (DatumGetFloat8(in->nodeLabels[i]) <= val) + best = i; + else + break; + } + + out->resultType = spgMatchNode; + out->result.matchNode.nodeN = best; + out->result.matchNode.levelAdd = 1; + out->result.matchNode.restDatum = in->leafDatum; + } + + PG_RETURN_VOID(); +} + + +/* + * spgist_tle_picksplit -- split a leaf page using equal-population strategy + * + * Sort by the current level's element (SMA at L0, inclination at L1), + * divide into floor(sqrt(nTuples)) bins clamped to [2, 16]. + * At level >= 2, emit a single allTheSame node. + */ +Datum +spgist_tle_picksplit(PG_FUNCTION_ARGS) +{ + spgPickSplitIn *in = (spgPickSplitIn *) PG_GETARG_POINTER(0); + spgPickSplitOut *out = (spgPickSplitOut *) PG_GETARG_POINTER(1); + + int level = in->level; + int nTuples = in->nTuples; + int nBins, perBin, remainder; + int i, bin, pos; + PicksplitEntry *entries; + + /* + * At level >= 2 we have no further partitioning dimension. + * Emit a single allTheSame node that accepts everything. + */ + if (level >= SPGIST_TLE_MAX_LEVEL) + { + out->nNodes = 1; + out->hasPrefix = false; + out->prefixDatum = (Datum) 0; + out->nodeLabels = (Datum *) palloc(sizeof(Datum)); + out->nodeLabels[0] = Float8GetDatum(0.0); + out->mapTuplesToNodes = (int *) palloc(sizeof(int) * nTuples); + out->leafTupleDatums = (Datum *) palloc(sizeof(Datum) * nTuples); + + for (i = 0; i < nTuples; i++) + { + out->mapTuplesToNodes[i] = 0; + out->leafTupleDatums[i] = in->datums[i]; + } + + PG_RETURN_VOID(); + } + + /* Extract and sort by the current level's element */ + entries = (PicksplitEntry *) palloc(sizeof(PicksplitEntry) * nTuples); + + for (i = 0; i < nTuples; i++) + { + pg_tle *tle = (pg_tle *) DatumGetPointer(in->datums[i]); + + entries[i].orig_index = i; + if (level == 0) + entries[i].key = tle_sma_km(tle); + else + entries[i].key = tle->inclination; + } + + qsort(entries, nTuples, sizeof(PicksplitEntry), picksplit_entry_cmp); + + /* Equal-population split: floor(sqrt(n)) bins, clamped [2, 16] */ + nBins = (int) floor(sqrt((double) nTuples)); + if (nBins < 2) + nBins = 2; + if (nBins > 16) + nBins = 16; + + perBin = nTuples / nBins; + remainder = nTuples % nBins; + + out->nNodes = nBins; + out->hasPrefix = false; + out->prefixDatum = (Datum) 0; + out->nodeLabels = (Datum *) palloc(sizeof(Datum) * nBins); + out->mapTuplesToNodes = (int *) palloc(sizeof(int) * nTuples); + out->leafTupleDatums = (Datum *) palloc(sizeof(Datum) * nTuples); + + pos = 0; + for (bin = 0; bin < nBins; bin++) + { + int count = perBin + (bin < remainder ? 1 : 0); + + /* Node label = key value of the first entry in this bin */ + out->nodeLabels[bin] = Float8GetDatum(entries[pos].key); + + for (i = 0; i < count; i++) + { + int orig = entries[pos + i].orig_index; + + out->mapTuplesToNodes[orig] = bin; + out->leafTupleDatums[orig] = in->datums[orig]; + } + + pos += count; + } + + pfree(entries); + + PG_RETURN_VOID(); +} + + +/* + * spgist_tle_inner_consistent -- prune child nodes during index scan + * + * L0: skip bins whose perigee altitude exceeds max visible altitude. + * The bin's SMA range is [label[i], label[i+1]). + * L1: skip bins whose inclination is too low to reach observer latitude. + * A satellite with inclination i has ground track bounded by [-i, +i]. + * Observer at latitude phi needs i + footprint >= |phi|. + * + * Propagates OrbitalTraversal state to surviving children via + * traversalMemoryContext for the RAAN filter at leaf level. + */ +Datum +spgist_tle_inner_consistent(PG_FUNCTION_ARGS) +{ + spgInnerConsistentIn *in = (spgInnerConsistentIn *) PG_GETARG_POINTER(0); + spgInnerConsistentOut *out = (spgInnerConsistentOut *) PG_GETARG_POINTER(1); + + int level = in->level; + int nkeys = in->nkeys; + int i; + ObserverWindow win; + bool have_query = false; + + /* Extract query from scankeys -- we need the &? operator's arg */ + for (i = 0; i < nkeys; i++) + { + if (in->scankeys[i].sk_strategy == 1) + { + HeapTupleHeader composite; + + composite = DatumGetHeapTupleHeader(in->scankeys[i].sk_argument); + extract_observer_window(composite, &win); + have_query = true; + break; + } + } + + /* Allocate output arrays */ + out->nodeNumbers = (int *) palloc(sizeof(int) * in->nNodes); + out->levelAdds = (int *) palloc(sizeof(int) * in->nNodes); + out->reconstructedValues = NULL; + out->traversalValues = (void **) palloc(sizeof(void *) * in->nNodes); + out->nNodes = 0; + + for (i = 0; i < in->nNodes; i++) + { + OrbitalTraversal *parent_trav; + OrbitalTraversal *child_trav; + double bin_low, bin_high; + bool dominated = false; + + /* Decode bin range from labels */ + bin_low = DatumGetFloat8(in->nodeLabels[i]); + if (i < in->nNodes - 1) + bin_high = DatumGetFloat8(in->nodeLabels[i + 1]); + else + bin_high = INFINITY; + + /* Inherit parent traversal state or initialize */ + if (in->traversalValue) + parent_trav = (OrbitalTraversal *) in->traversalValue; + else + parent_trav = NULL; + + /* Pruning logic per level */ + if (have_query && level == 0) + { + /* + * L0: SMA pruning. + * bin_low is the SMA of the lowest object in this bin. + * A satellite in this bin has perigee_alt >= bin_low - AE + * (approximately, for near-circular orbits in the bin). + * If the lowest possible altitude exceeds max_visible, skip. + * + * Conservative: use bin_low as a lower bound on SMA, + * compute perigee assuming e=0 (circular, worst case for + * "too high" pruning -- circular has highest perigee for + * a given SMA). + */ + double perigee_alt = bin_low - WGS72_AE; + double max_alt = max_visible_altitude_km(win.min_el_deg); + + if (perigee_alt > max_alt) + dominated = true; + } + else if (have_query && level == 1) + { + /* + * L1: Inclination pruning. + * bin_high is the upper bound on inclination in this bin. + * A satellite with inclination i has ground track [-i, +i]. + * The observer at latitude phi can see it if: + * i + footprint >= |phi| + * + * Use the parent SMA range to compute a conservative footprint. + * The largest footprint comes from the lowest altitude. + */ + double obs_lat = fabs(win.obs.lat); + double sma_for_footprint; + double footprint; + + if (parent_trav) + sma_for_footprint = parent_trav->sma_low; + else + sma_for_footprint = WGS72_AE + 200.0; /* conservative LEO */ + + footprint = ground_footprint_deg(sma_for_footprint, + win.min_el_deg) * DEG_TO_RAD; + + if (bin_high + footprint < obs_lat) + dominated = true; + } + + if (!dominated) + { + int idx = out->nNodes; + + /* Build child traversal state */ + child_trav = (OrbitalTraversal *) + MemoryContextAlloc(in->traversalMemoryContext, + sizeof(OrbitalTraversal)); + + if (parent_trav) + memcpy(child_trav, parent_trav, sizeof(OrbitalTraversal)); + else + { + child_trav->sma_low = 0.0; + child_trav->sma_high = INFINITY; + child_trav->inc_low = 0.0; + child_trav->inc_high = M_PI; + } + + /* Narrow bounds based on current level */ + if (level == 0) + { + child_trav->sma_low = bin_low; + child_trav->sma_high = bin_high; + } + else if (level == 1) + { + child_trav->inc_low = bin_low; + child_trav->inc_high = bin_high; + } + + out->nodeNumbers[idx] = i; + out->levelAdds[idx] = 1; + out->traversalValues[idx] = child_trav; + out->nNodes++; + } + } + + PG_RETURN_VOID(); +} + + +/* + * spgist_tle_leaf_consistent -- final check on a leaf TLE + * + * Applies three filters: + * 1. Perigee altitude check (with eccentricity) + * 2. Inclination + ground footprint vs observer latitude + * 3. RAAN query-time filter (J2 precession to query midpoint) + * + * recheck = false: the &? operator IS the superset filter. + * The user runs predict_passes() on survivors for SGP4 ground truth. + */ +Datum +spgist_tle_leaf_consistent(PG_FUNCTION_ARGS) +{ + spgLeafConsistentIn *in = (spgLeafConsistentIn *) PG_GETARG_POINTER(0); + spgLeafConsistentOut *out = (spgLeafConsistentOut *) PG_GETARG_POINTER(1); + + pg_tle *tle; + int i; + bool result = true; + + tle = (pg_tle *) DatumGetPointer(in->leafDatum); + out->leafValue = in->leafDatum; + out->recheck = false; + + for (i = 0; i < in->nkeys; i++) + { + if (in->scankeys[i].sk_strategy == 1) + { + ObserverWindow win; + HeapTupleHeader composite; + double sma, perigee_alt, max_alt; + double obs_lat_abs, footprint_rad; + double dt_days, raan_projected, lst; + double earth_rot_rad, raan_window_half, raan_diff; + + composite = DatumGetHeapTupleHeader(in->scankeys[i].sk_argument); + extract_observer_window(composite, &win); + + sma = tle_sma_km(tle); + + /* Filter 1: perigee altitude check */ + perigee_alt = sma * (1.0 - tle->eccentricity) - WGS72_AE; + max_alt = max_visible_altitude_km(win.min_el_deg); + if (perigee_alt > max_alt) + { + result = false; + break; + } + + /* Filter 2: inclination + footprint vs observer latitude */ + obs_lat_abs = fabs(win.obs.lat); + footprint_rad = ground_footprint_deg(sma, win.min_el_deg) + * DEG_TO_RAD; + + if (tle->inclination + footprint_rad < obs_lat_abs) + { + result = false; + break; + } + + /* + * Filter 3: RAAN query-time filter. + * + * Project RAAN to query midpoint via J2 precession rate. + * Check alignment with observer's local sidereal time. + * Earth rotates 360 deg/day, so the RAAN window widens + * with query duration. + */ + dt_days = win.jd_mid - tle->epoch; + raan_projected = tle->raan + + j2_raan_rate(sma, tle->inclination) * dt_days; + raan_projected = fmod(raan_projected, 2.0 * M_PI); + if (raan_projected < 0.0) + raan_projected += 2.0 * M_PI; + + /* Observer LST at query midpoint */ + lst = gmst_from_jd(win.jd_mid) + win.obs.lon; + lst = fmod(lst, 2.0 * M_PI); + if (lst < 0.0) + lst += 2.0 * M_PI; + + /* + * RAAN window: Earth rotation during query + footprint. + * Full rotation (24h+) means pass everything. + */ + earth_rot_rad = (win.jd_end - win.jd_start) * EARTH_ROT_RAD_PER_DAY; + raan_window_half = earth_rot_rad / 2.0 + + ground_footprint_deg(sma, win.min_el_deg) * DEG_TO_RAD; + + if (raan_window_half >= M_PI) + { + /* Full rotation or close: no RAAN filter for this key */ + continue; + } + + raan_diff = fabs(raan_projected - lst); + if (raan_diff > M_PI) + raan_diff = 2.0 * M_PI - raan_diff; + + if (raan_diff > raan_window_half) + { + result = false; + break; + } + } + } + + PG_RETURN_BOOL(result); +} + + +/* ================================================================ + * Operator function: &? (visibility cone check) + * ================================================================ + */ + + +/* + * tle_visibility_possible(observer_window, tle) -> bool + * + * Standalone operator: can the satellite possibly be visible from + * this observer during this time window? Combines altitude check, + * latitude/inclination check, and RAAN filter. + * + * This is the same logic as leaf_consistent, callable directly + * as a SQL operator for sequential scans or WHERE clauses. + */ +Datum +tle_visibility_possible(PG_FUNCTION_ARGS) +{ + HeapTupleHeader composite = PG_GETARG_HEAPTUPLEHEADER(0); + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(1); + ObserverWindow win; + double sma, perigee_alt, max_alt; + double obs_lat_abs, footprint_rad; + double dt_days, raan_projected, lst; + double earth_rot_rad, raan_window_half, raan_diff; + + extract_observer_window(composite, &win); + + sma = tle_sma_km(tle); + + /* Filter 1: perigee altitude */ + perigee_alt = sma * (1.0 - tle->eccentricity) - WGS72_AE; + max_alt = max_visible_altitude_km(win.min_el_deg); + if (perigee_alt > max_alt) + PG_RETURN_BOOL(false); + + /* Filter 2: inclination + footprint */ + obs_lat_abs = fabs(win.obs.lat); + footprint_rad = ground_footprint_deg(sma, win.min_el_deg) * DEG_TO_RAD; + + if (tle->inclination + footprint_rad < obs_lat_abs) + PG_RETURN_BOOL(false); + + /* Filter 3: RAAN alignment */ + dt_days = win.jd_mid - tle->epoch; + raan_projected = tle->raan + + j2_raan_rate(sma, tle->inclination) * dt_days; + raan_projected = fmod(raan_projected, 2.0 * M_PI); + if (raan_projected < 0.0) + raan_projected += 2.0 * M_PI; + + lst = gmst_from_jd(win.jd_mid) + win.obs.lon; + lst = fmod(lst, 2.0 * M_PI); + if (lst < 0.0) + lst += 2.0 * M_PI; + + earth_rot_rad = (win.jd_end - win.jd_start) * EARTH_ROT_RAD_PER_DAY; + raan_window_half = earth_rot_rad / 2.0 + + ground_footprint_deg(sma, win.min_el_deg) * DEG_TO_RAD; + + if (raan_window_half >= M_PI) + PG_RETURN_BOOL(true); /* full rotation -- pass everything */ + + raan_diff = fabs(raan_projected - lst); + if (raan_diff > M_PI) + raan_diff = 2.0 * M_PI - raan_diff; + + PG_RETURN_BOOL(raan_diff <= raan_window_half); +} diff --git a/test/expected/spgist_tle.out b/test/expected/spgist_tle.out new file mode 100644 index 0000000..9e0e714 --- /dev/null +++ b/test/expected/spgist_tle.out @@ -0,0 +1,213 @@ +-- Test SP-GiST orbital trie index and &? visibility cone operator +SET client_min_messages = warning; +CREATE EXTENSION IF NOT EXISTS pg_orrery; +RESET client_min_messages; +-- ============================================================ +-- Test table with mixed orbital regimes +-- ============================================================ +CREATE TABLE test_spgist ( + id serial, + name text, + tle tle +); +-- ISS (LEO, ~400km, 51.64 deg) +INSERT INTO test_spgist (name, tle) VALUES ('ISS', + '1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025 +2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'); +-- Hubble (LEO, ~540km, 28.47 deg) +INSERT INTO test_spgist (name, tle) VALUES ('Hubble', + '1 20580U 90037B 24001.50000000 .00000790 00000+0 39573-4 0 9992 +2 20580 28.4705 61.4398 0002797 317.3115 42.7577 15.09395228 00008'); +-- GPS IIR-M (MEO, ~20200km, 55.44 deg) +INSERT INTO test_spgist (name, tle) VALUES ('GPS-IIR', + '1 28874U 05038A 24001.50000000 .00000012 00000+0 00000+0 0 9993 +2 28874 55.4408 300.3467 0117034 51.6543 309.5420 2.00557079 00006'); +-- Equatorial-LEO (same altitude as ISS, 5 deg inclination) +INSERT INTO test_spgist (name, tle) VALUES ('Equatorial-LEO', + '1 99901U 24999A 24001.50000000 .00016717 00000-0 10270-3 0 9990 +2 99901 5.0000 208.9163 0006703 30.1694 61.7520 15.50100486 00001'); +-- SSO-800 (Sun-synchronous, ~800km, 98.7 deg) +INSERT INTO test_spgist (name, tle) VALUES ('SSO-800', + '1 99902U 24999B 24001.50000000 .00000100 00000+0 50000-4 0 9991 +2 99902 98.7000 120.0000 0001000 90.0000 270.0000 14.19553000 00001'); +-- GEO-SAT (Geostationary, ~35786km, 0.04 deg) +INSERT INTO test_spgist (name, tle) VALUES ('GEO-SAT', + '1 99903U 24999C 24001.50000000 .00000000 00000+0 00000+0 0 9992 +2 99903 0.0400 270.0000 0003000 0.0000 180.0000 1.00273791 00001'); +-- ============================================================ +-- Test 1: Operator standalone — ISS from Eagle Idaho (2h window) +-- Eagle Idaho: 43.6977N 116.3535W, 760m elevation +-- ISS passes altitude and inclination checks, but RAAN filter +-- rejects it — the orbital plane isn't overhead during this +-- specific 2-hour window (correct physics, see Test 5 for 24h). +-- ============================================================ +SELECT name, + ROW( + observer('43.6977N 116.3535W 760m'), + '2024-01-01 00:00:00+00'::timestamptz, + '2024-01-01 02:00:00+00'::timestamptz, + 10.0 + )::observer_window &? tle AS visible +FROM test_spgist +WHERE name = 'ISS'; + name | visible +------+--------- + ISS | f +(1 row) + +-- ============================================================ +-- Test 2: Equatorial-LEO NOT visible from Eagle Idaho +-- 5 deg inc + ~12 deg footprint = 17 deg < 43.7 deg latitude +-- ============================================================ +SELECT name, + ROW( + observer('43.6977N 116.3535W 760m'), + '2024-01-01 00:00:00+00'::timestamptz, + '2024-01-01 02:00:00+00'::timestamptz, + 10.0 + )::observer_window &? tle AS visible +FROM test_spgist +WHERE name = 'Equatorial-LEO'; + name | visible +----------------+--------- + Equatorial-LEO | f +(1 row) + +-- ============================================================ +-- Test 3: Create SP-GiST index, verify index scan with positive +-- results. Equatorial observer at 0E — SSO-800 RAAN (120 deg) +-- aligns with LST near 0E at this epoch, so it passes. +-- ============================================================ +CREATE INDEX test_spgist_idx ON test_spgist USING spgist (tle tle_spgist_ops); +SET enable_seqscan = off; +SELECT name +FROM test_spgist +WHERE ROW( + observer('0.0N 0.0E 0m'), + '2024-01-01 00:00:00+00'::timestamptz, + '2024-01-01 02:00:00+00'::timestamptz, + 10.0 +)::observer_window &? tle +ORDER BY name; + name +--------- + SSO-800 +(1 row) + +RESET enable_seqscan; +-- ============================================================ +-- Test 4: Seqscan vs index scan consistency — same query must +-- return identical results regardless of scan method. +-- ============================================================ +SET enable_indexscan = off; +SET enable_bitmapscan = off; +SELECT name +FROM test_spgist +WHERE ROW( + observer('0.0N 0.0E 0m'), + '2024-01-01 00:00:00+00'::timestamptz, + '2024-01-01 02:00:00+00'::timestamptz, + 10.0 +)::observer_window &? tle +ORDER BY name; + name +--------- + SSO-800 +(1 row) + +RESET enable_indexscan; +RESET enable_bitmapscan; +-- ============================================================ +-- Test 5: 24-hour window — RAAN filter bypassed (full Earth +-- rotation). Only ISS and SSO-800 pass inclination from Eagle +-- Idaho (43.7 deg). Hubble (28.5+14.8=43.3 deg) barely fails. +-- GPS-IIR and GEO-SAT filtered by altitude. +-- ============================================================ +SELECT name +FROM test_spgist +WHERE ROW( + observer('43.6977N 116.3535W 760m'), + '2024-01-01 00:00:00+00'::timestamptz, + '2024-01-02 00:00:00+00'::timestamptz, + 10.0 +)::observer_window &? tle +ORDER BY name; + name +--------- + ISS + SSO-800 +(2 rows) + +-- ============================================================ +-- Test 6: High min_el (45 deg) changes footprint — wider +-- footprint lets more inclinations through. Same 24h window. +-- ============================================================ +SELECT name +FROM test_spgist +WHERE ROW( + observer('43.6977N 116.3535W 760m'), + '2024-01-01 00:00:00+00'::timestamptz, + '2024-01-02 00:00:00+00'::timestamptz, + 45.0 +)::observer_window &? tle +ORDER BY name; + name +--------- + ISS + SSO-800 +(2 rows) + +-- ============================================================ +-- Test 7: GiST coexistence — both index types on same table +-- ============================================================ +CREATE INDEX test_gist_idx ON test_spgist USING gist (tle); +-- GiST overlap query still works +SELECT a.name AS sat_a, b.name AS sat_b, a.tle && b.tle AS overlaps +FROM test_spgist a, test_spgist b +WHERE a.name = 'ISS' AND b.name = 'Hubble'; + sat_a | sat_b | overlaps +-------+--------+---------- + ISS | Hubble | f +(1 row) + +-- SP-GiST query still works alongside GiST +SET enable_seqscan = off; +SELECT name +FROM test_spgist +WHERE ROW( + observer('43.6977N 116.3535W 760m'), + '2024-01-01 00:00:00+00'::timestamptz, + '2024-01-02 00:00:00+00'::timestamptz, + 10.0 +)::observer_window &? tle +ORDER BY name; + name +--------- + ISS + SSO-800 +(2 rows) + +RESET enable_seqscan; +-- ============================================================ +-- Test 8: NULL TLE handling — NULLs should be excluded +-- ============================================================ +INSERT INTO test_spgist (name, tle) VALUES ('NULL-SAT', NULL); +SELECT name +FROM test_spgist +WHERE ROW( + observer('43.6977N 116.3535W 760m'), + '2024-01-01 00:00:00+00'::timestamptz, + '2024-01-02 00:00:00+00'::timestamptz, + 10.0 +)::observer_window &? tle +ORDER BY name; + name +--------- + ISS + SSO-800 +(2 rows) + +-- ============================================================ +-- Cleanup +-- ============================================================ +DROP TABLE test_spgist; diff --git a/test/sql/spgist_tle.sql b/test/sql/spgist_tle.sql new file mode 100644 index 0000000..bb6f7c9 --- /dev/null +++ b/test/sql/spgist_tle.sql @@ -0,0 +1,197 @@ +-- Test SP-GiST orbital trie index and &? visibility cone operator +SET client_min_messages = warning; +CREATE EXTENSION IF NOT EXISTS pg_orrery; +RESET client_min_messages; + +-- ============================================================ +-- Test table with mixed orbital regimes +-- ============================================================ +CREATE TABLE test_spgist ( + id serial, + name text, + tle tle +); + +-- ISS (LEO, ~400km, 51.64 deg) +INSERT INTO test_spgist (name, tle) VALUES ('ISS', + '1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025 +2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'); + +-- Hubble (LEO, ~540km, 28.47 deg) +INSERT INTO test_spgist (name, tle) VALUES ('Hubble', + '1 20580U 90037B 24001.50000000 .00000790 00000+0 39573-4 0 9992 +2 20580 28.4705 61.4398 0002797 317.3115 42.7577 15.09395228 00008'); + +-- GPS IIR-M (MEO, ~20200km, 55.44 deg) +INSERT INTO test_spgist (name, tle) VALUES ('GPS-IIR', + '1 28874U 05038A 24001.50000000 .00000012 00000+0 00000+0 0 9993 +2 28874 55.4408 300.3467 0117034 51.6543 309.5420 2.00557079 00006'); + +-- Equatorial-LEO (same altitude as ISS, 5 deg inclination) +INSERT INTO test_spgist (name, tle) VALUES ('Equatorial-LEO', + '1 99901U 24999A 24001.50000000 .00016717 00000-0 10270-3 0 9990 +2 99901 5.0000 208.9163 0006703 30.1694 61.7520 15.50100486 00001'); + +-- SSO-800 (Sun-synchronous, ~800km, 98.7 deg) +INSERT INTO test_spgist (name, tle) VALUES ('SSO-800', + '1 99902U 24999B 24001.50000000 .00000100 00000+0 50000-4 0 9991 +2 99902 98.7000 120.0000 0001000 90.0000 270.0000 14.19553000 00001'); + +-- GEO-SAT (Geostationary, ~35786km, 0.04 deg) +INSERT INTO test_spgist (name, tle) VALUES ('GEO-SAT', + '1 99903U 24999C 24001.50000000 .00000000 00000+0 00000+0 0 9992 +2 99903 0.0400 270.0000 0003000 0.0000 180.0000 1.00273791 00001'); + + +-- ============================================================ +-- Test 1: Operator standalone — ISS from Eagle Idaho (2h window) +-- Eagle Idaho: 43.6977N 116.3535W, 760m elevation +-- ISS passes altitude and inclination checks, but RAAN filter +-- rejects it — the orbital plane isn't overhead during this +-- specific 2-hour window (correct physics, see Test 5 for 24h). +-- ============================================================ +SELECT name, + ROW( + observer('43.6977N 116.3535W 760m'), + '2024-01-01 00:00:00+00'::timestamptz, + '2024-01-01 02:00:00+00'::timestamptz, + 10.0 + )::observer_window &? tle AS visible +FROM test_spgist +WHERE name = 'ISS'; + + +-- ============================================================ +-- Test 2: Equatorial-LEO NOT visible from Eagle Idaho +-- 5 deg inc + ~12 deg footprint = 17 deg < 43.7 deg latitude +-- ============================================================ +SELECT name, + ROW( + observer('43.6977N 116.3535W 760m'), + '2024-01-01 00:00:00+00'::timestamptz, + '2024-01-01 02:00:00+00'::timestamptz, + 10.0 + )::observer_window &? tle AS visible +FROM test_spgist +WHERE name = 'Equatorial-LEO'; + + +-- ============================================================ +-- Test 3: Create SP-GiST index, verify index scan with positive +-- results. Equatorial observer at 0E — SSO-800 RAAN (120 deg) +-- aligns with LST near 0E at this epoch, so it passes. +-- ============================================================ +CREATE INDEX test_spgist_idx ON test_spgist USING spgist (tle tle_spgist_ops); + +SET enable_seqscan = off; + +SELECT name +FROM test_spgist +WHERE ROW( + observer('0.0N 0.0E 0m'), + '2024-01-01 00:00:00+00'::timestamptz, + '2024-01-01 02:00:00+00'::timestamptz, + 10.0 +)::observer_window &? tle +ORDER BY name; + +RESET enable_seqscan; + + +-- ============================================================ +-- Test 4: Seqscan vs index scan consistency — same query must +-- return identical results regardless of scan method. +-- ============================================================ +SET enable_indexscan = off; +SET enable_bitmapscan = off; + +SELECT name +FROM test_spgist +WHERE ROW( + observer('0.0N 0.0E 0m'), + '2024-01-01 00:00:00+00'::timestamptz, + '2024-01-01 02:00:00+00'::timestamptz, + 10.0 +)::observer_window &? tle +ORDER BY name; + +RESET enable_indexscan; +RESET enable_bitmapscan; + + +-- ============================================================ +-- Test 5: 24-hour window — RAAN filter bypassed (full Earth +-- rotation). Only ISS and SSO-800 pass inclination from Eagle +-- Idaho (43.7 deg). Hubble (28.5+14.8=43.3 deg) barely fails. +-- GPS-IIR and GEO-SAT filtered by altitude. +-- ============================================================ +SELECT name +FROM test_spgist +WHERE ROW( + observer('43.6977N 116.3535W 760m'), + '2024-01-01 00:00:00+00'::timestamptz, + '2024-01-02 00:00:00+00'::timestamptz, + 10.0 +)::observer_window &? tle +ORDER BY name; + + +-- ============================================================ +-- Test 6: High min_el (45 deg) changes footprint — wider +-- footprint lets more inclinations through. Same 24h window. +-- ============================================================ +SELECT name +FROM test_spgist +WHERE ROW( + observer('43.6977N 116.3535W 760m'), + '2024-01-01 00:00:00+00'::timestamptz, + '2024-01-02 00:00:00+00'::timestamptz, + 45.0 +)::observer_window &? tle +ORDER BY name; + + +-- ============================================================ +-- Test 7: GiST coexistence — both index types on same table +-- ============================================================ +CREATE INDEX test_gist_idx ON test_spgist USING gist (tle); + +-- GiST overlap query still works +SELECT a.name AS sat_a, b.name AS sat_b, a.tle && b.tle AS overlaps +FROM test_spgist a, test_spgist b +WHERE a.name = 'ISS' AND b.name = 'Hubble'; + +-- SP-GiST query still works alongside GiST +SET enable_seqscan = off; +SELECT name +FROM test_spgist +WHERE ROW( + observer('43.6977N 116.3535W 760m'), + '2024-01-01 00:00:00+00'::timestamptz, + '2024-01-02 00:00:00+00'::timestamptz, + 10.0 +)::observer_window &? tle +ORDER BY name; +RESET enable_seqscan; + + +-- ============================================================ +-- Test 8: NULL TLE handling — NULLs should be excluded +-- ============================================================ +INSERT INTO test_spgist (name, tle) VALUES ('NULL-SAT', NULL); + +SELECT name +FROM test_spgist +WHERE ROW( + observer('43.6977N 116.3535W 760m'), + '2024-01-01 00:00:00+00'::timestamptz, + '2024-01-02 00:00:00+00'::timestamptz, + 10.0 +)::observer_window &? tle +ORDER BY name; + + +-- ============================================================ +-- Cleanup +-- ============================================================ +DROP TABLE test_spgist;