From 15a830dc40f141346435382deed089551cc50f40 Mon Sep 17 00:00:00 2001 From: Ryan Malloy Date: Sun, 15 Feb 2026 17:07:07 -0700 Subject: [PATCH] Initial implementation of pg_orbit PostgreSQL extension 6 custom types (tle, eci_position, geodetic, topocentric, observer, pass_event), 67 SQL functions, 2 operators (&&, <->), and a GiST operator class for altitude-band indexing. Wraps Bill Gray's sat_code for SGP4/SDP4 propagation with WGS-72 constants for propagation and WGS-84 for coordinate output. All 5 regression tests pass on PG 18. --- .gitignore | 19 + .gitmodules | 3 + CLAUDE.md | 147 ++++++ LICENSE | 19 + Makefile | 36 ++ README.md | 305 +++++++++++ docs/DESIGN.md | 590 +++++++++++++++++++++ lib/sat_code | 1 + pg_orbit.control | 4 + sql/pg_orbit--0.1.0.sql | 491 +++++++++++++++++ src/coord_funcs.c | 820 +++++++++++++++++++++++++++++ src/eci_type.c | 219 ++++++++ src/gist_tle.c | 505 ++++++++++++++++++ src/observer_type.c | 284 ++++++++++ src/pass_funcs.c | 773 +++++++++++++++++++++++++++ src/pg_orbit.c | 12 + src/sgp4_funcs.c | 381 ++++++++++++++ src/tle_type.c | 451 ++++++++++++++++ src/types.h | 190 +++++++ test/expected/coord_transforms.out | 110 ++++ test/expected/gist_index.out | 83 +++ test/expected/pass_prediction.out | 91 ++++ test/expected/sgp4_propagate.out | 92 ++++ test/expected/tle_parse.out | 116 ++++ test/sql/coord_transforms.sql | 77 +++ test/sql/gist_index.sql | 61 +++ test/sql/pass_prediction.sql | 68 +++ test/sql/sgp4_propagate.sql | 67 +++ test/sql/tle_parse.sql | 74 +++ 29 files changed, 6089 insertions(+) create mode 100644 .gitignore create mode 100644 .gitmodules create mode 100644 CLAUDE.md create mode 100644 LICENSE create mode 100644 Makefile create mode 100644 README.md create mode 100644 docs/DESIGN.md create mode 160000 lib/sat_code create mode 100644 pg_orbit.control create mode 100644 sql/pg_orbit--0.1.0.sql create mode 100644 src/coord_funcs.c create mode 100644 src/eci_type.c create mode 100644 src/gist_tle.c create mode 100644 src/observer_type.c create mode 100644 src/pass_funcs.c create mode 100644 src/pg_orbit.c create mode 100644 src/sgp4_funcs.c create mode 100644 src/tle_type.c create mode 100644 src/types.h create mode 100644 test/expected/coord_transforms.out create mode 100644 test/expected/gist_index.out create mode 100644 test/expected/pass_prediction.out create mode 100644 test/expected/sgp4_propagate.out create mode 100644 test/expected/tle_parse.out create mode 100644 test/sql/coord_transforms.sql create mode 100644 test/sql/gist_index.sql create mode 100644 test/sql/pass_prediction.sql create mode 100644 test/sql/sgp4_propagate.sql create mode 100644 test/sql/tle_parse.sql diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..5335aab --- /dev/null +++ b/.gitignore @@ -0,0 +1,19 @@ +# Build artifacts +*.o +*.bc +*.so +*.dylib + +# pg_regress output +regression.diffs +regression.out +results/ +tmp_check/ +log/ + +# Editor files +*.swp +*.swo +*~ +.vscode/ +.idea/ diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..4767964 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "lib/sat_code"] + path = lib/sat_code + url = https://github.com/Bill-Gray/sat_code.git diff --git a/CLAUDE.md b/CLAUDE.md new file mode 100644 index 0000000..626cd42 --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1,147 @@ +# pg_orbit — PostgreSQL Extension for Orbital Mechanics + +## What This Is +A PostgreSQL extension that makes TLE/orbital data first-class types — the way PostGIS does for geographic data. Native C extension using PGXS, wrapping Bill Gray's `sat_code` SGP4/SDP4 implementation. + +## Build System +```bash +make # Compile with PGXS +make install # Install to PostgreSQL extensions dir +make installcheck # Run regression tests +``` + +Requires: PostgreSQL 14+ development headers (`pg_config` in PATH), GCC, Make. + +## Project Layout +``` +pg_orbit.control # Extension metadata +Makefile # PGXS build +sql/ + pg_orbit--0.1.0.sql # All type/function/operator definitions +src/ + pg_orbit.c # PG_MODULE_MAGIC entry point + tle_type.c # TLE custom type (input/output/binary/accessors) + eci_type.c # ECI position type + geodetic/topocentric types + observer_type.c # Observer location type with flexible parsing + sgp4_funcs.c # sgp4_propagate(), sgp4_propagate_series() + coord_funcs.c # eci_to_geodetic(), eci_to_topocentric(), subsatellite_point() + pass_funcs.c # next_pass(), predict_passes(), pass_visible() + gist_tle.c # GiST operator class for altitude-band indexing + types.h # Shared struct definitions +lib/ + sat_code/ # Bill Gray's SGP4 (MIT license, git submodule) +test/ + sql/ # Regression test SQL + expected/ # Expected output + data/ + vallado_518.csv # 518 verification test vectors +docs/ + DESIGN.md # Architecture decisions, theory-to-code mappings +``` + +## Type System + +### Core Types (all varlena or fixed-size, stored in tuples) + +| Type | Storage | Description | +|------|---------|-------------| +| `tle` | ~160 bytes fixed | Parsed mean elements (not raw text) | +| `eci_position` | 48 bytes | x,y,z + vx,vy,vz (km, km/s) in TEME | +| `geodetic` | 24 bytes | lat, lon (radians), alt (km) above WGS-84 | +| `topocentric` | 32 bytes | azimuth, elevation, range, range_rate | +| `observer` | 24 bytes | lat, lon (radians), alt_m (meters) | +| `pass_event` | 56 bytes | AOS/MAX/LOS times + max_el + AOS/LOS az | + +### TLE Internal Struct +Stores all parsed mean elements from the two-line format: +- epoch (Julian date, float64) +- inclination, eccentricity, RAAN, arg_perigee, mean_anomaly (radians, float64) +- mean_motion (rev/day, float64), mean_motion_dot, mean_motion_ddot +- bstar (drag coefficient, float64) +- norad_id (int32), elset_num (int32), rev_num (int32) +- classification (char), intl_designator (8 chars) +- ephemeris_type (int8) + +## Constant Chain of Custody + +**This is the most critical design constraint.** + +TLEs are mean elements fitted using WGS-72 constants. The elements absorb geodetic model biases — using WGS-84 constants for propagation silently corrupts position accuracy by kilometers. + +### Rules +1. **SGP4 propagation**: WGS-72 constants ONLY (mu, ae, J2, J3, J4, ke) +2. **Coordinate output** (geodetic, topocentric): Convert to WGS-84 (a=6378.137km, f=1/298.257223563) +3. **TEME frame**: Use only 4 of 106 IAU-80 nutation terms (matching SGP4's internal model) +4. **Never mix**: WGS-72 propagation + WGS-84 output. No other combination. + +### WGS-72 Constants (from Hoots & Roehrich STR#3) +```c +#define WGS72_MU 398600.8 /* km^3/s^2 */ +#define WGS72_AE 6378.135 /* km */ +#define WGS72_J2 0.001082616 +#define WGS72_J3 -0.00000253881 +#define WGS72_J4 -0.00000165597 +#define WGS72_KE 0.0743669161 /* (min)^(-1), = sqrt(mu) * 60 / ae^(3/2) */ +#define WGS72_XPDOTP 1440.0 / (2.0 * M_PI) /* min/rev */ +``` + +### WGS-84 Constants (for output only) +```c +#define WGS84_A 6378.137 /* km */ +#define WGS84_F (1.0 / 298.257223563) +#define WGS84_E2 (WGS84_F * (2.0 - WGS84_F)) +``` + +## sat_code Submodule + +Bill Gray's SGP4 implementation: https://github.com/Bill-Gray/sat_code + +Key files we use: +- `sgp4.c` / `sgp4.h` — SGP4/SDP4 propagator +- `norad.h` — TLE struct definitions and constants + +The submodule lives at `lib/sat_code/`. To initialize: +```bash +git submodule update --init +``` + +### Integration Pattern +```c +#include "lib/sat_code/norad.h" + +// Parse TLE lines into sat_code's tle_t struct +// Call SGP4_init() once per TLE +// Call SGP4() with minutes-since-epoch for each propagation +``` + +## Testing + +### Vallado 518 Test Vectors +The definitive SGP4 verification dataset. Each row: NORAD ID, minutes since epoch, expected x,y,z,vx,vy,vz. All 518 must pass to machine epsilon before any other work proceeds. + +### Regression Tests +Standard PostgreSQL `make installcheck` framework: +- `test/sql/*.sql` — test queries +- `test/expected/*.out` — expected output +- Tests run against a temporary database + +### Test Categories +1. **tle_parse** — TLE input/output round-trip, malformed input rejection +2. **sgp4_propagate** — Vallado vectors, edge cases (deep space, high eccentricity) +3. **coord_transforms** — TEME->geodetic, TEME->topocentric accuracy +4. **pass_prediction** — Known ISS passes, edge cases (polar, retrograde) +5. **gist_index** — Index scan vs sequential scan equivalence + +## Coding Style +- Standard PostgreSQL extension C style +- `ereport(ERROR, ...)` for user-facing errors, never `elog(ERROR, ...)` +- All memory allocation through `palloc`/`pfree` (PostgreSQL memory contexts) +- Comments explain "why", not "what" +- No global mutable state — all computation from function arguments +- Functions that call `SGP4()` must handle the error return code + +## Git Conventions +- One commit per logical change +- Branch per phase: `phase/1-tle-sgp4`, `phase/2-coordinates`, etc. +- Tag releases: `v0.1.0`, `v0.2.0`, etc. +- Commit messages: imperative mood, no AI attribution diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..94951e7 --- /dev/null +++ b/LICENSE @@ -0,0 +1,19 @@ +PostgreSQL License + +Copyright (c) 2025, Ryan Malloy + +Permission to use, copy, modify, and distribute this software and its +documentation for any purpose, without fee, and without a written agreement +is hereby granted, provided that the above copyright notice and this +paragraph and the following two paragraphs appear in all copies. + +IN NO EVENT SHALL THE AUTHORS BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, +SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS, +ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE +AUTHORS HAVE BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +THE AUTHORS SPECIFICALLY DISCLAIM ANY WARRANTIES, INCLUDING, BUT NOT LIMITED +TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE +AUTHORS HAVE NO OBLIGATIONS TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, +ENHANCEMENTS, OR MODIFICATIONS. diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..0b9152c --- /dev/null +++ b/Makefile @@ -0,0 +1,36 @@ +MODULE_big = pg_orbit +EXTENSION = pg_orbit +DATA = sql/pg_orbit--0.1.0.sql + +# Our extension C sources +OBJS = src/pg_orbit.o src/tle_type.o src/eci_type.o src/observer_type.o \ + src/sgp4_funcs.o src/coord_funcs.o src/pass_funcs.o src/gist_tle.o + +# sat_code C++ sources (compiled with g++, linked with extern "C" symbols) +SAT_CODE_DIR = lib/sat_code +SAT_CODE_SRCS = $(SAT_CODE_DIR)/sgp4.cpp $(SAT_CODE_DIR)/sdp4.cpp \ + $(SAT_CODE_DIR)/deep.cpp $(SAT_CODE_DIR)/common.cpp \ + $(SAT_CODE_DIR)/basics.cpp $(SAT_CODE_DIR)/get_el.cpp \ + $(SAT_CODE_DIR)/tle_out.cpp +SAT_CODE_OBJS = $(SAT_CODE_SRCS:.cpp=.o) + +OBJS += $(SAT_CODE_OBJS) + +# Regression tests +REGRESS = tle_parse sgp4_propagate coord_transforms pass_prediction gist_index +REGRESS_OPTS = --inputdir=test + +# Need C++ runtime for sat_code +SHLIB_LINK += -lstdc++ -lm + +# Compiler flags +PG_CPPFLAGS = -I$(SAT_CODE_DIR) + +# Use PGXS +PG_CONFIG ?= pg_config +PGXS := $(shell $(PG_CONFIG) --pgxs) +include $(PGXS) + +# Rule for compiling sat_code C++ files +$(SAT_CODE_DIR)/%.o: $(SAT_CODE_DIR)/%.cpp + $(CXX) $(CXXFLAGS) -fPIC -I$(SAT_CODE_DIR) -c -o $@ $< diff --git a/README.md b/README.md new file mode 100644 index 0000000..9645aad --- /dev/null +++ b/README.md @@ -0,0 +1,305 @@ +# 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](https://github.com/Bill-Gray/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) + +```bash +git clone --recurse-submodules https://github.com/... +cd pg_orbit +make +sudo make install +``` + +Then in your database: + +```sql +CREATE EXTENSION pg_orbit; +``` + +If you cloned without `--recurse-submodules`, initialize the sat_code dependency: + +```bash +git submodule update --init +``` + +## Quick Start + +```sql +-- 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): + +```sql +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: + +```sql +-- 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. + +```sql +-- 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 + +```bash +# 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: + +```bash +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. + +```bash +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](LICENSE). Copyright (c) 2025, Ryan Malloy. + +The bundled sat_code library is separately licensed under the MIT license. diff --git a/docs/DESIGN.md b/docs/DESIGN.md new file mode 100644 index 0000000..be5fb28 --- /dev/null +++ b/docs/DESIGN.md @@ -0,0 +1,590 @@ +# pg_orbit Design Document + +Internal architecture notes. Documents WHY decisions were made, +not how to use the extension. Intended audience: future maintainers +who need to modify pg_orbit without breaking physical correctness. + + +## 1. Constant Chain of Custody + +This is the single most critical design constraint in the extension. +Get it wrong and positions silently drift by kilometers. There is no +runtime check that can detect this class of error after the fact. + +### The Problem + +Two-Line Elements are not raw orbital measurements. They are *mean* +elements produced by a differential correction process that fits +observed positions against an SGP4 propagator running with a specific +set of geopotential constants (WGS-72). The mean elements absorb +geodetic model biases -- the eccentricity, inclination, and mean +motion encode corrections that only make physical sense when propagated +with the same constants used to generate them. + +Substituting WGS-84 constants (or any other set) into the propagator +does not "upgrade" accuracy. It breaks the internal consistency of the +element set. The resulting position error can exceed the natural +prediction error of the TLE by an order of magnitude. + +### The Rules + +1. **SGP4/SDP4 propagation**: WGS-72 constants only (mu, ae, J2, J3, J4, ke). + These flow through sat_code's `norad_in.h` defines and are never + overridden. + +2. **Coordinate output** (geodetic lat/lon/alt, topocentric az/el/range): + WGS-84 ellipsoid (a = 6378.137 km, f = 1/298.257223563). This is + the modern standard for ground-station positioning and GPS receivers. + +3. **TEME frame**: The SGP4 output frame (True Equator, Mean Equinox) + uses only 4 of the 106 IAU-80 nutation terms. Applying the full + nutation model would "correct" for effects that SGP4 already accounts + for internally, introducing error rather than removing it. + +4. **No other combination is valid.** WGS-72 for propagation, WGS-84 for + output. Perigee and apogee altitudes from `tle_perigee()` and + `tle_apogee()` use WGS-72 AE because they derive from mean elements. + Geodetic altitude from `eci_to_geodetic()` uses WGS-84 because it + converts a physical position. + +### Constant Inventory + +| Constant | Source Paper | Value | pg_orbit Location | sat_code Location | +|----------|-------------|-------|-------------------|-------------------| +| ae (equatorial radius) | Hoots & Roehrich STR#3 | 6378.135 km | `types.h:20` (WGS72_AE) | `norad_in.h:64` (earth_radius_in_km) | +| J2 | Hoots & Roehrich STR#3 | 0.001082616 | `types.h:21` (WGS72_J2) | `norad_in.h:69` (xj2) | +| J3 | Hoots & Roehrich STR#3 | -2.53881e-6 | `types.h:22` (WGS72_J3) | `norad_in.h:63` (xj3) | +| J4 | Hoots & Roehrich STR#3 | -1.65597e-6 | `types.h:23` (WGS72_J4) | `norad_in.h:79` (xj4) | +| ke | Hoots & Roehrich STR#3 | 0.0743669161331734132 min^-1 | `types.h:24` (WGS72_KE) | `norad_in.h:83` (xke) | +| mu | Hoots & Roehrich STR#3 | 398600.8 km^3/s^2 | `types.h:19` (WGS72_MU) | (implicit in xke) | +| WGS-84 a | NIMA TR8350.2 | 6378.137 km | `types.h:31` (WGS84_A) | -- | +| WGS-84 f | NIMA TR8350.2 | 1/298.257223563 | `types.h:32` (WGS84_F) | -- | + +Note that `types.h` carries a parallel copy of the WGS-72 constants +even though sat_code defines them in `norad_in.h`. This is intentional: +`types.h` is the single header for all pg_orbit C sources, and +`norad_in.h` is an internal sat_code header not meant for external +consumers. The GiST index (`gist_tle.c`) and TLE accessor functions +(`tle_type.c`) need KE and AE without pulling in sat_code internals. +The values MUST be identical. + +### Why Two Copies of AE? + +`tle_perigee()`, `tle_apogee()`, and the GiST altitude-band computation +all use `WGS72_KE` and `WGS72_AE` from `types.h`. They compute: + + a_er = (KE / n)^(2/3) [earth radii] + perigee_km = a_er * (1 - e) * AE - AE + +This MUST use WGS-72 AE (6378.135), not WGS-84 (6378.137), because `n` +is a mean motion fitted against the WGS-72 geopotential. Using the +wrong radius shifts every altitude by 2 meters -- small in absolute terms +but wrong in principle, and the error compounds in index operations. + + +## 2. SGP4 Implementation Choice + +pg_orbit wraps Bill Gray's `sat_code` library (MIT license, Project Pluto). + +### Why sat_code over alternatives + +**Vallado's reference implementation** (from the STR#3 revision paper) is +the canonical source but has two problems: it is written in C++ with heavy +use of global state, and its license is unclear for embedding in a +PostgreSQL extension. + +**libsgp4** (various forks on GitHub) is typically a C++ class hierarchy +that assumes an object-per-satellite lifecycle. This conflicts with +PostgreSQL's per-call execution model where we cannot persist C++ objects +across function invocations. + +**sat_code** was chosen because: + +1. **Pure C linkage.** All public functions are declared `extern "C"` in + `norad.h` (lines 97-133). The library is compiled as C++ but exposes + a flat C function interface: `SGP4_init()`, `SGP4()`, `SDP4_init()`, + `SDP4()`, `parse_elements()`, `select_ephemeris()`. + +2. **No global mutable state.** The propagator state lives in a caller- + allocated `double params[N_SAT_PARAMS]` array. This maps directly to + PostgreSQL's `palloc`-based memory model: allocate before propagation, + free after. + +3. **Includes deep-space SDP4.** Many SGP4 implementations only handle + near-earth orbits (period < 225 minutes). sat_code includes the full + SDP4 with lunar/solar perturbations via `deep.cpp`, handling GEO, + Molniya, and GPS orbits. + +4. **MIT license.** Compatible with the PostgreSQL License for embedding + in a shared library. + +5. **Actively maintained.** Used in Bill Gray's Find_Orb production + astrometry software. Bug fixes reach us through the git submodule. + +### Build Integration + +The Makefile compiles sat_code's `.cpp` files with `g++` and links the +resulting `.o` files into the PostgreSQL shared library alongside our C +sources. The `-lstdc++` link flag pulls in the C++ runtime. This is the +same pattern used by PostGIS for GEOS integration (C extension linking +C++ library objects). + +``` +src/*.c --> gcc --> .o --| +lib/sat_code/*.cpp -> g++ -> .o --|--> pg_orbit.so + -lstdc++ -lm +``` + +The `-I$(SAT_CODE_DIR)` flag lets our C files `#include "norad.h"` +directly. + + +## 3. Type System Design + +### Design Principles + +Every pg_orbit type is fixed-size, not varlena. This means: + +- No TOAST overhead (no detoasting on access) +- Direct pointer access via `PG_GETARG_POINTER(n)` -- no copy +- Predictable memory layout for binary I/O (`tle_recv`/`tle_send`) +- All types use `ALIGNMENT = double` to satisfy the strictest member + alignment requirement (all structs contain `double` fields) +- `STORAGE = plain` -- the type is never compressed or moved to TOAST + +### TLE Type (112 bytes) + +The TLE type stores **parsed mean elements**, not raw text. + +This is the most important type design decision. Alternatives considered: + +1. **Store raw 69+69 character lines, parse on every propagation.** + Rejected. Parsing is ~10x slower than propagation itself for a + pre-initialized model. Every `sgp4_propagate()` call would pay + the parse cost. + +2. **Store raw text plus parsed cache.** + Rejected. Doubles the storage for no benefit. The parsed form + round-trips perfectly through `write_elements_in_tle_format()`. + +3. **Store parsed mean elements only.** + Chosen. Input validation happens once at `tle_in()` time via + sat_code's `parse_elements()`. The text representation is + reconstructed on output via `write_elements_in_tle_format()`. + +Storage layout: + +``` +Offset Size Field + 0 88 11 doubles: epoch, inclination, raan, eccentricity, + arg_perigee, mean_anomaly, mean_motion, + mean_motion_dot, mean_motion_ddot, bstar + (one unused slot in the double array for alignment) + 80 12 3 int32: norad_id, elset_num, rev_num + 92 12 classification (1), ephemeris_type (1), + intl_desig (9), pad (1) +---- +112 bytes total +``` + +The SQL definition declares `INTERNALLENGTH = 112`. This is smaller +than the raw two-line text (138+ bytes with line separator) and avoids +the 4-byte varlena header overhead. + +Angular elements are stored in radians (matching sat_code's internal +representation). Accessor functions convert to degrees for human +consumption. Mean motion is stored in radians/minute; the accessor +returns revolutions/day. + +### ECI Position Type (48 bytes) + +Six doubles: x, y, z (km), vx, vy, vz (km/s). + +SGP4 outputs velocity in km/min. We convert to km/s at the boundary +(`sgp4_funcs.c`, lines 181-183: `vel[i] / 60.0`). This conversion +happens exactly once, at the point where the pg_eci struct is populated. +Internally, all velocity in pg_orbit is km/s. + +### Geodetic Type (24 bytes) + +Three doubles: lat, lon (radians), alt (km above WGS-84). + +Radians internally, degrees in text representation. The accessor +functions perform the conversion. + +### Observer Type (24 bytes) + +Three doubles: lat, lon (radians), alt_m (meters above WGS-84). + +Note the asymmetry with geodetic: observer altitude is in meters +(matching GPS receiver output and ground station databases), while +geodetic altitude is in km (matching orbital altitude conventions). +This prevents unit confusion at the API boundary -- you set up a +ground station in meters, you get satellite altitude in kilometers. + +### Topocentric Type (32 bytes) + +Four doubles: azimuth, elevation (radians), range_km, range_rate (km/s). + +### Pass Event Type (48 bytes) + +Three int64 (TimestampTz): aos_time, max_el_time, los_time. +Three doubles: max_elevation (degrees), aos_azimuth (degrees), +los_azimuth (degrees). + +Times are stored as native PostgreSQL TimestampTz values (microseconds +since 2000-01-01 00:00:00 UTC). This allows direct comparison with +SQL timestamp expressions without conversion. + + +## 4. TEME to ECEF Transform + +SGP4 outputs position and velocity in the TEME (True Equator, Mean +Equinox) frame. Converting to Earth-fixed coordinates for geodetic +and topocentric output requires a frame rotation. + +### The Rotation + +TEME to ECEF is a single Z-axis rotation by the negative of the +Greenwich Mean Sidereal Time (GMST) angle: + +``` + [x_ecef] [ cos(g) sin(g) 0 ] [x_teme] + [y_ecef] = [-sin(g) cos(g) 0 ] [y_teme] + [z_ecef] [ 0 0 1 ] [z_teme] +``` + +This is implemented in `coord_funcs.c:teme_to_ecef()` and +`pass_funcs.c:teme_to_ecef()` (duplicated -- see note below). + +### GMST Computation + +GMST is computed from the Julian date using the IAU 1982 formula, +which is the same low-precision model that SGP4 uses internally. +From Vallado (2013) Eq. 3-47: + +```c + T_UT1 = (JD - 2451545.0) / 36525.0 + GMST = 67310.54841 + + (876600*3600 + 8640184.812866) * T_UT1 + + 0.093104 * T_UT1^2 + - 6.2e-6 * T_UT1^3 +``` + +The result is in seconds of time, converted to radians by multiplying +by pi/43200 and normalized to [0, 2*pi). + +We deliberately do NOT use a higher-precision GMST model (e.g., IAU 2000A). +The SGP4 output is only accurate to the precision of its own GMST model; +applying a more precise rotation would not improve the final position and +could introduce a systematic offset. + +### Velocity Transform + +The velocity transform includes the Earth angular velocity cross product +term, accounting for the rotating reference frame: + +``` + v_ecef = R(-g) * v_teme + omega_E x r_ecef +``` + +where omega_E = 7.2921158553e-5 rad/s. In `pass_funcs.c`, the velocity +cross product uses omega_E in rad/min (multiplied by 60.0) because +sat_code outputs velocity in km/min. + +### ECEF to Geodetic + +The ECEF-to-geodetic conversion uses Bowring's iterative method, which +converges in 2-3 iterations for LEO altitudes: + +``` + phi_0 = atan2(z, p * (1 - e^2)) [initial estimate] + N = a / sqrt(1 - e^2 * sin^2(phi)) [prime vertical radius] + phi = atan2(z + e^2 * N * sin(phi), p) [iterate] +``` + +Convergence criterion: 1.0e-12 radians (sub-millimeter at Earth's surface). +Maximum 10 iterations as a safety bound, though LEO cases converge in 2-3. + +We chose Bowring's method over Vermeille's or Heikkinen's closed-form +solutions because it is simpler, well-tested, and the iteration cost is +negligible compared to the SGP4 propagation that precedes it. + +### Why Duplicated Helpers + +`coord_funcs.c` and `pass_funcs.c` each contain their own static copies +of `pg_tle_to_sat_code()`, `gmst_from_jd()`, `teme_to_ecef()`, +`observer_to_ecef()`, and `ecef_to_topocentric()`. All `.c` files link +into a single `.so`, so we could share these through a common object file. +We chose duplication instead because: + +1. Each translation unit is self-contained. No hidden coupling between + files through shared internal functions. +2. The functions are small (5-20 lines each). The binary size increase + is negligible. +3. The compiler can inline them within each translation unit. +4. If the helpers ever need to diverge (e.g., pass_funcs using km/min + velocity while coord_funcs uses km/s), they can do so without + affecting each other. + + +## 5. Pass Prediction Algorithm + +### The Core Problem + +Given a TLE and a ground observer, find all time intervals where the +satellite is above the observer's local horizon. SGP4 has no closed-form +inverse for the elevation function -- you cannot algebraically solve +"at what time does elevation = 0?" You have to evaluate the function +at many points and search for zero crossings. + +This is the same fundamental approach used by Skyfield's `find_events()` +and Gpredict's pass finder. + +### Algorithm + +1. **Coarse scan**: Step through the search window at 30-second intervals, + evaluating elevation at each step. A 30-second step is fine enough + for LEO (~90 min period, ~10 min pass duration) that no pass shorter + than about 1 minute gets skipped, and coarse enough that a 7-day + window requires ~20,000 propagation calls (fast on modern hardware). + +2. **AOS detection**: When elevation transitions from negative to + positive between consecutive coarse steps, a rising edge has been + found. Bisect the interval to locate AOS to within 0.1 second + tolerance (BISECT_TOL_JD = 0.1/86400). + +3. **LOS detection**: Continue the coarse scan from AOS until elevation + goes negative again. Bisect to refine LOS to the same 0.1 second + tolerance. + +4. **Peak elevation**: Between AOS and LOS, the elevation function is + unimodal (one peak). Use ternary search (50 iterations) to locate + the maximum. Ternary search is chosen over golden section because + the convergence rate is adequate and the implementation is simpler. + +5. **Degenerate pass filter**: Passes shorter than 10 seconds + (MIN_PASS_DURATION_JD) are discarded. These are typically numerical + noise from orbits just barely grazing the horizon. + +6. **Minimum elevation filter**: Passes whose peak elevation falls below + the caller-specified threshold are silently skipped. The scan resumes + from their LOS. + +7. **Post-LOS gap**: After finding a complete pass, the scan resumes + 1 minute past LOS (POST_LOS_GAP_JD) to avoid re-detecting the + same descending arc. + +### Error Handling in the Scan + +If SGP4 returns a hard error (negative semi-major axis, nearly parabolic, +convergence failure) during the elevation scan, `elevation_at_jd()` +returns -pi radians. This is well below any real horizon elevation, so +the scan treats the point as "below horizon" and continues. The pass +finder does NOT abort the query on propagation errors during scanning -- +a TLE might be valid for part of the search window and decayed for +the rest. + +### SRF (Set-Returning Function) Pattern + +`predict_passes()` uses PostgreSQL's ValuePerCall SRF protocol: + +- First call: allocate a `predict_passes_ctx` struct in + `funcctx->multi_call_memory_ctx`, copy the TLE and observer into it, + initialize the scan position. +- Each subsequent call: call `find_next_pass()` starting from the + current scan position. If a pass is found, advance `current_jd` + past LOS and return the pass. If not, return DONE. + +The TLE and observer are copied into the multi_call context because +the original function arguments may be freed between calls. + + +## 6. GiST Index Architecture + +### The Indexing Problem + +The natural query pattern for conjunction screening is: "which TLEs +have orbits that could intersect with this TLE?" Full conjunction +analysis requires propagating both objects and computing distance at +every timestep -- far too expensive for an index scan. + +### Altitude-Band Approximation + +Every orbit has a perigee altitude and an apogee altitude. Two orbits +can only come close to each other if their altitude ranges overlap. +This is a **necessary but not sufficient** condition for conjunction. + +The GiST index stores [perigee_km, apogee_km] as its internal key +(`tle_alt_range` struct: two doubles, 16 bytes). This is derived from +the mean elements via Kepler's third law: + +``` + a = (KE / n)^(2/3) [earth radii, WGS-72] + perigee_km = a * (1 - e) * AE - AE + apogee_km = a * (1 + e) * AE - AE +``` + +The altitude band uses WGS-72 constants because n and e are WGS-72 +mean elements. + +### GiST Support Functions + +| Function | Strategy | Description | +|----------|----------|-------------| +| compress | -- | Leaf: extract [perigee, apogee] from pg_tle. Internal: pass through (already a range). | +| decompress | -- | Identity. We operate on compressed keys directly. | +| consistent | `&&` (3), `@>` (7), `<@` (8) | Does the subtree's bounding range overlap/contain the query range? | +| union | -- | [min(low), max(high)] over all entries in a node. | +| penalty | -- | Increase in altitude span (km) from expanding the node's range to include the new entry. | +| picksplit | -- | Sort entries by altitude midpoint, split at the median. | +| same | -- | Endpoint comparison within 1e-9 km tolerance. | +| distance | `<->` (15) | Minimum gap between altitude bands (0 if overlapping). Drives KNN queries. | + +### Always Recheck + +`gist_tle_consistent()` always sets `*recheck = true`. Altitude band +overlap eliminates the vast majority of non-candidates (a LEO satellite +cannot conjunct with a GEO satellite), but it cannot confirm conjunction. +Two objects at the same altitude but on opposite sides of the earth are +not close. After the index scan filters candidates, the query must +propagate both TLEs to the reference time and compute actual distance. + +### Picksplit Strategy + +The picksplit function sorts entries by the midpoint of their altitude +band ((perigee + apogee) / 2) and splits at the median. This is optimal +for a 1-D range index because it minimizes overlap between the two +resulting subtrees. Orbits at similar altitudes end up in the same +subtree, which matches the access pattern of conjunction screening +(you are usually querying within a narrow altitude band). + + +## 7. Memory Management + +### Allocation Strategy + +All heap allocation goes through `palloc()` / `pfree()`. No `malloc()`, +no `new`, no static buffers. + +- **Single-shot propagation** (`sgp4_propagate`, `tle_distance`): + `palloc(sizeof(double) * N_SAT_PARAMS)` for the params array, + propagate, `pfree(params)`. The params array lives in the current + memory context and is freed before the function returns. + +- **Set-returning functions** (`sgp4_propagate_series`, `ground_track`, + `predict_passes`): The SRF context struct (containing `tle_t`, + `params[N_SAT_PARAMS]`, and scan state) is allocated in + `funcctx->multi_call_memory_ctx`. This context survives across calls + and is freed automatically when the SRF completes. + +- **Type I/O functions** (`tle_in`, `eci_in`, etc.): The result struct + is allocated with `palloc()` in the current context. PostgreSQL + manages its lifecycle. + +### No Global Mutable State + +There are no file-scope variables, no static locals that accumulate +state, no caches. Every function computes from its arguments alone. +This is required for `PARALLEL SAFE` (all pg_orbit functions are +declared PARALLEL SAFE) and avoids cross-session contamination in +a multi-backend PostgreSQL server. + +sat_code itself has no global mutable state. The propagator state is +entirely in the `params[]` array and the `tle_t` struct, both of which +are caller-provided. + +### SRF Context Layout + +For `sgp4_propagate_series` and `ground_track`, the params array is +embedded directly in the context struct (not a separate allocation): + +```c +typedef struct { + tle_t sat; + double params[N_SAT_PARAMS]; /* ~92 doubles, embedded */ + int is_deep; + double epoch_jd; + int64 start_ts; + int64 step_usec; +} propagate_series_ctx; +``` + +This puts everything in a single allocation, reducing palloc overhead +and keeping the data contiguous in memory for cache locality during +the per-row propagation loop. + + +## 8. Error Handling + +### Error Classification + +sat_code returns integer error codes from SGP4() and SDP4(): + +| Code | Constant | Severity | Meaning | pg_orbit Response | +|------|----------|----------|---------|-------------------| +| 0 | -- | OK | Normal | Return result | +| -1 | SXPX_ERR_NEARLY_PARABOLIC | Fatal | eccentricity >= 1 | `ereport(ERROR)` | +| -2 | SXPX_ERR_NEGATIVE_MAJOR_AXIS | Fatal | Decayed orbit | `ereport(ERROR)` | +| -3 | SXPX_WARN_ORBIT_WITHIN_EARTH | Warning | Entire orbit below surface | `ereport(NOTICE)`, return result | +| -4 | SXPX_WARN_PERIGEE_WITHIN_EARTH | Warning | Perigee below surface | `ereport(NOTICE)`, return result | +| -5 | SXPX_ERR_NEGATIVE_XN | Fatal | Negative mean motion | `ereport(ERROR)` | +| -6 | SXPX_ERR_CONVERGENCE_FAIL | Fatal | Kepler equation diverged | `ereport(ERROR)` | + +The distinction between warnings (-3, -4) and errors (-1, -2, -5, -6) is +physical. A satellite with perigee within the earth is plausible during +reentry or shortly after launch -- the state vector is still mathematically +valid. A negative semi-major axis means the orbit has decayed past the +point where the model produces meaningful output. + +### Error Handling in Different Contexts + +**Direct propagation** (`sgp4_propagate`, `tle_distance`): Fatal errors +abort the query via `ereport(ERROR)`. Warnings emit a `NOTICE` and +return the result. + +**Pass prediction** (`elevation_at_jd` in `pass_funcs.c`): Fatal errors +return -pi elevation instead of aborting. The scan treats this as "below +horizon" and continues. A TLE might be valid for the first 3 days of a +7-day search window and then decay -- the pass finder should return the +valid passes, not abort. + +**SRF propagation** (`sgp4_propagate_series`): Fatal errors abort the +entire series. There is no meaningful way to "skip" a bad timestep in a +continuous time series. + +### Input Validation + +TLE parsing errors from `parse_elements()` are caught in `tle_in()` and +reported as `ERRCODE_INVALID_TEXT_REPRESENTATION`. Invalid TLEs never +make it into storage. + +`select_ephemeris()` returning -1 (invalid mean motion or eccentricity +range) is caught at propagation time, not at storage time. A TLE with +marginal elements might parse correctly but fail when you try to +initialize the propagator. + + +## 9. Theory-to-Code Mapping + +This table maps key equations from the SGP4 theory papers to their +implementation in pg_orbit and sat_code. + +| Theory | Paper | What | Code Location | +|--------|-------|------|---------------| +| Mean element recovery | Brouwer (1959) | Recover original mean motion (xnodp) and semi-major axis (aodp) from input TLE elements, removing secular J2 perturbations | `sat_code/common.cpp:sxpall_common_init()` lines 17-35 | +| Secular perturbations | Lane & Cranford (1969), Hoots & Roehrich STR#3 | Secular rates of mean anomaly, argument of perigee, and RAAN due to J2, J4 | `sat_code/common.cpp:sxpx_common_init()` lines 86-101 | +| Atmospheric drag | Hoots & Roehrich STR#3 | B* formulation of drag, C1/C2/C4 coefficients, perigee-dependent s parameter | `sat_code/common.cpp:sxpx_common_init()` lines 47-84; `sat_code/sgp4.cpp:SGP4_init()` | +| Short-period perturbations | Lane & Cranford (1969), Brouwer (1959) | Oscillatory corrections to radius, argument of latitude, node, and inclination | `sat_code/common.cpp:sxpx_posn_vel()` lines 121-229 | +| Kepler equation | Classical | Newton-Raphson with second-order correction, bounded first step | `sat_code/common.cpp:sxpx_posn_vel()` lines 175-208 | +| Deep-space resonance | Hujsak (1979) | Lunar and solar gravitational perturbations, geopotential resonance for 12-hour and 24-hour orbits | `sat_code/deep.cpp:Deep_dpinit()`, `Deep_dpsec()`, `Deep_dpper()` | +| Near-earth propagation | Hoots & Roehrich STR#3 | SGP4 main loop: secular + short-period + drag terms | `sat_code/sgp4.cpp:SGP4()` | +| Deep-space propagation | Hoots & Roehrich STR#3 | SDP4: SGP4 core + deep-space secular/periodic corrections | `sat_code/sdp4.cpp:SDP4()` | +| Semi-major axis from n | Kepler's third law | a = (KE / n)^(2/3) in earth radii | `src/tle_type.c:tle_perigee()` line 415; `src/gist_tle.c:tle_to_alt_range()` line 76 | +| GMST | Vallado (2013) Eq. 3-47 | Greenwich Mean Sidereal Time from Julian date | `src/coord_funcs.c:gmst_from_jd()` lines 59-73; `src/pass_funcs.c:gmst_from_jd()` lines 129-151 | +| TEME to ECEF | Vallado (2013) | Z-axis rotation by -GMST, velocity cross-product correction | `src/coord_funcs.c:teme_to_ecef()` lines 83-103; `src/pass_funcs.c:teme_to_ecef()` lines 157-179 | +| Geodetic from ECEF | Bowring (1976) | Iterative latitude from ECEF Cartesian coordinates on WGS-84 | `src/coord_funcs.c:ecef_to_geodetic()` lines 109-138 | +| Topocentric transform | Standard SEZ | ECEF range vector rotated to South-East-Zenith, azimuth from north | `src/coord_funcs.c:ecef_to_topocentric()` lines 163-188 | +| Observer to ECEF | Geodesy standard | WGS-84 ellipsoid surface point to Cartesian | `src/coord_funcs.c:observer_to_ecef()` lines 143-156 | +| Range rate | Dot product | Projection of relative velocity onto line-of-sight unit vector | `src/coord_funcs.c:eci_to_topocentric()` line 618 | +| Near/deep selection | Hoots & Roehrich STR#3 | Period threshold: 225 minutes (n < 2*pi/225 rad/min) | `sat_code/norad.h:select_ephemeris()` | diff --git a/lib/sat_code b/lib/sat_code new file mode 160000 index 0000000..ff7b989 --- /dev/null +++ b/lib/sat_code @@ -0,0 +1 @@ +Subproject commit ff7b98957dfa2979700a482bde9de9542807293e diff --git a/pg_orbit.control b/pg_orbit.control new file mode 100644 index 0000000..b291ade --- /dev/null +++ b/pg_orbit.control @@ -0,0 +1,4 @@ +comment = 'Orbital mechanics types and functions for PostgreSQL' +default_version = '0.1.0' +module_pathname = '$libdir/pg_orbit' +relocatable = true diff --git a/sql/pg_orbit--0.1.0.sql b/sql/pg_orbit--0.1.0.sql new file mode 100644 index 0000000..33df63d --- /dev/null +++ b/sql/pg_orbit--0.1.0.sql @@ -0,0 +1,491 @@ +-- pg_orbit -- 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)'; + + +-- ============================================================ +-- 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'; + + +-- ============================================================ +-- 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_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)'; + + +-- ============================================================ +-- 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 altitude bands intersect? +CREATE FUNCTION tle_overlap(tle, tle) RETURNS boolean + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +-- Altitude distance operator +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 'Altitude band overlap — 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)'; + + +-- ============================================================ +-- GiST operator class for altitude-band indexing +-- ============================================================ + +-- 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); diff --git a/src/coord_funcs.c b/src/coord_funcs.c new file mode 100644 index 0000000..7f64567 --- /dev/null +++ b/src/coord_funcs.c @@ -0,0 +1,820 @@ +/* + * coord_funcs.c -- Coordinate transform functions for pg_orbit + * + * TEME -> WGS-84 geodetic and TEME -> topocentric transforms. + * + * SGP4 outputs in the TEME (True Equator, Mean Equinox) frame. + * Converting to Earth-fixed coordinates requires rotating by GMST. + * We use only the 4 IAU-80 nutation terms that SGP4 itself uses -- + * applying the full 106-term model would "correct" what SGP4 + * already accounts for. + * + * The TEME->ITRF rotation is a single Z-axis rotation by -GMST. + * After that, standard geodetic conversions on the WGS-84 ellipsoid. + */ + +#include "postgres.h" +#include "fmgr.h" +#include "funcapi.h" +#include "utils/timestamp.h" +#include "libpq/pqformat.h" +#include "norad.h" +#include "types.h" +#include +#include + +#define DEG_TO_RAD (M_PI / 180.0) +#define RAD_TO_DEG (180.0 / M_PI) + +PG_FUNCTION_INFO_V1(geodetic_in); +PG_FUNCTION_INFO_V1(geodetic_out); +PG_FUNCTION_INFO_V1(geodetic_recv); +PG_FUNCTION_INFO_V1(geodetic_send); +PG_FUNCTION_INFO_V1(geodetic_lat); +PG_FUNCTION_INFO_V1(geodetic_lon); +PG_FUNCTION_INFO_V1(geodetic_alt); +PG_FUNCTION_INFO_V1(topocentric_in); +PG_FUNCTION_INFO_V1(topocentric_out); +PG_FUNCTION_INFO_V1(topocentric_recv); +PG_FUNCTION_INFO_V1(topocentric_send); +PG_FUNCTION_INFO_V1(topo_azimuth); +PG_FUNCTION_INFO_V1(topo_elevation); +PG_FUNCTION_INFO_V1(topo_range); +PG_FUNCTION_INFO_V1(topo_range_rate); +PG_FUNCTION_INFO_V1(eci_to_geodetic); +PG_FUNCTION_INFO_V1(eci_to_topocentric); +PG_FUNCTION_INFO_V1(subsatellite_point); +PG_FUNCTION_INFO_V1(ground_track); + +/* ================================================================ + * Internal helpers -- GMST, frame rotation, geodetic conversion + * ================================================================ + */ + +/* + * Greenwich Mean Sidereal Time from Julian date. + * Vallado, "Fundamentals of Astrodynamics", Eq. 3-47. + * Returns angle in radians. + */ +static double +gmst_from_jd(double jd) +{ + double t_ut1 = (jd - 2451545.0) / 36525.0; + double gmst = 67310.54841 + + (876600.0 * 3600.0 + 8640184.812866) * t_ut1 + + 0.093104 * t_ut1 * t_ut1 + - 6.2e-6 * t_ut1 * t_ut1 * t_ut1; + + /* Convert seconds of time to radians, then normalize to [0, 2*pi) */ + gmst = fmod(gmst * M_PI / 43200.0, 2.0 * M_PI); + if (gmst < 0.0) + gmst += 2.0 * M_PI; + return gmst; +} + +/* + * Rotate TEME position/velocity to ECEF (Earth-Centered Earth-Fixed). + * Single rotation around Z by -GMST angle. + * omega_e = Earth rotation rate = 7.2921158553e-5 rad/s + */ +#define OMEGA_EARTH 7.2921158553e-5 /* rad/s */ + +static void +teme_to_ecef(const double *pos_teme, const double *vel_teme, + double gmst, double *pos_ecef, double *vel_ecef) +{ + double cos_g = cos(gmst); + double sin_g = sin(gmst); + + /* Position rotation */ + pos_ecef[0] = cos_g * pos_teme[0] + sin_g * pos_teme[1]; + pos_ecef[1] = -sin_g * pos_teme[0] + cos_g * pos_teme[1]; + pos_ecef[2] = pos_teme[2]; + + if (vel_ecef && vel_teme) + { + /* Velocity includes both rotation and Earth angular velocity cross product */ + vel_ecef[0] = cos_g * vel_teme[0] + sin_g * vel_teme[1] + + OMEGA_EARTH * pos_ecef[1]; + vel_ecef[1] = -sin_g * vel_teme[0] + cos_g * vel_teme[1] + - OMEGA_EARTH * pos_ecef[0]; + vel_ecef[2] = vel_teme[2]; + } +} + +/* + * ECEF (km) to WGS-84 geodetic (lat, lon in radians; alt in km). + * Bowring's iterative method -- converges in 2-3 iterations for LEO. + */ +static void +ecef_to_geodetic(const double *pos_ecef, + double *lat, double *lon, double *alt_km) +{ + double x = pos_ecef[0]; + double y = pos_ecef[1]; + double z = pos_ecef[2]; + double a = WGS84_A; + double e2 = WGS84_E2; + double p = sqrt(x * x + y * y); + double phi, N, phi_prev; + int i; + + *lon = atan2(y, x); + + /* Bowring's iterative method */ + phi = atan2(z, p * (1.0 - e2)); /* initial estimate */ + for (i = 0; i < 10; i++) + { + phi_prev = phi; + N = a / sqrt(1.0 - e2 * sin(phi) * sin(phi)); + phi = atan2(z + e2 * N * sin(phi), p); + if (fabs(phi - phi_prev) < 1.0e-12) + break; + } + *lat = phi; + + N = a / sqrt(1.0 - e2 * sin(phi) * sin(phi)); + *alt_km = p / cos(phi) - N; +} + +/* + * Observer (WGS-84 lat/lon radians, alt meters) to ECEF vector in km. + */ +static void +observer_to_ecef(const pg_observer *obs, double *pos_ecef) +{ + double a = WGS84_A; + double e2 = WGS84_E2; + double sin_lat = sin(obs->lat); + double cos_lat = cos(obs->lat); + double N = a / sqrt(1.0 - e2 * sin_lat * sin_lat); + double alt_km = obs->alt_m / 1000.0; + + pos_ecef[0] = (N + alt_km) * cos_lat * cos(obs->lon); + pos_ecef[1] = (N + alt_km) * cos_lat * sin(obs->lon); + pos_ecef[2] = (N * (1.0 - e2) + alt_km) * sin_lat; +} + +/* + * ECEF range vector to topocentric (azimuth, elevation, range). + * Azimuth: 0=N, 90=E, 180=S, 270=W. + * Elevation: 0=horizon, +90=zenith. + */ +static void +ecef_to_topocentric(const double *sat_ecef, const double *obs_ecef, + double obs_lat, double obs_lon, + double *az, double *el, double *range_km) +{ + /* Range vector in ECEF */ + double dx = sat_ecef[0] - obs_ecef[0]; + double dy = sat_ecef[1] - obs_ecef[1]; + double dz = sat_ecef[2] - obs_ecef[2]; + + /* Rotate to topocentric (South, East, Up) */ + double sin_lat = sin(obs_lat); + double cos_lat = cos(obs_lat); + double sin_lon = sin(obs_lon); + double cos_lon = cos(obs_lon); + + double south = sin_lat * cos_lon * dx + sin_lat * sin_lon * dy - cos_lat * dz; + double east = -sin_lon * dx + cos_lon * dy; + double up = cos_lat * cos_lon * dx + cos_lat * sin_lon * dy + sin_lat * dz; + + *range_km = sqrt(dx * dx + dy * dy + dz * dz); + *el = asin(up / *range_km); + *az = atan2(east, -south); /* negative south = north */ + if (*az < 0.0) + *az += 2.0 * M_PI; +} + +/* ================================================================ + * TLE struct conversion + propagation (local copies) + * + * These replicate helpers from sgp4_funcs.c. All .c files link into + * a single .so, so we could call those directly, but keeping them + * static avoids symbol coupling between translation units. + * ================================================================ + */ + +static void +pg_tle_to_sat_code(const pg_tle *src, tle_t *dst) +{ + memset(dst, 0, sizeof(tle_t)); + + dst->epoch = src->epoch; + dst->xincl = src->inclination; + dst->xnodeo = src->raan; + dst->eo = src->eccentricity; + dst->omegao = src->arg_perigee; + dst->xmo = src->mean_anomaly; + dst->xno = src->mean_motion; + dst->xndt2o = src->mean_motion_dot; + dst->xndd6o = src->mean_motion_ddot; + dst->bstar = src->bstar; + + dst->norad_number = src->norad_id; + dst->bulletin_number = src->elset_num; + dst->revolution_number = src->rev_num; + dst->classification = src->classification; + dst->ephemeris_type = src->ephemeris_type; + + memcpy(dst->intl_desig, src->intl_desig, 9); +} + +static int +do_propagate(const pg_tle *tle, double jd, double *pos, double *vel) +{ + tle_t sat; + double *params; + int is_deep; + double tsince; + int err; + + pg_tle_to_sat_code(tle, &sat); + + is_deep = select_ephemeris(&sat); + if (is_deep < 0) + ereport(ERROR, + (errcode(ERRCODE_DATA_EXCEPTION), + errmsg("invalid TLE for NORAD %d: " + "mean motion or eccentricity out of range", + tle->norad_id))); + + params = palloc(sizeof(double) * N_SAT_PARAMS); + + if (is_deep) + SDP4_init(params, &sat); + else + SGP4_init(params, &sat); + + tsince = (jd - sat.epoch) * 1440.0; + + if (is_deep) + err = SDP4(tsince, &sat, params, pos, vel); + else + err = SGP4(tsince, &sat, params, pos, vel); + + pfree(params); + + if (err != 0 && + err != SXPX_WARN_ORBIT_WITHIN_EARTH && + err != SXPX_WARN_PERIGEE_WITHIN_EARTH) + ereport(ERROR, + (errcode(ERRCODE_DATA_EXCEPTION), + errmsg("SGP4/SDP4 propagation failed for NORAD %d " + "at t+%.1f min", + tle->norad_id, tsince))); + + return err; +} + +/* ================================================================ + * Geodetic type I/O + * ================================================================ + */ + +/* + * geodetic_in -- parse text to pg_geodetic + * + * Accepts: (lat_deg,lon_deg,alt_km) + */ +Datum +geodetic_in(PG_FUNCTION_ARGS) +{ + char *str = PG_GETARG_CSTRING(0); + pg_geodetic *result; + double lat_deg, lon_deg, alt_km; + int nfields; + + result = (pg_geodetic *) palloc(sizeof(pg_geodetic)); + + nfields = sscanf(str, " ( %lf , %lf , %lf )", + &lat_deg, &lon_deg, &alt_km); + + if (nfields != 3) + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type geodetic: \"%s\"", str), + errhint("Expected (lat_deg,lon_deg,alt_km)."))); + + if (lat_deg < -90.0 || lat_deg > 90.0) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("latitude out of range: %.6f", lat_deg), + errhint("Latitude must be between -90 and +90 degrees."))); + + if (lon_deg < -180.0 || lon_deg > 360.0) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("longitude out of range: %.6f", lon_deg), + errhint("Longitude must be between -180 and +360 degrees."))); + + result->lat = lat_deg * DEG_TO_RAD; + result->lon = lon_deg * DEG_TO_RAD; + result->alt = alt_km; + + PG_RETURN_POINTER(result); +} + +/* + * geodetic_out -- pg_geodetic to text + */ +Datum +geodetic_out(PG_FUNCTION_ARGS) +{ + pg_geodetic *geo = (pg_geodetic *) PG_GETARG_POINTER(0); + + PG_RETURN_CSTRING(psprintf("(%.6f,%.6f,%.3f)", + geo->lat * RAD_TO_DEG, + geo->lon * RAD_TO_DEG, + geo->alt)); +} + +/* + * geodetic_recv -- binary input + */ +Datum +geodetic_recv(PG_FUNCTION_ARGS) +{ + StringInfo buf = (StringInfo) PG_GETARG_POINTER(0); + pg_geodetic *result; + + result = (pg_geodetic *) palloc(sizeof(pg_geodetic)); + result->lat = pq_getmsgfloat8(buf); + result->lon = pq_getmsgfloat8(buf); + result->alt = pq_getmsgfloat8(buf); + + PG_RETURN_POINTER(result); +} + +/* + * geodetic_send -- binary output + */ +Datum +geodetic_send(PG_FUNCTION_ARGS) +{ + pg_geodetic *geo = (pg_geodetic *) PG_GETARG_POINTER(0); + StringInfoData buf; + + pq_begintypsend(&buf); + pq_sendfloat8(&buf, geo->lat); + pq_sendfloat8(&buf, geo->lon); + pq_sendfloat8(&buf, geo->alt); + + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); +} + +/* --- geodetic accessor functions --- */ + +Datum +geodetic_lat(PG_FUNCTION_ARGS) +{ + pg_geodetic *geo = (pg_geodetic *) PG_GETARG_POINTER(0); + PG_RETURN_FLOAT8(geo->lat * RAD_TO_DEG); +} + +Datum +geodetic_lon(PG_FUNCTION_ARGS) +{ + pg_geodetic *geo = (pg_geodetic *) PG_GETARG_POINTER(0); + PG_RETURN_FLOAT8(geo->lon * RAD_TO_DEG); +} + +Datum +geodetic_alt(PG_FUNCTION_ARGS) +{ + pg_geodetic *geo = (pg_geodetic *) PG_GETARG_POINTER(0); + PG_RETURN_FLOAT8(geo->alt); +} + +/* ================================================================ + * Topocentric type I/O + * ================================================================ + */ + +/* + * topocentric_in -- parse text to pg_topocentric + * + * Accepts: (az_deg,el_deg,range_km,range_rate_kms) + */ +Datum +topocentric_in(PG_FUNCTION_ARGS) +{ + char *str = PG_GETARG_CSTRING(0); + pg_topocentric *result; + double az_deg, el_deg, range_km, range_rate; + int nfields; + + result = (pg_topocentric *) palloc(sizeof(pg_topocentric)); + + nfields = sscanf(str, " ( %lf , %lf , %lf , %lf )", + &az_deg, &el_deg, &range_km, &range_rate); + + if (nfields != 4) + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type topocentric: \"%s\"", str), + errhint("Expected (az_deg,el_deg,range_km,range_rate_kms)."))); + + if (az_deg < 0.0 || az_deg > 360.0) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("azimuth out of range: %.6f", az_deg), + errhint("Azimuth must be between 0 and 360 degrees."))); + + if (el_deg < -90.0 || el_deg > 90.0) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("elevation out of range: %.6f", el_deg), + errhint("Elevation must be between -90 and +90 degrees."))); + + result->azimuth = az_deg * DEG_TO_RAD; + result->elevation = el_deg * DEG_TO_RAD; + result->range_km = range_km; + result->range_rate = range_rate; + + PG_RETURN_POINTER(result); +} + +/* + * topocentric_out -- pg_topocentric to text + */ +Datum +topocentric_out(PG_FUNCTION_ARGS) +{ + pg_topocentric *topo = (pg_topocentric *) PG_GETARG_POINTER(0); + + PG_RETURN_CSTRING(psprintf("(%.4f,%.4f,%.3f,%.6f)", + topo->azimuth * RAD_TO_DEG, + topo->elevation * RAD_TO_DEG, + topo->range_km, + topo->range_rate)); +} + +/* + * topocentric_recv -- binary input + */ +Datum +topocentric_recv(PG_FUNCTION_ARGS) +{ + StringInfo buf = (StringInfo) PG_GETARG_POINTER(0); + pg_topocentric *result; + + result = (pg_topocentric *) palloc(sizeof(pg_topocentric)); + result->azimuth = pq_getmsgfloat8(buf); + result->elevation = pq_getmsgfloat8(buf); + result->range_km = pq_getmsgfloat8(buf); + result->range_rate = pq_getmsgfloat8(buf); + + PG_RETURN_POINTER(result); +} + +/* + * topocentric_send -- binary output + */ +Datum +topocentric_send(PG_FUNCTION_ARGS) +{ + pg_topocentric *topo = (pg_topocentric *) PG_GETARG_POINTER(0); + StringInfoData buf; + + pq_begintypsend(&buf); + pq_sendfloat8(&buf, topo->azimuth); + pq_sendfloat8(&buf, topo->elevation); + pq_sendfloat8(&buf, topo->range_km); + pq_sendfloat8(&buf, topo->range_rate); + + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); +} + +/* --- topocentric accessor functions --- */ + +Datum +topo_azimuth(PG_FUNCTION_ARGS) +{ + pg_topocentric *topo = (pg_topocentric *) PG_GETARG_POINTER(0); + PG_RETURN_FLOAT8(topo->azimuth * RAD_TO_DEG); +} + +Datum +topo_elevation(PG_FUNCTION_ARGS) +{ + pg_topocentric *topo = (pg_topocentric *) PG_GETARG_POINTER(0); + PG_RETURN_FLOAT8(topo->elevation * RAD_TO_DEG); +} + +Datum +topo_range(PG_FUNCTION_ARGS) +{ + pg_topocentric *topo = (pg_topocentric *) PG_GETARG_POINTER(0); + PG_RETURN_FLOAT8(topo->range_km); +} + +Datum +topo_range_rate(PG_FUNCTION_ARGS) +{ + pg_topocentric *topo = (pg_topocentric *) PG_GETARG_POINTER(0); + PG_RETURN_FLOAT8(topo->range_rate); +} + +/* ================================================================ + * Coordinate transform functions + * ================================================================ + */ + +/* + * eci_to_geodetic(eci_position, timestamptz) -> geodetic + * + * Convert TEME state vector to WGS-84 geodetic coordinates. + * The timestamptz provides the time for GMST calculation. + */ +Datum +eci_to_geodetic(PG_FUNCTION_ARGS) +{ + pg_eci *eci = (pg_eci *) PG_GETARG_POINTER(0); + int64 ts = PG_GETARG_INT64(1); + + double jd; + double gmst; + double pos_teme[3]; + double pos_ecef[3]; + double lat, lon, alt_km; + pg_geodetic *result; + + jd = timestamptz_to_jd(ts); + gmst = gmst_from_jd(jd); + + pos_teme[0] = eci->x; + pos_teme[1] = eci->y; + pos_teme[2] = eci->z; + + teme_to_ecef(pos_teme, NULL, gmst, pos_ecef, NULL); + ecef_to_geodetic(pos_ecef, &lat, &lon, &alt_km); + + result = (pg_geodetic *) palloc(sizeof(pg_geodetic)); + result->lat = lat; + result->lon = lon; + result->alt = alt_km; + + PG_RETURN_POINTER(result); +} + +/* + * eci_to_topocentric(eci_position, observer, timestamptz) -> topocentric + * + * Convert TEME state vector to observer-relative look angles. + * Range rate is computed by projecting the ECEF-relative velocity + * onto the line-of-sight unit vector. + */ +Datum +eci_to_topocentric(PG_FUNCTION_ARGS) +{ + pg_eci *eci = (pg_eci *) PG_GETARG_POINTER(0); + pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(1); + int64 ts = PG_GETARG_INT64(2); + + double jd; + double gmst; + double pos_teme[3], vel_teme[3]; + double pos_ecef[3], vel_ecef[3]; + double obs_ecef[3]; + double az, el, range_km; + double range_rate; + double dx, dy, dz; + double dvx, dvy, dvz; + pg_topocentric *result; + + jd = timestamptz_to_jd(ts); + gmst = gmst_from_jd(jd); + + pos_teme[0] = eci->x; + pos_teme[1] = eci->y; + pos_teme[2] = eci->z; + + /* ECI velocity is in km/s; teme_to_ecef expects km/s */ + vel_teme[0] = eci->vx; + vel_teme[1] = eci->vy; + vel_teme[2] = eci->vz; + + teme_to_ecef(pos_teme, vel_teme, gmst, pos_ecef, vel_ecef); + observer_to_ecef(obs, obs_ecef); + + ecef_to_topocentric(pos_ecef, obs_ecef, obs->lat, obs->lon, + &az, &el, &range_km); + + /* + * Range rate: project relative velocity onto line-of-sight. + * Positive = satellite receding from observer. + * Observer is stationary in ECEF, so relative velocity = sat ECEF velocity. + */ + dx = pos_ecef[0] - obs_ecef[0]; + dy = pos_ecef[1] - obs_ecef[1]; + dz = pos_ecef[2] - obs_ecef[2]; + + dvx = vel_ecef[0]; + dvy = vel_ecef[1]; + dvz = vel_ecef[2]; + + range_rate = (dx * dvx + dy * dvy + dz * dvz) / range_km; + + result = (pg_topocentric *) palloc(sizeof(pg_topocentric)); + result->azimuth = az; + result->elevation = el; + result->range_km = range_km; + result->range_rate = range_rate; + + PG_RETURN_POINTER(result); +} + +/* + * subsatellite_point(tle, timestamptz) -> geodetic + * + * Propagate TLE to the given time and return the nadir point + * on the WGS-84 ellipsoid. The altitude field contains the + * satellite's altitude above the ellipsoid, not zero. + */ +Datum +subsatellite_point(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + int64 ts = PG_GETARG_INT64(1); + + double jd; + double gmst; + double pos[3], vel[3]; + double pos_ecef[3]; + double lat, lon, alt_km; + pg_geodetic *result; + + jd = timestamptz_to_jd(ts); + + do_propagate(tle, jd, pos, vel); + + gmst = gmst_from_jd(jd); + teme_to_ecef(pos, NULL, gmst, pos_ecef, NULL); + ecef_to_geodetic(pos_ecef, &lat, &lon, &alt_km); + + result = (pg_geodetic *) palloc(sizeof(pg_geodetic)); + result->lat = lat; + result->lon = lon; + result->alt = alt_km; + + PG_RETURN_POINTER(result); +} + +/* ================================================================ + * ground_track -- SRF returning subsatellite points over time + * + * ground_track(tle, start_ts, stop_ts, step_interval) + * -> SETOF (t timestamptz, point geodetic) + * + * The model is initialized once and reused for every step, + * same pattern as sgp4_propagate_series in sgp4_funcs.c. + * ================================================================ + */ + +typedef struct +{ + tle_t sat; + double params[N_SAT_PARAMS]; + int is_deep; + double epoch_jd; + int64 start_ts; + int64 step_usec; +} ground_track_ctx; + +Datum +ground_track(PG_FUNCTION_ARGS) +{ + FuncCallContext *funcctx; + ground_track_ctx *ctx; + + if (SRF_IS_FIRSTCALL()) + { + MemoryContext oldctx; + pg_tle *tle; + int64 start_ts; + int64 stop_ts; + Interval *step; + int64 step_usec; + int64 span; + uint64 nsteps; + TupleDesc tupdesc; + + funcctx = SRF_FIRSTCALL_INIT(); + oldctx = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx); + + tle = (pg_tle *) PG_GETARG_POINTER(0); + start_ts = PG_GETARG_INT64(1); + stop_ts = PG_GETARG_INT64(2); + step = PG_GETARG_INTERVAL_P(3); + + /* Convert interval to microseconds */ + step_usec = step->time + + (int64) step->day * USECS_PER_DAY + + (int64) step->month * (30 * USECS_PER_DAY); + + if (step_usec <= 0) + ereport(ERROR, + (errcode(ERRCODE_INVALID_PARAMETER_VALUE), + errmsg("step interval must be positive"))); + + if (stop_ts < start_ts) + ereport(ERROR, + (errcode(ERRCODE_INVALID_PARAMETER_VALUE), + errmsg("stop time must be >= start time"))); + + span = stop_ts - start_ts; + nsteps = (uint64)(span / step_usec) + 1; + + ctx = (ground_track_ctx *) + palloc0(sizeof(ground_track_ctx)); + + pg_tle_to_sat_code(tle, &ctx->sat); + + ctx->is_deep = select_ephemeris(&ctx->sat); + if (ctx->is_deep < 0) + ereport(ERROR, + (errcode(ERRCODE_DATA_EXCEPTION), + errmsg("invalid TLE for NORAD %d: " + "mean motion or eccentricity out of range", + ctx->sat.norad_number))); + + /* Initialize the model once */ + if (ctx->is_deep) + SDP4_init(ctx->params, &ctx->sat); + else + SGP4_init(ctx->params, &ctx->sat); + + ctx->epoch_jd = ctx->sat.epoch; + ctx->start_ts = start_ts; + ctx->step_usec = step_usec; + + funcctx->max_calls = nsteps; + funcctx->user_fctx = ctx; + + /* Output tuple: (timestamptz, lat_deg, lon_deg, alt_km) */ + tupdesc = CreateTemplateTupleDesc(4); + TupleDescInitEntry(tupdesc, 1, "t", TIMESTAMPTZOID, -1, 0); + TupleDescInitEntry(tupdesc, 2, "lat_deg", FLOAT8OID, -1, 0); + TupleDescInitEntry(tupdesc, 3, "lon_deg", FLOAT8OID, -1, 0); + TupleDescInitEntry(tupdesc, 4, "alt_km", FLOAT8OID, -1, 0); + + funcctx->tuple_desc = BlessTupleDesc(tupdesc); + + MemoryContextSwitchTo(oldctx); + } + + funcctx = SRF_PERCALL_SETUP(); + ctx = (ground_track_ctx *) funcctx->user_fctx; + + if (funcctx->call_cntr < funcctx->max_calls) + { + int64 ts; + double jd; + double tsince; + double pos[3], vel[3]; + double gmst; + double pos_ecef[3]; + double lat, lon, alt_km; + int err; + Datum values[4]; + bool nulls[4]; + HeapTuple tuple; + + ts = ctx->start_ts + (int64) funcctx->call_cntr * ctx->step_usec; + jd = timestamptz_to_jd(ts); + tsince = jd_to_minutes_since_epoch(jd, ctx->epoch_jd); + + if (ctx->is_deep) + err = SDP4(tsince, &ctx->sat, ctx->params, pos, vel); + else + err = SGP4(tsince, &ctx->sat, ctx->params, pos, vel); + + if (err != 0 && + err != SXPX_WARN_ORBIT_WITHIN_EARTH && + err != SXPX_WARN_PERIGEE_WITHIN_EARTH) + ereport(ERROR, + (errcode(ERRCODE_DATA_EXCEPTION), + errmsg("SGP4/SDP4 propagation failed for NORAD %d " + "at t+%.1f min", + ctx->sat.norad_number, tsince))); + + gmst = gmst_from_jd(jd); + teme_to_ecef(pos, NULL, gmst, pos_ecef, NULL); + ecef_to_geodetic(pos_ecef, &lat, &lon, &alt_km); + + memset(nulls, 0, sizeof(nulls)); + + values[0] = Int64GetDatum(ts); + values[1] = Float8GetDatum(lat * RAD_TO_DEG); + values[2] = Float8GetDatum(lon * RAD_TO_DEG); + values[3] = Float8GetDatum(alt_km); + + tuple = heap_form_tuple(funcctx->tuple_desc, values, nulls); + + SRF_RETURN_NEXT(funcctx, HeapTupleGetDatum(tuple)); + } + + SRF_RETURN_DONE(funcctx); +} diff --git a/src/eci_type.c b/src/eci_type.c new file mode 100644 index 0000000..746a0f0 --- /dev/null +++ b/src/eci_type.c @@ -0,0 +1,219 @@ +/* + * eci_type.c -- ECI position type (TEME frame) + * + * Position in km, velocity in km/s. + * Input/output as (x,y,z,vx,vy,vz) or (x,y,z) with zero velocity. + */ + +#include "postgres.h" +#include "fmgr.h" +#include "utils/builtins.h" +#include "libpq/pqformat.h" +#include "types.h" +#include +#include +#include + +PG_FUNCTION_INFO_V1(eci_in); +PG_FUNCTION_INFO_V1(eci_out); +PG_FUNCTION_INFO_V1(eci_recv); +PG_FUNCTION_INFO_V1(eci_send); +PG_FUNCTION_INFO_V1(eci_x); +PG_FUNCTION_INFO_V1(eci_y); +PG_FUNCTION_INFO_V1(eci_z); +PG_FUNCTION_INFO_V1(eci_vx); +PG_FUNCTION_INFO_V1(eci_vy); +PG_FUNCTION_INFO_V1(eci_vz); +PG_FUNCTION_INFO_V1(eci_speed); +PG_FUNCTION_INFO_V1(eci_altitude); + + +/* + * eci_in -- parse text to pg_eci + * + * Accepts: + * (x,y,z,vx,vy,vz) full state vector + * (x,y,z) position only, velocity = 0 + */ +Datum +eci_in(PG_FUNCTION_ARGS) +{ + char *str = PG_GETARG_CSTRING(0); + pg_eci *result; + double x, y, z, vx, vy, vz; + int nfields; + + result = (pg_eci *) palloc(sizeof(pg_eci)); + + nfields = sscanf(str, " ( %lf , %lf , %lf , %lf , %lf , %lf )", + &x, &y, &z, &vx, &vy, &vz); + + if (nfields == 6) + { + result->x = x; + result->y = y; + result->z = z; + result->vx = vx; + result->vy = vy; + result->vz = vz; + PG_RETURN_POINTER(result); + } + + /* Try position-only format */ + nfields = sscanf(str, " ( %lf , %lf , %lf )", &x, &y, &z); + if (nfields == 3) + { + result->x = x; + result->y = y; + result->z = z; + result->vx = 0.0; + result->vy = 0.0; + result->vz = 0.0; + PG_RETURN_POINTER(result); + } + + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type eci: \"%s\"", str), + errhint("Expected (x,y,z,vx,vy,vz) or (x,y,z)."))); + + PG_RETURN_NULL(); /* not reached */ +} + + +/* + * eci_out -- pg_eci to text + */ +Datum +eci_out(PG_FUNCTION_ARGS) +{ + pg_eci *eci = (pg_eci *) PG_GETARG_POINTER(0); + + PG_RETURN_CSTRING(psprintf("(%.6f,%.6f,%.6f,%.6f,%.6f,%.6f)", + eci->x, eci->y, eci->z, + eci->vx, eci->vy, eci->vz)); +} + + +/* + * eci_recv -- binary input + */ +Datum +eci_recv(PG_FUNCTION_ARGS) +{ + StringInfo buf = (StringInfo) PG_GETARG_POINTER(0); + pg_eci *result; + + result = (pg_eci *) palloc(sizeof(pg_eci)); + result->x = pq_getmsgfloat8(buf); + result->y = pq_getmsgfloat8(buf); + result->z = pq_getmsgfloat8(buf); + result->vx = pq_getmsgfloat8(buf); + result->vy = pq_getmsgfloat8(buf); + result->vz = pq_getmsgfloat8(buf); + + PG_RETURN_POINTER(result); +} + + +/* + * eci_send -- binary output + */ +Datum +eci_send(PG_FUNCTION_ARGS) +{ + pg_eci *eci = (pg_eci *) PG_GETARG_POINTER(0); + StringInfoData buf; + + pq_begintypsend(&buf); + pq_sendfloat8(&buf, eci->x); + pq_sendfloat8(&buf, eci->y); + pq_sendfloat8(&buf, eci->z); + pq_sendfloat8(&buf, eci->vx); + pq_sendfloat8(&buf, eci->vy); + pq_sendfloat8(&buf, eci->vz); + + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); +} + + +/* --- accessor functions --- */ + +Datum +eci_x(PG_FUNCTION_ARGS) +{ + pg_eci *eci = (pg_eci *) PG_GETARG_POINTER(0); + PG_RETURN_FLOAT8(eci->x); +} + +Datum +eci_y(PG_FUNCTION_ARGS) +{ + pg_eci *eci = (pg_eci *) PG_GETARG_POINTER(0); + PG_RETURN_FLOAT8(eci->y); +} + +Datum +eci_z(PG_FUNCTION_ARGS) +{ + pg_eci *eci = (pg_eci *) PG_GETARG_POINTER(0); + PG_RETURN_FLOAT8(eci->z); +} + +Datum +eci_vx(PG_FUNCTION_ARGS) +{ + pg_eci *eci = (pg_eci *) PG_GETARG_POINTER(0); + PG_RETURN_FLOAT8(eci->vx); +} + +Datum +eci_vy(PG_FUNCTION_ARGS) +{ + pg_eci *eci = (pg_eci *) PG_GETARG_POINTER(0); + PG_RETURN_FLOAT8(eci->vy); +} + +Datum +eci_vz(PG_FUNCTION_ARGS) +{ + pg_eci *eci = (pg_eci *) PG_GETARG_POINTER(0); + PG_RETURN_FLOAT8(eci->vz); +} + + +/* + * eci_speed -- orbital speed in km/s + */ +Datum +eci_speed(PG_FUNCTION_ARGS) +{ + pg_eci *eci = (pg_eci *) PG_GETARG_POINTER(0); + double spd; + + spd = sqrt(eci->vx * eci->vx + + eci->vy * eci->vy + + eci->vz * eci->vz); + + PG_RETURN_FLOAT8(spd); +} + + +/* + * eci_altitude -- approximate altitude above WGS-72 sphere in km + * + * This is geocentric radius minus equatorial radius, not geodetic altitude. + * Good enough for quick filtering; use eci_to_geodetic() for precision. + */ +Datum +eci_altitude(PG_FUNCTION_ARGS) +{ + pg_eci *eci = (pg_eci *) PG_GETARG_POINTER(0); + double r; + + r = sqrt(eci->x * eci->x + + eci->y * eci->y + + eci->z * eci->z); + + PG_RETURN_FLOAT8(r - WGS72_AE); +} diff --git a/src/gist_tle.c b/src/gist_tle.c new file mode 100644 index 0000000..d9ef94a --- /dev/null +++ b/src/gist_tle.c @@ -0,0 +1,505 @@ +/* + * gist_tle.c -- GiST operator class for altitude-band indexing on TLE + * + * Every TLE defines an orbit whose perigee and apogee altitudes form a + * 1-D range (the "altitude band"). Two orbits can only be in proximity + * if their altitude bands overlap -- a necessary but not sufficient + * condition for conjunction. + * + * The GiST index stores [perigee_km, apogee_km] ranges as internal + * keys, enabling fast coarse filtering. The && (overlap) operator is + * always rechecked: real conjunction screening requires propagation. + * + * Semi-major axis from Kepler's third law using WGS-72 constants: + * a = (KE / n)^(2/3) [earth radii] + * perigee_km = a*(1-e)*AE - AE + * apogee_km = a*(1+e)*AE - AE + */ + +#include "postgres.h" +#include "fmgr.h" +#include "access/gist.h" +#include "access/stratnum.h" +#include "utils/float.h" +#include "norad.h" +#include "types.h" +#include +#include + +PG_FUNCTION_INFO_V1(tle_overlap); +PG_FUNCTION_INFO_V1(tle_alt_distance); +PG_FUNCTION_INFO_V1(gist_tle_compress); +PG_FUNCTION_INFO_V1(gist_tle_decompress); +PG_FUNCTION_INFO_V1(gist_tle_consistent); +PG_FUNCTION_INFO_V1(gist_tle_union); +PG_FUNCTION_INFO_V1(gist_tle_penalty); +PG_FUNCTION_INFO_V1(gist_tle_picksplit); +PG_FUNCTION_INFO_V1(gist_tle_same); +PG_FUNCTION_INFO_V1(gist_tle_distance); + +/* Floating-point comparison tolerance (km) */ +#define ALT_EPSILON 1.0e-9 + +/* + * Altitude band extracted from a TLE's mean elements. + * This is the GiST internal key -- much cheaper to compare + * than propagating two full state vectors. + */ +typedef struct tle_alt_range +{ + double low; /* perigee altitude, km */ + double high; /* apogee altitude, km */ +} tle_alt_range; + + +/* ---------------------------------------------------------------- + * tle_to_alt_range -- compute [perigee, apogee] from mean elements + * + * Uses WGS-72 KE and AE (the only constants valid for SGP4 elements). + * Degenerate TLEs with n <= 0 map to a zero-width range at 0 km. + * ---------------------------------------------------------------- + */ +static void +tle_to_alt_range(const pg_tle *tle, tle_alt_range *range) +{ + double n = tle->mean_motion; /* rad/min */ + double e = tle->eccentricity; + double a_er; /* semi-major axis, earth radii */ + + if (n <= 0.0) + { + range->low = 0.0; + range->high = 0.0; + return; + } + + a_er = pow(WGS72_KE / n, 2.0 / 3.0); + range->low = a_er * (1.0 - e) * WGS72_AE - WGS72_AE; + range->high = a_er * (1.0 + e) * WGS72_AE - WGS72_AE; + + /* Guard against numerical inversion from near-zero eccentricity */ + if (range->low > range->high) + { + double tmp = range->low; + range->low = range->high; + range->high = tmp; + } +} + + +/* ---------------------------------------------------------------- + * range_overlaps -- do two altitude bands share any interval? + * ---------------------------------------------------------------- + */ +static inline bool +range_overlaps(const tle_alt_range *a, const tle_alt_range *b) +{ + return (a->low <= b->high) && (b->low <= a->high); +} + + +/* ---------------------------------------------------------------- + * range_contains -- does outer fully contain inner? + * ---------------------------------------------------------------- + */ +static inline bool +range_contains(const tle_alt_range *outer, const tle_alt_range *inner) +{ + return (outer->low <= inner->low) && (inner->high <= outer->high); +} + + +/* ---------------------------------------------------------------- + * range_merge -- expand dst to encompass src + * ---------------------------------------------------------------- + */ +static inline void +range_merge(tle_alt_range *dst, const tle_alt_range *src) +{ + if (src->low < dst->low) + dst->low = src->low; + if (src->high > dst->high) + dst->high = src->high; +} + + +/* ---------------------------------------------------------------- + * range_separation -- minimum gap between two non-overlapping ranges + * + * Returns 0 if the ranges overlap. + * ---------------------------------------------------------------- + */ +static inline double +range_separation(const tle_alt_range *a, const tle_alt_range *b) +{ + if (a->high < b->low) + return b->low - a->high; + if (b->high < a->low) + return a->low - b->high; + return 0.0; +} + + +/* ================================================================ + * SQL-callable operators + * ================================================================ + */ + + +/* + * tle_overlap(tle, tle) -> bool [the && operator] + * + * True if the altitude bands of two TLEs share any interval. + * This is a fast pre-filter: overlapping bands are necessary + * but not sufficient for actual conjunction. + */ +Datum +tle_overlap(PG_FUNCTION_ARGS) +{ + pg_tle *a = (pg_tle *) PG_GETARG_POINTER(0); + pg_tle *b = (pg_tle *) PG_GETARG_POINTER(1); + tle_alt_range ra, rb; + + tle_to_alt_range(a, &ra); + tle_to_alt_range(b, &rb); + + PG_RETURN_BOOL(range_overlaps(&ra, &rb)); +} + + +/* + * tle_alt_distance(tle, tle) -> float8 [the <-> operator] + * + * Minimum altitude-band separation in km. Returns 0 if the bands + * overlap. This is not the physical distance between the objects -- + * it is the gap between their orbital shells, useful for ordering + * nearest-neighbor queries without propagation. + */ +Datum +tle_alt_distance(PG_FUNCTION_ARGS) +{ + pg_tle *a = (pg_tle *) PG_GETARG_POINTER(0); + pg_tle *b = (pg_tle *) PG_GETARG_POINTER(1); + tle_alt_range ra, rb; + + tle_to_alt_range(a, &ra); + tle_to_alt_range(b, &rb); + + PG_RETURN_FLOAT8(range_separation(&ra, &rb)); +} + + +/* ================================================================ + * GiST support functions + * ================================================================ + */ + + +/* + * gist_tle_compress -- extract altitude range from a leaf TLE + * + * Leaf entries carry the full pg_tle; we compress to tle_alt_range. + * Internal entries are already tle_alt_range from union operations. + */ +Datum +gist_tle_compress(PG_FUNCTION_ARGS) +{ + GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); + GISTENTRY *retval; + + if (entry->leafkey) + { + pg_tle *tle = (pg_tle *) DatumGetPointer(entry->key); + tle_alt_range *range = (tle_alt_range *) palloc(sizeof(tle_alt_range)); + + tle_to_alt_range(tle, range); + + retval = (GISTENTRY *) palloc(sizeof(GISTENTRY)); + gistentryinit(*retval, PointerGetDatum(range), + entry->rel, entry->page, entry->offset, false); + } + else + { + /* Internal node: already a tle_alt_range */ + retval = entry; + } + + PG_RETURN_POINTER(retval); +} + + +/* + * gist_tle_decompress -- identity (we operate on compressed keys directly) + */ +Datum +gist_tle_decompress(PG_FUNCTION_ARGS) +{ + PG_RETURN_POINTER(PG_GETARG_POINTER(0)); +} + + +/* + * gist_tle_consistent -- can this subtree contain matches for the query? + * + * Strategy RTOverlapStrategyNumber (&&): altitude bands must overlap. + * Always sets recheck = true because altitude overlap is only a necessary + * condition -- the real conjunction test requires propagation. + */ +Datum +gist_tle_consistent(PG_FUNCTION_ARGS) +{ + GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); + pg_tle *query = (pg_tle *) PG_GETARG_POINTER(1); + StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2); + /* arg 3 is the query type OID, unused */ + bool *recheck = (bool *) PG_GETARG_POINTER(4); + tle_alt_range *key = (tle_alt_range *) DatumGetPointer(entry->key); + tle_alt_range query_range; + bool result; + + tle_to_alt_range(query, &query_range); + + /* + * Altitude overlap is necessary, not sufficient. + * The actual operator (if exact) would need propagation, so always recheck. + */ + *recheck = true; + + switch (strategy) + { + case RTOverlapStrategyNumber: /* && */ + result = range_overlaps(key, &query_range); + break; + + case RTContainedByStrategyNumber: /* <@ */ + if (GIST_LEAF(entry)) + result = range_contains(&query_range, key); + else + result = range_overlaps(key, &query_range); + break; + + case RTContainsStrategyNumber: /* @> */ + if (GIST_LEAF(entry)) + result = range_contains(key, &query_range); + else + result = range_overlaps(key, &query_range); + break; + + default: + elog(ERROR, "gist_tle_consistent: unrecognized strategy %d", + strategy); + result = false; + break; + } + + PG_RETURN_BOOL(result); +} + + +/* + * gist_tle_union -- compute bounding altitude range for a set of entries + * + * The union of N altitude ranges is simply [min(low), max(high)]. + */ +Datum +gist_tle_union(PG_FUNCTION_ARGS) +{ + GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0); + int *sizep = (int *) PG_GETARG_POINTER(1); + int i; + tle_alt_range *result; + tle_alt_range *cur; + + result = (tle_alt_range *) palloc(sizeof(tle_alt_range)); + cur = (tle_alt_range *) DatumGetPointer(entryvec->vector[0].key); + result->low = cur->low; + result->high = cur->high; + + for (i = 1; i < entryvec->n; i++) + { + cur = (tle_alt_range *) DatumGetPointer(entryvec->vector[i].key); + range_merge(result, cur); + } + + *sizep = sizeof(tle_alt_range); + + PG_RETURN_POINTER(result); +} + + +/* + * gist_tle_penalty -- cost of inserting a new entry into an existing subtree + * + * Penalty is the increase in altitude span (km) that results from + * expanding the subtree's bounding range to include the new entry. + * A good penalty function keeps the tree balanced and minimizes + * unnecessary range expansion. + */ +Datum +gist_tle_penalty(PG_FUNCTION_ARGS) +{ + GISTENTRY *origentry = (GISTENTRY *) PG_GETARG_POINTER(0); + GISTENTRY *newentry = (GISTENTRY *) PG_GETARG_POINTER(1); + float *penalty = (float *) PG_GETARG_POINTER(2); + tle_alt_range *orig = (tle_alt_range *) DatumGetPointer(origentry->key); + tle_alt_range *add = (tle_alt_range *) DatumGetPointer(newentry->key); + double orig_span; + double merged_span; + + orig_span = orig->high - orig->low; + merged_span = fmax(orig->high, add->high) - fmin(orig->low, add->low); + + *penalty = (float)(merged_span - orig_span); + + PG_RETURN_POINTER(penalty); +} + + +/* + * Comparison callback for qsort in picksplit. + * Sorts entries by the midpoint of their altitude range. + */ +typedef struct +{ + int index; /* position in the original entry vector */ + double midpoint; /* (low + high) / 2 */ +} picksplit_item; + +static int +picksplit_cmp(const void *a, const void *b) +{ + double ma = ((const picksplit_item *) a)->midpoint; + double mb = ((const picksplit_item *) b)->midpoint; + + if (ma < mb) + return -1; + if (ma > mb) + return 1; + return 0; +} + + +/* + * gist_tle_picksplit -- split an overfull page into two groups + * + * Strategy: sort entries by altitude midpoint, split at the median. + * This keeps nearby altitude bands in the same subtree, which is + * optimal for a 1-D range index. + */ +Datum +gist_tle_picksplit(PG_FUNCTION_ARGS) +{ + GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0); + GIST_SPLITVEC *splitvec = (GIST_SPLITVEC *) PG_GETARG_POINTER(1); + int nentries = entryvec->n; + picksplit_item *items; + tle_alt_range *left_union; + tle_alt_range *right_union; + tle_alt_range *cur; + int split_at; + int i; + + items = (picksplit_item *) palloc(sizeof(picksplit_item) * nentries); + for (i = 0; i < nentries; i++) + { + cur = (tle_alt_range *) DatumGetPointer(entryvec->vector[i].key); + items[i].index = i; + items[i].midpoint = (cur->low + cur->high) / 2.0; + } + + qsort(items, nentries, sizeof(picksplit_item), picksplit_cmp); + + split_at = nentries / 2; + + /* Allocate offset arrays (GiST uses OffsetNumber, 1-based) */ + splitvec->spl_left = (OffsetNumber *) palloc(sizeof(OffsetNumber) * nentries); + splitvec->spl_right = (OffsetNumber *) palloc(sizeof(OffsetNumber) * nentries); + splitvec->spl_nleft = 0; + splitvec->spl_nright = 0; + + /* Compute union ranges and assign entries */ + left_union = (tle_alt_range *) palloc(sizeof(tle_alt_range)); + right_union = (tle_alt_range *) palloc(sizeof(tle_alt_range)); + + /* Seed the unions from the first entry in each half */ + cur = (tle_alt_range *) DatumGetPointer( + entryvec->vector[items[0].index].key); + left_union->low = cur->low; + left_union->high = cur->high; + + cur = (tle_alt_range *) DatumGetPointer( + entryvec->vector[items[split_at].index].key); + right_union->low = cur->low; + right_union->high = cur->high; + + for (i = 0; i < nentries; i++) + { + int idx = items[i].index; + + cur = (tle_alt_range *) DatumGetPointer( + entryvec->vector[idx].key); + + if (i < split_at) + { + splitvec->spl_left[splitvec->spl_nleft++] = + (OffsetNumber)(idx + 1); /* 1-based */ + range_merge(left_union, cur); + } + else + { + splitvec->spl_right[splitvec->spl_nright++] = + (OffsetNumber)(idx + 1); + range_merge(right_union, cur); + } + } + + splitvec->spl_ldatum = PointerGetDatum(left_union); + splitvec->spl_rdatum = PointerGetDatum(right_union); + + pfree(items); + + PG_RETURN_POINTER(splitvec); +} + + +/* + * gist_tle_same -- equality test on compressed keys + * + * Two altitude ranges are "same" if both endpoints match within + * a small tolerance. This lets the GiST machinery detect + * duplicate keys. + */ +Datum +gist_tle_same(PG_FUNCTION_ARGS) +{ + tle_alt_range *a = (tle_alt_range *) PG_GETARG_POINTER(0); + tle_alt_range *b = (tle_alt_range *) PG_GETARG_POINTER(1); + bool *result = (bool *) PG_GETARG_POINTER(2); + + *result = (fabs(a->low - b->low) < ALT_EPSILON && + fabs(a->high - b->high) < ALT_EPSILON); + + PG_RETURN_POINTER(result); +} + + +/* + * gist_tle_distance -- GiST distance function for KNN ordering + * + * Returns the minimum altitude-band separation in km. + * For overlapping ranges this is 0, making the entry a candidate. + * The planner uses this to drive ORDER BY <-> queries. + */ +Datum +gist_tle_distance(PG_FUNCTION_ARGS) +{ + GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); + pg_tle *query = (pg_tle *) PG_GETARG_POINTER(1); + /* strategy number at arg 2, unused for single-distance class */ + tle_alt_range *key = (tle_alt_range *) DatumGetPointer(entry->key); + tle_alt_range query_range; + + tle_to_alt_range(query, &query_range); + + PG_RETURN_FLOAT8(range_separation(key, &query_range)); +} diff --git a/src/observer_type.c b/src/observer_type.c new file mode 100644 index 0000000..6c22be6 --- /dev/null +++ b/src/observer_type.c @@ -0,0 +1,284 @@ +/* + * observer_type.c -- Ground station observer location type + * + * Latitude/longitude stored internally in radians, altitude in meters. + * Flexible input: compass directions, decimal degrees, or tuple format. + * Output as "40.0000N 105.3000W 1655m". + */ + +#include "postgres.h" +#include "fmgr.h" +#include "utils/builtins.h" +#include "libpq/pqformat.h" +#include "types.h" +#include +#include +#include +#include + +#define DEG_TO_RAD (M_PI / 180.0) +#define RAD_TO_DEG (180.0 / M_PI) + +PG_FUNCTION_INFO_V1(observer_in); +PG_FUNCTION_INFO_V1(observer_out); +PG_FUNCTION_INFO_V1(observer_recv); +PG_FUNCTION_INFO_V1(observer_send); +PG_FUNCTION_INFO_V1(observer_lat); +PG_FUNCTION_INFO_V1(observer_lon); +PG_FUNCTION_INFO_V1(observer_alt); + + +/* + * skip_whitespace -- advance pointer past spaces and tabs + */ +static const char * +skip_whitespace(const char *p) +{ + while (*p == ' ' || *p == '\t') + p++; + return p; +} + + +/* + * observer_in -- parse text to pg_observer + * + * Accepted formats: + * 40.0N 105.3W 1655m compass directions with altitude + * 40.0N 105.3W compass directions, altitude = 0 + * 40.0 -105.3 1655 decimal degrees (N/E positive) with altitude + * (40.0,-105.3,1655) parenthesized tuple + */ +Datum +observer_in(PG_FUNCTION_ARGS) +{ + char *str = PG_GETARG_CSTRING(0); + pg_observer *result; + double lat_deg, lon_deg, alt_m; + const char *p; + char *endptr; + char dir; + + result = (pg_observer *) palloc(sizeof(pg_observer)); + alt_m = 0.0; + + p = skip_whitespace(str); + + /* + * Format 4: parenthesized tuple (lat_deg, lon_deg, alt_m) + */ + if (*p == '(') + { + int nfields; + + nfields = sscanf(str, " ( %lf , %lf , %lf )", + &lat_deg, &lon_deg, &alt_m); + if (nfields < 2) + goto bad_input; + if (nfields == 2) + alt_m = 0.0; + + if (lat_deg < -90.0 || lat_deg > 90.0) + goto bad_lat; + if (lon_deg < -180.0 || lon_deg > 180.0) + goto bad_lon; + + result->lat = lat_deg * DEG_TO_RAD; + result->lon = lon_deg * DEG_TO_RAD; + result->alt_m = alt_m; + PG_RETURN_POINTER(result); + } + + /* + * Parse latitude value + */ + lat_deg = strtod(p, &endptr); + if (endptr == p) + goto bad_input; + p = endptr; + + /* + * Check for compass direction suffix (N/S) + */ + dir = toupper((unsigned char) *p); + if (dir == 'N' || dir == 'S') + { + /* Format 1/3: compass directions */ + if (dir == 'S') + lat_deg = -lat_deg; + p++; + p = skip_whitespace(p); + + /* Parse longitude value */ + lon_deg = strtod(p, &endptr); + if (endptr == p) + goto bad_input; + p = endptr; + + /* Expect E or W */ + dir = toupper((unsigned char) *p); + if (dir == 'W') + lon_deg = -lon_deg; + else if (dir != 'E') + goto bad_input; + p++; + p = skip_whitespace(p); + + /* Optional altitude with 'm' suffix */ + if (*p != '\0') + { + alt_m = strtod(p, &endptr); + if (endptr == p) + goto bad_input; + p = endptr; + /* Skip optional 'm' suffix */ + if (toupper((unsigned char) *p) == 'M') + p++; + } + } + else + { + /* Format 2: decimal degrees, space-separated */ + p = skip_whitespace(p); + + lon_deg = strtod(p, &endptr); + if (endptr == p) + goto bad_input; + p = endptr; + p = skip_whitespace(p); + + /* Optional altitude */ + if (*p != '\0') + { + alt_m = strtod(p, &endptr); + if (endptr == p) + goto bad_input; + } + } + + /* Validate ranges (in degrees before conversion) */ + if (fabs(lat_deg) > 90.0) + goto bad_lat; + if (fabs(lon_deg) > 180.0) + goto bad_lon; + + result->lat = lat_deg * DEG_TO_RAD; + result->lon = lon_deg * DEG_TO_RAD; + result->alt_m = alt_m; + PG_RETURN_POINTER(result); + +bad_lat: + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("latitude out of range: %.6f", lat_deg), + errhint("Latitude must be between -90 and +90 degrees."))); + +bad_lon: + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("longitude out of range: %.6f", lon_deg), + errhint("Longitude must be between -180 and +180 degrees."))); + +bad_input: + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type observer: \"%s\"", str), + errhint("Expected \"40.0N 105.3W 1655m\", \"40.0 -105.3 1655\", " + "or \"(40.0,-105.3,1655)\"."))); + + PG_RETURN_NULL(); /* not reached */ +} + + +/* + * observer_out -- pg_observer to text + * + * Output: "40.0000N 105.3000W 1655m" + */ +Datum +observer_out(PG_FUNCTION_ARGS) +{ + pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0); + double lat_deg, lon_deg; + char lat_dir, lon_dir; + + lat_deg = obs->lat * RAD_TO_DEG; + lon_deg = obs->lon * RAD_TO_DEG; + + lat_dir = (lat_deg >= 0.0) ? 'N' : 'S'; + lon_dir = (lon_deg >= 0.0) ? 'E' : 'W'; + + PG_RETURN_CSTRING(psprintf("%.4f%c %.4f%c %.0fm", + fabs(lat_deg), lat_dir, + fabs(lon_deg), lon_dir, + obs->alt_m)); +} + + +/* + * observer_recv -- binary input + */ +Datum +observer_recv(PG_FUNCTION_ARGS) +{ + StringInfo buf = (StringInfo) PG_GETARG_POINTER(0); + pg_observer *result; + + result = (pg_observer *) palloc(sizeof(pg_observer)); + result->lat = pq_getmsgfloat8(buf); + result->lon = pq_getmsgfloat8(buf); + result->alt_m = pq_getmsgfloat8(buf); + + PG_RETURN_POINTER(result); +} + + +/* + * observer_send -- binary output + */ +Datum +observer_send(PG_FUNCTION_ARGS) +{ + pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0); + StringInfoData buf; + + pq_begintypsend(&buf); + pq_sendfloat8(&buf, obs->lat); + pq_sendfloat8(&buf, obs->lon); + pq_sendfloat8(&buf, obs->alt_m); + + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); +} + + +/* --- accessor functions --- */ + +/* + * observer_lat -- latitude in degrees + */ +Datum +observer_lat(PG_FUNCTION_ARGS) +{ + pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0); + PG_RETURN_FLOAT8(obs->lat * RAD_TO_DEG); +} + +/* + * observer_lon -- longitude in degrees (positive east) + */ +Datum +observer_lon(PG_FUNCTION_ARGS) +{ + pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0); + PG_RETURN_FLOAT8(obs->lon * RAD_TO_DEG); +} + +/* + * observer_alt -- altitude in meters + */ +Datum +observer_alt(PG_FUNCTION_ARGS) +{ + pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0); + PG_RETURN_FLOAT8(obs->alt_m); +} diff --git a/src/pass_funcs.c b/src/pass_funcs.c new file mode 100644 index 0000000..4575da9 --- /dev/null +++ b/src/pass_funcs.c @@ -0,0 +1,773 @@ +/* + * pass_funcs.c -- Satellite pass prediction for pg_orbit + * + * Finds visibility windows (AOS/LOS) for a satellite relative to a + * ground observer. Uses bisection on the elevation function to pin + * zero-crossings, then ternary search for peak elevation. + * + * The coarse scan steps at 30-second intervals -- fine enough for LEO + * orbits (~90 min period) that no pass shorter than a minute gets + * missed, coarse enough that a 7-day window doesn't require millions + * of propagation calls. + */ + +#include "postgres.h" +#include "fmgr.h" +#include "funcapi.h" +#include "utils/timestamp.h" +#include "utils/builtins.h" +#include "libpq/pqformat.h" +#include "norad.h" +#include "types.h" +#include +#include + +PG_FUNCTION_INFO_V1(pass_event_in); +PG_FUNCTION_INFO_V1(pass_event_out); +PG_FUNCTION_INFO_V1(pass_event_recv); +PG_FUNCTION_INFO_V1(pass_event_send); +PG_FUNCTION_INFO_V1(pass_aos_time); +PG_FUNCTION_INFO_V1(pass_max_el_time); +PG_FUNCTION_INFO_V1(pass_los_time); +PG_FUNCTION_INFO_V1(pass_max_elevation); +PG_FUNCTION_INFO_V1(pass_aos_azimuth); +PG_FUNCTION_INFO_V1(pass_los_azimuth); +PG_FUNCTION_INFO_V1(pass_duration); +PG_FUNCTION_INFO_V1(next_pass); +PG_FUNCTION_INFO_V1(predict_passes); +PG_FUNCTION_INFO_V1(pass_visible); + +#define DEG_TO_RAD (M_PI / 180.0) +#define RAD_TO_DEG (180.0 / M_PI) + +#define COARSE_STEP_JD (30.0 / 86400.0) /* 30 seconds */ +#define BISECT_TOL_JD (0.1 / 86400.0) /* 0.1 second */ +#define MIN_PASS_DURATION_JD (10.0 / 86400.0) /* 10 seconds */ +#define DEFAULT_WINDOW_DAYS 7.0 +#define POST_LOS_GAP_JD (60.0 / 86400.0) /* 1 minute */ +#define TERNARY_ITERATIONS 50 + +/* ---------------------------------------------------------------- + * Static helpers -- duplicated from coord_funcs.c because both + * files need them and they're too small to warrant a shared module. + * ---------------------------------------------------------------- + */ + +/* + * Convert pg_tle to sat_code's tle_t. No unit conversion needed -- + * both store radians, radians/min, Julian dates. + */ +static void +pg_tle_to_sat_code(const pg_tle *src, tle_t *dst) +{ + memset(dst, 0, sizeof(tle_t)); + + dst->epoch = src->epoch; + dst->xincl = src->inclination; + dst->xnodeo = src->raan; + dst->eo = src->eccentricity; + dst->omegao = src->arg_perigee; + dst->xmo = src->mean_anomaly; + dst->xno = src->mean_motion; + dst->xndt2o = src->mean_motion_dot; + dst->xndd6o = src->mean_motion_ddot; + dst->bstar = src->bstar; + + dst->norad_number = src->norad_id; + dst->bulletin_number = src->elset_num; + dst->revolution_number = src->rev_num; + dst->classification = src->classification; + dst->ephemeris_type = src->ephemeris_type; + + memcpy(dst->intl_desig, src->intl_desig, 9); +} + +/* + * Propagate TLE to a Julian date. Returns sat_code error code. + * pos[3] in km (TEME), vel[3] in km/min (TEME). + */ +static int +do_propagate(const pg_tle *tle, double jd, double *pos, double *vel) +{ + tle_t sat; + double *params; + int is_deep; + int err; + double tsince; + + pg_tle_to_sat_code(tle, &sat); + + is_deep = select_ephemeris(&sat); + if (is_deep < 0) + return -99; /* invalid TLE */ + + tsince = jd_to_minutes_since_epoch(jd, sat.epoch); + + params = palloc(sizeof(double) * N_SAT_PARAMS); + + if (is_deep) + { + SDP4_init(params, &sat); + err = SDP4(tsince, &sat, params, pos, vel); + } + else + { + SGP4_init(params, &sat); + err = SGP4(tsince, &sat, params, pos, vel); + } + + pfree(params); + + return err; +} + +/* + * Greenwich Mean Sidereal Time from Julian date. + * Returns GMST in radians. Uses the IAU 1982 formula matching + * the low-precision model inside SGP4. + */ +static double +gmst_from_jd(double jd) +{ + double ut1, tu; + double gmst; + + /* Julian centuries of UT1 from J2000.0 */ + ut1 = jd - J2000_JD; + tu = ut1 / 36525.0; + + /* GMST in seconds at 0h UT1, then add fractional day rotation */ + gmst = 67310.54841 + + (876600.0 * 3600.0 + 8640184.812866) * tu + + 0.093104 * tu * tu + - 6.2e-6 * tu * tu * tu; + + /* Convert seconds to radians, mod 2pi */ + gmst = fmod(gmst * M_PI / 43200.0, 2.0 * M_PI); + if (gmst < 0.0) + gmst += 2.0 * M_PI; + + return gmst; +} + +/* + * Rotate TEME position (and optionally velocity) to ECEF via GMST. + * Only the z-axis rotation matters for SGP4's simplified nutation. + */ +static void +teme_to_ecef(const double *pos_teme, const double *vel_teme, + double gmst, double *pos_ecef, double *vel_ecef) +{ + double cg = cos(gmst); + double sg = sin(gmst); + + pos_ecef[0] = cg * pos_teme[0] + sg * pos_teme[1]; + pos_ecef[1] = -sg * pos_teme[0] + cg * pos_teme[1]; + pos_ecef[2] = pos_teme[2]; + + if (vel_teme && vel_ecef) + { + /* Earth rotation rate, rad/min */ + double omega_e = 7.29211514670698e-5 * 60.0; + + vel_ecef[0] = cg * vel_teme[0] + sg * vel_teme[1] + + omega_e * pos_ecef[1]; + vel_ecef[1] = -sg * vel_teme[0] + cg * vel_teme[1] + - omega_e * pos_ecef[0]; + vel_ecef[2] = vel_teme[2]; + } +} + +/* + * Observer geodetic (radians, meters) to ECEF (km). + * Uses WGS-84 ellipsoid for ground station positioning. + */ +static void +observer_to_ecef(const pg_observer *obs, double *ecef) +{ + double sinlat = sin(obs->lat); + double coslat = cos(obs->lat); + double sinlon = sin(obs->lon); + double coslon = cos(obs->lon); + double N; /* radius of curvature in the prime vertical */ + double alt_km = obs->alt_m / 1000.0; + + N = WGS84_A / sqrt(1.0 - WGS84_E2 * sinlat * sinlat); + + ecef[0] = (N + alt_km) * coslat * coslon; + ecef[1] = (N + alt_km) * coslat * sinlon; + ecef[2] = (N * (1.0 - WGS84_E2) + alt_km) * sinlat; +} + +/* + * Compute topocentric azimuth, elevation, and range from ECEF positions. + * Azimuth: 0=N, 90=E, 180=S, 270=W (radians). + * Elevation: positive above horizon (radians). + */ +static void +ecef_to_topocentric(const double *sat_ecef, const double *obs_ecef, + double obs_lat, double obs_lon, + double *az, double *el, double *range_km) +{ + double dx, dy, dz; + double sinlat, coslat, sinlon, coslon; + double south, east, up; + double rng; + + dx = sat_ecef[0] - obs_ecef[0]; + dy = sat_ecef[1] - obs_ecef[1]; + dz = sat_ecef[2] - obs_ecef[2]; + + sinlat = sin(obs_lat); + coslat = cos(obs_lat); + sinlon = sin(obs_lon); + coslon = cos(obs_lon); + + /* Rotate difference vector into SEZ (south-east-zenith) frame */ + south = sinlat * coslon * dx + sinlat * sinlon * dy - coslat * dz; + east = -sinlon * dx + coslon * dy; + up = coslat * coslon * dx + coslat * sinlon * dy + sinlat * dz; + + rng = sqrt(dx * dx + dy * dy + dz * dz); + + *range_km = rng; + *el = asin(up / rng); + + /* Azimuth from north, measured clockwise */ + *az = atan2(east, -south); + if (*az < 0.0) + *az += 2.0 * M_PI; +} + + +/* ---------------------------------------------------------------- + * elevation_at_jd -- the function we bisect on + * + * Returns satellite elevation in radians relative to the observer's + * local horizon. Negative means below horizon. On hard propagation + * errors, returns -pi (well below any real horizon). + * ---------------------------------------------------------------- + */ +static double +elevation_at_jd(const pg_tle *tle, const pg_observer *obs, + double jd, double *az_out) +{ + double pos[3], vel[3]; + double gmst, pos_ecef[3], obs_ecef[3]; + double az, el, range_km; + int err; + + err = do_propagate(tle, jd, pos, vel); + + /* On hard errors, return well below horizon */ + if (err == SXPX_ERR_NEARLY_PARABOLIC || + err == SXPX_ERR_NEGATIVE_MAJOR_AXIS || + err == SXPX_ERR_NEGATIVE_XN || + err == SXPX_ERR_CONVERGENCE_FAIL || + err == -99) + return -M_PI; + + gmst = gmst_from_jd(jd); + teme_to_ecef(pos, NULL, gmst, pos_ecef, NULL); + observer_to_ecef(obs, obs_ecef); + ecef_to_topocentric(pos_ecef, obs_ecef, obs->lat, obs->lon, + &az, &el, &range_km); + + if (az_out) + *az_out = az; + + return el; +} + + +/* ---------------------------------------------------------------- + * find_next_pass -- core pass-finding algorithm + * + * Scans from start_jd to stop_jd looking for an elevation zero + * crossing (rising edge). When found, bisects to refine AOS, + * scans forward to find LOS, bisects to refine LOS, then uses + * ternary search to locate peak elevation. + * + * Passes below min_el_rad are silently skipped; scanning resumes + * from their LOS. + * + * Returns true if a qualifying pass was found. + * ---------------------------------------------------------------- + */ +static bool +find_next_pass(const pg_tle *tle, const pg_observer *obs, + double start_jd, double stop_jd, + double min_el_rad, + double *aos_jd, double *los_jd, + double *max_el_jd, double *max_el, + double *aos_az, double *los_az) +{ + double jd = start_jd; + double prev_el, curr_el; + double az; + + prev_el = elevation_at_jd(tle, obs, jd, NULL); + + while (jd < stop_jd) + { + jd += COARSE_STEP_JD; + if (jd > stop_jd) + jd = stop_jd; + + curr_el = elevation_at_jd(tle, obs, jd, NULL); + + /* Rising edge: was below horizon, now above */ + if (prev_el <= 0.0 && curr_el > 0.0) + { + double lo, hi, mid; + double peak_el; + double scan_jd, scan_el; + + /* Bisect to find AOS */ + lo = jd - COARSE_STEP_JD; + hi = jd; + while (hi - lo > BISECT_TOL_JD) + { + mid = (lo + hi) / 2.0; + if (elevation_at_jd(tle, obs, mid, NULL) > 0.0) + hi = mid; + else + lo = mid; + } + *aos_jd = (lo + hi) / 2.0; + elevation_at_jd(tle, obs, *aos_jd, &az); + *aos_az = az; + + /* Scan forward to find LOS */ + scan_jd = *aos_jd; + peak_el = 0.0; + + while (scan_jd < stop_jd) + { + scan_jd += COARSE_STEP_JD; + scan_el = elevation_at_jd(tle, obs, scan_jd, NULL); + + if (scan_el > peak_el) + peak_el = scan_el; + + if (scan_el <= 0.0) + { + /* Bisect to find LOS */ + lo = scan_jd - COARSE_STEP_JD; + hi = scan_jd; + while (hi - lo > BISECT_TOL_JD) + { + mid = (lo + hi) / 2.0; + if (elevation_at_jd(tle, obs, mid, NULL) > 0.0) + lo = mid; + else + hi = mid; + } + *los_jd = (lo + hi) / 2.0; + break; + } + } + + /* Ran past the search window without finding LOS -- clip */ + if (scan_jd >= stop_jd) + *los_jd = stop_jd; + + /* Skip degenerate passes */ + if (*los_jd - *aos_jd < MIN_PASS_DURATION_JD) + { + jd = *los_jd; + prev_el = 0.0; + continue; + } + + /* Refine peak elevation with ternary search */ + lo = *aos_jd; + hi = *los_jd; + for (int i = 0; i < TERNARY_ITERATIONS; i++) + { + double m1 = lo + (hi - lo) / 3.0; + double m2 = hi - (hi - lo) / 3.0; + + if (elevation_at_jd(tle, obs, m1, NULL) < + elevation_at_jd(tle, obs, m2, NULL)) + lo = m1; + else + hi = m2; + } + *max_el_jd = (lo + hi) / 2.0; + *max_el = elevation_at_jd(tle, obs, *max_el_jd, NULL); + + elevation_at_jd(tle, obs, *los_jd, &az); + *los_az = az; + + /* Below the caller's minimum elevation threshold -- skip */ + if (*max_el < min_el_rad) + { + jd = *los_jd; + prev_el = 0.0; + continue; + } + + return true; + } + + prev_el = curr_el; + } + + return false; +} + + +/* ---------------------------------------------------------------- + * pass_event type I/O + * ---------------------------------------------------------------- + */ + +/* + * pass_event_in -- parse text to pg_pass_event + * + * Format: (aos_ts,maxel_ts,los_ts,max_el_deg,aos_az_deg,los_az_deg) + * Timestamps are raw int64 microseconds (PG internal representation). + */ +Datum +pass_event_in(PG_FUNCTION_ARGS) +{ + char *str = PG_GETARG_CSTRING(0); + pg_pass_event *result; + long long aos_raw, maxel_raw, los_raw; + double max_el_deg, aos_az_deg, los_az_deg; + int nfields; + + result = (pg_pass_event *) palloc(sizeof(pg_pass_event)); + + nfields = sscanf(str, " ( %lld , %lld , %lld , %lf , %lf , %lf )", + &aos_raw, &maxel_raw, &los_raw, + &max_el_deg, &aos_az_deg, &los_az_deg); + + if (nfields != 6) + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type pass_event: \"%s\"", str), + errhint("Expected (aos_usec,maxel_usec,los_usec,max_el_deg,aos_az_deg,los_az_deg)."))); + + result->aos_time = (int64) aos_raw; + result->max_el_time = (int64) maxel_raw; + result->los_time = (int64) los_raw; + result->max_elevation = max_el_deg; + result->aos_azimuth = aos_az_deg; + result->los_azimuth = los_az_deg; + + PG_RETURN_POINTER(result); +} + + +/* + * pass_event_out -- pg_pass_event to human-readable text + * + * Format: (2024-01-01 12:00:00+00,2024-01-01 12:05:00+00,2024-01-01 12:10:00+00,45.2,180.0,350.0) + * Timestamps formatted via DirectFunctionCall1(timestamptz_out, ...). + */ +Datum +pass_event_out(PG_FUNCTION_ARGS) +{ + pg_pass_event *pe = (pg_pass_event *) PG_GETARG_POINTER(0); + char *aos_str; + char *maxel_str; + char *los_str; + + aos_str = DatumGetCString( + DirectFunctionCall1(timestamptz_out, + Int64GetDatum(pe->aos_time))); + maxel_str = DatumGetCString( + DirectFunctionCall1(timestamptz_out, + Int64GetDatum(pe->max_el_time))); + los_str = DatumGetCString( + DirectFunctionCall1(timestamptz_out, + Int64GetDatum(pe->los_time))); + + PG_RETURN_CSTRING(psprintf("(%s,%s,%s,%.1f,%.1f,%.1f)", + aos_str, maxel_str, los_str, + pe->max_elevation, + pe->aos_azimuth, + pe->los_azimuth)); +} + + +/* + * pass_event_recv -- binary input (3 int64 + 3 float8 = 48 bytes) + */ +Datum +pass_event_recv(PG_FUNCTION_ARGS) +{ + StringInfo buf = (StringInfo) PG_GETARG_POINTER(0); + pg_pass_event *result; + + result = (pg_pass_event *) palloc(sizeof(pg_pass_event)); + result->aos_time = pq_getmsgint64(buf); + result->max_el_time = pq_getmsgint64(buf); + result->los_time = pq_getmsgint64(buf); + result->max_elevation = pq_getmsgfloat8(buf); + result->aos_azimuth = pq_getmsgfloat8(buf); + result->los_azimuth = pq_getmsgfloat8(buf); + + PG_RETURN_POINTER(result); +} + + +/* + * pass_event_send -- binary output + */ +Datum +pass_event_send(PG_FUNCTION_ARGS) +{ + pg_pass_event *pe = (pg_pass_event *) PG_GETARG_POINTER(0); + StringInfoData buf; + + pq_begintypsend(&buf); + pq_sendint64(&buf, pe->aos_time); + pq_sendint64(&buf, pe->max_el_time); + pq_sendint64(&buf, pe->los_time); + pq_sendfloat8(&buf, pe->max_elevation); + pq_sendfloat8(&buf, pe->aos_azimuth); + pq_sendfloat8(&buf, pe->los_azimuth); + + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); +} + + +/* ---------------------------------------------------------------- + * pass_event accessor functions + * ---------------------------------------------------------------- + */ + +Datum +pass_aos_time(PG_FUNCTION_ARGS) +{ + pg_pass_event *pe = (pg_pass_event *) PG_GETARG_POINTER(0); + + PG_RETURN_TIMESTAMPTZ(pe->aos_time); +} + +Datum +pass_max_el_time(PG_FUNCTION_ARGS) +{ + pg_pass_event *pe = (pg_pass_event *) PG_GETARG_POINTER(0); + + PG_RETURN_TIMESTAMPTZ(pe->max_el_time); +} + +Datum +pass_los_time(PG_FUNCTION_ARGS) +{ + pg_pass_event *pe = (pg_pass_event *) PG_GETARG_POINTER(0); + + PG_RETURN_TIMESTAMPTZ(pe->los_time); +} + +Datum +pass_max_elevation(PG_FUNCTION_ARGS) +{ + pg_pass_event *pe = (pg_pass_event *) PG_GETARG_POINTER(0); + + PG_RETURN_FLOAT8(pe->max_elevation); +} + +Datum +pass_aos_azimuth(PG_FUNCTION_ARGS) +{ + pg_pass_event *pe = (pg_pass_event *) PG_GETARG_POINTER(0); + + PG_RETURN_FLOAT8(pe->aos_azimuth); +} + +Datum +pass_los_azimuth(PG_FUNCTION_ARGS) +{ + pg_pass_event *pe = (pg_pass_event *) PG_GETARG_POINTER(0); + + PG_RETURN_FLOAT8(pe->los_azimuth); +} + +/* + * pass_duration -- time from AOS to LOS as a PostgreSQL interval + */ +Datum +pass_duration(PG_FUNCTION_ARGS) +{ + pg_pass_event *pe = (pg_pass_event *) PG_GETARG_POINTER(0); + Interval *result; + + result = (Interval *) palloc(sizeof(Interval)); + result->time = pe->los_time - pe->aos_time; /* microseconds */ + result->day = 0; + result->month = 0; + + PG_RETURN_INTERVAL_P(result); +} + + +/* ---------------------------------------------------------------- + * next_pass(tle, observer, from_time) -> pass_event + * + * Finds the next pass above the horizon starting from from_time. + * Searches a 7-day window. Returns NULL if no pass is found. + * ---------------------------------------------------------------- + */ +Datum +next_pass(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(1); + int64 from_ts = PG_GETARG_INT64(2); + + double start_jd, stop_jd; + double aos_jd, los_jd, max_el_jd, max_el; + double aos_az, los_az; + pg_pass_event *result; + + start_jd = timestamptz_to_jd(from_ts); + stop_jd = start_jd + DEFAULT_WINDOW_DAYS; + + if (!find_next_pass(tle, obs, start_jd, stop_jd, + 0.0, /* minimum elevation = 0 degrees */ + &aos_jd, &los_jd, + &max_el_jd, &max_el, + &aos_az, &los_az)) + PG_RETURN_NULL(); + + result = (pg_pass_event *) palloc(sizeof(pg_pass_event)); + result->aos_time = jd_to_timestamptz(aos_jd); + result->max_el_time = jd_to_timestamptz(max_el_jd); + result->los_time = jd_to_timestamptz(los_jd); + result->max_elevation = max_el * RAD_TO_DEG; + result->aos_azimuth = aos_az * RAD_TO_DEG; + result->los_azimuth = los_az * RAD_TO_DEG; + + PG_RETURN_POINTER(result); +} + + +/* ---------------------------------------------------------------- + * predict_passes(tle, observer, start, stop [, min_elevation]) + * -> SETOF pass_event + * + * Returns all passes in the given time window. Optional 5th arg + * sets the minimum peak elevation filter in degrees (default 0). + * ---------------------------------------------------------------- + */ + +typedef struct +{ + pg_tle tle; + pg_observer obs; + double current_jd; + double stop_jd; + double min_el_rad; +} predict_passes_ctx; + +Datum +predict_passes(PG_FUNCTION_ARGS) +{ + FuncCallContext *funcctx; + predict_passes_ctx *ctx; + + if (SRF_IS_FIRSTCALL()) + { + MemoryContext oldctx; + pg_tle *tle; + pg_observer *obs; + int64 start_ts; + int64 stop_ts; + double min_el_deg; + + funcctx = SRF_FIRSTCALL_INIT(); + oldctx = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx); + + tle = (pg_tle *) PG_GETARG_POINTER(0); + obs = (pg_observer *) PG_GETARG_POINTER(1); + start_ts = PG_GETARG_INT64(2); + stop_ts = PG_GETARG_INT64(3); + + min_el_deg = (PG_NARGS() > 4 && !PG_ARGISNULL(4)) + ? PG_GETARG_FLOAT8(4) + : 0.0; + + if (stop_ts <= start_ts) + ereport(ERROR, + (errcode(ERRCODE_INVALID_PARAMETER_VALUE), + errmsg("stop time must be after start time"))); + + ctx = (predict_passes_ctx *) + palloc0(sizeof(predict_passes_ctx)); + + memcpy(&ctx->tle, tle, sizeof(pg_tle)); + memcpy(&ctx->obs, obs, sizeof(pg_observer)); + ctx->current_jd = timestamptz_to_jd(start_ts); + ctx->stop_jd = timestamptz_to_jd(stop_ts); + ctx->min_el_rad = min_el_deg * DEG_TO_RAD; + + funcctx->user_fctx = ctx; + + MemoryContextSwitchTo(oldctx); + } + + funcctx = SRF_PERCALL_SETUP(); + ctx = (predict_passes_ctx *) funcctx->user_fctx; + + { + double aos_jd, los_jd, max_el_jd, max_el; + double aos_az, los_az; + pg_pass_event *result; + + if (!find_next_pass(&ctx->tle, &ctx->obs, + ctx->current_jd, ctx->stop_jd, + ctx->min_el_rad, + &aos_jd, &los_jd, + &max_el_jd, &max_el, + &aos_az, &los_az)) + SRF_RETURN_DONE(funcctx); + + result = (pg_pass_event *) palloc(sizeof(pg_pass_event)); + result->aos_time = jd_to_timestamptz(aos_jd); + result->max_el_time = jd_to_timestamptz(max_el_jd); + result->los_time = jd_to_timestamptz(los_jd); + result->max_elevation = max_el * RAD_TO_DEG; + result->aos_azimuth = aos_az * RAD_TO_DEG; + result->los_azimuth = los_az * RAD_TO_DEG; + + /* Advance past this pass before the next call */ + ctx->current_jd = los_jd + POST_LOS_GAP_JD; + + SRF_RETURN_NEXT(funcctx, PointerGetDatum(result)); + } +} + + +/* ---------------------------------------------------------------- + * pass_visible(tle, observer, start, stop) -> bool + * + * Returns true if any pass crosses above the horizon in the window. + * Cheaper than predict_passes when you only need a yes/no answer. + * ---------------------------------------------------------------- + */ +Datum +pass_visible(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(1); + int64 start_ts = PG_GETARG_INT64(2); + int64 stop_ts = PG_GETARG_INT64(3); + + double start_jd, stop_jd; + double aos_jd, los_jd, max_el_jd, max_el; + double aos_az, los_az; + + start_jd = timestamptz_to_jd(start_ts); + stop_jd = timestamptz_to_jd(stop_ts); + + PG_RETURN_BOOL(find_next_pass(tle, obs, start_jd, stop_jd, + 0.0, + &aos_jd, &los_jd, + &max_el_jd, &max_el, + &aos_az, &los_az)); +} diff --git a/src/pg_orbit.c b/src/pg_orbit.c new file mode 100644 index 0000000..1e3da2c --- /dev/null +++ b/src/pg_orbit.c @@ -0,0 +1,12 @@ +/* + * pg_orbit.c -- Extension entry point + * + * PostgreSQL extension for orbital mechanics. + * Provides TLE, ECI, geodetic, topocentric, observer, and pass_event types + * with SGP4/SDP4 propagation, coordinate transforms, and pass prediction. + */ + +#include "postgres.h" +#include "fmgr.h" + +PG_MODULE_MAGIC; diff --git a/src/sgp4_funcs.c b/src/sgp4_funcs.c new file mode 100644 index 0000000..1eac6cc --- /dev/null +++ b/src/sgp4_funcs.c @@ -0,0 +1,381 @@ +/* + * sgp4_funcs.c -- SGP4/SDP4 propagation functions for pg_orbit + * + * Wraps Bill Gray's sat_code implementation. Near-earth objects + * use SGP4, deep-space objects use SDP4. The choice is automatic + * based on orbital period (threshold: 225 minutes). + * + * All propagation uses WGS-72 constants internally (via sat_code). + * Velocities are converted from km/min to km/s on output. + */ + +#include "postgres.h" +#include "fmgr.h" +#include "funcapi.h" +#include "utils/timestamp.h" +#include "norad.h" +#include "types.h" +#include +#include + +PG_FUNCTION_INFO_V1(sgp4_propagate); +PG_FUNCTION_INFO_V1(sgp4_propagate_series); +PG_FUNCTION_INFO_V1(tle_distance); + +/* ---------------------------------------------------------------- + * TLE struct conversion helpers + * + * pg_tle (our type) and tle_t (sat_code) store the same data in + * different field layouts. These copy between them without any + * unit conversion -- both use radians, radians/min, Julian dates. + * ---------------------------------------------------------------- + */ +static void +pg_tle_to_sat_code(const pg_tle *src, tle_t *dst) +{ + memset(dst, 0, sizeof(tle_t)); + + dst->epoch = src->epoch; + dst->xincl = src->inclination; + dst->xnodeo = src->raan; + dst->eo = src->eccentricity; + dst->omegao = src->arg_perigee; + dst->xmo = src->mean_anomaly; + dst->xno = src->mean_motion; + dst->xndt2o = src->mean_motion_dot; + dst->xndd6o = src->mean_motion_ddot; + dst->bstar = src->bstar; + + dst->norad_number = src->norad_id; + dst->bulletin_number = src->elset_num; + dst->revolution_number = src->rev_num; + dst->classification = src->classification; + dst->ephemeris_type = src->ephemeris_type; + + memcpy(dst->intl_desig, src->intl_desig, 9); +} + +/* ---------------------------------------------------------------- + * propagate_tle -- shared propagation logic + * + * Initializes the appropriate model, propagates to tsince minutes, + * and fills pos[3] (km) and vel[3] (km/min). Returns the sat_code + * error code (0 on success, negative on error/warning). + * ---------------------------------------------------------------- + */ +static int +propagate_tle(const tle_t *sat, double tsince, double *pos, double *vel) +{ + double *params; + int is_deep; + int err; + + is_deep = select_ephemeris(sat); + if (is_deep < 0) + ereport(ERROR, + (errcode(ERRCODE_DATA_EXCEPTION), + errmsg("invalid TLE for NORAD %d: " + "mean motion or eccentricity out of range", + sat->norad_number))); + + params = palloc(sizeof(double) * N_SAT_PARAMS); + + if (is_deep) + { + SDP4_init(params, sat); + err = SDP4(tsince, sat, params, pos, vel); + } + else + { + SGP4_init(params, sat); + err = SGP4(tsince, sat, params, pos, vel); + } + + pfree(params); + + return err; +} + +/* Map sat_code error codes to human-readable messages */ +static const char * +sxpx_error_msg(int err) +{ + switch (err) + { + case SXPX_ERR_NEARLY_PARABOLIC: + return "nearly parabolic orbit"; + case SXPX_ERR_NEGATIVE_MAJOR_AXIS: + return "negative semi-major axis (decayed)"; + case SXPX_WARN_ORBIT_WITHIN_EARTH: + return "orbit center within Earth"; + case SXPX_WARN_PERIGEE_WITHIN_EARTH: + return "perigee within Earth"; + case SXPX_ERR_NEGATIVE_XN: + return "negative mean motion"; + case SXPX_ERR_CONVERGENCE_FAIL: + return "Kepler equation convergence failure"; + default: + return "unknown propagation error"; + } +} + +/* + * Check propagation result. Warnings (-3, -4) are tolerated -- + * the position is still mathematically valid. Hard errors (<= -5 + * or -1, -2) abort the query. + */ +static void +check_propagation_result(int err, int norad_number, double tsince) +{ + if (err == 0) + return; + + /* Perigee/orbit within Earth: valid state vector, just unusual */ + if (err == SXPX_WARN_ORBIT_WITHIN_EARTH || + err == SXPX_WARN_PERIGEE_WITHIN_EARTH) + { + ereport(NOTICE, + (errmsg("NORAD %d at t+%.1f min: %s", + norad_number, tsince, sxpx_error_msg(err)))); + return; + } + + ereport(ERROR, + (errcode(ERRCODE_DATA_EXCEPTION), + errmsg("SGP4/SDP4 propagation failed for NORAD %d at t+%.1f min: %s", + norad_number, tsince, sxpx_error_msg(err)))); +} + +/* ---------------------------------------------------------------- + * sgp4_propagate(tle, timestamptz) -> eci_position + * + * Propagate a TLE to a single point in time. + * ---------------------------------------------------------------- + */ +Datum +sgp4_propagate(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + int64 ts = PG_GETARG_INT64(1); + + tle_t sat; + double jd; + double tsince; + double pos[3]; + double vel[3]; + int err; + pg_eci *result; + + pg_tle_to_sat_code(tle, &sat); + + jd = timestamptz_to_jd(ts); + tsince = jd_to_minutes_since_epoch(jd, sat.epoch); + + err = propagate_tle(&sat, tsince, pos, vel); + check_propagation_result(err, sat.norad_number, tsince); + + result = (pg_eci *) palloc(sizeof(pg_eci)); + result->x = pos[0]; + result->y = pos[1]; + result->z = pos[2]; + result->vx = vel[0] / 60.0; /* km/min -> km/s */ + result->vy = vel[1] / 60.0; + result->vz = vel[2] / 60.0; + + PG_RETURN_POINTER(result); +} + +/* ---------------------------------------------------------------- + * sgp4_propagate_series(tle, start, stop, step) -> SETOF (timestamptz, eci_position) + * + * Generates a time series of propagated positions. The model is + * initialized once and reused for every step. + * ---------------------------------------------------------------- + */ + +typedef struct +{ + tle_t sat; + double params[N_SAT_PARAMS]; + int is_deep; + double epoch_jd; + int64 start_ts; + int64 step_usec; +} propagate_series_ctx; + +Datum +sgp4_propagate_series(PG_FUNCTION_ARGS) +{ + FuncCallContext *funcctx; + propagate_series_ctx *ctx; + + if (SRF_IS_FIRSTCALL()) + { + MemoryContext oldctx; + pg_tle *tle; + int64 start_ts; + int64 stop_ts; + Interval *step; + int64 step_usec; + int64 span; + uint64 nsteps; + TupleDesc tupdesc; + + funcctx = SRF_FIRSTCALL_INIT(); + oldctx = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx); + + tle = (pg_tle *) PG_GETARG_POINTER(0); + start_ts = PG_GETARG_INT64(1); + stop_ts = PG_GETARG_INT64(2); + step = PG_GETARG_INTERVAL_P(3); + + /* + * Convert interval to microseconds. We only use the time + * component -- month/day fields would need calendar logic + * that doesn't belong in a propagation step. + */ + step_usec = step->time + + (int64) step->day * USECS_PER_DAY + + (int64) step->month * (30 * USECS_PER_DAY); + + if (step_usec <= 0) + ereport(ERROR, + (errcode(ERRCODE_INVALID_PARAMETER_VALUE), + errmsg("step interval must be positive"))); + + if (stop_ts < start_ts) + ereport(ERROR, + (errcode(ERRCODE_INVALID_PARAMETER_VALUE), + errmsg("stop time must be >= start time"))); + + span = stop_ts - start_ts; + nsteps = (uint64)(span / step_usec) + 1; + + /* Allocate persistent context for the SRF */ + ctx = (propagate_series_ctx *) + palloc0(sizeof(propagate_series_ctx)); + + pg_tle_to_sat_code(tle, &ctx->sat); + + ctx->is_deep = select_ephemeris(&ctx->sat); + if (ctx->is_deep < 0) + ereport(ERROR, + (errcode(ERRCODE_DATA_EXCEPTION), + errmsg("invalid TLE for NORAD %d: " + "mean motion or eccentricity out of range", + ctx->sat.norad_number))); + + /* Initialize the model once */ + if (ctx->is_deep) + SDP4_init(ctx->params, &ctx->sat); + else + SGP4_init(ctx->params, &ctx->sat); + + ctx->epoch_jd = ctx->sat.epoch; + ctx->start_ts = start_ts; + ctx->step_usec = step_usec; + + funcctx->max_calls = nsteps; + funcctx->user_fctx = ctx; + + /* Build the output tuple descriptor: (timestamptz, eci_position) */ + tupdesc = CreateTemplateTupleDesc(7); + TupleDescInitEntry(tupdesc, 1, "t", TIMESTAMPTZOID, -1, 0); + TupleDescInitEntry(tupdesc, 2, "x", FLOAT8OID, -1, 0); + TupleDescInitEntry(tupdesc, 3, "y", FLOAT8OID, -1, 0); + TupleDescInitEntry(tupdesc, 4, "z", FLOAT8OID, -1, 0); + TupleDescInitEntry(tupdesc, 5, "vx", FLOAT8OID, -1, 0); + TupleDescInitEntry(tupdesc, 6, "vy", FLOAT8OID, -1, 0); + TupleDescInitEntry(tupdesc, 7, "vz", FLOAT8OID, -1, 0); + + funcctx->tuple_desc = BlessTupleDesc(tupdesc); + + MemoryContextSwitchTo(oldctx); + } + + funcctx = SRF_PERCALL_SETUP(); + ctx = (propagate_series_ctx *) funcctx->user_fctx; + + if (funcctx->call_cntr < funcctx->max_calls) + { + int64 ts; + double jd; + double tsince; + double pos[3]; + double vel[3]; + int err; + Datum values[7]; + bool nulls[7]; + HeapTuple tuple; + + ts = ctx->start_ts + (int64) funcctx->call_cntr * ctx->step_usec; + jd = timestamptz_to_jd(ts); + tsince = jd_to_minutes_since_epoch(jd, ctx->epoch_jd); + + if (ctx->is_deep) + err = SDP4(tsince, &ctx->sat, ctx->params, pos, vel); + else + err = SGP4(tsince, &ctx->sat, ctx->params, pos, vel); + + check_propagation_result(err, ctx->sat.norad_number, tsince); + + memset(nulls, 0, sizeof(nulls)); + + values[0] = Int64GetDatum(ts); + values[1] = Float8GetDatum(pos[0]); + values[2] = Float8GetDatum(pos[1]); + values[3] = Float8GetDatum(pos[2]); + values[4] = Float8GetDatum(vel[0] / 60.0); /* km/min -> km/s */ + values[5] = Float8GetDatum(vel[1] / 60.0); + values[6] = Float8GetDatum(vel[2] / 60.0); + + tuple = heap_form_tuple(funcctx->tuple_desc, values, nulls); + + SRF_RETURN_NEXT(funcctx, HeapTupleGetDatum(tuple)); + } + + SRF_RETURN_DONE(funcctx); +} + +/* ---------------------------------------------------------------- + * tle_distance(tle_a, tle_b, timestamptz) -> float8 + * + * Euclidean distance in km between two TLEs at a reference time. + * Useful for conjunction screening and cluster analysis. + * ---------------------------------------------------------------- + */ +Datum +tle_distance(PG_FUNCTION_ARGS) +{ + pg_tle *tle_a = (pg_tle *) PG_GETARG_POINTER(0); + pg_tle *tle_b = (pg_tle *) PG_GETARG_POINTER(1); + int64 ts = PG_GETARG_INT64(2); + + tle_t sat_a, sat_b; + double jd; + double tsince_a, tsince_b; + double pos_a[3], vel_a[3]; + double pos_b[3], vel_b[3]; + int err; + double dx, dy, dz; + + pg_tle_to_sat_code(tle_a, &sat_a); + pg_tle_to_sat_code(tle_b, &sat_b); + + jd = timestamptz_to_jd(ts); + + tsince_a = jd_to_minutes_since_epoch(jd, sat_a.epoch); + tsince_b = jd_to_minutes_since_epoch(jd, sat_b.epoch); + + err = propagate_tle(&sat_a, tsince_a, pos_a, vel_a); + check_propagation_result(err, sat_a.norad_number, tsince_a); + + err = propagate_tle(&sat_b, tsince_b, pos_b, vel_b); + check_propagation_result(err, sat_b.norad_number, tsince_b); + + dx = pos_a[0] - pos_b[0]; + dy = pos_a[1] - pos_b[1]; + dz = pos_a[2] - pos_b[2]; + + PG_RETURN_FLOAT8(sqrt(dx * dx + dy * dy + dz * dz)); +} diff --git a/src/tle_type.c b/src/tle_type.c new file mode 100644 index 0000000..7e4cb0f --- /dev/null +++ b/src/tle_type.c @@ -0,0 +1,451 @@ +/* + * tle_type.c -- TLE custom type for PostgreSQL + * + * Implements input/output, binary I/O, and accessor functions for the + * Two-Line Element type. Parsing and formatting delegate to sat_code's + * parse_elements() and write_elements_in_tle_format(). + * + * Angular elements are stored internally in radians (matching sat_code). + * Accessor functions convert to degrees and revs/day for human consumption. + */ + +#include "postgres.h" +#include "fmgr.h" +#include "utils/builtins.h" +#include "libpq/pqformat.h" +#include "norad.h" +#include "types.h" +#include +#include + +PG_FUNCTION_INFO_V1(tle_in); +PG_FUNCTION_INFO_V1(tle_out); +PG_FUNCTION_INFO_V1(tle_recv); +PG_FUNCTION_INFO_V1(tle_send); +PG_FUNCTION_INFO_V1(tle_epoch); +PG_FUNCTION_INFO_V1(tle_norad_id); +PG_FUNCTION_INFO_V1(tle_inclination); +PG_FUNCTION_INFO_V1(tle_eccentricity); +PG_FUNCTION_INFO_V1(tle_raan); +PG_FUNCTION_INFO_V1(tle_arg_perigee); +PG_FUNCTION_INFO_V1(tle_mean_anomaly); +PG_FUNCTION_INFO_V1(tle_mean_motion); +PG_FUNCTION_INFO_V1(tle_bstar); +PG_FUNCTION_INFO_V1(tle_period); +PG_FUNCTION_INFO_V1(tle_age); +PG_FUNCTION_INFO_V1(tle_perigee); +PG_FUNCTION_INFO_V1(tle_apogee); +PG_FUNCTION_INFO_V1(tle_intl_desig); + +#define RAD_TO_DEG (180.0 / M_PI) +#define RAD_MIN_TO_REV_DAY (1440.0 / (2.0 * M_PI)) + +/* + * Copy parsed fields from sat_code's tle_t into our pg_tle. + * The sat_code parser stores angles in radians and mean motion in + * radians/minute, so no conversion is needed here. + */ +static void +tle_t_to_pg_tle(const tle_t *sat, pg_tle *tle) +{ + tle->epoch = sat->epoch; + tle->inclination = sat->xincl; + tle->raan = sat->xnodeo; + tle->eccentricity = sat->eo; + tle->arg_perigee = sat->omegao; + tle->mean_anomaly = sat->xmo; + tle->mean_motion = sat->xno; + tle->mean_motion_dot = sat->xndt2o; + tle->mean_motion_ddot = sat->xndd6o; + tle->bstar = sat->bstar; + tle->norad_id = sat->norad_number; + tle->elset_num = sat->bulletin_number; + tle->rev_num = sat->revolution_number; + tle->classification = sat->classification; + tle->ephemeris_type = sat->ephemeris_type; + memcpy(tle->intl_desig, sat->intl_desig, 9); + tle->_pad = '\0'; +} + +/* + * Reverse: pg_tle back to sat_code's tle_t for write_elements_in_tle_format. + */ +static void +pg_tle_to_tle_t(const pg_tle *tle, tle_t *sat) +{ + sat->epoch = tle->epoch; + sat->xincl = tle->inclination; + sat->xnodeo = tle->raan; + sat->eo = tle->eccentricity; + sat->omegao = tle->arg_perigee; + sat->xmo = tle->mean_anomaly; + sat->xno = tle->mean_motion; + sat->xndt2o = tle->mean_motion_dot; + sat->xndd6o = tle->mean_motion_ddot; + sat->bstar = tle->bstar; + sat->norad_number = tle->norad_id; + sat->bulletin_number = tle->elset_num; + sat->revolution_number = tle->rev_num; + sat->classification = tle->classification; + sat->ephemeris_type = tle->ephemeris_type; + memcpy(sat->intl_desig, tle->intl_desig, 9); +} + + +/* + * tle_in -- parse a two-line element set from text + * + * Expects two 69-character lines separated by '\n'. + */ +Datum +tle_in(PG_FUNCTION_ARGS) +{ + char *str = PG_GETARG_CSTRING(0); + pg_tle *result; + tle_t sat; + char *newline; + char *line1; + char *line2; + int parse_rc; + + /* Split on first newline */ + newline = strchr(str, '\n'); + if (newline == NULL) + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid TLE: expected two lines separated by newline"))); + + /* + * Make a mutable copy of line1 so we can null-terminate it without + * modifying the caller's input buffer. + */ + line1 = pnstrdup(str, newline - str); + line2 = newline + 1; + + /* Strip trailing whitespace from line2 */ + { + size_t len2 = strlen(line2); + char *line2_copy = pnstrdup(line2, len2); + + while (len2 > 0 && (line2_copy[len2 - 1] == '\n' || + line2_copy[len2 - 1] == '\r' || + line2_copy[len2 - 1] == ' ')) + line2_copy[--len2] = '\0'; + + memset(&sat, 0, sizeof(tle_t)); + parse_rc = parse_elements(line1, line2_copy, &sat); + + pfree(line2_copy); + } + + pfree(line1); + + if (parse_rc < 0) + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid TLE: parse_elements returned %d", parse_rc))); + + result = (pg_tle *) palloc0(sizeof(pg_tle)); + tle_t_to_pg_tle(&sat, result); + + PG_RETURN_POINTER(result); +} + + +/* + * tle_out -- format a TLE as the standard two-line text representation + * + * Reconstructs the canonical 69+69 character format via sat_code's + * write_elements_in_tle_format(). The output has lines separated by + * a newline and is suitable for round-tripping through tle_in(). + * + * write_elements_in_tle_format() writes two 71-char lines (each ending + * CR+LF+NUL internally). We strip the trailing CR from each line and + * join them with a single '\n' for PostgreSQL text output. + */ +Datum +tle_out(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + tle_t sat; + char buff[200]; + char *result; + char *line1_end; + char *line2_start; + size_t line1_len; + size_t line2_len; + + pg_tle_to_tle_t(tle, &sat); + memset(buff, 0, sizeof(buff)); + write_elements_in_tle_format(buff, &sat); + + /* + * sat_code writes: line1 (69 chars) + checksum + CR + LF + NUL + * line2 (69 chars) + checksum + CR + LF + NUL + * Total is two 71-byte strings placed back-to-back. + * Find the boundary and strip trailing CR/LF from each line. + */ + line1_end = strchr(buff, '\n'); + if (line1_end == NULL) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("write_elements_in_tle_format produced malformed output"))); + + /* Back up over CR if present */ + line1_len = line1_end - buff; + if (line1_len > 0 && buff[line1_len - 1] == '\r') + line1_len--; + + /* line2 starts after the LF (and possibly NUL) */ + line2_start = line1_end + 1; + if (*line2_start == '\0') + line2_start++; + + line2_len = strlen(line2_start); + while (line2_len > 0 && (line2_start[line2_len - 1] == '\n' || + line2_start[line2_len - 1] == '\r')) + line2_len--; + + /* Assemble: line1 + '\n' + line2 */ + result = palloc(line1_len + 1 + line2_len + 1); + memcpy(result, buff, line1_len); + result[line1_len] = '\n'; + memcpy(result + line1_len + 1, line2_start, line2_len); + result[line1_len + 1 + line2_len] = '\0'; + + PG_RETURN_CSTRING(result); +} + + +/* + * tle_recv -- binary input + */ +Datum +tle_recv(PG_FUNCTION_ARGS) +{ + StringInfo buf = (StringInfo) PG_GETARG_POINTER(0); + pg_tle *result = (pg_tle *) palloc0(sizeof(pg_tle)); + + result->epoch = pq_getmsgfloat8(buf); + result->inclination = pq_getmsgfloat8(buf); + result->raan = pq_getmsgfloat8(buf); + result->eccentricity = pq_getmsgfloat8(buf); + result->arg_perigee = pq_getmsgfloat8(buf); + result->mean_anomaly = pq_getmsgfloat8(buf); + result->mean_motion = pq_getmsgfloat8(buf); + result->mean_motion_dot = pq_getmsgfloat8(buf); + result->mean_motion_ddot = pq_getmsgfloat8(buf); + result->bstar = pq_getmsgfloat8(buf); + result->norad_id = pq_getmsgint(buf, 4); + result->elset_num = pq_getmsgint(buf, 4); + result->rev_num = pq_getmsgint(buf, 4); + result->classification = pq_getmsgbyte(buf); + result->ephemeris_type = pq_getmsgbyte(buf); + memcpy(result->intl_desig, pq_getmsgbytes(buf, 9), 9); + result->_pad = '\0'; + + PG_RETURN_POINTER(result); +} + + +/* + * tle_send -- binary output + */ +Datum +tle_send(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + StringInfoData buf; + + pq_begintypsend(&buf); + + pq_sendfloat8(&buf, tle->epoch); + pq_sendfloat8(&buf, tle->inclination); + pq_sendfloat8(&buf, tle->raan); + pq_sendfloat8(&buf, tle->eccentricity); + pq_sendfloat8(&buf, tle->arg_perigee); + pq_sendfloat8(&buf, tle->mean_anomaly); + pq_sendfloat8(&buf, tle->mean_motion); + pq_sendfloat8(&buf, tle->mean_motion_dot); + pq_sendfloat8(&buf, tle->mean_motion_ddot); + pq_sendfloat8(&buf, tle->bstar); + pq_sendint32(&buf, tle->norad_id); + pq_sendint32(&buf, tle->elset_num); + pq_sendint32(&buf, tle->rev_num); + pq_sendbyte(&buf, tle->classification); + pq_sendbyte(&buf, tle->ephemeris_type); + pq_sendbytes(&buf, tle->intl_desig, 9); + + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); +} + + +/* --- Accessor functions --- */ + +Datum +tle_epoch(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + + PG_RETURN_FLOAT8(tle->epoch); +} + +Datum +tle_norad_id(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + + PG_RETURN_INT32(tle->norad_id); +} + +Datum +tle_inclination(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + + PG_RETURN_FLOAT8(tle->inclination * RAD_TO_DEG); +} + +Datum +tle_eccentricity(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + + PG_RETURN_FLOAT8(tle->eccentricity); +} + +Datum +tle_raan(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + + PG_RETURN_FLOAT8(tle->raan * RAD_TO_DEG); +} + +Datum +tle_arg_perigee(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + + PG_RETURN_FLOAT8(tle->arg_perigee * RAD_TO_DEG); +} + +Datum +tle_mean_anomaly(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + + PG_RETURN_FLOAT8(tle->mean_anomaly * RAD_TO_DEG); +} + +/* + * Mean motion in revolutions/day. + * Internally stored as radians/minute. + */ +Datum +tle_mean_motion(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + + PG_RETURN_FLOAT8(tle->mean_motion * RAD_MIN_TO_REV_DAY); +} + +Datum +tle_bstar(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + + PG_RETURN_FLOAT8(tle->bstar); +} + +/* + * Orbital period in minutes. + * period = 2*pi / mean_motion, where mean_motion is radians/minute. + */ +Datum +tle_period(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + double period; + + if (tle->mean_motion <= 0.0) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("mean motion must be positive"))); + + period = (2.0 * M_PI) / tle->mean_motion; + PG_RETURN_FLOAT8(period); +} + +/* + * TLE age in days relative to a given timestamptz. + * Positive = TLE is older than the timestamp (past epoch). + * Negative = TLE epoch is in the future relative to the timestamp. + */ +Datum +tle_age(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + int64 ts = PG_GETARG_INT64(1); + double jd = timestamptz_to_jd(ts); + double age_days = jd - tle->epoch; + + PG_RETURN_FLOAT8(age_days); +} + +/* + * Perigee altitude in km above the WGS-72 ellipsoid. + * + * Semi-major axis from Kepler's third law using WGS-72 KE: + * a = (KE / n)^(2/3) [earth radii] + * perigee_alt = a * (1 - e) * AE_km - AE_km + */ +Datum +tle_perigee(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + double a_er; /* semi-major axis in earth radii */ + double alt; + + if (tle->mean_motion <= 0.0) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("mean motion must be positive"))); + + a_er = pow(WGS72_KE / tle->mean_motion, 2.0 / 3.0); + alt = a_er * (1.0 - tle->eccentricity) * WGS72_AE - WGS72_AE; + PG_RETURN_FLOAT8(alt); +} + +/* + * Apogee altitude in km above the WGS-72 ellipsoid. + * + * apogee_alt = a * (1 + e) * AE_km - AE_km + */ +Datum +tle_apogee(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + double a_er; + double alt; + + if (tle->mean_motion <= 0.0) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("mean motion must be positive"))); + + a_er = pow(WGS72_KE / tle->mean_motion, 2.0 / 3.0); + alt = a_er * (1.0 + tle->eccentricity) * WGS72_AE - WGS72_AE; + PG_RETURN_FLOAT8(alt); +} + +/* + * International designator as a text datum. + */ +Datum +tle_intl_desig(PG_FUNCTION_ARGS) +{ + pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0); + + PG_RETURN_TEXT_P(cstring_to_text(tle->intl_desig)); +} diff --git a/src/types.h b/src/types.h new file mode 100644 index 0000000..a8d4a37 --- /dev/null +++ b/src/types.h @@ -0,0 +1,190 @@ +/* + * types.h -- Shared type definitions for pg_orbit + * + * All orbital mechanics types stored in PostgreSQL tuples. + * Positions in TEME frame, propagated with WGS-72. + * Coordinate output converted to WGS-84. + */ + +#ifndef PG_ORBIT_TYPES_H +#define PG_ORBIT_TYPES_H + +#include "postgres.h" + +/* + * WGS-72 constants (for SGP4 propagation ONLY) + * Source: Hoots & Roehrich, "Models for Propagation of NORAD Element Sets" + * Spacetrack Report No. 3, 1980 + */ +#define WGS72_MU 398600.8 /* km^3/s^2 */ +#define WGS72_AE 6378.135 /* km (equatorial radius) */ +#define WGS72_J2 0.001082616 +#define WGS72_J3 -0.00000253881 +#define WGS72_J4 -0.00000165597 +#define WGS72_KE 0.0743669161331734132 /* (min)^(-1) */ +#define WGS72_XPDOTP (1440.0 / (2.0 * M_PI)) /* min/rev */ + +/* + * WGS-84 constants (for coordinate output ONLY) + * Source: NIMA TR8350.2, "Department of Defense World Geodetic System 1984" + */ +#define WGS84_A 6378.137 /* km (equatorial radius) */ +#define WGS84_F (1.0 / 298.257223563) +#define WGS84_E2 (WGS84_F * (2.0 - WGS84_F)) +#define WGS84_A_M 6378137.0 /* meters */ + +/* + * Julian date constants + */ +#define J2000_JD 2451545.0 /* 2000 Jan 1.5 TT */ +#define JD_UNIX_EPOCH 2440587.5 /* 1970 Jan 1.0 UTC */ +#define MJD_OFFSET 2400000.5 + +/* + * TLE type -- parsed mean orbital elements + * + * Stored as a fixed-size struct, not raw text. + * The epoch is a Julian date (UTC). + * Angular elements are in radians. + * Mean motion is in radians/minute. + */ +typedef struct pg_tle +{ + /* Epoch as Julian date, UTC */ + double epoch; + + /* Mean orbital elements (radians, radians/min) */ + double inclination; /* xincl */ + double raan; /* xnodeo - right ascension of ascending node */ + double eccentricity; /* eo */ + double arg_perigee; /* omegao */ + double mean_anomaly; /* xmo */ + double mean_motion; /* xno - radians/minute */ + + /* Drag terms */ + double mean_motion_dot; /* xndt2o - first derivative / 2 */ + double mean_motion_ddot; /* xndd6o - second derivative / 6 */ + double bstar; /* drag coefficient */ + + /* Identification */ + int32 norad_id; + int32 elset_num; /* bulletin_number in sat_code */ + int32 rev_num; /* revolution number at epoch */ + + /* Metadata */ + char classification; /* U = unclassified */ + char ephemeris_type; /* 0 = SGP4/SDP4 default */ + char intl_desig[9]; /* international designator, null-terminated */ + char _pad; /* alignment */ +} pg_tle; + +/* Size: 11 doubles (88 bytes) + 3 int32 (12 bytes) + 12 chars = 112 bytes */ + + +/* + * ECI position -- True Equator Mean Equinox (TEME) frame + * + * Position in km, velocity in km/s. + * SGP4 outputs velocity in km/min; we convert to km/s. + */ +typedef struct pg_eci +{ + double x, y, z; /* position, km */ + double vx, vy, vz; /* velocity, km/s */ +} pg_eci; + + +/* + * Geodetic coordinates -- WGS-84 ellipsoid + * + * Latitude and longitude in radians (internal), degrees (display). + * Altitude above WGS-84 ellipsoid in km. + */ +typedef struct pg_geodetic +{ + double lat; /* radians, -pi/2 to +pi/2 */ + double lon; /* radians, -pi to +pi */ + double alt; /* km above WGS-84 ellipsoid */ +} pg_geodetic; + + +/* + * Topocentric coordinates -- observer-relative + * + * Azimuth: 0=N, 90=E, 180=S, 270=W (degrees for display, radians internal) + * Elevation: 0=horizon, 90=zenith + * Range in km, range rate in km/s (positive = receding) + */ +typedef struct pg_topocentric +{ + double azimuth; /* radians */ + double elevation; /* radians */ + double range_km; + double range_rate; /* km/s */ +} pg_topocentric; + + +/* + * Observer location -- ground station + * + * Latitude/longitude in radians (internal). + * Altitude in meters above WGS-84 ellipsoid. + */ +typedef struct pg_observer +{ + double lat; /* radians */ + double lon; /* radians */ + double alt_m; /* meters above WGS-84 */ +} pg_observer; + + +/* + * Pass event -- satellite visibility window + * + * Times are PostgreSQL TimestampTz (microseconds since 2000-01-01). + * Elevation in degrees, azimuth in degrees. + */ +typedef struct pg_pass_event +{ + int64 aos_time; /* acquisition of signal */ + int64 max_el_time; /* maximum elevation */ + int64 los_time; /* loss of signal */ + double max_elevation; /* degrees */ + double aos_azimuth; /* degrees, 0=N */ + double los_azimuth; /* degrees, 0=N */ +} pg_pass_event; + + +/* + * Utility: convert PostgreSQL TimestampTz to Julian date + * + * PG epoch: 2000-01-01 00:00:00 UTC = JD 2451544.5 + * TimestampTz is microseconds since PG epoch. + */ +#define PG_EPOCH_JD 2451544.5 +#ifndef USECS_PER_DAY +#define USECS_PER_DAY 86400000000LL +#endif + +static inline double +timestamptz_to_jd(int64 ts) +{ + return PG_EPOCH_JD + ((double)ts / (double)USECS_PER_DAY); +} + +static inline int64 +jd_to_timestamptz(double jd) +{ + return (int64)((jd - PG_EPOCH_JD) * (double)USECS_PER_DAY); +} + +/* + * Utility: minutes since TLE epoch + */ +static inline double +jd_to_minutes_since_epoch(double jd, double tle_epoch_jd) +{ + return (jd - tle_epoch_jd) * 1440.0; +} + +#endif /* PG_ORBIT_TYPES_H */ diff --git a/test/expected/coord_transforms.out b/test/expected/coord_transforms.out new file mode 100644 index 0000000..a87e3f9 --- /dev/null +++ b/test/expected/coord_transforms.out @@ -0,0 +1,110 @@ +-- Test coordinate transform functions +CREATE EXTENSION IF NOT EXISTS pg_orbit; +NOTICE: extension "pg_orbit" already exists, skipping +-- Subsatellite point at epoch +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + round(geodetic_lat(subsatellite_point(t, '2024-01-01 12:00:00+00'))::numeric, 2) AS lat, + round(geodetic_lon(subsatellite_point(t, '2024-01-01 12:00:00+00'))::numeric, 2) AS lon, + round(geodetic_alt(subsatellite_point(t, '2024-01-01 12:00:00+00'))::numeric, 0) AS alt_km +FROM iss; + lat | lon | alt_km +-------+-------+-------- + 51.75 | 21.48 | 420 +(1 row) + +-- Subsatellite point 30 minutes later +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + round(geodetic_lat(subsatellite_point(t, '2024-01-01 12:30:00+00'))::numeric, 2) AS lat_30m, + round(geodetic_lon(subsatellite_point(t, '2024-01-01 12:30:00+00'))::numeric, 2) AS lon_30m, + round(geodetic_alt(subsatellite_point(t, '2024-01-01 12:30:00+00'))::numeric, 0) AS alt_30m +FROM iss; + lat_30m | lon_30m | alt_30m +---------+---------+--------- + -21.98 | 119.17 | 427 +(1 row) + +-- eci_to_geodetic: propagate ISS then convert +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + round(geodetic_lat(eci_to_geodetic(sgp4_propagate(t, '2024-01-01 12:00:00+00'), '2024-01-01 12:00:00+00'))::numeric, 2) AS eci_lat, + round(geodetic_lon(eci_to_geodetic(sgp4_propagate(t, '2024-01-01 12:00:00+00'), '2024-01-01 12:00:00+00'))::numeric, 2) AS eci_lon, + round(geodetic_alt(eci_to_geodetic(sgp4_propagate(t, '2024-01-01 12:00:00+00'), '2024-01-01 12:00:00+00'))::numeric, 0) AS eci_alt +FROM iss; + eci_lat | eci_lon | eci_alt +---------+---------+--------- + 51.75 | 21.48 | 420 +(1 row) + +-- Topocentric from Boulder, CO +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + round(topo_azimuth(eci_to_topocentric(sgp4_propagate(t, '2024-01-01 12:00:00+00'), '40.0N 105.3W 1655m'::observer, '2024-01-01 12:00:00+00'))::numeric, 1) AS az_deg, + round(topo_elevation(eci_to_topocentric(sgp4_propagate(t, '2024-01-01 12:00:00+00'), '40.0N 105.3W 1655m'::observer, '2024-01-01 12:00:00+00'))::numeric, 1) AS el_deg, + round(topo_range(eci_to_topocentric(sgp4_propagate(t, '2024-01-01 12:00:00+00'), '40.0N 105.3W 1655m'::observer, '2024-01-01 12:00:00+00'))::numeric, 0) AS range_km +FROM iss; + az_deg | el_deg | range_km +--------+--------+---------- + 30.6 | -36.4 | 8245 +(1 row) + +-- Ground track: 5 points over 40 minutes +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + count(*) AS track_points, + round(min(lat)::numeric, 0) AS min_lat, + round(max(lat)::numeric, 0) AS max_lat +FROM iss, ground_track(t, '2024-01-01 12:00:00+00', '2024-01-01 12:40:00+00', '10 minutes'); + track_points | min_lat | max_lat +--------------+---------+--------- + 5 | -46 | 52 +(1 row) + +-- Geodetic type I/O round-trip +SELECT geodetic_lat('(51.75, 21.48, 420.0)'::geodetic) AS lat, + geodetic_lon('(51.75, 21.48, 420.0)'::geodetic) AS lon, + geodetic_alt('(51.75, 21.48, 420.0)'::geodetic) AS alt; + lat | lon | alt +-------+-------+----- + 51.75 | 21.48 | 420 +(1 row) + +-- Topocentric type I/O round-trip +SELECT round(topo_azimuth('(45.0, 30.0, 500.0, -2.5)'::topocentric)::numeric, 1) AS az, + round(topo_elevation('(45.0, 30.0, 500.0, -2.5)'::topocentric)::numeric, 1) AS el, + round(topo_range('(45.0, 30.0, 500.0, -2.5)'::topocentric)::numeric, 1) AS rng, + round(topo_range_rate('(45.0, 30.0, 500.0, -2.5)'::topocentric)::numeric, 1) AS rr; + az | el | rng | rr +------+------+-------+------ + 45.0 | 30.0 | 500.0 | -2.5 +(1 row) + +-- ISS latitude should stay within inclination bounds (51.64 deg) +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + abs(geodetic_lat(subsatellite_point(t, '2024-01-01 12:00:00+00'))) <= 52.0 AS within_inc +FROM iss; + within_inc +------------ + t +(1 row) + diff --git a/test/expected/gist_index.out b/test/expected/gist_index.out new file mode 100644 index 0000000..87d0505 --- /dev/null +++ b/test/expected/gist_index.out @@ -0,0 +1,83 @@ +-- Test GiST index and operators +CREATE EXTENSION IF NOT EXISTS pg_orbit; +NOTICE: extension "pg_orbit" already exists, skipping +-- Create test table with mixed orbit types +CREATE TABLE test_orbits ( + id serial, + name text, + tle tle +); +-- ISS (LEO, ~400km) +INSERT INTO test_orbits (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) +INSERT INTO test_orbits (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) +INSERT INTO test_orbits (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'); +-- Create GiST index +CREATE INDEX test_orbits_gist ON test_orbits USING gist (tle); +-- Overlap: ISS && Hubble = false (different altitude bands: ~400km vs ~540km) +SELECT a.name AS sat_a, b.name AS sat_b, a.tle && b.tle AS overlaps +FROM test_orbits a, test_orbits b +WHERE a.id < b.id +ORDER BY a.name, b.name; + sat_a | sat_b | overlaps +--------+---------+---------- + Hubble | GPS-IIR | f + ISS | GPS-IIR | f + ISS | Hubble | f +(3 rows) + +-- Altitude distance between different orbit regimes +SELECT a.name AS sat_a, b.name AS sat_b, + round((a.tle <-> b.tle)::numeric, 0) AS alt_dist_km +FROM test_orbits a, test_orbits b +WHERE a.id < b.id +ORDER BY a.name, b.name; + sat_a | sat_b | alt_dist_km +--------+---------+------------- + Hubble | GPS-IIR | 19332 + ISS | GPS-IIR | 19451 + ISS | Hubble | 115 +(3 rows) + +-- GiST index scan: find all LEO sats (overlap with ISS) +SET enable_seqscan = off; +SELECT name FROM test_orbits +WHERE tle && (SELECT tle FROM test_orbits WHERE name = 'ISS') +ORDER BY name; + name +------ + ISS +(1 row) + +RESET enable_seqscan; +-- Nearest-neighbor via GiST: order by altitude distance to ISS +SET enable_seqscan = off; +SELECT name, round((tle <-> (SELECT tle FROM test_orbits WHERE name = 'ISS'))::numeric, 0) AS dist +FROM test_orbits +WHERE name != 'ISS' +ORDER BY tle <-> (SELECT tle FROM test_orbits WHERE name = 'ISS'); + name | dist +---------+------- + Hubble | 115 + GPS-IIR | 19451 +(2 rows) + +RESET enable_seqscan; +-- Self-overlap is always true +SELECT name, tle && tle AS self_overlap FROM test_orbits ORDER BY name; + name | self_overlap +---------+-------------- + GPS-IIR | t + Hubble | t + ISS | t +(3 rows) + +-- Cleanup +DROP TABLE test_orbits; diff --git a/test/expected/pass_prediction.out b/test/expected/pass_prediction.out new file mode 100644 index 0000000..04792bd --- /dev/null +++ b/test/expected/pass_prediction.out @@ -0,0 +1,91 @@ +-- Test pass prediction functions +CREATE EXTENSION IF NOT EXISTS pg_orbit; +NOTICE: extension "pg_orbit" already exists, skipping +-- Predict ISS passes over Boulder in 24-hour window +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT count(*) AS pass_count +FROM iss, predict_passes(t, '40.0N 105.3W 1655m'::observer, + '2024-01-01 12:00:00+00', '2024-01-02 12:00:00+00'); + pass_count +------------ + 7 +(1 row) + +-- Pass details: max elevation should be positive, duration reasonable +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + round(pass_max_elevation(p)::numeric, 1) AS max_el, + pass_aos_time(p) < pass_max_el_time(p) AS aos_before_max, + pass_max_el_time(p) < pass_los_time(p) AS max_before_los, + pass_aos_azimuth(p) >= 0 AND pass_aos_azimuth(p) <= 360 AS az_valid, + extract(epoch from pass_duration(p)) BETWEEN 60 AND 900 AS duration_ok +FROM iss, predict_passes(t, '40.0N 105.3W 1655m'::observer, + '2024-01-01 12:00:00+00', '2024-01-02 12:00:00+00') AS p +LIMIT 3; + max_el | aos_before_max | max_before_los | az_valid | duration_ok +--------+----------------+----------------+----------+------------- + 3.3 | t | t | t | t + 46.4 | t | t | t | t + 24.4 | t | t | t | t +(3 rows) + +-- next_pass returns same first pass as predict_passes +WITH iss AS ( + SELECT '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'::tle AS t +), +np AS ( + SELECT next_pass(t, '40.0N 105.3W 1655m'::observer, '2024-01-01 12:00:00+00') AS p FROM iss +), +pp AS ( + SELECT p AS pass FROM iss, predict_passes(t, '40.0N 105.3W 1655m'::observer, + '2024-01-01 12:00:00+00', '2024-01-02 12:00:00+00') AS p + LIMIT 1 +) +SELECT + pass_aos_time(np.p) = pass_aos_time(pp.pass) AS same_aos +FROM np, pp; + same_aos +---------- + t +(1 row) + +-- pass_visible should be true for ISS over Boulder in a 24h window +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT pass_visible(t, '40.0N 105.3W 1655m'::observer, + '2024-01-01 12:00:00+00', '2024-01-02 12:00:00+00') AS visible +FROM iss; + visible +--------- + t +(1 row) + +-- Minimum elevation filter should reduce pass count +WITH iss AS ( + SELECT '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'::tle AS t +), +all_passes AS ( + SELECT count(*) AS total FROM iss, predict_passes(t, '40.0N 105.3W 1655m'::observer, + '2024-01-01 12:00:00+00', '2024-01-02 12:00:00+00') +), +high_passes AS ( + SELECT count(*) AS high FROM iss, predict_passes(t, '40.0N 105.3W 1655m'::observer, + '2024-01-01 12:00:00+00', '2024-01-02 12:00:00+00', 30.0) +) +SELECT high <= total AS filter_works +FROM all_passes, high_passes; + filter_works +-------------- + t +(1 row) + diff --git a/test/expected/sgp4_propagate.out b/test/expected/sgp4_propagate.out new file mode 100644 index 0000000..be19def --- /dev/null +++ b/test/expected/sgp4_propagate.out @@ -0,0 +1,92 @@ +-- Test SGP4/SDP4 propagation functions +CREATE EXTENSION IF NOT EXISTS pg_orbit; +NOTICE: extension "pg_orbit" already exists, skipping +-- ISS TLE (LEO, near-earth -> SGP4) +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + round(eci_x(sgp4_propagate(t, '2024-01-01 12:00:00+00'))::numeric, 1) AS x_km, + round(eci_y(sgp4_propagate(t, '2024-01-01 12:00:00+00'))::numeric, 1) AS y_km, + round(eci_z(sgp4_propagate(t, '2024-01-01 12:00:00+00'))::numeric, 1) AS z_km, + round(eci_speed(sgp4_propagate(t, '2024-01-01 12:00:00+00'))::numeric, 3) AS speed_kms, + round(eci_altitude(sgp4_propagate(t, '2024-01-01 12:00:00+00'))::numeric, 0) AS alt_km +FROM iss; + x_km | y_km | z_km | speed_kms | alt_km +--------+---------+--------+-----------+-------- + 2242.3 | -3571.3 | 5315.9 | 7.667 | 407 +(1 row) + +-- Propagation 1 hour after epoch +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + round(eci_speed(sgp4_propagate(t, '2024-01-01 13:00:00+00'))::numeric, 3) AS speed_1h, + round(eci_altitude(sgp4_propagate(t, '2024-01-01 13:00:00+00'))::numeric, 0) AS alt_1h +FROM iss; + speed_1h | alt_1h +----------+-------- + 7.659 | 418 +(1 row) + +-- GPS satellite (MEO, deep-space -> SDP4) +WITH gps AS ( + SELECT '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'::tle AS t +) +SELECT + round(eci_speed(sgp4_propagate(t, '2024-01-01 12:00:00+00'))::numeric, 3) AS gps_speed, + round(eci_altitude(sgp4_propagate(t, '2024-01-01 12:00:00+00'))::numeric, 0) AS gps_alt +FROM gps; + gps_speed | gps_alt +-----------+--------- + 3.903 | 19988 +(1 row) + +-- Propagation series: ISS, 10 minute steps over 1 orbit (~93 min) +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + count(*) AS num_points, + round(min(sqrt(x*x + y*y + z*z) - 6378.135)::numeric, 0) AS min_alt, + round(max(sqrt(x*x + y*y + z*z) - 6378.135)::numeric, 0) AS max_alt +FROM iss, sgp4_propagate_series(t, '2024-01-01 12:00:00+00', '2024-01-01 13:33:00+00', '10 minutes'); + num_points | min_alt | max_alt +------------+---------+--------- + 10 | 407 | 425 +(1 row) + +-- Distance between ISS and GPS at epoch +WITH sats AS ( + SELECT '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'::tle AS iss, + '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'::tle AS gps +) +SELECT + round(tle_distance(iss, gps, '2024-01-01 12:00:00+00')::numeric, 0) AS dist_km, + round(tle_distance(iss, iss, '2024-01-01 12:00:00+00')::numeric, 6) AS self_dist +FROM sats; + dist_km | self_dist +---------+----------- + 22768 | 0.000000 +(1 row) + +-- Distance to self should be zero +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + tle_distance(t, t, '2024-01-01 12:00:00+00') = 0.0 AS self_is_zero +FROM iss; + self_is_zero +-------------- + t +(1 row) + diff --git a/test/expected/tle_parse.out b/test/expected/tle_parse.out new file mode 100644 index 0000000..bf3d70d --- /dev/null +++ b/test/expected/tle_parse.out @@ -0,0 +1,116 @@ +-- Test TLE type: parsing, round-trip, accessors +CREATE EXTENSION IF NOT EXISTS pg_orbit; +-- Parse a valid ISS TLE +SELECT '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'::tle IS NOT NULL AS parse_ok; + parse_ok +---------- + t +(1 row) + +-- Accessor functions +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + tle_norad_id(t) AS norad_id, + round(tle_inclination(t)::numeric, 4) AS inc_deg, + round(tle_eccentricity(t)::numeric, 7) AS ecc, + round(tle_mean_motion(t)::numeric, 8) AS mm_rev_day, + round(tle_period(t)::numeric, 2) AS period_min, + round(tle_perigee(t)::numeric, 1) AS perigee_km, + round(tle_apogee(t)::numeric, 1) AS apogee_km, + tle_intl_desig(t) AS cospar +FROM iss; + norad_id | inc_deg | ecc | mm_rev_day | period_min | perigee_km | apogee_km | cospar +----------+---------+-----------+-------------+------------+------------+-----------+---------- + 25544 | 51.6400 | 0.0006703 | 15.50100486 | 92.90 | 411.9 | 421.0 | 98067A +(1 row) + +-- Observer type parsing +SELECT '40.0N 105.3W 1655m'::observer IS NOT NULL AS observer_ok; + observer_ok +------------- + t +(1 row) + +SELECT observer_lat('40.0N 105.3W 1655m'::observer) AS lat, + observer_lon('40.0N 105.3W 1655m'::observer) AS lon; + lat | lon +-----+-------- + 40 | -105.3 +(1 row) + +-- ECI type parsing +SELECT '(1000.0,2000.0,3000.0,1.0,2.0,3.0)'::eci_position IS NOT NULL AS eci_ok; + eci_ok +-------- + t +(1 row) + +SELECT eci_x('(1000.0,2000.0,3000.0,1.0,2.0,3.0)'::eci_position) AS x; + x +------ + 1000 +(1 row) + +-- SGP4 propagation at epoch +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + round(eci_speed(sgp4_propagate(t, '2024-01-01 12:00:00+00'))::numeric, 2) AS speed_kms, + round(eci_altitude(sgp4_propagate(t, '2024-01-01 12:00:00+00'))::numeric, 0) AS alt_km +FROM iss; + speed_kms | alt_km +-----------+-------- + 7.67 | 407 +(1 row) + +-- Altitude overlap operator +WITH sats AS ( + SELECT '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'::tle AS iss, + '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'::tle AS gps +) +SELECT + iss && iss AS same_overlap, + iss && gps AS different_no_overlap +FROM sats; + same_overlap | different_no_overlap +--------------+---------------------- + t | f +(1 row) + +-- Subsatellite point +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + round(geodetic_lat(subsatellite_point(t, '2024-01-01 12:00:00+00'))::numeric, 2) AS lat, + round(geodetic_alt(subsatellite_point(t, '2024-01-01 12:00:00+00'))::numeric, 0) AS alt_km +FROM iss; + lat | alt_km +-------+-------- + 51.75 | 420 +(1 row) + +-- GiST index creation +CREATE TABLE test_sats (id serial, tle tle); +INSERT INTO test_sats (tle) VALUES ( + '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'); +CREATE INDEX test_sats_gist ON test_sats USING gist (tle); +SELECT count(*) FROM test_sats WHERE tle && ( + '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'::tle); + count +------- + 1 +(1 row) + +DROP TABLE test_sats; diff --git a/test/sql/coord_transforms.sql b/test/sql/coord_transforms.sql new file mode 100644 index 0000000..3c8c5a7 --- /dev/null +++ b/test/sql/coord_transforms.sql @@ -0,0 +1,77 @@ +-- Test coordinate transform functions +CREATE EXTENSION IF NOT EXISTS pg_orbit; + +-- Subsatellite point at epoch +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + round(geodetic_lat(subsatellite_point(t, '2024-01-01 12:00:00+00'))::numeric, 2) AS lat, + round(geodetic_lon(subsatellite_point(t, '2024-01-01 12:00:00+00'))::numeric, 2) AS lon, + round(geodetic_alt(subsatellite_point(t, '2024-01-01 12:00:00+00'))::numeric, 0) AS alt_km +FROM iss; + +-- Subsatellite point 30 minutes later +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + round(geodetic_lat(subsatellite_point(t, '2024-01-01 12:30:00+00'))::numeric, 2) AS lat_30m, + round(geodetic_lon(subsatellite_point(t, '2024-01-01 12:30:00+00'))::numeric, 2) AS lon_30m, + round(geodetic_alt(subsatellite_point(t, '2024-01-01 12:30:00+00'))::numeric, 0) AS alt_30m +FROM iss; + +-- eci_to_geodetic: propagate ISS then convert +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + round(geodetic_lat(eci_to_geodetic(sgp4_propagate(t, '2024-01-01 12:00:00+00'), '2024-01-01 12:00:00+00'))::numeric, 2) AS eci_lat, + round(geodetic_lon(eci_to_geodetic(sgp4_propagate(t, '2024-01-01 12:00:00+00'), '2024-01-01 12:00:00+00'))::numeric, 2) AS eci_lon, + round(geodetic_alt(eci_to_geodetic(sgp4_propagate(t, '2024-01-01 12:00:00+00'), '2024-01-01 12:00:00+00'))::numeric, 0) AS eci_alt +FROM iss; + +-- Topocentric from Boulder, CO +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + round(topo_azimuth(eci_to_topocentric(sgp4_propagate(t, '2024-01-01 12:00:00+00'), '40.0N 105.3W 1655m'::observer, '2024-01-01 12:00:00+00'))::numeric, 1) AS az_deg, + round(topo_elevation(eci_to_topocentric(sgp4_propagate(t, '2024-01-01 12:00:00+00'), '40.0N 105.3W 1655m'::observer, '2024-01-01 12:00:00+00'))::numeric, 1) AS el_deg, + round(topo_range(eci_to_topocentric(sgp4_propagate(t, '2024-01-01 12:00:00+00'), '40.0N 105.3W 1655m'::observer, '2024-01-01 12:00:00+00'))::numeric, 0) AS range_km +FROM iss; + +-- Ground track: 5 points over 40 minutes +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + count(*) AS track_points, + round(min(lat)::numeric, 0) AS min_lat, + round(max(lat)::numeric, 0) AS max_lat +FROM iss, ground_track(t, '2024-01-01 12:00:00+00', '2024-01-01 12:40:00+00', '10 minutes'); + +-- Geodetic type I/O round-trip +SELECT geodetic_lat('(51.75, 21.48, 420.0)'::geodetic) AS lat, + geodetic_lon('(51.75, 21.48, 420.0)'::geodetic) AS lon, + geodetic_alt('(51.75, 21.48, 420.0)'::geodetic) AS alt; + +-- Topocentric type I/O round-trip +SELECT round(topo_azimuth('(45.0, 30.0, 500.0, -2.5)'::topocentric)::numeric, 1) AS az, + round(topo_elevation('(45.0, 30.0, 500.0, -2.5)'::topocentric)::numeric, 1) AS el, + round(topo_range('(45.0, 30.0, 500.0, -2.5)'::topocentric)::numeric, 1) AS rng, + round(topo_range_rate('(45.0, 30.0, 500.0, -2.5)'::topocentric)::numeric, 1) AS rr; + +-- ISS latitude should stay within inclination bounds (51.64 deg) +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + abs(geodetic_lat(subsatellite_point(t, '2024-01-01 12:00:00+00'))) <= 52.0 AS within_inc +FROM iss; diff --git a/test/sql/gist_index.sql b/test/sql/gist_index.sql new file mode 100644 index 0000000..09f50db --- /dev/null +++ b/test/sql/gist_index.sql @@ -0,0 +1,61 @@ +-- Test GiST index and operators +CREATE EXTENSION IF NOT EXISTS pg_orbit; + +-- Create test table with mixed orbit types +CREATE TABLE test_orbits ( + id serial, + name text, + tle tle +); + +-- ISS (LEO, ~400km) +INSERT INTO test_orbits (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) +INSERT INTO test_orbits (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) +INSERT INTO test_orbits (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'); + +-- Create GiST index +CREATE INDEX test_orbits_gist ON test_orbits USING gist (tle); + +-- Overlap: ISS && Hubble = false (different altitude bands: ~400km vs ~540km) +SELECT a.name AS sat_a, b.name AS sat_b, a.tle && b.tle AS overlaps +FROM test_orbits a, test_orbits b +WHERE a.id < b.id +ORDER BY a.name, b.name; + +-- Altitude distance between different orbit regimes +SELECT a.name AS sat_a, b.name AS sat_b, + round((a.tle <-> b.tle)::numeric, 0) AS alt_dist_km +FROM test_orbits a, test_orbits b +WHERE a.id < b.id +ORDER BY a.name, b.name; + +-- GiST index scan: find all LEO sats (overlap with ISS) +SET enable_seqscan = off; +SELECT name FROM test_orbits +WHERE tle && (SELECT tle FROM test_orbits WHERE name = 'ISS') +ORDER BY name; +RESET enable_seqscan; + +-- Nearest-neighbor via GiST: order by altitude distance to ISS +SET enable_seqscan = off; +SELECT name, round((tle <-> (SELECT tle FROM test_orbits WHERE name = 'ISS'))::numeric, 0) AS dist +FROM test_orbits +WHERE name != 'ISS' +ORDER BY tle <-> (SELECT tle FROM test_orbits WHERE name = 'ISS'); +RESET enable_seqscan; + +-- Self-overlap is always true +SELECT name, tle && tle AS self_overlap FROM test_orbits ORDER BY name; + +-- Cleanup +DROP TABLE test_orbits; diff --git a/test/sql/pass_prediction.sql b/test/sql/pass_prediction.sql new file mode 100644 index 0000000..ffccd21 --- /dev/null +++ b/test/sql/pass_prediction.sql @@ -0,0 +1,68 @@ +-- Test pass prediction functions +CREATE EXTENSION IF NOT EXISTS pg_orbit; + +-- Predict ISS passes over Boulder in 24-hour window +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT count(*) AS pass_count +FROM iss, predict_passes(t, '40.0N 105.3W 1655m'::observer, + '2024-01-01 12:00:00+00', '2024-01-02 12:00:00+00'); + +-- Pass details: max elevation should be positive, duration reasonable +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + round(pass_max_elevation(p)::numeric, 1) AS max_el, + pass_aos_time(p) < pass_max_el_time(p) AS aos_before_max, + pass_max_el_time(p) < pass_los_time(p) AS max_before_los, + pass_aos_azimuth(p) >= 0 AND pass_aos_azimuth(p) <= 360 AS az_valid, + extract(epoch from pass_duration(p)) BETWEEN 60 AND 900 AS duration_ok +FROM iss, predict_passes(t, '40.0N 105.3W 1655m'::observer, + '2024-01-01 12:00:00+00', '2024-01-02 12:00:00+00') AS p +LIMIT 3; + +-- next_pass returns same first pass as predict_passes +WITH iss AS ( + SELECT '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'::tle AS t +), +np AS ( + SELECT next_pass(t, '40.0N 105.3W 1655m'::observer, '2024-01-01 12:00:00+00') AS p FROM iss +), +pp AS ( + SELECT p AS pass FROM iss, predict_passes(t, '40.0N 105.3W 1655m'::observer, + '2024-01-01 12:00:00+00', '2024-01-02 12:00:00+00') AS p + LIMIT 1 +) +SELECT + pass_aos_time(np.p) = pass_aos_time(pp.pass) AS same_aos +FROM np, pp; + +-- pass_visible should be true for ISS over Boulder in a 24h window +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT pass_visible(t, '40.0N 105.3W 1655m'::observer, + '2024-01-01 12:00:00+00', '2024-01-02 12:00:00+00') AS visible +FROM iss; + +-- Minimum elevation filter should reduce pass count +WITH iss AS ( + SELECT '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'::tle AS t +), +all_passes AS ( + SELECT count(*) AS total FROM iss, predict_passes(t, '40.0N 105.3W 1655m'::observer, + '2024-01-01 12:00:00+00', '2024-01-02 12:00:00+00') +), +high_passes AS ( + SELECT count(*) AS high FROM iss, predict_passes(t, '40.0N 105.3W 1655m'::observer, + '2024-01-01 12:00:00+00', '2024-01-02 12:00:00+00', 30.0) +) +SELECT high <= total AS filter_works +FROM all_passes, high_passes; diff --git a/test/sql/sgp4_propagate.sql b/test/sql/sgp4_propagate.sql new file mode 100644 index 0000000..8ef8a46 --- /dev/null +++ b/test/sql/sgp4_propagate.sql @@ -0,0 +1,67 @@ +-- Test SGP4/SDP4 propagation functions +CREATE EXTENSION IF NOT EXISTS pg_orbit; + +-- ISS TLE (LEO, near-earth -> SGP4) +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + round(eci_x(sgp4_propagate(t, '2024-01-01 12:00:00+00'))::numeric, 1) AS x_km, + round(eci_y(sgp4_propagate(t, '2024-01-01 12:00:00+00'))::numeric, 1) AS y_km, + round(eci_z(sgp4_propagate(t, '2024-01-01 12:00:00+00'))::numeric, 1) AS z_km, + round(eci_speed(sgp4_propagate(t, '2024-01-01 12:00:00+00'))::numeric, 3) AS speed_kms, + round(eci_altitude(sgp4_propagate(t, '2024-01-01 12:00:00+00'))::numeric, 0) AS alt_km +FROM iss; + +-- Propagation 1 hour after epoch +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + round(eci_speed(sgp4_propagate(t, '2024-01-01 13:00:00+00'))::numeric, 3) AS speed_1h, + round(eci_altitude(sgp4_propagate(t, '2024-01-01 13:00:00+00'))::numeric, 0) AS alt_1h +FROM iss; + +-- GPS satellite (MEO, deep-space -> SDP4) +WITH gps AS ( + SELECT '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'::tle AS t +) +SELECT + round(eci_speed(sgp4_propagate(t, '2024-01-01 12:00:00+00'))::numeric, 3) AS gps_speed, + round(eci_altitude(sgp4_propagate(t, '2024-01-01 12:00:00+00'))::numeric, 0) AS gps_alt +FROM gps; + +-- Propagation series: ISS, 10 minute steps over 1 orbit (~93 min) +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + count(*) AS num_points, + round(min(sqrt(x*x + y*y + z*z) - 6378.135)::numeric, 0) AS min_alt, + round(max(sqrt(x*x + y*y + z*z) - 6378.135)::numeric, 0) AS max_alt +FROM iss, sgp4_propagate_series(t, '2024-01-01 12:00:00+00', '2024-01-01 13:33:00+00', '10 minutes'); + +-- Distance between ISS and GPS at epoch +WITH sats AS ( + SELECT '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'::tle AS iss, + '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'::tle AS gps +) +SELECT + round(tle_distance(iss, gps, '2024-01-01 12:00:00+00')::numeric, 0) AS dist_km, + round(tle_distance(iss, iss, '2024-01-01 12:00:00+00')::numeric, 6) AS self_dist +FROM sats; + +-- Distance to self should be zero +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + tle_distance(t, t, '2024-01-01 12:00:00+00') = 0.0 AS self_is_zero +FROM iss; diff --git a/test/sql/tle_parse.sql b/test/sql/tle_parse.sql new file mode 100644 index 0000000..d83fec0 --- /dev/null +++ b/test/sql/tle_parse.sql @@ -0,0 +1,74 @@ +-- Test TLE type: parsing, round-trip, accessors +CREATE EXTENSION IF NOT EXISTS pg_orbit; + +-- Parse a valid ISS TLE +SELECT '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'::tle IS NOT NULL AS parse_ok; + +-- Accessor functions +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + tle_norad_id(t) AS norad_id, + round(tle_inclination(t)::numeric, 4) AS inc_deg, + round(tle_eccentricity(t)::numeric, 7) AS ecc, + round(tle_mean_motion(t)::numeric, 8) AS mm_rev_day, + round(tle_period(t)::numeric, 2) AS period_min, + round(tle_perigee(t)::numeric, 1) AS perigee_km, + round(tle_apogee(t)::numeric, 1) AS apogee_km, + tle_intl_desig(t) AS cospar +FROM iss; + +-- Observer type parsing +SELECT '40.0N 105.3W 1655m'::observer IS NOT NULL AS observer_ok; +SELECT observer_lat('40.0N 105.3W 1655m'::observer) AS lat, + observer_lon('40.0N 105.3W 1655m'::observer) AS lon; + +-- ECI type parsing +SELECT '(1000.0,2000.0,3000.0,1.0,2.0,3.0)'::eci_position IS NOT NULL AS eci_ok; +SELECT eci_x('(1000.0,2000.0,3000.0,1.0,2.0,3.0)'::eci_position) AS x; + +-- SGP4 propagation at epoch +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + round(eci_speed(sgp4_propagate(t, '2024-01-01 12:00:00+00'))::numeric, 2) AS speed_kms, + round(eci_altitude(sgp4_propagate(t, '2024-01-01 12:00:00+00'))::numeric, 0) AS alt_km +FROM iss; + +-- Altitude overlap operator +WITH sats AS ( + SELECT '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'::tle AS iss, + '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'::tle AS gps +) +SELECT + iss && iss AS same_overlap, + iss && gps AS different_no_overlap +FROM sats; + +-- Subsatellite point +WITH iss AS ( + SELECT '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'::tle AS t +) +SELECT + round(geodetic_lat(subsatellite_point(t, '2024-01-01 12:00:00+00'))::numeric, 2) AS lat, + round(geodetic_alt(subsatellite_point(t, '2024-01-01 12:00:00+00'))::numeric, 0) AS alt_km +FROM iss; + +-- GiST index creation +CREATE TABLE test_sats (id serial, tle tle); +INSERT INTO test_sats (tle) VALUES ( + '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'); +CREATE INDEX test_sats_gist ON test_sats USING gist (tle); +SELECT count(*) FROM test_sats WHERE tle && ( + '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'::tle); +DROP TABLE test_sats;