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

Initial implementation of pg_orbit PostgreSQL extension

Browse Source 
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.
This commit is contained in:
Ryan Malloy 2026-02-15 17:07:07 -07:00
commit 15a830dc40
29 changed files with 6089 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
.gitignore vendored Normal file
View File

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

3
.gitmodules vendored Normal file
View File

@ -0,0 +1,3 @@
[submodule "lib/sat_code"]
path = lib/sat_code
url = https://github.com/Bill-Gray/sat_code.git

147
CLAUDE.md Normal file
View File

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

19
LICENSE Normal file
View File

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

36
Makefile Normal file
View File

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

305
README.md Normal file
View File

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

590
docs/DESIGN.md Normal file
View File

@ -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()` |

1
lib/sat_code Submodule

@ -0,0 +1 @@
Subproject commit ff7b98957dfa2979700a482bde9de9542807293e

4
pg_orbit.control Normal file
View File

@ -0,0 +1,4 @@
comment = 'Orbital mechanics types and functions for PostgreSQL'
default_version = '0.1.0'
module_pathname = '$libdir/pg_orbit'
relocatable = true

491
sql/pg_orbit--0.1.0.sql Normal file
View File

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

820
src/coord_funcs.c Normal file
View File

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

219
src/eci_type.c Normal file
View File

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

505
src/gist_tle.c Normal file
View File

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

284
src/observer_type.c Normal file
View File

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

773
src/pass_funcs.c Normal file
View File

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

12
src/pg_orbit.c Normal file
View File

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

381
src/sgp4_funcs.c Normal file
View File

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

451
src/tle_type.c Normal file
View File

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

190
src/types.h Normal file
View File

@ -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 */

110
test/expected/coord_transforms.out Normal file
View File

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

83
test/expected/gist_index.out Normal file
View File

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

91
test/expected/pass_prediction.out Normal file
View File

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

92
test/expected/sgp4_propagate.out Normal file
View File

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

116
test/expected/tle_parse.out Normal file
View File

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

77
test/sql/coord_transforms.sql Normal file
View File

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

61
test/sql/gist_index.sql Normal file
View File

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

68
test/sql/pass_prediction.sql Normal file
View File

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

67
test/sql/sgp4_propagate.sql Normal file
View File

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

74
test/sql/tle_parse.sql Normal file
View File

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