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

Harden SP-GiST orbital trie after Hamilton review

Browse Source 
Fix M2: clamp picksplit nBins to nTuples to prevent out-of-bounds
read on the entries array when called with a single tuple.

Fix H2: use WGS72_AE as effective bin_low for the first bin in L0
inner_consistent, preventing false negatives when objects with lower
SMA than the first label are inserted after index creation.

Fix H3: reject degenerate TLEs (mean_motion <= 0) early in the
visibility filter rather than propagating nonsensical values.

Fix L1: extract shared tle_passes_visibility_filter() to eliminate
duplicated 3-stage filter logic between leaf_consistent and the
standalone tle_visibility_possible operator.

Add boundary tests: degenerate TLE, polar observer, zero-duration
window, and post-insert index-vs-seqscan consistency check.
This commit is contained in:
Ryan Malloy 2026-02-17 20:36:47 -07:00
parent 6f9be428f9
commit 2a7240e739
3 changed files with 269 additions and 122 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

214
src/spgist_tle.c
View File

@ -116,8 +116,13 @@ max_visible_altitude_km(double min_el_deg)
double el_rad, sin_el; double el_rad, sin_el;
double rho, h_max; double rho, h_max;
/*
* Below 5 deg elevation the horizon geometry allows visibility to
* GEO+ altitudes. Disable the altitude filter rather than compute
* an impractically large bound. 50000 km exceeds GEO (35786 km).
*/
if (min_el_deg < 5.0) if (min_el_deg < 5.0)
return 50000.0; /* effectively no altitude filter */ return 50000.0;
el_rad = min_el_deg * DEG_TO_RAD; el_rad = min_el_deg * DEG_TO_RAD;
sin_el = sin(el_rad); sin_el = sin(el_rad);
@ -131,7 +136,7 @@ max_visible_altitude_km(double min_el_deg)
* *
* Simpler conservative bound: h_max = rho_max / sin(el) for el > 0. * Simpler conservative bound: h_max = rho_max / sin(el) for el > 0.
*/ */
rho = 5000.0; rho = 5000.0; /* max practical slant range; well beyond LEO/MEO */
h_max = sqrt(rho * rho + 2.0 * WGS72_AE * rho * sin_el h_max = sqrt(rho * rho + 2.0 * WGS72_AE * rho * sin_el
+ WGS72_AE * WGS72_AE) - WGS72_AE; + WGS72_AE * WGS72_AE) - WGS72_AE;
@ -176,6 +181,12 @@ ground_footprint_deg(double sma_km, double min_el_deg)
* *
* where n is mean motion in rad/s, J2 and Re are WGS-72. * where n is mean motion in rad/s, J2 and Re are WGS-72.
* Result in rad/day (multiply rad/s by 86400). * Result in rad/day (multiply rad/s by 86400).
*
* Uses only the J2 zonal harmonic (no J3/J4 short-period terms).
* For LEO this can accumulate ~0.5 deg/day error from J3 short-
* period oscillations. Acceptable because (a) the RAAN window
* includes footprint + Earth rotation padding, and (b) any query
* window >= ~12 hours bypasses the RAAN filter entirely.
* ---------------------------------------------------------------- * ----------------------------------------------------------------
*/ */
static inline double static inline double
@ -465,6 +476,9 @@ spgist_tle_picksplit(PG_FUNCTION_ARGS)
nBins = 2; nBins = 2;
if (nBins > 16) if (nBins > 16)
nBins = 16; nBins = 16;
/* Prevent over-read: never more bins than tuples */
if (nBins > nTuples)
nBins = nTuples;
perBin = nTuples / nBins; perBin = nTuples / nBins;
remainder = nTuples % nBins; remainder = nTuples % nBins;
@ -571,17 +585,17 @@ spgist_tle_inner_consistent(PG_FUNCTION_ARGS)
{ {
/* /*
* L0: SMA pruning. * L0: SMA pruning.
* bin_low is the SMA of the lowest object in this bin. * bin_low is the SMA of the lowest object in this bin at
* A satellite in this bin has perigee_alt >= bin_low - AE * index build time. Later inserts with lower SMA route to
* (approximately, for near-circular orbits in the bin). * bin 0, so for the first bin we use AE (physical minimum
* If the lowest possible altitude exceeds max_visible, skip. * SMA) to avoid false negatives.
* *
* Conservative: use bin_low as a lower bound on SMA, * Conservative: assume e=0 (circular, worst case for
* compute perigee assuming e=0 (circular, worst case for
* "too high" pruning -- circular has highest perigee for * "too high" pruning -- circular has highest perigee for
* a given SMA). * a given SMA).
*/ */
double perigee_alt = bin_low - WGS72_AE; double effective_low = (i == 0) ? WGS72_AE : bin_low;
double perigee_alt = effective_low - WGS72_AE;
double max_alt = max_visible_altitude_km(win.min_el_deg); double max_alt = max_visible_altitude_km(win.min_el_deg);
if (perigee_alt > max_alt) if (perigee_alt > max_alt)
@ -657,14 +671,78 @@ spgist_tle_inner_consistent(PG_FUNCTION_ARGS)
} }
/* /* ----------------------------------------------------------------
* spgist_tle_leaf_consistent -- final check on a leaf TLE * Shared filter: three-stage visibility check on a single TLE.
* *
* Applies three filters:
* 1. Perigee altitude check (with eccentricity) * 1. Perigee altitude check (with eccentricity)
* 2. Inclination + ground footprint vs observer latitude * 2. Inclination + ground footprint vs observer latitude
* 3. RAAN query-time filter (J2 precession to query midpoint) * 3. RAAN query-time filter (J2 precession to query midpoint)
* *
* Called from both leaf_consistent (index scan) and
* tle_visibility_possible (sequential scan / standalone operator).
* ----------------------------------------------------------------
*/
static bool
tle_passes_visibility_filter(const pg_tle *tle, const ObserverWindow *win)
{
double sma, perigee_alt, max_alt;
double obs_lat_abs, footprint_rad;
double dt_days, raan_projected, lst;
double earth_rot_rad, raan_window_half, raan_diff;
/* Reject degenerate TLEs (decay, error data) */
if (tle->mean_motion <= 0.0)
return false;
sma = tle_sma_km(tle);
/* Filter 1: perigee altitude */
perigee_alt = sma * (1.0 - tle->eccentricity) - WGS72_AE;
max_alt = max_visible_altitude_km(win->min_el_deg);
if (perigee_alt > max_alt)
return false;
/* Filter 2: inclination + footprint vs observer latitude */
obs_lat_abs = fabs(win->obs.lat);
footprint_rad = ground_footprint_deg(sma, win->min_el_deg) * DEG_TO_RAD;
if (tle->inclination + footprint_rad < obs_lat_abs)
return false;
/* Filter 3: RAAN alignment via J2 secular precession */
dt_days = win->jd_mid - tle->epoch;
raan_projected = tle->raan
+ j2_raan_rate(sma, tle->inclination) * dt_days;
raan_projected = fmod(raan_projected, 2.0 * M_PI);
if (raan_projected < 0.0)
raan_projected += 2.0 * M_PI;
/* Observer LST at query midpoint */
lst = gmst_from_jd(win->jd_mid) + win->obs.lon;
lst = fmod(lst, 2.0 * M_PI);
if (lst < 0.0)
lst += 2.0 * M_PI;
/* RAAN window: Earth rotation during query + footprint pad */
earth_rot_rad = (win->jd_end - win->jd_start) * EARTH_ROT_RAD_PER_DAY;
raan_window_half = earth_rot_rad / 2.0
+ ground_footprint_deg(sma, win->min_el_deg) * DEG_TO_RAD;
if (raan_window_half >= M_PI)
return true; /* full rotation -- pass everything */
raan_diff = fabs(raan_projected - lst);
if (raan_diff > M_PI)
raan_diff = 2.0 * M_PI - raan_diff;
return (raan_diff <= raan_window_half);
}
/*
* spgist_tle_leaf_consistent -- final check on a leaf TLE
*
* Delegates to tle_passes_visibility_filter() for the &? operator.
* recheck = false: the &? operator IS the superset filter. * recheck = false: the &? operator IS the superset filter.
* The user runs predict_passes() on survivors for SGP4 ground truth. * The user runs predict_passes() on survivors for SGP4 ground truth.
*/ */
@ -688,76 +766,11 @@ spgist_tle_leaf_consistent(PG_FUNCTION_ARGS)
{ {
ObserverWindow win; ObserverWindow win;
HeapTupleHeader composite; HeapTupleHeader composite;
double sma, perigee_alt, max_alt;
double obs_lat_abs, footprint_rad;
double dt_days, raan_projected, lst;
double earth_rot_rad, raan_window_half, raan_diff;
composite = DatumGetHeapTupleHeader(in->scankeys[i].sk_argument); composite = DatumGetHeapTupleHeader(in->scankeys[i].sk_argument);
extract_observer_window(composite, &win); extract_observer_window(composite, &win);
sma = tle_sma_km(tle); if (!tle_passes_visibility_filter(tle, &win))
/* Filter 1: perigee altitude check */
perigee_alt = sma * (1.0 - tle->eccentricity) - WGS72_AE;
max_alt = max_visible_altitude_km(win.min_el_deg);
if (perigee_alt > max_alt)
{
result = false;
break;
}
/* Filter 2: inclination + footprint vs observer latitude */
obs_lat_abs = fabs(win.obs.lat);
footprint_rad = ground_footprint_deg(sma, win.min_el_deg)
* DEG_TO_RAD;
if (tle->inclination + footprint_rad < obs_lat_abs)
{
result = false;
break;
}
/*
* Filter 3: RAAN query-time filter.
*
* Project RAAN to query midpoint via J2 precession rate.
* Check alignment with observer's local sidereal time.
* Earth rotates 360 deg/day, so the RAAN window widens
* with query duration.
*/
dt_days = win.jd_mid - tle->epoch;
raan_projected = tle->raan
+ j2_raan_rate(sma, tle->inclination) * dt_days;
raan_projected = fmod(raan_projected, 2.0 * M_PI);
if (raan_projected < 0.0)
raan_projected += 2.0 * M_PI;
/* Observer LST at query midpoint */
lst = gmst_from_jd(win.jd_mid) + win.obs.lon;
lst = fmod(lst, 2.0 * M_PI);
if (lst < 0.0)
lst += 2.0 * M_PI;
/*
* RAAN window: Earth rotation during query + footprint.
* Full rotation (24h+) means pass everything.
*/
earth_rot_rad = (win.jd_end - win.jd_start) * EARTH_ROT_RAD_PER_DAY;
raan_window_half = earth_rot_rad / 2.0
+ ground_footprint_deg(sma, win.min_el_deg) * DEG_TO_RAD;
if (raan_window_half >= M_PI)
{
/* Full rotation or close: no RAAN filter for this key */
continue;
}
raan_diff = fabs(raan_projected - lst);
if (raan_diff > M_PI)
raan_diff = 2.0 * M_PI - raan_diff;
if (raan_diff > raan_window_half)
{ {
result = false; result = false;
break; break;
@ -791,51 +804,8 @@ tle_visibility_possible(PG_FUNCTION_ARGS)
HeapTupleHeader composite = PG_GETARG_HEAPTUPLEHEADER(0); HeapTupleHeader composite = PG_GETARG_HEAPTUPLEHEADER(0);
pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(1); pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(1);
ObserverWindow win; ObserverWindow win;
double sma, perigee_alt, max_alt;
double obs_lat_abs, footprint_rad;
double dt_days, raan_projected, lst;
double earth_rot_rad, raan_window_half, raan_diff;
extract_observer_window(composite, &win); extract_observer_window(composite, &win);
sma = tle_sma_km(tle); PG_RETURN_BOOL(tle_passes_visibility_filter(tle, &win));
/* Filter 1: perigee altitude */
perigee_alt = sma * (1.0 - tle->eccentricity) - WGS72_AE;
max_alt = max_visible_altitude_km(win.min_el_deg);
if (perigee_alt > max_alt)
PG_RETURN_BOOL(false);
/* Filter 2: inclination + footprint */
obs_lat_abs = fabs(win.obs.lat);
footprint_rad = ground_footprint_deg(sma, win.min_el_deg) * DEG_TO_RAD;
if (tle->inclination + footprint_rad < obs_lat_abs)
PG_RETURN_BOOL(false);
/* Filter 3: RAAN alignment */
dt_days = win.jd_mid - tle->epoch;
raan_projected = tle->raan
+ j2_raan_rate(sma, tle->inclination) * dt_days;
raan_projected = fmod(raan_projected, 2.0 * M_PI);
if (raan_projected < 0.0)
raan_projected += 2.0 * M_PI;
lst = gmst_from_jd(win.jd_mid) + win.obs.lon;
lst = fmod(lst, 2.0 * M_PI);
if (lst < 0.0)
lst += 2.0 * M_PI;
earth_rot_rad = (win.jd_end - win.jd_start) * EARTH_ROT_RAD_PER_DAY;
raan_window_half = earth_rot_rad / 2.0
+ ground_footprint_deg(sma, win.min_el_deg) * DEG_TO_RAD;
if (raan_window_half >= M_PI)
PG_RETURN_BOOL(true); /* full rotation -- pass everything */
raan_diff = fabs(raan_projected - lst);
if (raan_diff > M_PI)
raan_diff = 2.0 * M_PI - raan_diff;
PG_RETURN_BOOL(raan_diff <= raan_window_half);
} }

97
test/expected/spgist_tle.out
View File

@ -207,6 +207,103 @@ ORDER BY name;
SSO-800 SSO-800
(2 rows) (2 rows)
-- ============================================================
-- Test 9: Degenerate TLE (mean_motion = 0) — rejected by filter
-- ============================================================
INSERT INTO test_spgist (name, tle) VALUES ('DECAYED',
'1 99904U 24999D 24001.50000000 .00000000 00000+0 00000+0 0 9993
2 99904 0.0000 0.0000 0000000 0.0000 0.0000 0.00000000 00001');
SELECT name,
ROW(
observer('0.0N 0.0E 0m'),
'2024-01-01 00:00:00+00'::timestamptz,
'2024-01-02 00:00:00+00'::timestamptz,
10.0
)::observer_window &? tle AS visible
FROM test_spgist
WHERE name = 'DECAYED';
name | visible
---------+---------
DECAYED | f
(1 row)
-- ============================================================
-- Test 10: Polar observer (90N) — only ISS and SSO-800 reach
-- the pole. ISS (51.6 + footprint) < 90, so only SSO-800
-- (retrograde, 98.7 deg inc > 90 deg) passes. 24h window.
-- ============================================================
SELECT name
FROM test_spgist
WHERE ROW(
observer('90.0N 0.0E 0m'),
'2024-01-01 00:00:00+00'::timestamptz,
'2024-01-02 00:00:00+00'::timestamptz,
10.0
)::observer_window &? tle
ORDER BY name;
name
---------
SSO-800
(1 row)
-- ============================================================
-- Test 11: Zero-duration window — sees only what is directly
-- overhead at the instant. RAAN window = footprint only.
-- ============================================================
SELECT name,
ROW(
observer('0.0N 0.0E 0m'),
'2024-01-01 00:00:00+00'::timestamptz,
'2024-01-01 00:00:00+00'::timestamptz,
10.0
)::observer_window &? tle AS visible
FROM test_spgist
WHERE name = 'ISS';
name | visible
------+---------
ISS | f
(1 row)
-- ============================================================
-- Test 12: Index-vs-seqscan consistency on 24h Eagle Idaho
-- (the primary correctness test, now after all inserts)
-- ============================================================
SET enable_seqscan = off;
SELECT name
FROM test_spgist
WHERE ROW(
observer('43.6977N 116.3535W 760m'),
'2024-01-01 00:00:00+00'::timestamptz,
'2024-01-02 00:00:00+00'::timestamptz,
10.0
)::observer_window &? tle
ORDER BY name;
name
---------
ISS
SSO-800
(2 rows)
RESET enable_seqscan;
SET enable_indexscan = off;
SET enable_bitmapscan = off;
SELECT name
FROM test_spgist
WHERE ROW(
observer('43.6977N 116.3535W 760m'),
'2024-01-01 00:00:00+00'::timestamptz,
'2024-01-02 00:00:00+00'::timestamptz,
10.0
)::observer_window &? tle
ORDER BY name;
name
---------
ISS
SSO-800
(2 rows)
RESET enable_indexscan;
RESET enable_bitmapscan;
-- ============================================================ -- ============================================================
-- Cleanup -- Cleanup
-- ============================================================ -- ============================================================

80
test/sql/spgist_tle.sql
View File

@ -191,6 +191,86 @@ WHERE ROW(
ORDER BY name; ORDER BY name;
-- ============================================================
-- Test 9: Degenerate TLE (mean_motion = 0) — rejected by filter
-- ============================================================
INSERT INTO test_spgist (name, tle) VALUES ('DECAYED',
'1 99904U 24999D 24001.50000000 .00000000 00000+0 00000+0 0 9993
2 99904 0.0000 0.0000 0000000 0.0000 0.0000 0.00000000 00001');
SELECT name,
ROW(
observer('0.0N 0.0E 0m'),
'2024-01-01 00:00:00+00'::timestamptz,
'2024-01-02 00:00:00+00'::timestamptz,
10.0
)::observer_window &? tle AS visible
FROM test_spgist
WHERE name = 'DECAYED';
-- ============================================================
-- Test 10: Polar observer (90N) — only ISS and SSO-800 reach
-- the pole. ISS (51.6 + footprint) < 90, so only SSO-800
-- (retrograde, 98.7 deg inc > 90 deg) passes. 24h window.
-- ============================================================
SELECT name
FROM test_spgist
WHERE ROW(
observer('90.0N 0.0E 0m'),
'2024-01-01 00:00:00+00'::timestamptz,
'2024-01-02 00:00:00+00'::timestamptz,
10.0
)::observer_window &? tle
ORDER BY name;
-- ============================================================
-- Test 11: Zero-duration window — sees only what is directly
-- overhead at the instant. RAAN window = footprint only.
-- ============================================================
SELECT name,
ROW(
observer('0.0N 0.0E 0m'),
'2024-01-01 00:00:00+00'::timestamptz,
'2024-01-01 00:00:00+00'::timestamptz,
10.0
)::observer_window &? tle AS visible
FROM test_spgist
WHERE name = 'ISS';
-- ============================================================
-- Test 12: Index-vs-seqscan consistency on 24h Eagle Idaho
-- (the primary correctness test, now after all inserts)
-- ============================================================
SET enable_seqscan = off;
SELECT name
FROM test_spgist
WHERE ROW(
observer('43.6977N 116.3535W 760m'),
'2024-01-01 00:00:00+00'::timestamptz,
'2024-01-02 00:00:00+00'::timestamptz,
10.0
)::observer_window &? tle
ORDER BY name;
RESET enable_seqscan;
SET enable_indexscan = off;
SET enable_bitmapscan = off;
SELECT name
FROM test_spgist
WHERE ROW(
observer('43.6977N 116.3535W 760m'),
'2024-01-01 00:00:00+00'::timestamptz,
'2024-01-02 00:00:00+00'::timestamptz,
10.0
)::observer_window &? tle
ORDER BY name;
RESET enable_indexscan;
RESET enable_bitmapscan;
-- ============================================================ -- ============================================================
-- Cleanup -- Cleanup
-- ============================================================ -- ============================================================