warehack.ing/pg_orrery
2
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add SP-GiST orbital trie index for satellite pass prediction (v0.7.0)

Browse Source 
2-level SP-GiST index on TLE data: SMA at L0, inclination at L1, with
query-time RAAN filter via J2 secular precession.  New &? operator
(observer_window &? tle) returns true when a satellite might be visible
from a ground observer during a time window.

Index prunes by altitude band, inclination+footprint vs observer
latitude, and RAAN alignment against local sidereal time.  Operator
class tle_spgist_ops is opt-in (not default), coexists with existing
GiST tle_ops.  Equal-population picksplit with sqrt(n) bins.
This commit is contained in:
Ryan Malloy 2026-02-17 20:27:54 -07:00
parent bb235f51fa
commit 6f9be428f9
7 changed files with 2294 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
Makefile
View File

@ -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.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.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.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 # Our extension C sources
OBJS = src/pg_orrery.o src/tle_type.o src/eci_type.o src/observer_type.o \ 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/moon_funcs.o src/radio_funcs.o \
src/lambert.o src/transfer_funcs.o \ src/lambert.o src/transfer_funcs.o \
src/de_reader.o src/eph_provider.o src/de_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) # Vendored SGP4/SDP4 sources (pure C, from Bill Gray's sat_code, MIT license)
SGP4_DIR = src/sgp4 SGP4_DIR = src/sgp4
@ -31,7 +33,7 @@ OBJS += $(SGP4_OBJS)
# Regression tests # Regression tests
REGRESS = tle_parse sgp4_propagate coord_transforms pass_prediction gist_index convenience \ REGRESS = tle_parse sgp4_propagate coord_transforms pass_prediction gist_index convenience \
star_observe kepler_comet planet_observe moon_observe lambert_transfer \ 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 REGRESS_OPTS = --inputdir=test
# Pure C — no C++ runtime needed. LAPACK for OD solver (dgelss_). # Pure C — no C++ runtime needed. LAPACK for OD solver (dgelss_).

2
pg_orrery.control
View File

@ -1,4 +1,4 @@
comment = 'A database orrery — celestial mechanics types and functions for PostgreSQL' 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' module_pathname = '$libdir/pg_orrery'
relocatable = true relocatable = true

76
sql/pg_orrery--0.6.0--0.7.0.sql Normal file
View File

@ -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);

961
sql/pg_orrery--0.7.0.sql Normal file
View File

@ -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);

841
src/spgist_tle.c Normal file
View File

@ -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 <math.h>
#include <float.h>
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);
}

213
test/expected/spgist_tle.out Normal file
View File

@ -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;

197
test/sql/spgist_tle.sql Normal file
View File

@ -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;