Ryan Malloy a792e7e083 Extend GiST index to 2-D: altitude + inclination
The 1-D altitude-band index only pruned ~25% of the 22k satellite
catalog (eliminates MEO/GEO/HEO but 75% is LEO).  Adding inclination
as a second indexed dimension prunes an additional ~40% of remaining
candidates — objects in equatorial or low-inclination orbits that
geometrically cannot pass over the observer's latitude.

Key changes:
- tle_alt_range (16 bytes) → tle_orbital_key (32 bytes) with
  inc_low/inc_high fields
- All 8 GiST support functions updated for 2-D bounding boxes
- Penalty uses margin (half-perimeter) not area to avoid degeneracy
  when leaf entries have zero-width inclination ranges
- Picksplit selects split dimension by normalized spread
- && operator now checks altitude AND inclination overlap
- <-> operator remains altitude-only (conjunction screening is
  altitude-dominant)
- SQL operator comments updated for 2-D semantics
- Test adds Equatorial-LEO satellite at ISS altitude but 5° inclination
  to validate inclination-based pruning
2026-02-15 18:10:19 -07:00

pg_orbit

Orbital mechanics types and functions for PostgreSQL.

pg_orbit adds native SQL types for TLEs, orbital positions, ground stations, and satellite passes. It wraps Bill Gray's sat_code (MIT) for SGP4/SDP4 propagation, provides coordinate transforms between inertial and ground-referenced frames, predicts passes over observer locations, and supports GiST-indexed conjunction screening on altitude bands.

Think PostGIS, but for objects in orbit.

Installation

Requirements:

  • PostgreSQL 14+ development headers (pg_config in PATH)
  • GCC and Make
  • C++ compiler (for sat_code)
git clone --recurse-submodules https://github.com/...
cd pg_orbit
make
sudo make install

Then in your database:

CREATE EXTENSION pg_orbit;

If you cloned without --recurse-submodules, initialize the sat_code dependency:

git submodule update --init

Quick Start

-- Create a table with a TLE column
CREATE TABLE satellites (
    norad_id  int PRIMARY KEY,
    name      text NOT NULL,
    tle       tle NOT NULL
);

