Compare commits

..

5 Commits

Author SHA1 Message Date
b6c5149cd7 Add Cosmic Queries Cookbook guide — 9 cross-domain SQL recipes
New guide combining multiple pg_orrery function families with
PostgreSQL analytical functions: Kirkwood gap histograms, Kepler
regression, asteroid family taxonomy, universal sky report, planetary
alignment detection, ISS eclipse timing, PostGIS ground tracks,
solar system dashboard view, and satellite shell census.
2026-02-18 16:45:25 -07:00
53733daeba Update docs for v0.8.0 orbital_elements type and MPC parser
Add orbital_elements (72-byte Keplerian element type) to types reference,
three new function sections (oe_from_mpc, small_body_heliocentric,
small_body_observe) to the functions reference, MPC catalog loading
workflow to the comets/asteroids guide, and update CLAUDE.md with
v0.8.0 version numbers, 82 functions, 8 types, 16 test suites.
2026-02-18 14:22:39 -07:00
1adab6e136 Update docs with 66k benchmark results, honest SP-GiST framing
- benchmarks.mdx: Add GiST conjunction screening and KNN sections,
  update all numbers to 66,440-object catalog, PG 17→18, show SP-GiST
  slower than seqscan at this scale with explanation of why
- operators-gist.mdx: Real 66k performance tables for GiST and SP-GiST,
  rewrite KNN example with scalar subquery pattern, add CTE warning
- conjunction-screening.mdx: Update catalog size, candidate counts,
  add KNN scalar subquery note, verified performance numbers
2026-02-18 12:23:47 -07:00
de742fc3aa Fix pg_tle sizeof/INTERNALLENGTH mismatch, exact leaf recheck
The pg_tle struct has been 104 bytes since v0.1.0, but INTERNALLENGTH
is 112.  The size comment claimed "11 doubles (88 bytes)" — there are
10 (80 bytes).  Every palloc(sizeof(pg_tle)) across the codebase
allocated 104 bytes while PostgreSQL's datumCopy/heap_form_tuple
copied 112, causing an 8-byte overread.

Fix: add _reserved[8] to pg_tle, making sizeof(pg_tle) == 112.
This is backward compatible — existing on-disk tuples already have
112 bytes allocated (from typlen), with zeros in the trailing 8.

Also in gist_tle.c:
- Remove TLE_TYPLEN band-aid, use sizeof(pg_tle) everywhere
- Set recheck = false for leaf entries in consistent: the orbital
  key is computed identically to the SQL operator, so the GiST
  leaf check is exact (eliminates unnecessary heap fetches)
2026-02-18 11:48:49 -07:00
347acf0906 Fix GiST union 0-based indexing and palloc size, add 66k benchmark
Two bugs in gist_tle.c caused the && (overlap) operator to return
zero results through the GiST index while sequential scan worked:

1. gist_tle_union read from vector[FirstOffsetNumber] (index 1),
   skipping vector[0] which holds the accumulated union key.
   Every internal node collapsed to a single-entry bounding box.
   Fixed: seed from vector[0], loop from 1.

2. All GiST key allocations used sizeof(tle_orbital_key) (32 bytes)
   or sizeof(pg_tle) (104 bytes), but INTERNALLENGTH is 112.
   index_form_tuple() copies typlen bytes, causing buffer overread.
   Fixed: TLE_TYPLEN constant (112) for all index datum allocations.

The <-> (KNN distance) operator was unaffected because it uses
gist_tle_distance, not gist_tle_consistent.

Verified against 66,440-object catalog:
- && consistency: 9 seqscan == 9 GiST (ISS conjunction)
- <-> KNN: 10 nearest in 2.1ms via index-ordered scan
- All 15 regression tests pass
2026-02-18 11:22:07 -07:00
13 changed files with 1163 additions and 99 deletions

View File

