From 2a7240e73903f6cc0de783f5fed8d1c94001c817 Mon Sep 17 00:00:00 2001 From: Ryan Malloy Date: Tue, 17 Feb 2026 20:36:47 -0700 Subject: [PATCH] Harden SP-GiST orbital trie after Hamilton review 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. --- src/spgist_tle.c | 214 +++++++++++++++-------------------- test/expected/spgist_tle.out | 97 ++++++++++++++++ test/sql/spgist_tle.sql | 80 +++++++++++++ 3 files changed, 269 insertions(+), 122 deletions(-) diff --git a/src/spgist_tle.c b/src/spgist_tle.c index 18ef774..0fcfacc 100644 --- a/src/spgist_tle.c +++ b/src/spgist_tle.c @@ -116,8 +116,13 @@ max_visible_altitude_km(double min_el_deg) double el_rad, sin_el; 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) - return 50000.0; /* effectively no altitude filter */ + return 50000.0; el_rad = min_el_deg * DEG_TO_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. */ - 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 + 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. * 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 @@ -465,6 +476,9 @@ spgist_tle_picksplit(PG_FUNCTION_ARGS) nBins = 2; if (nBins > 16) nBins = 16; + /* Prevent over-read: never more bins than tuples */ + if (nBins > nTuples) + nBins = nTuples; perBin = nTuples / nBins; remainder = nTuples % nBins; @@ -571,17 +585,17 @@ spgist_tle_inner_consistent(PG_FUNCTION_ARGS) { /* * L0: SMA pruning. - * bin_low is the SMA of the lowest object in this bin. - * A satellite in this bin has perigee_alt >= bin_low - AE - * (approximately, for near-circular orbits in the bin). - * If the lowest possible altitude exceeds max_visible, skip. + * bin_low is the SMA of the lowest object in this bin at + * index build time. Later inserts with lower SMA route to + * bin 0, so for the first bin we use AE (physical minimum + * SMA) to avoid false negatives. * - * Conservative: use bin_low as a lower bound on SMA, - * compute perigee assuming e=0 (circular, worst case for + * Conservative: assume e=0 (circular, worst case for * "too high" pruning -- circular has highest perigee for * 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); 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) * 2. Inclination + ground footprint vs observer latitude * 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. * The user runs predict_passes() on survivors for SGP4 ground truth. */ @@ -688,76 +766,11 @@ spgist_tle_leaf_consistent(PG_FUNCTION_ARGS) { ObserverWindow win; 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); extract_observer_window(composite, &win); - sma = tle_sma_km(tle); - - /* 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) + if (!tle_passes_visibility_filter(tle, &win)) { result = false; break; @@ -791,51 +804,8 @@ tle_visibility_possible(PG_FUNCTION_ARGS) HeapTupleHeader composite = PG_GETARG_HEAPTUPLEHEADER(0); pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(1); 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); - 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) - 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); + PG_RETURN_BOOL(tle_passes_visibility_filter(tle, &win)); } diff --git a/test/expected/spgist_tle.out b/test/expected/spgist_tle.out index 9e0e714..8e45528 100644 --- a/test/expected/spgist_tle.out +++ b/test/expected/spgist_tle.out @@ -207,6 +207,103 @@ ORDER BY name; SSO-800 (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 -- ============================================================ diff --git a/test/sql/spgist_tle.sql b/test/sql/spgist_tle.sql index bb6f7c9..c3fc4d8 100644 --- a/test/sql/spgist_tle.sql +++ b/test/sql/spgist_tle.sql @@ -191,6 +191,86 @@ WHERE ROW( 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 -- ============================================================