-- Insert a TLE (standard two-line format, concatenated with newline)
INSERT INTO satellites VALUES (25544, 'ISS',
    '1 25544U 98067A   24001.50000000  .00016717  00000-0  10270-3 0  9002
2 25544  51.6400 208.5000 0006000  30.0000 330.0000 15.50000000400000');

-- Propagate to a point in time
SELECT sgp4_propagate(tle, now()) FROM satellites WHERE norad_id = 25544;

-- Subsatellite point (nadir)
SELECT subsatellite_point(tle, now()) FROM satellites WHERE norad_id = 25544;

-- All passes over Boulder, CO in the next 24 hours (min 10 deg elevation)
SELECT sat.name, p.*
FROM satellites sat,
     LATERAL predict_passes(sat.tle, '40.0N 105.3W 1655m'::observer,
                            now(), now() + '24h', 10.0) p
ORDER BY pass_aos_time(p);

-- Conjunction screening with GiST index
CREATE INDEX ON satellites USING gist (tle);

SELECT a.name, b.name, tle_distance(a.tle, b.tle, now())
FROM satellites a, satellites b
WHERE a.tle && b.tle AND a.norad_id < b.norad_id
  AND tle_distance(a.tle, b.tle, now()) < 10.0;

Types

Type Size Description
tle 112 bytes Parsed mean orbital elements (epoch, Keplerian elements, drag terms, identifiers). Stored as a fixed-size struct, not raw text.
eci_position 48 bytes Position (km) and velocity (km/s) in the True Equator Mean Equinox (TEME) frame.
geodetic 24 bytes Latitude/longitude (degrees) and altitude (km) on the WGS-84 ellipsoid.
topocentric 32 bytes Azimuth, elevation (degrees), slant range (km), and range rate (km/s) relative to an observer.
observer 24 bytes Ground station location. Accepts human-readable input: '40.0N 105.3W 1655m' or decimal degrees.
pass_event 48 bytes Satellite pass with AOS/MAX/LOS times, max elevation, and AOS/LOS azimuths.

Input Formats

tle -- Standard two-line format (lines joined by newline):

SELECT '1 25544U 98067A   24001.50000000  .00016717  00000-0  10270-3 0  9002
2 25544  51.6400 208.5000 0006000  30.0000 330.0000 15.50000000400000'::tle;

observer -- Flexible ground station input:

-- Compass notation with altitude
SELECT '40.0N 105.3W 1655m'::observer;

-- Decimal degrees (positive East, altitude in meters)
SELECT '40.0 -105.3 1655'::observer;

Functions

TLE Accessors

Function Returns Description
tle_norad_id(tle) int4 NORAD catalog number
tle_epoch(tle) float8 Epoch as Julian date (UTC)
tle_inclination(tle) float8 Inclination in degrees
tle_eccentricity(tle) float8 Eccentricity (dimensionless)
tle_raan(tle) float8 Right ascension of ascending node (degrees)
tle_arg_perigee(tle) float8 Argument of perigee (degrees)
tle_mean_anomaly(tle) float8 Mean anomaly (degrees)
tle_mean_motion(tle) float8 Mean motion (rev/day)
tle_bstar(tle) float8 B* drag coefficient (1/earth-radii)
tle_period(tle) float8 Orbital period (minutes)
tle_perigee(tle) float8 Perigee altitude (km above WGS-72 ellipsoid)
tle_apogee(tle) float8 Apogee altitude (km above WGS-72 ellipsoid)
tle_age(tle, timestamptz) float8 TLE age in days (positive = past epoch)
tle_intl_desig(tle) text International designator (COSPAR ID)

Propagation

Function Returns Description
sgp4_propagate(tle, timestamptz) eci_position Propagate to a point in time. Uses SGP4 for near-earth, SDP4 for deep-space.
sgp4_propagate_series(tle, start, stop, step) SETOF (t, x, y, z, vx, vy, vz) Time series of TEME positions at regular intervals.
tle_distance(tle, tle, timestamptz) float8 Euclidean distance (km) between two objects at a reference time.

Coordinate Transforms

Function Returns Description
eci_to_geodetic(eci_position, timestamptz) geodetic TEME to WGS-84 geodetic (lat/lon/alt). Requires time for Earth rotation.
eci_to_topocentric(eci_position, observer, timestamptz) topocentric TEME to observer-relative az/el/range.
subsatellite_point(tle, timestamptz) geodetic Nadir point on WGS-84 ellipsoid. Propagates internally.
ground_track(tle, start, stop, step) SETOF (t, lat, lon, alt) Ground track as time series of subsatellite points.

ECI Accessors

Function Returns Description
eci_x(eci_position) float8 X position (km, TEME)
eci_y(eci_position) float8 Y position (km, TEME)
eci_z(eci_position) float8 Z position (km, TEME)
eci_vx(eci_position) float8 X velocity (km/s)
eci_vy(eci_position) float8 Y velocity (km/s)
eci_vz(eci_position) float8 Z velocity (km/s)
eci_speed(eci_position) float8 Velocity magnitude (km/s)
eci_altitude(eci_position) float8 Geocentric altitude (km)

Topocentric Accessors

Function Returns Description
topo_azimuth(topocentric) float8 Azimuth in degrees (0=N, 90=E, 180=S, 270=W)
topo_elevation(topocentric) float8 Elevation in degrees (0=horizon, 90=zenith)
topo_range(topocentric) float8 Slant range (km)
topo_range_rate(topocentric) float8 Range rate (km/s, positive = receding)

Pass Prediction

Function Returns Description
next_pass(tle, observer, timestamptz) pass_event Next pass over observer. Searches up to 7 days.
predict_passes(tle, observer, start, stop [, min_el]) SETOF pass_event All passes in a time window. Optional minimum elevation (degrees).
pass_visible(tle, observer, start, stop) boolean True if any pass occurs in the time window.

Pass Event Accessors

Function Returns Description
pass_aos_time(pass_event) timestamptz Acquisition of signal time
pass_max_el_time(pass_event) timestamptz Maximum elevation time
pass_los_time(pass_event) timestamptz Loss of signal time
pass_max_elevation(pass_event) float8 Maximum elevation (degrees)
pass_aos_azimuth(pass_event) float8 AOS azimuth (degrees, 0=N)
pass_los_azimuth(pass_event) float8 LOS azimuth (degrees, 0=N)
pass_duration(pass_event) interval Pass duration (LOS - AOS)

Operators

Operator Operands Description
&& tle, tle Altitude band overlap. Necessary (not sufficient) condition for conjunction.
<-> tle, tle Minimum altitude-band separation in km. Returns 0 if bands overlap.

Both operators are supported by the GiST tle_ops operator class for indexed scans.

GiST Indexing

The tle_ops operator class indexes TLEs by their altitude band (perigee to apogee). This provides fast filtering for conjunction screening: only pairs whose altitude bands overlap can possibly be close to each other.

-- Create the index
CREATE INDEX idx_tle_alt ON satellites USING gist (tle);

-- The && operator triggers index scans
EXPLAIN SELECT a.name, b.name
FROM satellites a, satellites b
WHERE a.tle && b.tle AND a.norad_id < b.norad_id;

-- KNN ordering by altitude-band distance
SELECT name, tle <-> (SELECT tle FROM satellites WHERE norad_id = 25544) AS sep
FROM satellites
ORDER BY tle <-> (SELECT tle FROM satellites WHERE norad_id = 25544)
LIMIT 20;

The index reduces conjunction candidate pairs from O(n^2) to the set of objects with intersecting altitude bands, which is then refined by computing actual tle_distance() at a specific time.

Geodetic Constants

TLEs are mean elements fitted using WGS-72 constants. Using WGS-84 constants for propagation introduces kilometer-scale position errors because the elements absorb geodetic model biases during the fitting process.

pg_orbit enforces this:

  • Propagation (SGP4/SDP4): WGS-72 constants only (mu, ae, J2, J3, J4, ke)
  • Coordinate output (geodetic, topocentric): WGS-84 (a=6378.137 km, f=1/298.257223563)

These are not interchangeable. Mixing them is a silent accuracy loss.

Building from Source

# Build (requires pg_config in PATH)
make

# Install to PostgreSQL extension directory
sudo make install

# Run regression tests against a live database
make installcheck

Override pg_config location if needed:

make PG_CONFIG=/usr/lib/postgresql/16/bin/pg_config
sudo make PG_CONFIG=/usr/lib/postgresql/16/bin/pg_config install

Project Layout

pg_orbit.control           Extension metadata (version 0.1.0)
Makefile                   PGXS build
sql/
  pg_orbit--0.1.0.sql      Type, function, operator, and GiST definitions
src/
  pg_orbit.c               PG_MODULE_MAGIC entry point
  tle_type.c               TLE input/output/binary/accessors
  eci_type.c               ECI position type
  observer_type.c          Observer type with flexible parsing
  sgp4_funcs.c             SGP4 propagation and distance
  coord_funcs.c            Coordinate transforms (TEME/geodetic/topocentric)
  pass_funcs.c             Pass prediction (next_pass, predict_passes)
  gist_tle.c               GiST operator class for altitude-band indexing
  types.h                  Shared struct definitions and constants
lib/
  sat_code/                Bill Gray's SGP4/SDP4 (MIT, git submodule)
test/
  sql/                     Regression test SQL
  expected/                Expected output
  data/
    vallado_518.csv        518 verification test vectors (Vallado et al.)

Testing

pg_orbit uses the standard PostgreSQL regression test framework.

make installcheck

Test categories:

Suite Coverage
tle_parse TLE input/output round-trip, malformed input rejection
sgp4_propagate Vallado 518 test vectors, deep-space and high-eccentricity edge cases
coord_transforms TEME to geodetic, TEME to topocentric accuracy
pass_prediction Known ISS passes, polar and retrograde orbits
gist_index Index scan vs. sequential scan result equivalence

The Vallado 518 test vectors are the standard SGP4 verification dataset. Each row contains a NORAD ID, minutes since epoch, and expected position/velocity. All 518 must pass to machine epsilon.

License

PostgreSQL License. Copyright (c) 2025, Ryan Malloy.

The bundled sat_code library is separately licensed under the MIT license.

Description
No description provided
https://pg-orrery.warehack.ing
Readme 49 MiB
Languages
C 97.8%
PLpgSQL 1.9%
Shell 0.1%