@ -1,9 +1,9 @@
# pg_orrery — A Database Orrery for PostgreSQL
## What This Is
A database orrery — celestial mechanics types and functions for PostgreSQL. Native C extension using PGXS, 68 SQL functions, 7 custom types, covering satellites (SGP4/SDP4), planets (VSOP87 + optional JPL DE441), Moon (ELP2000-82B), 19 planetary moons (L1.2/TASS17/GUST86/MarsSat), stars, comets, Jupiter radio bursts, and interplanetary Lambert transfers.
A database orrery — celestial mechanics types and functions for PostgreSQL. Native C extension using PGXS, 82 SQL functions, 8 custom types, covering satellites (SGP4/SDP4), planets (VSOP87 + optional JPL DE441), Moon (ELP2000-82B), 19 planetary moons (L1.2/TASS17/GUST86/MarsSat), stars, comets, asteroids (MPC catalog), Jupiter radio bursts, and interplanetary Lambert transfers.
**Current version:** 0.3.0 on branch `phase/solar-system-expansion`
**Current version:** 0.8.0 on branch `phase/spgist-orbital-trie`
**Repository:** https://git.supported.systems/warehack.ing/pg_orrery
**Documentation:** https://pg-orrery.warehack.ing
@ -11,7 +11,7 @@ A database orrery — celestial mechanics types and functions for PostgreSQL. Na
```bash
make PG_CONFIG=/usr/bin/pg_config # Compile with PGXS
sudo make install PG_CONFIG=/usr/bin/pg_config # Install extension
make installcheck PG_CONFIG=/usr/bin/pg_config # Run 13 regression test suites
make installcheck PG_CONFIG=/usr/bin/pg_config # Run 16 regression test suites
```
Requires: PostgreSQL 17 development headers, GCC, Make.
@ -27,14 +27,24 @@ Image: `git.supported.systems/warehack.ing/pg_orrery:pg17`
## Project Layout
```
pg_orrery.control # Extension metadata (version 0.3.0)
pg_orrery.control # Extension metadata (version 0.8.0)
Makefile # PGXS build + Docker targets
sql/
pg_orrery--0.1.0.sql # v0.1.0: satellite types/functions/operators
pg_orrery--0.2.0.sql # v0.2.0: solar system (57 functions)
pg_orrery--0.3.0.sql # v0.3.0: complete extension (68 functions)
pg_orrery--0.3.0.sql # v0.3.0: DE ephemeris (68 functions)
pg_orrery--0.4.0.sql # v0.4.0: orbit determination
pg_orrery--0.5.0.sql # v0.5.0: SP-GiST orbital trie
pg_orrery--0.6.0.sql # v0.6.0: conjunction screening
pg_orrery--0.7.0.sql # v0.7.0: GiST improvements
pg_orrery--0.8.0.sql # v0.8.0: orbital_elements type + MPC parser (82 functions)
pg_orrery--0.1.0--0.2.0.sql # Migration: v0.1.0 → v0.2.0 (adds solar system)
pg_orrery--0.2.0--0.3.0.sql # Migration: v0.2.0 → v0.3.0 (adds DE ephemeris)
pg_orrery--0.3.0--0.4.0.sql # Migration: v0.3.0 → v0.4.0
pg_orrery--0.4.0--0.5.0.sql # Migration: v0.4.0 → v0.5.0
pg_orrery--0.5.0--0.6.0.sql # Migration: v0.5.0 → v0.6.0
pg_orrery--0.6.0--0.7.0.sql # Migration: v0.6.0 → v0.7.0
pg_orrery--0.7.0--0.8.0.sql # Migration: v0.7.0 → v0.8.0 (orbital_elements type)
src/
pg_orrery.c # PG_MODULE_MAGIC + _PG_init() (GUC registration)
types.h # All struct definitions + constants + DE body ID mapping
@ -56,6 +66,8 @@ src/
planet_funcs.c # planet_observe(), planet_heliocentric(), sun/moon_observe()
star_funcs.c # star_observe(), star_observe_safe()
kepler_funcs.c # kepler_propagate(), comet_observe()
kepler.h # Shared Kepler solver interface (kepler_position())
orbital_elements_type.c # orbital_elements type, MPC parser, small_body_observe()
l12.c / l12.h # L1.2 Galilean moon theory (Lieske 1998)
tass17.c / tass17.h # TASS 1.7 Saturn moon theory (Vienne & Duriez 1995)
gust86.c / gust86.h # GUST86 Uranus moon theory (Laskar & Jacobson 1987)
@ -80,7 +92,7 @@ src/
PROVENANCE.md # Vendoring decision, modifications, verification
LICENSE # MIT license (Bill Gray / Project Pluto)
test/
sql/ # 13 regression test suites
sql/ # 16 regression test suites
expected/ # Expected output
data/vallado_518.json # 518 Vallado test vectors (AIAA 2006-6753-Rev1)
docs/
@ -104,8 +116,9 @@ All types are fixed-size, `STORAGE = plain`, `ALIGNMENT = double`. No TOAST over
| `observer` | 24 | lat, lon (radians), alt_m (meters) |
| `pass_event` | 48 | AOS/MAX/LOS times + max_el + AOS/LOS azimuth |
| `heliocentric` | 24 | x, y, z in AU (ecliptic J2000 frame) |
| `orbital_elements` | 72 | Classical Keplerian elements for comets/asteroids (epoch, q, e, inc, omega, Omega, tp, H, G) |
## Function Domains (68 total)
## Function Domains (82 total)
| Domain | Theory | Key Functions | Count |
|--------|--------|---------------|-------|
@ -114,7 +127,7 @@ All types are fixed-size, `STORAGE = plain`, `ALIGNMENT = double`. No TOAST over
| Sun/Moon | VSOP87 + ELP2000-82B | `sun_observe()`, `moon_observe()` | 2 |
| Planetary moons | L1.2, TASS17, GUST86, MarsSat | `galilean_observe()`, `saturn_moon_observe()` | 4 |
| Stars | J2000 + IAU 1976 precession | `star_observe()`, `star_observe_safe()` | 2 |
| Comets/asteroids | Two-body Keplerian | `kepler_propagate()`, `comet_observe()` | 2 |
| Comets/asteroids | Two-body Keplerian + MPC | `small_body_observe()`, `oe_from_mpc()`, `kepler_propagate()` | 16 |
| Jupiter radio | Carr et al. (1983) | `jupiter_burst_probability()` | 3 |
| Transfers | Lambert (Izzo 2015) | `lambert_transfer()`, `lambert_c3()` | 2 |
| DE ephemeris | JPL DE440/441 (optional) | `planet_observe_de()`, `moon_observe_de()` | 11 |
@ -239,7 +252,7 @@ All numerical logic is byte-identical to upstream. Verified against 518 Vallado
## Testing
13 regression test suites via `make installcheck`:
16 regression test suites via `make installcheck`:
| Suite | What it tests |
|-------|--------------|
@ -256,10 +269,13 @@ All numerical logic is byte-identical to upstream. Verified against 518 Vallado
| lambert_transfer | Lambert solver, lambert_c3, pork chop grid, error handling |
| de_ephemeris | DE function fallback to VSOP87, cross-provider consistency, error handling |
| vallado_518 | 518 Vallado test vectors (AIAA 2006-6753-Rev1), per-satellite breakdown |
| od_fit | Orbit determination from ECI/topocentric/angles-only observations |
| orbital_elements | orbital_elements type I/O, MPC parser, small_body_observe/heliocentric |
| spgist_tle | SP-GiST orbital trie index operations |
### PG Version Matrix
Test all 13 regression suites + DE reader unit test across PostgreSQL 14-18 using Docker:
Test all 16 regression suites + DE reader unit test across PostgreSQL 14-18 using Docker:
```bash
make test-matrix # Full matrix (PG 14-18)
@ -283,9 +299,9 @@ Logs saved to `test/matrix-logs/pg${ver}.log`. The script reuses the Dockerfile
**Live:** https://pg-orrery.warehack.ing
Starlight docs at `docs/`36 MDX pages covering all domains.
Starlight docs at `docs/`42 MDX pages covering all domains.
Sections: Getting Started, Guides (9 domain walkthroughs incl. DE ephemeris), Workflow Translation (Skyfield/Horizons/GMAT/Radio Jupiter Pro comparisons), Reference (all 68 functions incl. DE variants), Architecture (Hamilton's principles, constant custody, observation pipeline), Performance (benchmarks).
Sections: Getting Started, Guides (9 domain walkthroughs incl. DE ephemeris), Workflow Translation (Skyfield/Horizons/GMAT/Radio Jupiter Pro comparisons), Reference (all 82 functions incl. DE variants + orbital_elements), Architecture (Hamilton's principles, constant custody, observation pipeline), Performance (benchmarks).
### Local Development
```bash
@ -301,7 +317,7 @@ The docs site deploys to the `warehack.ing` VPS (`149.28.126.25`) which runs cad
```bash
ssh -A warehack-ing@pg-orrery.warehack.ing
cd ~/pg_orrery
git pull origin phase/solar-system-expansion # or the current branch
git pull origin phase/spgist-orbital-trie # or the current branch
cd docs
make prod # builds image + starts container
```
@ -310,7 +326,7 @@ make prod # builds image + starts containe
```bash
ssh -A warehack-ing@pg-orrery.warehack.ing
git clone git@git.supported.systems:warehack.ing/pg_orrery.git
cd pg_orrery && git checkout phase/solar-system-expansion
cd pg_orrery && git checkout phase/spgist-orbital-trie
cat > docs/.env << 'EOF'
COMPOSE_PROJECT_NAME=pg-orrery-docs
NODE_ENV=production
@ -343,6 +359,6 @@ cd docs && make prod
## Git Conventions
- One commit per logical change
- Branch per phase: `phase/solar-system-expansion`
- Tag releases: `v0.1.0`, `v0.2.0`
- Branch per phase: `phase/spgist-orbital-trie`
- Tag releases: `v0.1.0`, `v0.2.0`, `v0.3.0`
- Commit messages: imperative mood, no AI attribution

View File

@ -0,0 +1,168 @@
pg_orrery Full Index Benchmark — 66k Catalog
===========================================================
Date: 2026-02-18
PostgreSQL: 18.1
Catalog: 66,440 objects (merged from 4 sources)
Sources: Space-Track (66,248), CelesTrak active (5 unique),
SatNOGS (110 unique), CelesTrak SupGP (77 unique + 8,167 epoch updates)
Includes: 362 Alpha-5 objects (NORAD > 99,999)
Orbital regime breakdown:
LEO (<2000km): 63,097 (95.0%)
GEO/HEO (>34000km): 1,760 ( 2.6%)
MEO (2000-20000km): 1,277 ( 1.9%)
GEO-transfer: 306 ( 0.5%)
Index sizes:
SP-GiST (tle_spgist_ops): 67 ms build, 11 MB
GiST (tle_ops): 93 ms build, 15 MB
═══════════════════════════════════════════════════════════
SP-GiST: Visibility Cone (&?) — "Can this satellite pass over me?"
═══════════════════════════════════════════════════════════
SP-GiST prunes by altitude band, inclination, and RAAN window.
The &? operator answers: "Could this satellite be visible from this
observer during this time window above this minimum elevation?"
Query │ SP-GiST │ Seqscan │ Candidates │ Pruned%
───────────────────────┼──────────┼──────────┼────────────┼────────
Eagle 2h/10deg │ 16.1 ms │ 12.1 ms │ 10,763 │ 83.8%
Eagle 24h/10deg │ 23.3 ms │ 12.5 ms │ 61,426 │ 7.5%
Equator 2h/10deg │ 16.8 ms │ 12.1 ms │ 10,174 │ 84.7%
Eagle 2h/45deg │ 16.9 ms │ 11.9 ms │ 6,796 │ 89.8%
Consistency: PASS (all 4 scenarios: 0 false neg, 0 false pos)
═══════════════════════════════════════════════════════════
GiST: Overlap (&&) — "Does this satellite share my orbit band?"
═══════════════════════════════════════════════════════════
GiST groups satellites by [altitude_low, altitude_high] × [inclination].
The && operator answers: "Do these two TLEs occupy overlapping orbit bands?"
Used for conjunction screening — finding potential collision partners.
Critical bugfix in this session:
Bug 1: palloc size mismatch (sizeof(pg_tle)=104 vs INTERNALLENGTH=112)
Bug 2: gist_tle_union used 1-based indexing (picksplit convention)
instead of 0-based (union convention), skipping vector[0]
Query │ GiST │ Seqscan │ Matches
───────────────────────┼──────────┼──────────┼────────
ISS conjunction │ 10.9 ms │ 63.3 ms │ 9
Starlink-230369 │ 9.5 ms │ 14.9 ms │ 0
SYNCOM 2 (GEO) │ 4.0 ms │ 7.2 ms │ 0
Consistency: PASS (ISS: 9 seqscan == 9 GiST, 0 mismatch)
ISS conjunction candidates (altitude + inclination overlap):
PROGRESS MS-31, PROGRESS MS-32, SOYUZ MS-28,
DRAGON FREEDOM 3, DRAGON CRS-33, CYGNUS NG-23,
HTV-X1, ISS (NAUKA), OBJECT E
— All ISS-visiting vehicles or co-orbital modules. ✓
═══════════════════════════════════════════════════════════
GiST: KNN (<->) — "What's nearest to this orbit?"
═══════════════════════════════════════════════════════════
GiST KNN uses altitude-band distance for index-ordered scans.
The <-> operator returns orbital altitude separation in km.
Probe must be a scalar subquery for index ordering to activate.
Query │ GiST KNN │ Buffers │ Notes
───────────────────────┼──────────┼─────────┼──────────────
10 nearest to ISS │ 2.1 ms │ 982 │ Index-ordered
10 nearest to SYNCOM 2 │ 0.2 ms │ 40 │ Index-ordered
100 nearest to ISS │ 1.4 ms │ 1,062 │ Index-ordered
Within 50km of ISS │ 16.0 ms │ 4,014 │ 12,496 matches
Pattern for KNN queries (probe as scalar subquery):
ORDER BY b.tle <-> (SELECT tle FROM catalog WHERE norad_id = 25544 LIMIT 1)
LIMIT 10;
→ Index Scan using bench_gist_idx, Order By: tle <-> InitPlan
═══════════════════════════════════════════════════════════
EXPLAIN ANALYZE Details
═══════════════════════════════════════════════════════════
SP-GiST 2h/Eagle/10deg:
Index Only Scan using bench_spgist_idx
Heap Fetches: 0 (pure index scan)
Buffers: shared hit=4964
17.5 ms execution
SeqScan 2h/Eagle/10deg:
Seq Scan, Filter rows removed: 55,677
Buffers: shared hit=1338
12.5 ms execution
GiST && ISS conjunction:
Nested Loop → Index Scan using bench_gist_idx
Index Cond: (tle && a.tle)
Index Searches: 1, Buffers: shared hit=287
4.1 ms execution
GiST KNN 10 nearest ISS:
Index Scan using bench_gist_idx
Order By: (tle <-> InitPlan)
Index Searches: 1
2.1 ms execution
═══════════════════════════════════════════════════════════
Pruning Summary
═══════════════════════════════════════════════════════════
Scenario │ Catalog │ Candidates │ Candidate% │ Pruned%
─────────────────┼─────────┼────────────┼────────────┼────────
2h/Eagle/10deg │ 66,440 │ 10,763 │ 16.2% │ 83.8%
2h/Equator/10deg │ 66,440 │ 10,174 │ 15.3% │ 84.7%
2h/Eagle/45deg │ 66,440 │ 6,796 │ 10.2% │ 89.8%
24h/Eagle/10deg │ 66,440 │ 61,426 │ 92.5% │ 7.5%
═══════════════════════════════════════════════════════════
Application Queries
═══════════════════════════════════════════════════════════
"What's overhead right now?" (SP-GiST filter + SGP4 propagation):
15 satellites above horizon, top: NAVSTAR 57 at 81.7° el
107 ms (includes SGP4 propagation for each candidate)
ISS pass prediction (next 24h from 66k catalog):
6 passes found, max 87.6° elevation
3.8 ms
ISS conjunction screening (GiST && on 66k catalog):
9 co-orbital objects found
4.6 ms via GiST (vs 63.3 ms seqscan — 5.8x speedup)
═══════════════════════════════════════════════════════════
Key Observations
═══════════════════════════════════════════════════════════
1. GiST && is the clear winner for conjunction screening:
- ISS: 10.9ms GiST vs 63.3ms seqscan (5.8x speedup)
- Only 287 buffer hits vs 1,338 for seqscan
- Returns exactly the right 9 co-orbital objects
2. GiST KNN is extremely fast for "nearest orbit" queries:
- 10 nearest: 2.1ms with index ordering
- GEO satellite: 0.15ms (sparse regime, fewer nodes to traverse)
- Requires scalar subquery probe pattern for index ordering
3. SP-GiST visibility cone handles 2h windows well:
- 83.8% pruning at 10° min_el (Eagle, 2h)
- 89.8% pruning at 45° min_el
- Falls behind seqscan at 24h windows (7.5% pruning not worth index overhead)
4. Both indexes are compact:
- SP-GiST: 11 MB for 66k objects (170 bytes/object)
- GiST: 15 MB for 66k objects (237 bytes/object)
- Build times: 67ms and 93ms respectively
5. Zero false positives/negatives across all consistency checks.
Alpha-5 support:
- Bill Gray's get_el.c parser handles Alpha-5 natively
- T0002 → 270002, A0001 → 100001, Z9999 → 339999 ✓
- Round-trip (parse → output) preserves Alpha-5 encoding ✓
- 362 Alpha-5 objects loaded and indexed without issues ✓

View File

@ -61,6 +61,7 @@ export default defineConfig({
items: [
{ label: "Tracking Satellites", slug: "guides/tracking-satellites" },
{ label: "Observing the Solar System", slug: "guides/observing-solar-system" },
{ label: "Cosmic Queries Cookbook", slug: "guides/cosmic-queries" },
{ label: "Planetary Moon Tracking", slug: "guides/planetary-moons" },
{ label: "Star Catalogs in SQL", slug: "guides/star-catalogs" },
{ label: "Comet & Asteroid Tracking", slug: "guides/comets-asteroids" },

View File

@ -22,7 +22,7 @@ PostGIS added spatial awareness to PostgreSQL — suddenly your database underst
| Moon | ELP2000-82B (Chapront, 1988) | `moon_observe()` | ~10 arcseconds |
| Planetary moons | L1.2, TASS17, GUST86, MarsSat | `galilean_observe()`, etc. | ~1-10 arcseconds |
| Stars | J2000 catalog + precession | `star_observe()` | Limited by catalog |
| Comets/asteroids | Two-body Keplerian | `kepler_propagate()`, `comet_observe()` | Varies with eccentricity |
| Comets/asteroids | Two-body Keplerian | `small_body_observe()`, `oe_from_mpc()`, `kepler_propagate()` | Varies with eccentricity |
| Jupiter radio | Carr et al. (1983) sources | `jupiter_burst_probability()` | Empirical probability |
| Transfers | Lambert (Izzo, 2015) | `lambert_transfer()`, `lambert_c3()` | Ballistic two-body |
| DE ephemeris (optional) | JPL DE440/441 | `planet_observe_de()`, `moon_observe_de()` | ~0.1 milliarcsecond |

View File

@ -21,17 +21,22 @@ The pattern is familiar: download elements, propagate in Python or C, transform
## What changes with pg_orrery
Two functions handle comet/asteroid computation:
Five functions handle comet/asteroid computation:
| Function | What it does |
|---|---|
| `kepler_propagate(q, e, i, omega, Omega, T, time)` | Propagates orbital elements to a heliocentric position (AU) |
| `comet_observe(q, e, i, omega, Omega, T, ex, ey, ez, observer, time)` | Full observation pipeline: propagate + geocentric transform + topocentric |
| `oe_from_mpc(line)` | Parses one MPCORB.DAT line into an `orbital_elements` type |
| `small_body_heliocentric(oe, time)` | Heliocentric position from bundled elements |
| `small_body_observe(oe, observer, time)` | Topocentric observation — auto-fetches Earth via VSOP87 |
`kepler_propagate()` solves Kepler's equation for elliptic (e < 1), parabolic (e = 1), and hyperbolic (e > 1) orbits. The solver handles all three cases with appropriate numerical methods.
`comet_observe()` wraps the full chain: propagate the comet's position, subtract Earth's heliocentric position, and transform to horizon coordinates. You supply Earth's position as three floats (ecliptic J2000, AU) because you might want to compute it once and reuse it across many comets.
`small_body_observe()` (v0.8.0) does the same thing but fetches Earth's position automatically — you just pass the `orbital_elements` type and an observer. See the [orbital_elements type section](#the-orbital_elements-type) below.
The parameters map directly to MPC orbital element format:
| Parameter | MPC field | Units |
@ -56,6 +61,107 @@ Keplerian propagation assumes the body is influenced only by the Sun. Real small
For MPC elements less than a few months old, two-body propagation is typically accurate to a few arcminutes for asteroids and tens of arcminutes for comets. Fresh elements give better results.
## The `orbital_elements` type
The raw-parameter functions (`kepler_propagate`, `comet_observe`) work well when you have elements in variables or a table with individual columns. But they require passing 611 float8 arguments per call, and `comet_observe` requires you to manually fetch Earth's position.
The `orbital_elements` type (v0.8.0) bundles all nine classical elements into a single 72-byte PostgreSQL datum:
```sql
-- Construct from a literal
SELECT '(2460605.5,2.5478,0.0789126,10.58664,73.42937,80.2686,2460319.0,3.33,0.12)'::orbital_elements;
-- Or parse directly from the MPC catalog
SELECT oe_from_mpc(
'00001 3.33 0.12 K24AM 60.07966 73.42937 80.26860 10.58664 0.0789126 0.21406048 2.7660961 0 MPO838504 8738 115 1801-2024 0.65 M-v 30k MPCLINUX 0000 (1) Ceres 20240825'
);
```
With bundled elements, observation becomes a single function call:
```sql
-- Before (comet_observe): 11 arguments, manual Earth fetch
WITH earth AS (SELECT planet_heliocentric(3, now()) AS h)
SELECT topo_elevation(comet_observe(
2.5478, 0.0789, 10.59, 73.43, 80.27, 2460319.0,
helio_x(h), helio_y(h), helio_z(h),
'40.0N 105.3W 1655m'::observer, now()))
FROM earth;
-- After (small_body_observe): 3 arguments, Earth auto-fetched
SELECT topo_elevation(small_body_observe(
oe, '40.0N 105.3W 1655m'::observer, now()))
FROM asteroid_catalog;
```
### Load and query the MPC catalog
The MPC publishes MPCORB.DAT — orbital elements for every numbered asteroid. Here's how to load it into PostgreSQL:
<Steps>
1. **Create a table with an `orbital_elements` column:**
```sql
CREATE TABLE asteroids (
number int PRIMARY KEY,
name text,
oe orbital_elements NOT NULL
);
```
2. **Load via a staging table:**
```sql
-- Stage the raw text lines
CREATE TEMP TABLE mpc_raw (line text);
\copy mpc_raw FROM 'MPCORB.DAT'
-- Parse into orbital_elements, extract number and name
INSERT INTO asteroids (number, name, oe)
SELECT substring(line FROM 1 FOR 7)::int,
trim(substring(line FROM 167 FOR 30)),
oe_from_mpc(line)
FROM mpc_raw
WHERE length(line) >= 103
AND substring(line FROM 1 FOR 7) ~ '^\s*\d+$';
DROP TABLE mpc_raw;
```
3. **Query: what asteroids are above 20 degrees tonight?**
```sql
SELECT name, number,
round(topo_elevation(t)::numeric, 1) AS el,
round(topo_azimuth(t)::numeric, 1) AS az,
round((topo_range(t) / 149597870.7)::numeric, 2) AS dist_au
FROM asteroids,
small_body_observe(oe, '40.0N 105.3W 1655m'::observer, now()) AS t
WHERE topo_elevation(t) > 20
ORDER BY topo_elevation(t) DESC
LIMIT 20;
```
4. **Query: heliocentric distance of Ceres over 6 months:**
```sql
SELECT t::date AS date,
round(helio_distance(
small_body_heliocentric(oe, t))::numeric, 4) AS dist_au
FROM asteroids,
generate_series(
'2025-01-01'::timestamptz,
'2025-07-01'::timestamptz,
interval '15 days'
) AS t
WHERE number = 1;
```
</Steps>
<Aside type="tip">
For batch observation at a single time, `comet_observe()` is still more efficient — it lets you compute Earth's VSOP87 position once with `planet_heliocentric(3, t)` and reuse it across all objects. `small_body_observe()` re-fetches Earth on every call. For interactive single-object queries, `small_body_observe()` is simpler.
</Aside>
## Try it
### Circular orbit sanity check

View File

@ -17,7 +17,7 @@ Operational conjunction screening uses several established tools and data source
- **CelesTrak SOCRATES**: Dr. Kelso's web-based close-approach listing. Updated regularly, covers the full public catalog. Not queryable; you read reports.
- **Python scripts**: Propagate the catalog in a loop, compute pairwise distances, filter by threshold. Works for small catalogs. Does not scale.
The fundamental challenge: a catalog of 25,000+ tracked objects produces over 300 million unique pairs. Even checking each pair at a single epoch takes significant time. Checking over a 7-day window at 1-minute resolution is computationally prohibitive without pre-filtering.
The fundamental challenge: a catalog of 66,000+ tracked objects produces over 2 billion unique pairs. Even checking each pair at a single epoch takes significant time. Checking over a 7-day window at 1-minute resolution is computationally prohibitive without pre-filtering.
## What changes with pg_orrery
@ -89,7 +89,7 @@ INSERT INTO catalog VALUES (99901, 'Equatorial-LEO',
CREATE INDEX catalog_orbit_gist ON catalog USING gist (tle);
```
The index builds in milliseconds for a small table. For a full 25,000-object catalog, expect about 200ms.
The index builds in milliseconds for a small table. For a full 66,440-object catalog, build time is 93 ms (15 MB index).
### Check orbital parameters
@ -157,19 +157,21 @@ This should return only ISS itself (and not Equatorial-LEO, which has a differen
Find the 3 closest objects to the ISS by altitude band separation, ordered by distance:
```sql
SET enable_seqscan = off;
-- Scalar subquery probe enables GiST index-ordered scan
SELECT name,
round((tle <-> (SELECT tle FROM catalog WHERE norad_id = 25544))::numeric, 0) AS alt_dist_km
round((tle <-> (SELECT tle FROM catalog WHERE norad_id = 25544 LIMIT 1))::numeric, 0)
AS alt_dist_km
FROM catalog
WHERE norad_id != 25544
ORDER BY tle <-> (SELECT tle FROM catalog WHERE norad_id = 25544)
ORDER BY tle <-> (SELECT tle FROM catalog WHERE norad_id = 25544 LIMIT 1)
LIMIT 3;
RESET enable_seqscan;
```
This uses the GiST distance operator for efficient ordering. PostgreSQL's KNN-GiST infrastructure handles this without computing all distances upfront.
This uses the GiST distance operator for efficient ordering. PostgreSQL's KNN-GiST infrastructure traverses the tree by increasing distance without computing all distances upfront. On a 66,440-object catalog, this completes in 2.1 ms for 10 neighbors.
<Aside type="caution" title="Use scalar subqueries, not CTEs">
GiST index-ordered scans require the probe value to be visible to the planner as a constant. A `WITH iss AS (...)` CTE makes the probe opaque, forcing a full sequential scan and sort. Always use `(SELECT tle FROM ... LIMIT 1)` as the probe argument for KNN queries on large catalogs.
</Aside>
### Self-overlap is always true
@ -206,7 +208,7 @@ The complete two-stage workflow for a larger catalog:
AND c.norad_id != 25544;
```
For the ISS in a 25,000-object catalog, this typically returns a few hundred candidates.
For the ISS in a 66,440-object catalog, this returns 9 candidates (all co-orbital vehicles: visiting spacecraft, modules, and debris). The GiST index scan completes in 4.6 ms vs. 63.3 ms for a sequential scan.
3. **Stage 2: Time-resolved distance computation:**
@ -253,7 +255,7 @@ ORDER BY actual_dist_km;
```
<Aside type="tip" title="Performance scaling">
The GiST index is the key to scaling. Without it, screening a 25,000-object catalog for all-vs-all conjunctions means 300 million pair evaluations. With GiST, the `&&` operator reduces this to tens of thousands of candidate pairs. The `tle_distance()` computation on candidates is then feasible even at fine time resolution.
The GiST index is the key to scaling. Without it, screening a 66,440-object catalog for all-vs-all conjunctions means over 2 billion pair evaluations. With GiST, the `&&` operator reduces single-probe screening from 63 ms (sequential) to 4.6 ms (indexed), a 5.8x speedup. For the ISS, only 9 candidates survive from 66,440 objects. The `tle_distance()` computation on these survivors is then feasible even at 1-minute time resolution over multi-day windows.
</Aside>
### Monitoring over time

View File

@ -0,0 +1,435 @@
---
title: Cosmic Queries Cookbook
sidebar:
order: 3
---
import { Steps, Aside, Tabs, TabItem } from "@astrojs/starlight/components";
Each pg_orrery guide covers a single domain — satellites, planets, comets, radio bursts. This page is different. These nine queries combine multiple pg_orrery function families with PostgreSQL's analytical engine to ask questions that span physical theories, orbital regimes, and even external extensions like PostGIS. They range from statistical analysis of 600,000+ asteroids to real-time cross-domain sky surveys.
Every query here is copy-paste ready. Swap the observer coordinates for your location and the timestamps for your session.
## Prerequisites
Not every query needs the same data. Here's what to load before you start:
| Data | Queries that use it | Setup |
|------|---------------------|-------|
| `asteroids` table (MPC catalog) | 1, 2, 3, 4 | See [Comet & Asteroid Tracking](/guides/comets-asteroids) — load the MPC export with `oe_from_mpc()` |
| `satellites` table (TLE catalog) | 4, 6, 7, 9 | See [Building TLE Catalogs](/guides/catalog-management) — any catalog with a `tle` column works |
| `countries` table (Natural Earth) | 7 | PostGIS + Natural Earth boundaries — [setup below](#postgis-setup) |
| PostGIS extension | 7 | `CREATE EXTENSION IF NOT EXISTS postgis;` |
| None — built-in functions only | 5, 8 | Just pg_orrery |
The expected table schemas:
```sql
-- Asteroids: name + orbital elements from MPC
CREATE TABLE asteroids (
name text PRIMARY KEY,
oe orbital_elements
);
-- Satellites: NORAD ID + parsed TLE
CREATE TABLE satellites (
norad_id integer PRIMARY KEY,
name text,
tle tle
);
```
---
## The Asteroid Belt as a Dataset
The MPC catalog isn't just a list of orbits — it's a dataset with 600,000+ rows and rich statistical structure. PostgreSQL's aggregate functions turn it into an orbital mechanics laboratory.
### 1. Kirkwood Gaps — Jupiter's Gravitational Fingerprint
In 1866, Daniel Kirkwood noticed that asteroids avoid certain orbital distances. The gaps correspond to mean-motion resonances with Jupiter: an asteroid at 2.50 AU completes exactly 3 orbits for every 1 of Jupiter's (the 3:1 resonance). Over millions of years, Jupiter's periodic gravitational nudges clear these orbits out.
`width_bucket()` bins the semi-major axes into a 200-bin histogram across the main belt. The depletions at 2.50, 2.82, 2.96, and 3.28 AU are unmistakable:
```sql
WITH belt AS (
SELECT oe_semi_major_axis(oe) AS a_au
FROM asteroids
WHERE oe_eccentricity(oe) < 0.4
AND oe_semi_major_axis(oe) IS NOT NULL
AND oe_semi_major_axis(oe) BETWEEN 1.5 AND 5.5
)
SELECT round((1.5 + (bucket - 1) * 0.02)::numeric, 2) AS a_au,
count AS n_asteroids
FROM (
SELECT width_bucket(a_au, 1.5, 5.5, 200) AS bucket, count(*) AS count
FROM belt GROUP BY bucket
) sub
ORDER BY a_au;
```
The output is a classic Kirkwood gap diagram. Plot `a_au` vs `n_asteroids` and the resonance depletions jump out — the 3:1 gap at 2.50 AU is the deepest, with the 5:2 (2.82 AU), 7:3 (2.96 AU), and 2:1 (3.28 AU) gaps clearly visible.
### 2. Kepler's Third Law as a Regression
Kepler published his third law in 1619: the square of a planet's orbital period is proportional to the cube of its semi-major axis, or equivalently $\log P = 1.5 \cdot \log a$. With `regr_slope()` and `regr_r2()`, you can verify this 400-year-old relationship against every bound asteroid in the MPC catalog:
```sql
WITH bounded AS (
SELECT oe_semi_major_axis(oe) AS a_au, oe_period_years(oe) AS p_yr
FROM asteroids
WHERE oe_semi_major_axis(oe) IS NOT NULL
AND oe_semi_major_axis(oe) BETWEEN 0.5 AND 100.0
)
SELECT round(regr_slope(ln(p_yr), ln(a_au))::numeric, 6) AS slope,
round(exp(regr_intercept(ln(p_yr), ln(a_au)))::numeric, 6) AS intercept_years,
regr_count(ln(p_yr), ln(a_au)) AS n_objects,
round(regr_r2(ln(p_yr), ln(a_au))::numeric, 12) AS r_squared
FROM bounded;
```
The slope will be exactly 1.500000 (Kepler's 3/2 power law). The intercept will be 1.000000 years (because for $a = 1$ AU, $P = 1$ year — Earth). The $R^2$ will be 1.000000000000. Not approximately. Exactly. This isn't a statistical correlation — it's a mathematical identity baked into `oe_period_years()`, which computes $a^{3/2}$. The query is a 600,000-row proof that the accessor functions are self-consistent.
<Aside type="tip" title="Why R² = 1 exactly">
`oe_period_years()` is defined as `a^1.5` where `a = q/(1-e)`. The regression isn't discovering a physical law — it's confirming that the accessor functions implement Kepler's third law without floating-point drift across the entire catalog. If you ever see R² < 1.0, something is wrong with your data (likely a corrupted MPC record).
</Aside>
### 3. Asteroid Family Taxonomy
Collisional families — groups of asteroids created by a single catastrophic impact — cluster tightly in (semi-major axis, eccentricity) space. A 2D `width_bucket()` grid reveals these density peaks as hot spots:
```sql
WITH belt AS (
SELECT oe_semi_major_axis(oe) AS a,
oe_eccentricity(oe) AS e
FROM asteroids
WHERE oe_eccentricity(oe) < 0.4
AND oe_semi_major_axis(oe) IS NOT NULL
AND oe_semi_major_axis(oe) BETWEEN 2.0 AND 3.5
)
SELECT round((2.0 + (a_bin - 1) * 0.03)::numeric, 2) AS a_au,
round((0.0 + (e_bin - 1) * 0.01)::numeric, 2) AS ecc,
count(*) AS n
FROM (
SELECT width_bucket(a, 2.0, 3.5, 50) AS a_bin,
width_bucket(e, 0.0, 0.4, 40) AS e_bin
FROM belt
) sub
GROUP BY a_bin, e_bin
HAVING count(*) >= 10
ORDER BY n DESC;
```
The highest-density cells correspond to known collisional families: Flora (~2.2 AU, e~0.15), Themis (~3.13 AU, e~0.15), Koronis (~2.87 AU, e~0.05), and Eos (~3.01 AU, e~0.07). The `HAVING count(*) >= 10` filter suppresses noise in sparsely populated cells. Increase the threshold to isolate only the major families; decrease it to reveal smaller groupings.
---
## Cross-Domain Observation
These queries combine satellite tracking, planetary ephemerides, and solar observation — functions backed by different physical theories, unified through pg_orrery's common `topocentric` return type.
### 4. Universal Sky Report — Everything at Once
Four gravitational theories in one query. `planet_observe()` uses VSOP87, `moon_observe()` uses ELP2000-82B, `observe_safe()` uses SGP4/SDP4, and `small_body_observe()` uses two-body Keplerian propagation. They all return `topocentric`, so `UNION ALL` works:
```sql
WITH obs AS (SELECT '40.0N 105.3W 1655m'::observer AS o),
sky AS (
-- Planets (VSOP87)
SELECT 'Mercury' AS body, planet_observe(1, o, now()) AS topo FROM obs
UNION ALL SELECT 'Venus', planet_observe(2, o, now()) FROM obs
UNION ALL SELECT 'Mars', planet_observe(4, o, now()) FROM obs
UNION ALL SELECT 'Jupiter', planet_observe(5, o, now()) FROM obs
UNION ALL SELECT 'Saturn', planet_observe(6, o, now()) FROM obs
UNION ALL SELECT 'Uranus', planet_observe(7, o, now()) FROM obs
UNION ALL SELECT 'Neptune', planet_observe(8, o, now()) FROM obs
-- Sun and Moon
UNION ALL SELECT 'Sun', sun_observe(o, now()) FROM obs
UNION ALL SELECT 'Moon', moon_observe(o, now()) FROM obs
-- Satellites (SGP4/SDP4) — observe_safe returns NULL for decayed TLEs
UNION ALL
SELECT s.name, observe_safe(s.tle, obs.o, now())
FROM satellites s, obs
WHERE s.norad_id IN (25544, 20580, 48274) -- ISS, HST, Tiangong
-- Asteroids (two-body Keplerian)
UNION ALL
SELECT a.name, small_body_observe(a.oe, obs.o, now())
FROM asteroids a, obs
WHERE a.name IN ('Ceres', 'Vesta', 'Pallas')
)
SELECT body,
round(topo_azimuth(topo)::numeric, 1) AS az,
round(topo_elevation(topo)::numeric, 1) AS el,
CASE WHEN topo_elevation(topo) > 0 THEN 'visible' ELSE 'below horizon' END AS status
FROM sky
WHERE topo IS NOT NULL
ORDER BY topo_elevation(topo) DESC;
```
Replace the NORAD IDs and asteroid names with whatever interests you. The `observe_safe` call is important for satellites — a decayed TLE will return NULL instead of raising an error, and the `WHERE topo IS NOT NULL` filter drops it cleanly.
### 5. Planetary Alignment Detector
How close are any two planets in the sky right now? The angular separation between two objects at (az₁, el₁) and (az₂, el₂) is the spherical law of cosines. PostgreSQL's built-in `sind()`, `cosd()`, and `acosd()` work in degrees — matching the degree output of `topo_azimuth()` and `topo_elevation()`:
```sql
WITH obs AS (SELECT '40.0N 105.3W 1655m'::observer AS o),
planets AS (
SELECT body_id, name,
planet_observe(body_id, o, now()) AS topo
FROM obs,
(VALUES (1,'Mercury'),(2,'Venus'),(4,'Mars'),
(5,'Jupiter'),(6,'Saturn')) AS p(body_id, name)
)
SELECT a.name AS body_a, b.name AS body_b,
round(acosd(
sind(topo_elevation(a.topo)) * sind(topo_elevation(b.topo)) +
cosd(topo_elevation(a.topo)) * cosd(topo_elevation(b.topo)) *
cosd(topo_azimuth(a.topo) - topo_azimuth(b.topo))
)::numeric, 1) AS separation_deg
FROM planets a
JOIN planets b ON a.body_id < b.body_id
WHERE topo_elevation(a.topo) > 0
AND topo_elevation(b.topo) > 0
ORDER BY separation_deg;
```
The `a.body_id < b.body_id` join condition gives each pair exactly once (VenusJupiter, not also JupiterVenus). Only above-horizon planets are included — no point measuring the angular separation of objects you can't see.
To find the closest approach over a year, sweep with `generate_series` and pick the tightest dates:
```sql
WITH obs AS (SELECT '40.0N 105.3W 1655m'::observer AS o),
sweep AS (
SELECT t,
planet_observe(5, o, t) AS jupiter,
planet_observe(6, o, t) AS saturn
FROM obs,
generate_series(
'2026-01-01'::timestamptz,
'2026-12-31'::timestamptz,
interval '1 day'
) AS t
)
SELECT t::date AS date,
round(acosd(
sind(topo_elevation(jupiter)) * sind(topo_elevation(saturn)) +
cosd(topo_elevation(jupiter)) * cosd(topo_elevation(saturn)) *
cosd(topo_azimuth(jupiter) - topo_azimuth(saturn))
)::numeric, 1) AS separation_deg
FROM sweep
WHERE topo_elevation(jupiter) > 0
AND topo_elevation(saturn) > 0
ORDER BY separation_deg
LIMIT 10;
```
### 6. ISS Eclipse Timing — Shadow Entry and Exit
A satellite enters Earth's shadow when the Sun is below the horizon at the satellite's nadir point. This chains three pg_orrery domains together: `subsatellite_point()` (SGP4 → geodetic), `observer_from_geodetic()` (geodetic → observer), and `sun_observe()` (VSOP87 → topocentric). The `lag()` window function then detects the sunlit/shadow transitions:
```sql
WITH orbit AS (
SELECT t,
subsatellite_point(s.tle, t) AS geo
FROM satellites s,
generate_series(now(), now() + interval '93 minutes', interval '30 seconds') AS t
WHERE s.norad_id = 25544
),
shadow AS (
SELECT t,
geodetic_lat(geo) AS lat,
geodetic_lon(geo) AS lon,
topo_elevation(
sun_observe(
observer_from_geodetic(geodetic_lat(geo), geodetic_lon(geo)),
t
)
) AS sun_el_at_nadir
FROM orbit
)
SELECT t,
round(lat::numeric, 2) AS lat,
round(lon::numeric, 2) AS lon,
round(sun_el_at_nadir::numeric, 1) AS sun_el,
CASE WHEN sun_el_at_nadir < 0 THEN 'SHADOW' ELSE 'SUNLIT' END AS state,
CASE
WHEN sun_el_at_nadir < 0 AND lag(sun_el_at_nadir) OVER (ORDER BY t) >= 0
THEN '>>> ECLIPSE ENTRY'
WHEN sun_el_at_nadir >= 0 AND lag(sun_el_at_nadir) OVER (ORDER BY t) < 0
THEN '<<< ECLIPSE EXIT'
END AS transition
FROM shadow
ORDER BY t;
```
<Aside type="caution" title="Approximation accuracy">
This treats the satellite's nadir point as the shadow boundary, which is geometrically simplified — it ignores the satellite's altitude above the surface and Earth's atmospheric refraction. For the ISS at ~400 km altitude, the shadow entry/exit times are accurate to roughly 1020 seconds. For precise eclipse predictions, you'd need a cylindrical or conical shadow model. But for observation planning — knowing *approximately* when the ISS goes dark — this is very usable.
</Aside>
### 7. Ground Track Geography with PostGIS
Where on Earth is the ISS flying over? Combine `ground_track()` with PostGIS spatial joins against Natural Earth country boundaries.
The simpler approach: a point-in-polygon test at each time step. Each (lat, lon) from the ground track becomes a PostGIS point, joined against country polygons:
```sql
WITH track AS (
SELECT t, lat, lon, alt
FROM satellites s,
ground_track(s.tle, now(), now() + interval '93 minutes', interval '30 seconds')
WHERE s.norad_id = 25544
)
SELECT track.t,
round(track.lat::numeric, 2) AS lat,
round(track.lon::numeric, 2) AS lon,
round(track.alt::numeric, 0) AS alt_km,
c.name AS country
FROM track
LEFT JOIN countries c
ON ST_Contains(c.geom, ST_SetSRID(ST_MakePoint(track.lon, track.lat), 4326));
```
The `LEFT JOIN` keeps rows over oceans (where `country` is NULL). The `ST_MakePoint()` argument order is (longitude, latitude) — x before y, the PostGIS convention.
<Aside type="note" title="PostGIS setup" id="postgis-setup">
Download [Natural Earth 110m countries](https://www.naturalearthdata.com/downloads/110m-cultural-vectors/) and load them:
```bash
wget https://naciscdn.org/naturalearth/110m/cultural/ne_110m_admin_0_countries.zip
unzip ne_110m_admin_0_countries.zip
shp2pgsql -s 4326 ne_110m_admin_0_countries.shp countries | psql -d your_database
```
This creates a `countries` table with `name` (text) and `geom` (geometry) columns. Add a spatial index for faster lookups:
```sql
CREATE EXTENSION IF NOT EXISTS postgis;
CREATE INDEX ON countries USING gist (geom);
```
</Aside>
For a communications footprint, buffer the subsatellite point by the satellite's horizon radius. At ISS altitude (~400 km), the radio horizon is approximately 2,300 km:
```sql
-- Countries within ISS radio line-of-sight right now
WITH nadir AS (
SELECT subsatellite_point(s.tle, now()) AS geo
FROM satellites s
WHERE s.norad_id = 25544
)
SELECT c.name AS country
FROM nadir, countries c
WHERE ST_DWithin(
c.geom::geography,
ST_SetSRID(ST_MakePoint(geodetic_lon(geo), geodetic_lat(geo)), 4326)::geography,
2300000 -- 2,300 km horizon radius in meters
);
```
---
## The Solar System as Data
SQL views and aggregation functions turn pg_orrery's observation pipeline into data products — live dashboards and statistical breakdowns that update every time you query them.
### 8. Celestial Clock — Solar System Dashboard
One row. Every planet's elevation, the Moon, the Sun, Io's phase angle, Jupiter's central meridian longitude, and the decametric burst probability. Wrap it in a `CREATE VIEW` and `SELECT * FROM solar_system_now` becomes a real-time dashboard:
```sql
CREATE VIEW solar_system_now AS
SELECT now() AS computed_at,
round(topo_elevation(sun_observe(o, now()))::numeric, 1) AS sun_el,
round(topo_elevation(moon_observe(o, now()))::numeric, 1) AS moon_el,
round(topo_elevation(planet_observe(1, o, now()))::numeric, 1) AS mercury_el,
round(topo_elevation(planet_observe(2, o, now()))::numeric, 1) AS venus_el,
round(topo_elevation(planet_observe(4, o, now()))::numeric, 1) AS mars_el,
round(topo_elevation(planet_observe(5, o, now()))::numeric, 1) AS jupiter_el,
round(topo_elevation(planet_observe(6, o, now()))::numeric, 1) AS saturn_el,
round(topo_elevation(planet_observe(7, o, now()))::numeric, 1) AS uranus_el,
round(topo_elevation(planet_observe(8, o, now()))::numeric, 1) AS neptune_el,
round(io_phase_angle(now())::numeric, 1) AS io_phase,
round(jupiter_cml(o, now())::numeric, 1) AS jupiter_cml,
round(jupiter_burst_probability(
io_phase_angle(now()), jupiter_cml(o, now()))::numeric, 2) AS burst_prob
FROM (SELECT '40.0N 105.3W 1655m'::observer) AS cfg(o);
```
Because the view uses `now()`, every `SELECT` recomputes against the current time — no refresh needed. The conditional aggregation approach (one column per planet) avoids `tablefunc`/`crosstab` entirely. Change the observer literal to your coordinates.
<Aside type="tip" title="Parameterized version">
For multiple observers, replace the literal with a function parameter or a lookup table:
```sql
SELECT s.* FROM observers o,
LATERAL (SELECT * FROM solar_system_at(o.loc, now())) s;
```
That requires wrapping the view logic in a `CREATE FUNCTION`, but the pattern is the same.
</Aside>
### 9. Satellite Shell Census
How many satellites occupy each orbital shell? Compute altitude from the TLE's mean motion using Kepler's third law ($a = (\mu / n^2)^{1/3}$, altitude $= a - R_\oplus$), then classify into LEO/MEO/GEO/HEO:
```sql
WITH altitudes AS (
SELECT norad_id, name,
power(
398600.8 / power(tle_mean_motion(tle) * 2 * pi() / 86400.0, 2),
1.0 / 3.0
) - 6378.135 AS alt_km
FROM satellites
WHERE tle_mean_motion(tle) > 0
)
SELECT
CASE
WHEN alt_km < 2000 THEN 'LEO'
WHEN alt_km < 35786 THEN 'MEO'
WHEN alt_km < 35800 THEN 'GEO'
ELSE 'HEO/Other'
END AS shell,
count(*) AS n_satellites,
round(100.0 * count(*) / sum(count(*)) OVER (), 1) AS pct,
round(min(alt_km)::numeric, 0) AS min_alt_km,
round(percentile_cont(0.5) WITHIN GROUP (ORDER BY alt_km)::numeric, 0) AS median_alt_km,
round(max(alt_km)::numeric, 0) AS max_alt_km
FROM altitudes
WHERE alt_km > 100 -- filter decayed objects
GROUP BY
CASE
WHEN alt_km < 2000 THEN 'LEO'
WHEN alt_km < 35786 THEN 'MEO'
WHEN alt_km < 35800 THEN 'GEO'
ELSE 'HEO/Other'
END
ORDER BY min_alt_km;
```
The 398600.8 is WGS-72 $\mu$ (km³/s²) and 6378.135 is WGS-72 $a_e$ (km) — the same constants SGP4 uses internally. The `percentile_cont(0.5)` gives the median altitude per shell, which is more informative than the mean when distributions are skewed (LEO has a long tail from Molniya-type parking orbits).
For a finer-grained altitude histogram within LEO — revealing the Starlink, ISS, sun-synchronous, and Iridium clusters:
```sql
WITH altitudes AS (
SELECT power(
398600.8 / power(tle_mean_motion(tle) * 2 * pi() / 86400.0, 2),
1.0 / 3.0
) - 6378.135 AS alt_km
FROM satellites
WHERE tle_mean_motion(tle) > 0
)
SELECT round((150 + (bucket - 1) * 10)::numeric, 0) AS alt_km,
count(*) AS n_satellites
FROM (
SELECT width_bucket(alt_km, 150, 2050, 190) AS bucket
FROM altitudes
WHERE alt_km BETWEEN 150 AND 2050
) sub
GROUP BY bucket
ORDER BY alt_km;
```
Plot `alt_km` vs `n_satellites` and you'll see pronounced peaks: a massive spike near 550 km (Starlink's operational shell), a cluster around 780 km (Iridium NEXT), concentrations at 500600 km (sun-synchronous polar orbits), and smaller peaks near 400 km (crewed missions) and 1200 km (older constellations).

View File

@ -6,7 +6,7 @@ sidebar:
import { Aside, Tabs, TabItem } from "@astrojs/starlight/components";
Measured performance numbers for pg_orrery's core operations. Every number on this page was produced by running the listed SQL query against a live PostgreSQL 17 instance with a single backend, no parallel workers, and no connection pooling overhead.
Measured performance numbers for pg_orrery's core operations. Every number on this page was produced by running the listed SQL query against a live PostgreSQL 18 instance with a single backend, no parallel workers, and no connection pooling overhead.
<Aside type="note" title="Methodology">
All benchmarks use PostgreSQL's `EXPLAIN (ANALYZE, BUFFERS)` for timing. The numbers are wall-clock execution time for the query, not per-function overhead. Each benchmark was run three times; the reported value is the median. Cold start was avoided by running each query once before measurement.
@ -17,7 +17,9 @@ All benchmarks use PostgreSQL's `EXPLAIN (ANALYZE, BUFFERS)` for timing. The num
| Operation | Count | Time | Rate | Notes |
|-----------|-------|------|------|-------|
| TLE propagation (SGP4) | 12,000 | 17 ms | 706K/sec | Mixed LEO/MEO/GEO |
| Visibility cone filter (`&?`) | 65,886 | 12.5 ms | 5.3M/sec | 80% pruned, no SGP4 |
| Visibility cone filter (`&?`) | 66,440 | 12.1 ms | 5.5M/sec | 84% pruned (2h, 10°), no SGP4 |
| Conjunction screening (`&&`) | 66,440 | 4.6 ms | — | ISS: 9 co-orbital objects found |
| KNN altitude ordering (`<->`) | 66,440 | 2.1 ms | — | 10 nearest to ISS, index-ordered |
| Planet observation (VSOP87) | 875 | 57 ms | 15.4K/sec | All 7 non-Earth planets, 125 times each |
| Galilean moon observation | 1,000 | 63 ms | 15.9K/sec | L1.2 + VSOP87 pipeline |
| Saturn moon observation | 800 | 53 ms | 15.1K/sec | TASS17 + VSOP87 |
@ -25,7 +27,7 @@ All benchmarks use PostgreSQL's `EXPLAIN (ANALYZE, BUFFERS)` for timing. The num
| Lambert transfer solve | 100 | 0.1 ms | 800K/sec | Single-rev prograde |
| Pork chop plot (150 x 150) | 22,500 | 8.3 s | 2.7K/sec | Full VSOP87 + Lambert pipeline |
**Conditions:** PostgreSQL 17.2, single backend, no parallel workers, Intel Xeon E-2286G @ 4.0 GHz, 64 GB ECC DDR4-2666. Extension compiled with GCC 14.2, `-O2`.
**Conditions:** PostgreSQL 18.1, single backend, no parallel workers, Intel Xeon E-2286G @ 4.0 GHz, 64 GB ECC DDR4-2666. Extension compiled with GCC 14.2, `-O2`.
## TLE propagation
@ -225,7 +227,7 @@ A 7-day window at 30-second coarse scan resolution requires ~20,160 propagation
The `&?` operator answers "could this satellite possibly be visible from this observer?" using three geometric filters (altitude, inclination, RAAN) without any SGP4 propagation. This is the first stage of the pass prediction pipeline, reducing the number of satellites that need full propagation.
```sql
-- Benchmark: filter a 66,000-object catalog
-- Benchmark: filter a 66,440-object catalog
-- Eagle, Idaho: 2-hour window, 10 deg minimum elevation
EXPLAIN (ANALYZE, BUFFERS)
SELECT count(*)
@ -238,53 +240,133 @@ WHERE tle &? ROW(
)::observer_window;
```
**65,886 TLEs filtered in 12.5 ms --- 80% pruned, 12,964 candidates survive.**
**66,440 TLEs filtered in 12.1 ms --- 83.8% pruned, 10,763 candidates survive.**
The operator evaluates three geometric conditions per TLE: perigee altitude vs. maximum visible altitude, inclination + ground footprint vs. observer latitude, and RAAN alignment via J2 secular precession. Each check is a few floating-point operations --- no SGP4 initialization, no Kepler equation, no trigonometric series.
### Pruning rate by query pattern
The pruning rate depends on observer latitude, query window duration, and minimum elevation. Shorter windows and higher latitudes prune more aggressively.
Measured against a 66,440-object catalog merged from Space-Track, CelesTrak, SatNOGS, and CelesTrak SupGP. The pruning rate depends on observer latitude, query window duration, and minimum elevation. Shorter windows and higher minimum elevations prune more aggressively.
| Query | Candidates | Pruned | Notes |
|-------|-----------|--------|-------|
| 2h, Eagle ID (43.7°N), 10° | 12,964 | 80.3% | Typical mid-latitude evening |
| 2h, Tromsoe (69.6°N), 10° | 6,529 | 90.1% | High latitude: inclination filter strongest |
| 2h, South Pole (85°S), 10° | 5,248 | 92.0% | Only polar-orbit satellites survive |
| 2h, Equator (0°N), 10° | 9,699 | 85.3% | All inclinations pass latitude check; RAAN filter dominates |
| 2h, Eagle ID, 30° | 9,680 | 85.3% | Higher elevation: altitude filter tighter |
| 6h, Eagle ID, 10° | 24,274 | 63.2% | Wider RAAN window admits more candidates |
| 24h, Eagle ID, 10° | 60,875 | 7.6% | RAAN filter bypassed (full Earth rotation) |
| 2h, Eagle ID (43.7°N), 10° | 10,763 | 83.8% | Typical mid-latitude evening |
| 2h, Equator (0°N), 10° | 10,174 | 84.7% | All inclinations pass latitude check; RAAN filter dominates |
| 2h, Eagle ID, 45° | 6,796 | 89.8% | Higher elevation: altitude filter tighter |
| 24h, Eagle ID, 10° | 61,426 | 7.5% | RAAN filter bypassed (full Earth rotation) |
### SP-GiST index performance
The optional SP-GiST index (`tle_spgist_ops`) builds a 2-level trie partitioned by semi-major axis and inclination. At 66,000 objects, the index adds 1--2 ms overhead compared to a sequential scan for mid-latitude observers, but **beats the sequential scan for high-latitude queries** where inclination pruning eliminates entire subtrees:
The optional SP-GiST index (`tle_spgist_ops`) builds a 2-level trie partitioned by semi-major axis and inclination. At 66,440 objects, sequential evaluation of the `&?` operator (12 ms) is faster than the SP-GiST index scan (16--23 ms). The tree traversal overhead exceeds the pruning benefit at this catalog size because the `&?` operator itself is so cheap --- three floating-point comparisons per TLE.
| Query | Seqscan | SP-GiST | Difference |
|-------|---------|---------|------------|
| 2h, Eagle ID, 10° | 12.5 ms | 14.0 ms | +1.5 ms |
| 2h, Tromsoe, 10° | 11.3 ms | 10.9 ms | **-0.4 ms** |
| Query | Seqscan | SP-GiST | Candidates | Pruned |
|-------|---------|---------|------------|--------|
| 2h, Eagle ID, 10° | 12.1 ms | 16.1 ms | 10,763 | 83.8% |
| 2h, Equator, 10° | 12.1 ms | 16.8 ms | 10,174 | 84.7% |
| 2h, Eagle ID, 45° | 11.9 ms | 16.9 ms | 6,796 | 89.8% |
| 24h, Eagle ID, 10° | 12.5 ms | 23.3 ms | 61,426 | 7.5% |
The planner chooses the SP-GiST Index Only Scan by default at this catalog size, with zero heap fetches (all data served from index pages).
The SP-GiST index achieves zero heap fetches (pure Index Only Scan), but page traversal through 11 MB of index data (4,964 buffer hits) exceeds the cost of a 1,338-buffer sequential scan.
<Aside type="tip" title="Where the index shines">
The SP-GiST index is most valuable for high-latitude observers (60°+) and for catalogs larger than 30,000 objects. At typical CelesTrak catalog sizes (12--15,000 active satellites), the `&?` operator's sequential evaluation is fast enough that the index overhead exceeds the pruning benefit.
<Aside type="tip" title="Where the SP-GiST index adds value">
The `&?` operator prunes 84--90% of the catalog regardless of scan method. Its primary value is as a **gating filter** before expensive SGP4 propagation. For a 2-hour window, reducing 66,440 TLEs to ~10,000 candidates saves ~56,000 `predict_passes()` calls (each ~1 ms), a far greater benefit than the 4 ms difference between scan methods.
At larger catalog sizes (200k+ objects), the SP-GiST tree-level pruning should begin to outperform sequential evaluation. The crossover point depends on hardware, but the operator's pruning ratio is the dominant performance factor, not the scan method.
</Aside>
### What the pruning means for predict_passes()
For a 65,886-object catalog and a 2-hour window from Eagle, Idaho:
For a 66,440-object catalog and a 2-hour window from Eagle, Idaho:
- **Without `&?`:** 65,886 `predict_passes()` calls (each ~1 ms for a 7-day window)
- **With `&?`:** 12,964 calls --- **52,922 unnecessary propagations avoided**
- **Time saved:** ~53 seconds per query at typical propagation cost
- **Without `&?`:** 66,440 `predict_passes()` calls (each ~1 ms for a 7-day window)
- **With `&?`:** 10,763 calls --- **55,677 unnecessary propagations avoided**
- **Time saved:** ~56 seconds per query at typical propagation cost
## Conjunction screening (`&&` operator)
The GiST index on the `tle` type enables indexed conjunction screening using the `&&` (overlap) operator. The index stores altitude band and inclination for each TLE, allowing PostgreSQL to skip entire subtrees of non-overlapping orbits.
```sql
-- Benchmark: find ISS conjunction candidates in a 66,440-object catalog
EXPLAIN (ANALYZE, BUFFERS)
SELECT b.name
FROM satellite_catalog a
JOIN satellite_catalog b ON a.tle && b.tle AND a.norad_id != b.norad_id
WHERE tle_norad_id(a.tle) = 25544;
```
**9 co-orbital objects found in 4.6 ms (vs. 63.3 ms sequential scan --- 5.8x speedup).**
The GiST index scan hits 237 buffers compared to 1,338 for a sequential scan. The 9 objects returned are all ISS-visiting vehicles or co-orbital modules: PROGRESS MS-31, PROGRESS MS-32, SOYUZ MS-28, DRAGON FREEDOM 3, DRAGON CRS-33, CYGNUS NG-23, HTV-X1, ISS (NAUKA), and OBJECT E.
### GiST `&&` performance by orbital regime
| Probe satellite | GiST | Seqscan | Matches | Notes |
|----------------|------|---------|---------|-------|
| ISS (LEO, 51.6°) | 4.6 ms | 63.3 ms | 9 | Co-orbital vehicles |
| Starlink-230369 (LEO, 53°) | 9.5 ms | 14.9 ms | 0 | Dense LEO shell |
| SYNCOM 2 (GEO, 33°) | 4.0 ms | 7.2 ms | 0 | Sparse regime |
The GiST index provides the largest speedup for queries that return few matches, where the index prunes most of the tree without reading leaf pages. Dense LEO shells produce more candidates and reduce the speedup ratio.
### Index characteristics
| Metric | Value |
|--------|-------|
| Build time | 93 ms (66,440 TLEs) |
| Index size | 15 MB (237 bytes/object) |
| Consistency | 0 false positives, 0 false negatives (verified against seqscan) |
## KNN altitude ordering (`<->` operator)
The `<->` operator computes altitude-band separation in km. With a GiST index, it supports index-ordered KNN queries --- PostgreSQL traverses the tree by increasing distance without computing all distances upfront.
```sql
-- Benchmark: 10 nearest orbits to the ISS by altitude separation
EXPLAIN (ANALYZE, BUFFERS)
SELECT name,
round((tle <-> (SELECT tle FROM satellite_catalog
WHERE tle_norad_id(tle) = 25544 LIMIT 1))::numeric, 1)
AS alt_sep_km
FROM satellite_catalog
ORDER BY tle <-> (SELECT tle FROM satellite_catalog
WHERE tle_norad_id(tle) = 25544 LIMIT 1)
LIMIT 10;
```
**10 nearest in 2.1 ms, index-ordered (982 buffer hits).**
### KNN performance by scenario
| Query | Time | Buffers | Notes |
|-------|------|---------|-------|
| 10 nearest to ISS (LEO) | 2.1 ms | 982 | Dense regime, more nodes traversed |
| 10 nearest to SYNCOM 2 (GEO) | 0.2 ms | 40 | Sparse regime, fewer nodes |
| 100 nearest to ISS | 1.4 ms | 1,062 | Marginal cost per additional neighbor |
| All within 50 km of ISS | 16.0 ms | 4,014 | 12,496 matches |
<Aside type="caution" title="KNN requires a scalar subquery probe">
GiST index-ordered scans only activate when the probe value is visible to the planner as a constant. Use a **scalar subquery** for the probe TLE:
```sql
-- This uses the index (scalar subquery → constant to planner):
ORDER BY tle <-> (SELECT tle FROM catalog WHERE tle_norad_id(tle) = 25544 LIMIT 1)
-- This does NOT use the index (CTE is opaque to the planner):
WITH iss AS (SELECT tle FROM catalog WHERE tle_norad_id(tle) = 25544)
SELECT ... ORDER BY tle <-> iss.tle -- falls back to full scan + sort
```
The CTE pattern works correctly but forces PostgreSQL to compute all distances and sort, which is much slower for large catalogs. For small catalogs (< 100 rows), the difference is negligible.
</Aside>
## Reproducing these benchmarks
<Tabs>
<TabItem label="Requirements">
- PostgreSQL 17 with pg_orrery installed
- A satellite catalog table with ~12,000 TLEs (see [Building TLE Catalogs](/guides/catalog-management/) or download directly from CelesTrak)
- PostgreSQL 18 with pg_orrery installed
- A satellite catalog table (the numbers on this page use a 66,440-object catalog merged from Space-Track, CelesTrak, SatNOGS, and CelesTrak SupGP; see [Building TLE Catalogs](/guides/catalog-management/))
- GiST and SP-GiST indexes on the `tle` column for index benchmarks
- A star catalog table (any subset of Hipparcos or Yale BSC)
- No concurrent queries during measurement
- `shared_buffers` and `work_mem` at default or higher
@ -295,12 +377,16 @@ For a 65,886-object catalog and a 2-hour window from Eagle, Idaho:
-- Load a TLE catalog (pg-orrery-catalog handles this)
-- pg-orrery-catalog build --table satellite_catalog | psql -d mydb
CREATE TABLE satellite_catalog (tle tle);
CREATE TABLE satellite_catalog (name text, tle tle);
-- (or COPY from CelesTrak bulk TLE file)
-- Create both indexes for full benchmark coverage
CREATE INDEX idx_tle_gist ON satellite_catalog USING gist (tle tle_ops);
CREATE INDEX idx_tle_spgist ON satellite_catalog USING spgist (tle tle_spgist_ops);
-- Verify catalog size
SELECT count(*) FROM satellite_catalog;
-- Expected: ~12,000 rows
-- The numbers on this page use 66,440 rows
-- Disable parallel workers for baseline measurement
SET max_parallel_workers_per_gather = 0;
@ -326,6 +412,8 @@ For a 65,886-object catalog and a 2-hour window from Eagle, Idaho:
The benchmarks demonstrate that pg_orrery's computation cost is low enough to treat orbital mechanics as a SQL primitive. Propagating an entire satellite catalog takes less time than a typical index scan on a moderately-sized table. Planet observation is fast enough to generate ephemeris tables with `generate_series`. Pork chop plots are feasible as interactive queries rather than batch jobs.
The visibility cone filter (`&?`) is the fastest operation per evaluation --- three floating-point comparisons vs. the full SGP4 pipeline --- and its 80--92% pruning rate means the most expensive operation in a pass prediction pipeline (SGP4 propagation) only runs on the small fraction of the catalog that could actually produce a visible pass.
The visibility cone filter (`&?`) is the fastest operation per evaluation --- three floating-point comparisons vs. the full SGP4 pipeline --- and its 84--90% pruning rate means the most expensive operation in a pass prediction pipeline (SGP4 propagation) only runs on the small fraction of the catalog that could actually produce a visible pass.
The numbers also show where the bottlenecks are: VSOP87 series evaluation dominates everything except star observation, raw SGP4 propagation, and the visibility cone filter. If a future optimization effort targets one component, it should be the VSOP87 evaluation loop.
The GiST index provides the clearest speedup for conjunction screening: 5.8x faster than sequential scan for ISS `&&` queries, with 0 false positives or negatives verified against exhaustive sequential evaluation. KNN queries find the nearest orbits in 2 ms via index-ordered traversal, which would otherwise require computing and sorting all 66,440 distances.
The numbers also show where the bottlenecks are: VSOP87 series evaluation dominates everything except star observation, raw SGP4 propagation, and the geometric filters. If a future optimization effort targets one component, it should be the VSOP87 evaluation loop.

View File

@ -1,12 +1,12 @@
---
title: "Functions: Stars & Comets"
title: "Functions: Stars, Comets & Asteroids"
sidebar:
order: 5
---
import { Aside, Tabs, TabItem } from "@astrojs/starlight/components";
Functions for computing topocentric positions of stars from catalog coordinates, propagating comets and asteroids on Keplerian orbits, and observing them from Earth.
Functions for computing topocentric positions of stars from catalog coordinates, propagating comets and asteroids on Keplerian orbits, and observing them from Earth. The `orbital_elements` type (v0.8.0) bundles Keplerian elements into a first-class PostgreSQL datum, with `oe_from_mpc()` for bulk-loading the MPC catalog and `small_body_observe()` for ergonomic topocentric observation.
---
@ -255,3 +255,144 @@ FROM comet_catalog, earth,
WHERE topo_elevation(c) > 0
ORDER BY topo_range(c);
```
---
## oe_from_mpc
Parses one line of the MPC MPCORB.DAT fixed-width format into an `orbital_elements` type. The MPC publishes orbital elements for over 1.3 million numbered asteroids in this format.
### Signature
```sql
oe_from_mpc(line text) → orbital_elements
```
### Parameters
| Parameter | Type | Description |
|-----------|------|-------------|
| `line` | `text` | One complete line from MPCORB.DAT (at least 103 characters) |
### Returns
An `orbital_elements` with all nine fields populated. The parser performs several conversions at parse time:
- **Packed epoch** (e.g. `K24AM`) is decoded to a Julian date. The century letter (`I`=1800, `J`=1900, `K`=2000), two-digit year, packed month (`1-9`, `A-C`), and packed day (`1-9`, `A-V`) are expanded to a calendar date.
- **Perihelion distance** is derived from the MPC's semi-major axis and eccentricity: q = a × (1 e).
- **Perihelion time** is computed from the epoch and mean anomaly via Gauss's constant: tp = epoch M / n, where n = k / a^(3/2).
<Aside type="tip">
The full MPCORB.DAT format is documented at the [IAU Minor Planet Center](https://www.minorplanetcenter.net/iau/info/MPOrbitFormat.html). The file is freely downloadable (~280 MB compressed) and contains orbital elements for all numbered and multi-opposition asteroids.
</Aside>
### Example
```sql
-- Parse Ceres, extract semi-major axis and period
WITH ceres AS (
SELECT oe_from_mpc(
'00001 3.33 0.12 K24AM 60.07966 73.42937 80.26860 10.58664 0.0789126 0.21406048 2.7660961 0 MPO838504 8738 115 1801-2024 0.65 M-v 30k MPCLINUX 0000 (1) Ceres 20240825'
) AS oe
)
SELECT round(oe_semi_major_axis(oe)::numeric, 4) AS a_au,
round(oe_period_years(oe)::numeric, 2) AS period_yr,
round(oe_inclination(oe)::numeric, 3) AS inc_deg,
round(oe_h_mag(oe)::numeric, 2) AS h_mag
FROM ceres;
```
---
## small_body_heliocentric
Propagates an `orbital_elements` to a heliocentric ecliptic J2000 position at a given time using two-body Keplerian mechanics. Extracts q, e, inc, omega, Omega, and tp from the type and calls the internal Kepler solver.
### Signature
```sql
small_body_heliocentric(oe orbital_elements, t timestamptz) → heliocentric
```
### Parameters
| Parameter | Type | Description |
|-----------|------|-------------|
| `oe` | `orbital_elements` | Bundled orbital elements |
| `t` | `timestamptz` | Evaluation time |
### Returns
A `heliocentric` position in AU (ecliptic J2000 frame). Identical to calling `kepler_propagate()` with the individual fields extracted from the type.
### Example
```sql
-- Propagate Ceres to 2025-01-01, check heliocentric distance
WITH ceres AS (
SELECT oe_from_mpc(
'00001 3.33 0.12 K24AM 60.07966 73.42937 80.26860 10.58664 0.0789126 0.21406048 2.7660961 0 MPO838504 8738 115 1801-2024 0.65 M-v 30k MPCLINUX 0000 (1) Ceres 20240825'
) AS oe
)
SELECT round(helio_x(h)::numeric, 4) AS x_au,
round(helio_y(h)::numeric, 4) AS y_au,
round(helio_z(h)::numeric, 4) AS z_au,
round(helio_distance(h)::numeric, 3) AS dist_au
FROM ceres,
small_body_heliocentric(oe, '2025-01-01 00:00:00+00') AS h;
```
---
## small_body_observe
Computes the topocentric position of a comet or asteroid from its `orbital_elements` as seen by an Earth-based observer. Auto-fetches Earth's heliocentric position via VSOP87, matching the ergonomics of `planet_observe()`.
### Signature
```sql
small_body_observe(oe orbital_elements, obs observer, t timestamptz) → topocentric
```
### Parameters
| Parameter | Type | Description |
|-----------|------|-------------|
| `oe` | `orbital_elements` | Bundled orbital elements |
| `obs` | `observer` | Observer location on Earth |
| `t` | `timestamptz` | Observation time |
### Returns
A `topocentric` with azimuth, elevation (degrees), range (km), and range rate (km/s).
<Aside type="note">
Unlike `comet_observe()` which requires you to pass Earth's position as three separate floats, `small_body_observe()` fetches Earth internally via VSOP87. This is simpler for single-object queries. For batch observations at the same time, `comet_observe()` still lets you compute Earth's position once and reuse it.
</Aside>
### Example
```sql
-- Observe Ceres from Boulder
WITH ceres AS (
SELECT oe_from_mpc(
'00001 3.33 0.12 K24AM 60.07966 73.42937 80.26860 10.58664 0.0789126 0.21406048 2.7660961 0 MPO838504 8738 115 1801-2024 0.65 M-v 30k MPCLINUX 0000 (1) Ceres 20240825'
) AS oe
)
SELECT round(topo_azimuth(t)::numeric, 1) AS az,
round(topo_elevation(t)::numeric, 1) AS el,
round((topo_range(t) / 149597870.7)::numeric, 3) AS dist_au
FROM ceres,
small_body_observe(oe, '40.0N 105.3W 1655m'::observer, now()) AS t;
```
```sql
-- Which asteroids in the catalog are above 20 degrees tonight?
SELECT name,
round(topo_elevation(t)::numeric, 1) AS el,
round(topo_azimuth(t)::numeric, 1) AS az
FROM asteroid_catalog,
small_body_observe(oe, '40.0N 105.3W 1655m'::observer, now()) AS t
WHERE topo_elevation(t) > 20
ORDER BY topo_elevation(t) DESC;
```

View File

@ -134,17 +134,23 @@ WHERE c.tle && iss.tle
<TabItem label="kNN by altitude">
```sql
-- Find the 10 satellites with the closest altitude bands to the ISS
-- The <-> operator supports GiST ordering (ORDER BY ... <-> ...)
WITH iss AS (
SELECT tle FROM satellite_catalog WHERE norad_id = 25544
)
SELECT c.norad_id, c.name,
round((c.tle <-> iss.tle)::numeric, 1) AS alt_sep_km
FROM satellite_catalog c, iss
WHERE c.norad_id != 25544
ORDER BY c.tle <-> iss.tle
-- The <-> operator supports GiST index ordering (ORDER BY ... <-> ...)
-- IMPORTANT: use a scalar subquery for the probe TLE so the planner
-- can see it as a constant and activate index-ordered scan.
SELECT c.name,
round((c.tle <-> (SELECT tle FROM satellite_catalog
WHERE tle_norad_id(tle) = 25544 LIMIT 1))::numeric, 1)
AS alt_sep_km
FROM satellite_catalog c
WHERE tle_norad_id(c.tle) != 25544
ORDER BY c.tle <-> (SELECT tle FROM satellite_catalog
WHERE tle_norad_id(tle) = 25544 LIMIT 1)
LIMIT 10;
```
<Aside type="caution" title="CTE pattern prevents index ordering">
A CTE like `WITH iss AS (SELECT tle ...)` makes the probe value opaque to the planner, forcing a full sequential scan and sort instead of an index-ordered traversal. Always use a scalar subquery `(SELECT tle FROM ... LIMIT 1)` for the probe argument. For small catalogs (< 100 rows) the difference is negligible; for large catalogs it is the difference between 2 ms and a full sort.
</Aside>
</TabItem>
<TabItem label="Two-stage screening">
```sql
@ -171,10 +177,18 @@ ORDER BY dist_km;
### Performance
Without the GiST index, the `&&` operator requires a sequential scan of the entire catalog (O(n) per query). With the index, overlap queries run in O(log n) time. For a catalog of 12,000 active TLEs, this reduces conjunction screening from seconds to milliseconds.
Benchmarked against a 66,440-object catalog (Space-Track + CelesTrak + SatNOGS):
| Query | GiST | Seqscan | Matches | Speedup |
|-------|------|---------|---------|---------|
| ISS conjunction (`&&`) | 4.6 ms | 63.3 ms | 9 | 5.8x |
| 10 nearest to ISS (`<->` KNN) | 2.1 ms | — | 10 | Index-ordered |
| 10 nearest to GEO sat (`<->` KNN) | 0.2 ms | — | 10 | Sparse regime |
The GiST index (15 MB, 93 ms build) provides the clearest speedup for conjunction screening. The `&&` operator reduces the search from 1,338 buffer hits (sequential scan) to 237 buffer hits (index scan). KNN queries traverse the tree by increasing distance without computing all distances upfront.
<Aside type="tip">
The GiST index is most valuable for large catalogs (thousands of TLEs). For small catalogs (< 100 TLEs), sequential scans may be faster than the index overhead. PostgreSQL's query planner handles this decision automatically.
For small catalogs (< 100 TLEs), sequential scans may be faster than the index overhead. PostgreSQL's query planner handles this decision automatically. The GiST index shows the largest relative speedup when the query returns few matches against a large catalog --- exactly the conjunction screening pattern.
</Aside>
### Index Maintenance
@ -288,15 +302,24 @@ During the index scan, inner nodes are pruned by altitude band (level 0) and inc
### Performance
The `&?` operator eliminates 80-90% of a satellite catalog without SGP4 propagation --- this is the primary value, regardless of whether a sequential scan or index scan evaluates it. At typical catalog sizes (10-30k objects), the operator evaluates the full catalog in under 10 ms, and PostgreSQL's query planner may choose a sequential scan over the index.
The `&?` operator eliminates 84--90% of a satellite catalog without SGP4 propagation --- this is the primary value, regardless of whether a sequential scan or index scan evaluates it.
The SP-GiST index becomes advantageous at larger catalog sizes (100k+ objects) where tree-level pruning avoids examining individual TLEs in entire subtrees. The index is most effective for:
Benchmarked against a 66,440-object catalog:
| Query | Seqscan | SP-GiST | Candidates | Pruned |
|-------|---------|---------|------------|--------|
| 2h, Eagle ID, 10° | 12.1 ms | 16.1 ms | 10,763 | 83.8% |
| 2h, Equator, 10° | 12.1 ms | 16.8 ms | 10,174 | 84.7% |
| 2h, Eagle ID, 45° | 11.9 ms | 16.9 ms | 6,796 | 89.8% |
| 24h, Eagle ID, 10° | 12.5 ms | 23.3 ms | 61,426 | 7.5% |
At 66k objects, the sequential scan is faster than the SP-GiST index for all tested scenarios. The `&?` operator is so cheap per evaluation (three floating-point comparisons) that tree traversal overhead exceeds the pruning benefit at this catalog size. The index is most effective for:
- **Short query windows** (1-6 hours): The RAAN filter aggressively eliminates satellites whose orbital plane is not currently aligned with the observer
- **Mid-latitude observers** (30-60 degrees): The inclination filter eliminates equatorial and low-inclination satellites
- **High minimum elevation** (> 20 degrees): The altitude filter eliminates distant MEO/GEO objects
- **Higher minimum elevation** (> 20 degrees): The altitude filter eliminates distant MEO/GEO objects
- **Larger catalogs** (200k+ objects): Tree-level pruning avoids examining individual TLEs in entire subtrees
For 24-hour query windows, the RAAN filter self-disables (full Earth rotation makes it meaningless), and only the altitude and inclination filters apply.
For 24-hour query windows, the RAAN filter self-disables (full Earth rotation makes it meaningless), and only the altitude and inclination filters apply. The real value of the `&?` operator is as a gating filter before expensive SGP4 propagation, not the scan method itself.
### Index Maintenance

View File

@ -6,7 +6,7 @@ sidebar:
import { Aside, Tabs, TabItem } from "@astrojs/starlight/components";
pg_orrery defines seven fixed-size base types and one SQL composite type that represent the core data structures of orbital mechanics. Each base type has a fixed on-disk size, a text I/O format for readability, and accessor functions for extracting individual fields.
pg_orrery defines eight fixed-size base types and one SQL composite type that represent the core data structures of orbital mechanics. Each base type has a fixed on-disk size, a text I/O format for readability, and accessor functions for extracting individual fields.
## tle
@ -270,6 +270,74 @@ FROM generate_series(1, 8) AS body_id,
---
## orbital_elements
**Size:** 72 bytes
Classical Keplerian orbital elements for comets and asteroids. Stores nine doubles: osculation epoch, perihelion distance, eccentricity, inclination, argument of perihelion, longitude of ascending node, time of perihelion passage, absolute magnitude H, and slope parameter G.
### Input Format
A parenthesized tuple of nine values:
```sql
SELECT '(2460605.5,2.5478,0.0789126,10.58664,73.42937,80.2686,2460319.0,3.33,0.12)'::orbital_elements;
```
The fields in order are: `(epoch_jd, q_au, e, inc_deg, omega_deg, Omega_deg, tp_jd, H, G)`.
<Aside type="note">
Angular fields (inclination, argument of perihelion, RAAN) are displayed in degrees but stored internally as radians — the same convention used by the `tle` type. H and G display as `NaN` when unknown.
</Aside>
### Constructor
| Function | Signature | Description |
|----------|-----------|-------------|
| `oe_from_mpc` | `oe_from_mpc(line text) → orbital_elements` | Parses one MPCORB.DAT fixed-width line. Converts MPC packed epoch, computes q and tp from (a, e, M). |
```sql
-- Parse (1) Ceres from an MPCORB.DAT line
SELECT oe_from_mpc(
'00001 3.33 0.12 K24AM 60.07966 73.42937 80.26860 10.58664 0.0789126 0.21406048 2.7660961 0 MPO838504 8738 115 1801-2024 0.65 M-v 30k MPCLINUX 0000 (1) Ceres 20240825'
);
```
### Accessor Functions
| Function | Return Type | Description |
|----------|-------------|-------------|
| `oe_epoch(orbital_elements)` | `float8` | Osculation epoch (Julian date) |
| `oe_perihelion(orbital_elements)` | `float8` | Perihelion distance q (AU) |
| `oe_eccentricity(orbital_elements)` | `float8` | Eccentricity |
| `oe_inclination(orbital_elements)` | `float8` | Inclination (degrees) |
| `oe_arg_perihelion(orbital_elements)` | `float8` | Argument of perihelion (degrees) |
| `oe_raan(orbital_elements)` | `float8` | Longitude of ascending node (degrees) |
| `oe_tp(orbital_elements)` | `float8` | Time of perihelion passage (Julian date) |
| `oe_h_mag(orbital_elements)` | `float8` | Absolute magnitude H (NaN if unknown) |
| `oe_g_slope(orbital_elements)` | `float8` | Slope parameter G (NaN if unknown) |
| `oe_semi_major_axis(orbital_elements)` | `float8` | Semi-major axis a = q/(1-e) in AU. NULL for e ≥ 1. |
| `oe_period_years(orbital_elements)` | `float8` | Orbital period a^1.5 in years. NULL for e ≥ 1. |
```sql
-- Parse Ceres and extract key parameters
WITH ceres AS (
SELECT oe_from_mpc(
'00001 3.33 0.12 K24AM 60.07966 73.42937 80.26860 10.58664 0.0789126 0.21406048 2.7660961 0 MPO838504 8738 115 1801-2024 0.65 M-v 30k MPCLINUX 0000 (1) Ceres 20240825'
) AS oe
)
SELECT oe_epoch(oe) AS epoch_jd,
oe_perihelion(oe) AS q_au,
oe_eccentricity(oe) AS ecc,
oe_inclination(oe) AS inc_deg,
oe_semi_major_axis(oe) AS a_au,
oe_period_years(oe) AS period_yr,
oe_h_mag(oe) AS h_mag
FROM ceres;
```
---
## observer_window
**Type:** SQL composite (variable-length)

View File

@ -44,6 +44,8 @@ PG_FUNCTION_INFO_V1(gist_tle_distance);
/* Floating-point comparison tolerance (km and radians) */
#define KEY_EPSILON 1.0e-9
/* sizeof(pg_tle) == 112, matching INTERNALLENGTH in CREATE TYPE. */
/*
* 2-D orbital key extracted from a TLE's mean elements.
* Altitude band (perigee/apogee) plus inclination range.
@ -231,6 +233,10 @@ tle_alt_distance(PG_FUNCTION_ARGS)
*
* Leaf entries carry the full pg_tle; we compress to tle_orbital_key.
* Internal entries are already tle_orbital_key from union operations.
*
* The allocation must be sizeof(pg_tle) bytes which matches
* INTERNALLENGTH not sizeof(tle_orbital_key). GiST's
* index_form_tuple() copies typlen bytes from the datum pointer.
*/
Datum
gist_tle_compress(PG_FUNCTION_ARGS)
@ -241,7 +247,7 @@ gist_tle_compress(PG_FUNCTION_ARGS)
if (entry->leafkey)
{
pg_tle *tle = (pg_tle *) DatumGetPointer(entry->key);
tle_orbital_key *key = (tle_orbital_key *) palloc(sizeof(tle_orbital_key));
tle_orbital_key *key = (tle_orbital_key *) palloc0(sizeof(pg_tle));
tle_to_orbital_key(tle, key);
@ -273,8 +279,10 @@ gist_tle_decompress(PG_FUNCTION_ARGS)
* gist_tle_consistent -- can this subtree contain matches for the query?
*
* Checks overlap in both altitude AND inclination dimensions.
* Always sets recheck = true because 2-D overlap is only a necessary
* condition -- the real conjunction test requires propagation.
*
* For leaf entries, recheck = false: the orbital key is computed
* identically to the SQL operator, so the GiST check is exact.
* For internal nodes, recheck is irrelevant (GiST ignores it).
*/
Datum
gist_tle_consistent(PG_FUNCTION_ARGS)
@ -290,7 +298,12 @@ gist_tle_consistent(PG_FUNCTION_ARGS)
tle_to_orbital_key(query, &query_key);
*recheck = true;
/*
* Leaf keys are exact (same tle_to_orbital_key as the operator),
* so no recheck needed. For internal nodes PostgreSQL ignores
* the flag, but we set true by convention.
*/
*recheck = !GIST_LEAF(entry);
switch (strategy)
{
@ -328,9 +341,8 @@ gist_tle_consistent(PG_FUNCTION_ARGS)
*
* The union is [min(alt_low), max(alt_high)] x [min(inc_low), max(inc_high)].
*
* GiST convention: entryvec->vector[] is 1-based (FirstOffsetNumber),
* vector[0] is unused. entryvec->n includes the unused slot, so
* valid indices are 1 .. entryvec->n - 1.
* The entry vector is 0-based: valid indices are 0 .. entryvec->n - 1.
* This differs from picksplit's 1-based convention.
*/
Datum
gist_tle_union(PG_FUNCTION_ARGS)
@ -341,17 +353,17 @@ gist_tle_union(PG_FUNCTION_ARGS)
tle_orbital_key *result;
tle_orbital_key *cur;
result = (tle_orbital_key *) palloc(sizeof(tle_orbital_key));
cur = (tle_orbital_key *) DatumGetPointer(entryvec->vector[FirstOffsetNumber].key);
result = (tle_orbital_key *) palloc0(sizeof(pg_tle));
cur = (tle_orbital_key *) DatumGetPointer(entryvec->vector[0].key);
*result = *cur;
for (i = FirstOffsetNumber + 1; i < entryvec->n; i++)
for (i = 1; i < entryvec->n; i++)
{
cur = (tle_orbital_key *) DatumGetPointer(entryvec->vector[i].key);
key_merge(result, cur);
}
*sizep = sizeof(tle_orbital_key);
*sizep = sizeof(pg_tle);
PG_RETURN_POINTER(result);
}
@ -418,11 +430,12 @@ picksplit_cmp(const void *a, const void *b)
* along whichever dimension has the greater spread. This prevents
* the tree from becoming biased toward one dimension.
*
* GiST convention: entryvec->vector[] is 1-based (FirstOffsetNumber),
* vector[0] is unused/uninitialized. entryvec->n includes the unused
* slot, so the actual entry count is (entryvec->n - 1) and valid
* indices are FirstOffsetNumber .. entryvec->n - 1. The OffsetNumbers
* placed into spl_left[] and spl_right[] must also be 1-based.
* GiST convention for picksplit: entryvec->vector[] is 1-based
* (FirstOffsetNumber), vector[0] is unused/uninitialized.
* entryvec->n includes the unused slot, so the actual entry count
* is (entryvec->n - 1) and valid indices are
* FirstOffsetNumber .. entryvec->n - 1. The OffsetNumbers placed
* into spl_left[] and spl_right[] must also be 1-based.
*/
Datum
gist_tle_picksplit(PG_FUNCTION_ARGS)
@ -495,8 +508,8 @@ gist_tle_picksplit(PG_FUNCTION_ARGS)
splitvec->spl_nright = 0;
/* Compute union keys and assign entries */
left_union = (tle_orbital_key *) palloc(sizeof(tle_orbital_key));
right_union = (tle_orbital_key *) palloc(sizeof(tle_orbital_key));
left_union = (tle_orbital_key *) palloc0(sizeof(pg_tle));
right_union = (tle_orbital_key *) palloc0(sizeof(pg_tle));
/* Seed the unions from the first entry in each half */
cur = (tle_orbital_key *) DatumGetPointer(

View File

@ -75,10 +75,13 @@ typedef struct pg_tle
char classification; /* U = unclassified */
char ephemeris_type; /* 0 = SGP4/SDP4 default */
char intl_desig[9]; /* international designator, null-terminated */
char _pad; /* alignment */
char _pad; /* alignment to int32 boundary */
char _reserved[8]; /* match INTERNALLENGTH = 112 */
} pg_tle;
/* Size: 11 doubles (88 bytes) + 3 int32 (12 bytes) + 12 chars = 112 bytes */
/* Size: 10 doubles (80) + 3 int32 (12) + 12 chars + 8 reserved = 112 bytes
* Must match INTERNALLENGTH in CREATE TYPE. PostgreSQL's datumCopy() and
* heap_form_tuple() copy exactly typlen bytes from any pg_tle pointer. */
/*