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

Extend GiST index to 2-D: altitude + inclination

Browse Source 
The 1-D altitude-band index only pruned ~25% of the 22k satellite
catalog (eliminates MEO/GEO/HEO but 75% is LEO).  Adding inclination
as a second indexed dimension prunes an additional ~40% of remaining
candidates — objects in equatorial or low-inclination orbits that
geometrically cannot pass over the observer's latitude.

Key changes:
- tle_alt_range (16 bytes) → tle_orbital_key (32 bytes) with
  inc_low/inc_high fields
- All 8 GiST support functions updated for 2-D bounding boxes
- Penalty uses margin (half-perimeter) not area to avoid degeneracy
  when leaf entries have zero-width inclination ranges
- Picksplit selects split dimension by normalized spread
- && operator now checks altitude AND inclination overlap
- <-> operator remains altitude-only (conjunction screening is
  altitude-dominant)
- SQL operator comments updated for 2-D semantics
- Test adds Equatorial-LEO satellite at ISS altitude but 5° inclination
  to validate inclination-based pruning
This commit is contained in:
Ryan Malloy 2026-02-15 18:10:19 -07:00
parent 5552bf3280
commit a792e7e083
4 changed files with 284 additions and 194 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
sql/pg_orbit--0.1.0.sql
View File

@ -444,11 +444,11 @@ COMMENT ON FUNCTION pass_visible(tle, observer, timestamptz, timestamptz) IS
-- GiST operator support functions
-- ============================================================
-- Overlap operator: do altitude bands intersect?
-- Overlap operator: do orbital keys overlap in altitude AND inclination?
CREATE FUNCTION tle_overlap(tle, tle) RETURNS boolean
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
-- Altitude distance operator
-- Altitude distance operator (altitude-only, for KNN ordering)
CREATE FUNCTION tle_alt_distance(tle, tle) RETURNS float8
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
@ -461,7 +461,7 @@ CREATE OPERATOR && (
JOIN = areajoinsel
);
COMMENT ON OPERATOR && (tle, tle) IS 'Altitude band overlap — necessary condition for conjunction';
COMMENT ON OPERATOR && (tle, tle) IS 'Orbital key overlap (altitude band AND inclination range) — necessary condition for conjunction';
CREATE OPERATOR <-> (
LEFTARG = tle,
@ -470,11 +470,11 @@ CREATE OPERATOR <-> (
COMMUTATOR = <->
);
COMMENT ON OPERATOR <-> (tle, tle) IS 'Minimum altitude-band separation in km (0 if overlapping)';
COMMENT ON OPERATOR <-> (tle, tle) IS 'Minimum altitude-band separation in km (0 if overlapping). Altitude-only — does not account for inclination. Use && for 2-D filtering.';
-- ============================================================
-- GiST operator class for altitude-band indexing
-- GiST operator class for 2-D orbital indexing (altitude + inclination)
-- ============================================================
-- GiST internal support functions

333
src/gist_tle.c
View File

@ -1,14 +1,18 @@
/*
* gist_tle.c -- GiST operator class for altitude-band indexing on TLE
* gist_tle.c -- GiST operator class for 2-D orbital indexing on TLE
*
* Every TLE defines an orbit whose perigee and apogee altitudes form a
* 1-D range (the "altitude band"). Two orbits can only be in proximity
* if their altitude bands overlap -- a necessary but not sufficient
* condition for conjunction.
* Every TLE defines an orbit whose perigee/apogee altitudes and
* inclination form a 2-D bounding box in (altitude, inclination) space.
* Two orbits can only be in proximity if their altitude bands overlap
* AND their inclination ranges overlap -- both necessary but not
* sufficient conditions for conjunction.
*
* The GiST index stores [perigee_km, apogee_km] ranges as internal
* keys, enabling fast coarse filtering. The && (overlap) operator is
* always rechecked: real conjunction screening requires propagation.
* The GiST index stores [alt_low, alt_high] x [inc_low, inc_high] keys,
* enabling fast coarse filtering. For leaf entries inc_low == inc_high
* (exact value); internal nodes widen to the bounding box of children.
*
* The && (overlap) operator is always rechecked: real conjunction
* screening requires propagation.
*
* Semi-major axis from Kepler's third law using WGS-72 constants:
* a = (KE / n)^(2/3) [earth radii]
@ -37,30 +41,37 @@ PG_FUNCTION_INFO_V1(gist_tle_picksplit);
PG_FUNCTION_INFO_V1(gist_tle_same);
PG_FUNCTION_INFO_V1(gist_tle_distance);
/* Floating-point comparison tolerance (km) */
#define ALT_EPSILON 1.0e-9
/* Floating-point comparison tolerance (km and radians) */
#define KEY_EPSILON 1.0e-9
/*
* Altitude band extracted from a TLE's mean elements.
* 2-D orbital key extracted from a TLE's mean elements.
* Altitude band (perigee/apogee) plus inclination range.
* This is the GiST internal key -- much cheaper to compare
* than propagating two full state vectors.
*
* For leaf entries: inc_low == inc_high (exact inclination).
* For internal nodes: bounding box over all children.
*/
typedef struct tle_alt_range
typedef struct tle_orbital_key
{
double low; /* perigee altitude, km */
double high; /* apogee altitude, km */
} tle_alt_range;
double alt_low; /* perigee altitude, km */
double alt_high; /* apogee altitude, km */
double inc_low; /* inclination lower bound, radians */
double inc_high; /* inclination upper bound, radians */
} tle_orbital_key;
/* ----------------------------------------------------------------
* tle_to_alt_range -- compute [perigee, apogee] from mean elements
* tle_to_orbital_key -- compute [perigee, apogee] x [inc, inc]
*
* Uses WGS-72 KE and AE (the only constants valid for SGP4 elements).
* Degenerate TLEs with n <= 0 map to a zero-width range at 0 km.
* Degenerate TLEs with n <= 0 map to zero-width ranges at 0.
* Inclination is stored in radians (same as pg_tle).
* ----------------------------------------------------------------
*/
static void
tle_to_alt_range(const pg_tle *tle, tle_alt_range *range)
tle_to_orbital_key(const pg_tle *tle, tle_orbital_key *key)
{
double n = tle->mean_motion; /* rad/min */
double e = tle->eccentricity;
@ -68,74 +79,91 @@ tle_to_alt_range(const pg_tle *tle, tle_alt_range *range)
if (n <= 0.0)
{
range->low = 0.0;
range->high = 0.0;
key->alt_low = 0.0;
key->alt_high = 0.0;
key->inc_low = 0.0;
key->inc_high = 0.0;
return;
}
a_er = pow(WGS72_KE / n, 2.0 / 3.0);
range->low = a_er * (1.0 - e) * WGS72_AE - WGS72_AE;
range->high = a_er * (1.0 + e) * WGS72_AE - WGS72_AE;
key->alt_low = a_er * (1.0 - e) * WGS72_AE - WGS72_AE;
key->alt_high = a_er * (1.0 + e) * WGS72_AE - WGS72_AE;
/* Guard against numerical inversion from near-zero eccentricity */
if (range->low > range->high)
if (key->alt_low > key->alt_high)
{
double tmp = range->low;
range->low = range->high;
range->high = tmp;
double tmp = key->alt_low;
key->alt_low = key->alt_high;
key->alt_high = tmp;
}
/* Leaf entry: exact inclination (radians) */
key->inc_low = tle->inclination;
key->inc_high = tle->inclination;
}
/* ----------------------------------------------------------------
* range_overlaps -- do two altitude bands share any interval?
* key_overlaps -- do two orbital keys overlap in BOTH dimensions?
*
* Altitude bands AND inclination ranges must both overlap.
* ----------------------------------------------------------------
*/
static inline bool
range_overlaps(const tle_alt_range *a, const tle_alt_range *b)
key_overlaps(const tle_orbital_key *a, const tle_orbital_key *b)
{
return (a->low <= b->high) && (b->low <= a->high);
return (a->alt_low <= b->alt_high) && (b->alt_low <= a->alt_high)
&& (a->inc_low <= b->inc_high) && (b->inc_low <= a->inc_high);
}
/* ----------------------------------------------------------------
* range_contains -- does outer fully contain inner?
* key_contains -- does outer fully contain inner in both dimensions?
* ----------------------------------------------------------------
*/
static inline bool
range_contains(const tle_alt_range *outer, const tle_alt_range *inner)
key_contains(const tle_orbital_key *outer, const tle_orbital_key *inner)
{
return (outer->low <= inner->low) && (inner->high <= outer->high);
return (outer->alt_low <= inner->alt_low)
&& (inner->alt_high <= outer->alt_high)
&& (outer->inc_low <= inner->inc_low)
&& (inner->inc_high <= outer->inc_high);
}
/* ----------------------------------------------------------------
* range_merge -- expand dst to encompass src
* key_merge -- expand dst bounding box to encompass src
* ----------------------------------------------------------------
*/
static inline void
range_merge(tle_alt_range *dst, const tle_alt_range *src)
key_merge(tle_orbital_key *dst, const tle_orbital_key *src)
{
if (src->low < dst->low)
dst->low = src->low;
if (src->high > dst->high)
dst->high = src->high;
if (src->alt_low < dst->alt_low)
dst->alt_low = src->alt_low;
if (src->alt_high > dst->alt_high)
dst->alt_high = src->alt_high;
if (src->inc_low < dst->inc_low)
dst->inc_low = src->inc_low;
if (src->inc_high > dst->inc_high)
dst->inc_high = src->inc_high;
}
/* ----------------------------------------------------------------
* range_separation -- minimum gap between two non-overlapping ranges
* alt_separation -- minimum altitude gap between two keys
*
* Returns 0 if the ranges overlap.
* Returns 0 if the altitude bands overlap.
* Used for KNN distance (altitude-dominant ordering).
* ----------------------------------------------------------------
*/
static inline double
range_separation(const tle_alt_range *a, const tle_alt_range *b)
alt_separation(const tle_orbital_key *a, const tle_orbital_key *b)
{
if (a->high < b->low)
return b->low - a->high;
if (b->high < a->low)
return a->low - b->high;
if (a->alt_high < b->alt_low)
return b->alt_low - a->alt_high;
if (b->alt_high < a->alt_low)
return a->alt_low - b->alt_high;
return 0.0;
}
@ -149,8 +177,8 @@ range_separation(const tle_alt_range *a, const tle_alt_range *b)
/*
* tle_overlap(tle, tle) -> bool [the && operator]
*
* True if the altitude bands of two TLEs share any interval.
* This is a fast pre-filter: overlapping bands are necessary
* True if two TLEs overlap in both altitude AND inclination.
* This is a fast pre-filter: overlapping keys are necessary
* but not sufficient for actual conjunction.
*/
Datum
@ -158,12 +186,12 @@ tle_overlap(PG_FUNCTION_ARGS)
{
pg_tle *a = (pg_tle *) PG_GETARG_POINTER(0);
pg_tle *b = (pg_tle *) PG_GETARG_POINTER(1);
tle_alt_range ra, rb;
tle_orbital_key ka, kb;
tle_to_alt_range(a, &ra);
tle_to_alt_range(b, &rb);
tle_to_orbital_key(a, &ka);
tle_to_orbital_key(b, &kb);
PG_RETURN_BOOL(range_overlaps(&ra, &rb));
PG_RETURN_BOOL(key_overlaps(&ka, &kb));
}
@ -174,18 +202,21 @@ tle_overlap(PG_FUNCTION_ARGS)
* overlap. This is not the physical distance between the objects --
* it is the gap between their orbital shells, useful for ordering
* nearest-neighbor queries without propagation.
*
* Altitude-only: inclination weighting adds complexity without
* meaningful benefit for KNN conjunction screening.
*/
Datum
tle_alt_distance(PG_FUNCTION_ARGS)
{
pg_tle *a = (pg_tle *) PG_GETARG_POINTER(0);
pg_tle *b = (pg_tle *) PG_GETARG_POINTER(1);
tle_alt_range ra, rb;
tle_orbital_key ka, kb;
tle_to_alt_range(a, &ra);
tle_to_alt_range(b, &rb);
tle_to_orbital_key(a, &ka);
tle_to_orbital_key(b, &kb);
PG_RETURN_FLOAT8(range_separation(&ra, &rb));
PG_RETURN_FLOAT8(alt_separation(&ka, &kb));
}
@ -196,10 +227,10 @@ tle_alt_distance(PG_FUNCTION_ARGS)
/*
* gist_tle_compress -- extract altitude range from a leaf TLE
* gist_tle_compress -- extract orbital key from a leaf TLE
*
* Leaf entries carry the full pg_tle; we compress to tle_alt_range.
* Internal entries are already tle_alt_range from union operations.
* Leaf entries carry the full pg_tle; we compress to tle_orbital_key.
* Internal entries are already tle_orbital_key from union operations.
*/
Datum
gist_tle_compress(PG_FUNCTION_ARGS)
@ -210,17 +241,17 @@ gist_tle_compress(PG_FUNCTION_ARGS)
if (entry->leafkey)
{
pg_tle *tle = (pg_tle *) DatumGetPointer(entry->key);
tle_alt_range *range = (tle_alt_range *) palloc(sizeof(tle_alt_range));
tle_orbital_key *key = (tle_orbital_key *) palloc(sizeof(tle_orbital_key));
tle_to_alt_range(tle, range);
tle_to_orbital_key(tle, key);
retval = (GISTENTRY *) palloc(sizeof(GISTENTRY));
gistentryinit(*retval, PointerGetDatum(range),
gistentryinit(*retval, PointerGetDatum(key),
entry->rel, entry->page, entry->offset, false);
}
else
{
/* Internal node: already a tle_alt_range */
/* Internal node: already a tle_orbital_key */
retval = entry;
}
@ -241,8 +272,8 @@ gist_tle_decompress(PG_FUNCTION_ARGS)
/*
* gist_tle_consistent -- can this subtree contain matches for the query?
*
* Strategy RTOverlapStrategyNumber (&&): altitude bands must overlap.
* Always sets recheck = true because altitude overlap is only a necessary
* 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.
*/
Datum
@ -253,36 +284,32 @@ gist_tle_consistent(PG_FUNCTION_ARGS)
StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
/* arg 3 is the query type OID, unused */
bool *recheck = (bool *) PG_GETARG_POINTER(4);
tle_alt_range *key = (tle_alt_range *) DatumGetPointer(entry->key);
tle_alt_range query_range;
tle_orbital_key *key = (tle_orbital_key *) DatumGetPointer(entry->key);
tle_orbital_key query_key;
bool result;
tle_to_alt_range(query, &query_range);
tle_to_orbital_key(query, &query_key);
/*
* Altitude overlap is necessary, not sufficient.
* The actual operator (if exact) would need propagation, so always recheck.
*/
*recheck = true;
switch (strategy)
{
case RTOverlapStrategyNumber: /* && */
result = range_overlaps(key, &query_range);
result = key_overlaps(key, &query_key);
break;
case RTContainedByStrategyNumber: /* <@ */
if (GIST_LEAF(entry))
result = range_contains(&query_range, key);
result = key_contains(&query_key, key);
else
result = range_overlaps(key, &query_range);
result = key_overlaps(key, &query_key);
break;
case RTContainsStrategyNumber: /* @> */
if (GIST_LEAF(entry))
result = range_contains(key, &query_range);
result = key_contains(key, &query_key);
else
result = range_overlaps(key, &query_range);
result = key_overlaps(key, &query_key);
break;
default:
@ -297,9 +324,9 @@ gist_tle_consistent(PG_FUNCTION_ARGS)
/*
* gist_tle_union -- compute bounding altitude range for a set of entries
* gist_tle_union -- compute 2-D bounding box for a set of entries
*
* The union of N altitude ranges is simply [min(low), max(high)].
* The union is [min(alt_low), max(alt_high)] x [min(inc_low), max(inc_high)].
*/
Datum
gist_tle_union(PG_FUNCTION_ARGS)
@ -307,21 +334,20 @@ gist_tle_union(PG_FUNCTION_ARGS)
GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0);
int *sizep = (int *) PG_GETARG_POINTER(1);
int i;
tle_alt_range *result;
tle_alt_range *cur;
tle_orbital_key *result;
tle_orbital_key *cur;
result = (tle_alt_range *) palloc(sizeof(tle_alt_range));
cur = (tle_alt_range *) DatumGetPointer(entryvec->vector[0].key);
result->low = cur->low;
result->high = cur->high;
result = (tle_orbital_key *) palloc(sizeof(tle_orbital_key));
cur = (tle_orbital_key *) DatumGetPointer(entryvec->vector[0].key);
*result = *cur;
for (i = 1; i < entryvec->n; i++)
{
cur = (tle_alt_range *) DatumGetPointer(entryvec->vector[i].key);
range_merge(result, cur);
cur = (tle_orbital_key *) DatumGetPointer(entryvec->vector[i].key);
key_merge(result, cur);
}
*sizep = sizeof(tle_alt_range);
*sizep = sizeof(tle_orbital_key);
PG_RETURN_POINTER(result);
}
@ -330,10 +356,10 @@ gist_tle_union(PG_FUNCTION_ARGS)
/*
* gist_tle_penalty -- cost of inserting a new entry into an existing subtree
*
* Penalty is the increase in altitude span (km) that results from
* expanding the subtree's bounding range to include the new entry.
* A good penalty function keeps the tree balanced and minimizes
* unnecessary range expansion.
* Uses margin (half-perimeter) instead of area. Leaf entries have
* inc_low == inc_high, giving zero area -- area-based penalty would
* degenerate to 0 for every insertion, making the tree random.
* Margin remains non-zero for degenerate (zero-width) bounding boxes.
*/
Datum
gist_tle_penalty(PG_FUNCTION_ARGS)
@ -341,15 +367,17 @@ gist_tle_penalty(PG_FUNCTION_ARGS)
GISTENTRY *origentry = (GISTENTRY *) PG_GETARG_POINTER(0);
GISTENTRY *newentry = (GISTENTRY *) PG_GETARG_POINTER(1);
float *penalty = (float *) PG_GETARG_POINTER(2);
tle_alt_range *orig = (tle_alt_range *) DatumGetPointer(origentry->key);
tle_alt_range *add = (tle_alt_range *) DatumGetPointer(newentry->key);
double orig_span;
double merged_span;
tle_orbital_key *orig = (tle_orbital_key *) DatumGetPointer(origentry->key);
tle_orbital_key *add = (tle_orbital_key *) DatumGetPointer(newentry->key);
double orig_margin;
double merged_margin;
orig_span = orig->high - orig->low;
merged_span = fmax(orig->high, add->high) - fmin(orig->low, add->low);
orig_margin = (orig->alt_high - orig->alt_low)
+ (orig->inc_high - orig->inc_low);
merged_margin = (fmax(orig->alt_high, add->alt_high) - fmin(orig->alt_low, add->alt_low))
+ (fmax(orig->inc_high, add->inc_high) - fmin(orig->inc_low, add->inc_low));
*penalty = (float)(merged_span - orig_span);
*penalty = (float)(merged_margin - orig_margin);
PG_RETURN_POINTER(penalty);
}
@ -357,19 +385,19 @@ gist_tle_penalty(PG_FUNCTION_ARGS)
/*
* Comparison callback for qsort in picksplit.
* Sorts entries by the midpoint of their altitude range.
* Sorts entries by a sort key chosen at call time.
*/
typedef struct
{
int index; /* position in the original entry vector */
double midpoint; /* (low + high) / 2 */
double sortval; /* midpoint in the chosen dimension */
} picksplit_item;
static int
picksplit_cmp(const void *a, const void *b)
{
double ma = ((const picksplit_item *) a)->midpoint;
double mb = ((const picksplit_item *) b)->midpoint;
double ma = ((const picksplit_item *) a)->sortval;
double mb = ((const picksplit_item *) b)->sortval;
if (ma < mb)
return -1;
@ -382,9 +410,9 @@ picksplit_cmp(const void *a, const void *b)
/*
* gist_tle_picksplit -- split an overfull page into two groups
*
* Strategy: sort entries by altitude midpoint, split at the median.
* This keeps nearby altitude bands in the same subtree, which is
* optimal for a 1-D range index.
* Standard R-tree approach: compute spread in both dimensions, split
* along whichever dimension has the greater spread. This prevents
* the tree from becoming biased toward one dimension.
*/
Datum
gist_tle_picksplit(PG_FUNCTION_ARGS)
@ -393,18 +421,55 @@ gist_tle_picksplit(PG_FUNCTION_ARGS)
GIST_SPLITVEC *splitvec = (GIST_SPLITVEC *) PG_GETARG_POINTER(1);
int nentries = entryvec->n;
picksplit_item *items;
tle_alt_range *left_union;
tle_alt_range *right_union;
tle_alt_range *cur;
tle_orbital_key *left_union;
tle_orbital_key *right_union;
tle_orbital_key *cur;
int split_at;
int i;
double alt_min, alt_max, inc_min, inc_max;
double alt_spread, inc_spread;
bool split_on_alt;
/* First pass: compute spread in both dimensions */
cur = (tle_orbital_key *) DatumGetPointer(entryvec->vector[0].key);
alt_min = (cur->alt_low + cur->alt_high) / 2.0;
alt_max = alt_min;
inc_min = (cur->inc_low + cur->inc_high) / 2.0;
inc_max = inc_min;
for (i = 1; i < nentries; i++)
{
double alt_mid, inc_mid;
cur = (tle_orbital_key *) DatumGetPointer(entryvec->vector[i].key);
alt_mid = (cur->alt_low + cur->alt_high) / 2.0;
inc_mid = (cur->inc_low + cur->inc_high) / 2.0;
if (alt_mid < alt_min) alt_min = alt_mid;
if (alt_mid > alt_max) alt_max = alt_mid;
if (inc_mid < inc_min) inc_min = inc_mid;
if (inc_mid > inc_max) inc_max = inc_mid;
}
/*
* Normalize spreads so they're comparable. Altitude is km above
* surface (range ~0-35786), inclination in radians (range 0-pi).
* Divide each by its domain width to get a 0-1 fraction.
*/
alt_spread = (alt_max - alt_min) / 35786.0; /* GEO altitude above surface */
inc_spread = (inc_max - inc_min) / M_PI;
split_on_alt = (alt_spread >= inc_spread);
/* Second pass: compute sort values in the chosen dimension */
items = (picksplit_item *) palloc(sizeof(picksplit_item) * nentries);
for (i = 0; i < nentries; i++)
{
cur = (tle_alt_range *) DatumGetPointer(entryvec->vector[i].key);
cur = (tle_orbital_key *) DatumGetPointer(entryvec->vector[i].key);
items[i].index = i;
items[i].midpoint = (cur->low + cur->high) / 2.0;
if (split_on_alt)
items[i].sortval = (cur->alt_low + cur->alt_high) / 2.0;
else
items[i].sortval = (cur->inc_low + cur->inc_high) / 2.0;
}
qsort(items, nentries, sizeof(picksplit_item), picksplit_cmp);
@ -417,39 +482,37 @@ gist_tle_picksplit(PG_FUNCTION_ARGS)
splitvec->spl_nleft = 0;
splitvec->spl_nright = 0;
/* Compute union ranges and assign entries */
left_union = (tle_alt_range *) palloc(sizeof(tle_alt_range));
right_union = (tle_alt_range *) palloc(sizeof(tle_alt_range));
/* 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));
/* Seed the unions from the first entry in each half */
cur = (tle_alt_range *) DatumGetPointer(
cur = (tle_orbital_key *) DatumGetPointer(
entryvec->vector[items[0].index].key);
left_union->low = cur->low;
left_union->high = cur->high;
*left_union = *cur;
cur = (tle_alt_range *) DatumGetPointer(
cur = (tle_orbital_key *) DatumGetPointer(
entryvec->vector[items[split_at].index].key);
right_union->low = cur->low;
right_union->high = cur->high;
*right_union = *cur;
for (i = 0; i < nentries; i++)
{
int idx = items[i].index;
cur = (tle_alt_range *) DatumGetPointer(
cur = (tle_orbital_key *) DatumGetPointer(
entryvec->vector[idx].key);
if (i < split_at)
{
splitvec->spl_left[splitvec->spl_nleft++] =
(OffsetNumber)(idx + 1); /* 1-based */
range_merge(left_union, cur);
key_merge(left_union, cur);
}
else
{
splitvec->spl_right[splitvec->spl_nright++] =
(OffsetNumber)(idx + 1);
range_merge(right_union, cur);
key_merge(right_union, cur);
}
}
@ -465,19 +528,20 @@ gist_tle_picksplit(PG_FUNCTION_ARGS)
/*
* gist_tle_same -- equality test on compressed keys
*
* Two altitude ranges are "same" if both endpoints match within
* a small tolerance. This lets the GiST machinery detect
* duplicate keys.
* Two orbital keys are "same" if all four bounds match within
* tolerance. This lets the GiST machinery detect duplicate keys.
*/
Datum
gist_tle_same(PG_FUNCTION_ARGS)
{
tle_alt_range *a = (tle_alt_range *) PG_GETARG_POINTER(0);
tle_alt_range *b = (tle_alt_range *) PG_GETARG_POINTER(1);
tle_orbital_key *a = (tle_orbital_key *) PG_GETARG_POINTER(0);
tle_orbital_key *b = (tle_orbital_key *) PG_GETARG_POINTER(1);
bool *result = (bool *) PG_GETARG_POINTER(2);
*result = (fabs(a->low - b->low) < ALT_EPSILON &&
fabs(a->high - b->high) < ALT_EPSILON);
*result = (fabs(a->alt_low - b->alt_low) < KEY_EPSILON
&& fabs(a->alt_high - b->alt_high) < KEY_EPSILON
&& fabs(a->inc_low - b->inc_low) < KEY_EPSILON
&& fabs(a->inc_high - b->inc_high) < KEY_EPSILON);
PG_RETURN_POINTER(result);
}
@ -489,6 +553,9 @@ gist_tle_same(PG_FUNCTION_ARGS)
* Returns the minimum altitude-band separation in km.
* For overlapping ranges this is 0, making the entry a candidate.
* The planner uses this to drive ORDER BY <-> queries.
*
* Altitude-only: conjunction screening is altitude-dominant.
* Inclination weighting can be added later if needed.
*/
Datum
gist_tle_distance(PG_FUNCTION_ARGS)
@ -496,10 +563,10 @@ gist_tle_distance(PG_FUNCTION_ARGS)
GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
pg_tle *query = (pg_tle *) PG_GETARG_POINTER(1);
/* strategy number at arg 2, unused for single-distance class */
tle_alt_range *key = (tle_alt_range *) DatumGetPointer(entry->key);
tle_alt_range query_range;
tle_orbital_key *key = (tle_orbital_key *) DatumGetPointer(entry->key);
tle_orbital_key query_key;
tle_to_alt_range(query, &query_range);
tle_to_orbital_key(query, &query_key);
PG_RETURN_FLOAT8(range_separation(key, &query_range));
PG_RETURN_FLOAT8(alt_separation(key, &query_key));
}

45
test/expected/gist_index.out
View File

@ -1,4 +1,4 @@
-- Test GiST index and operators
-- Test GiST index and operators (2-D: altitude + inclination)
CREATE EXTENSION IF NOT EXISTS pg_orbit;
NOTICE: extension "pg_orbit" already exists, skipping
-- Create test table with mixed orbit types
@ -7,46 +7,59 @@ CREATE TABLE test_orbits (
name text,
tle tle
);
-- ISS (LEO, ~400km)
-- ISS (LEO, ~400km, inclination 51.64 deg)
INSERT INTO test_orbits (name, tle) VALUES ('ISS',
'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025
2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001');
-- Hubble (LEO, ~540km)
-- Hubble (LEO, ~540km, inclination 28.47 deg)
INSERT INTO test_orbits (name, tle) VALUES ('Hubble',
'1 20580U 90037B 24001.50000000 .00000790 00000+0 39573-4 0 9992
2 20580 28.4705 61.4398 0002797 317.3115 42.7577 15.09395228 00008');
-- GPS IIR-M (MEO, ~20200km)
-- GPS IIR-M (MEO, ~20200km, inclination 55.44 deg)
INSERT INTO test_orbits (name, tle) VALUES ('GPS-IIR',
'1 28874U 05038A 24001.50000000 .00000012 00000+0 00000+0 0 9993
2 28874 55.4408 300.3467 0117034 51.6543 309.5420 2.00557079 00006');
-- Equatorial-LEO: same altitude as ISS but near-equatorial inclination (5 deg).
-- Same mean motion (15.501 rev/day) and eccentricity as ISS → same altitude band.
-- Different inclination → should NOT overlap with ISS under 2-D key.
INSERT INTO test_orbits (name, tle) VALUES ('Equatorial-LEO',
'1 99901U 24999A 24001.50000000 .00016717 00000-0 10270-3 0 9990
2 99901 5.0000 208.9163 0006703 30.1694 61.7520 15.50100486 00001');
-- Create GiST index
CREATE INDEX test_orbits_gist ON test_orbits USING gist (tle);
-- Overlap: ISS && Hubble = false (different altitude bands: ~400km vs ~540km)
-- 2-D overlap: ISS && Equatorial-LEO should be false (same altitude, different inclination)
SELECT a.name AS sat_a, b.name AS sat_b, a.tle && b.tle AS overlaps
FROM test_orbits a, test_orbits b
WHERE a.id < b.id
ORDER BY a.name, b.name;
sat_a | sat_b | overlaps
--------+---------+----------
---------+----------------+----------
GPS-IIR | Equatorial-LEO | f
Hubble | Equatorial-LEO | f
Hubble | GPS-IIR | f
ISS | Equatorial-LEO | f
ISS | GPS-IIR | f
ISS | Hubble | f
(3 rows)
(6 rows)
-- Altitude distance between different orbit regimes
-- Altitude distance: ISS <-> Equatorial-LEO should be ~0 (same altitude shell)
SELECT a.name AS sat_a, b.name AS sat_b,
round((a.tle <-> b.tle)::numeric, 0) AS alt_dist_km
FROM test_orbits a, test_orbits b
WHERE a.id < b.id
ORDER BY a.name, b.name;
sat_a | sat_b | alt_dist_km
--------+---------+-------------
---------+----------------+-------------
GPS-IIR | Equatorial-LEO | 19451
Hubble | Equatorial-LEO | 115
Hubble | GPS-IIR | 19332
ISS | Equatorial-LEO | 0
ISS | GPS-IIR | 19451
ISS | Hubble | 115
(3 rows)
(6 rows)
-- GiST index scan: find all LEO sats (overlap with ISS)
-- GiST index scan: find all sats overlapping ISS (altitude AND inclination)
-- Should return only ISS itself (Equatorial-LEO is at same altitude but wrong inclination)
SET enable_seqscan = off;
SELECT name FROM test_orbits
WHERE tle && (SELECT tle FROM test_orbits WHERE name = 'ISS')
@ -64,20 +77,22 @@ FROM test_orbits
WHERE name != 'ISS'
ORDER BY tle <-> (SELECT tle FROM test_orbits WHERE name = 'ISS');
name | dist
---------+-------
----------------+-------
Equatorial-LEO | 0
Hubble | 115
GPS-IIR | 19451
(2 rows)
(3 rows)
RESET enable_seqscan;
-- Self-overlap is always true
SELECT name, tle && tle AS self_overlap FROM test_orbits ORDER BY name;
name | self_overlap
---------+--------------
----------------+--------------
Equatorial-LEO | t
GPS-IIR | t
Hubble | t
ISS | t
(3 rows)
(4 rows)
-- Cleanup
DROP TABLE test_orbits;

22
test/sql/gist_index.sql
View File

@ -1,4 +1,4 @@
-- Test GiST index and operators
-- Test GiST index and operators (2-D: altitude + inclination)
CREATE EXTENSION IF NOT EXISTS pg_orbit;
-- Create test table with mixed orbit types
@ -8,38 +8,46 @@ CREATE TABLE test_orbits (
tle tle
);
-- ISS (LEO, ~400km)
-- ISS (LEO, ~400km, inclination 51.64 deg)
INSERT INTO test_orbits (name, tle) VALUES ('ISS',
'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025
2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001');
-- Hubble (LEO, ~540km)
-- Hubble (LEO, ~540km, inclination 28.47 deg)
INSERT INTO test_orbits (name, tle) VALUES ('Hubble',
'1 20580U 90037B 24001.50000000 .00000790 00000+0 39573-4 0 9992
2 20580 28.4705 61.4398 0002797 317.3115 42.7577 15.09395228 00008');
-- GPS IIR-M (MEO, ~20200km)
-- GPS IIR-M (MEO, ~20200km, inclination 55.44 deg)
INSERT INTO test_orbits (name, tle) VALUES ('GPS-IIR',
'1 28874U 05038A 24001.50000000 .00000012 00000+0 00000+0 0 9993
2 28874 55.4408 300.3467 0117034 51.6543 309.5420 2.00557079 00006');
-- Equatorial-LEO: same altitude as ISS but near-equatorial inclination (5 deg).
-- Same mean motion (15.501 rev/day) and eccentricity as ISS → same altitude band.
-- Different inclination → should NOT overlap with ISS under 2-D key.
INSERT INTO test_orbits (name, tle) VALUES ('Equatorial-LEO',
'1 99901U 24999A 24001.50000000 .00016717 00000-0 10270-3 0 9990
2 99901 5.0000 208.9163 0006703 30.1694 61.7520 15.50100486 00001');
-- Create GiST index
CREATE INDEX test_orbits_gist ON test_orbits USING gist (tle);
-- Overlap: ISS && Hubble = false (different altitude bands: ~400km vs ~540km)
-- 2-D overlap: ISS && Equatorial-LEO should be false (same altitude, different inclination)
SELECT a.name AS sat_a, b.name AS sat_b, a.tle && b.tle AS overlaps
FROM test_orbits a, test_orbits b
WHERE a.id < b.id
ORDER BY a.name, b.name;
-- Altitude distance between different orbit regimes
-- Altitude distance: ISS <-> Equatorial-LEO should be ~0 (same altitude shell)
SELECT a.name AS sat_a, b.name AS sat_b,
round((a.tle <-> b.tle)::numeric, 0) AS alt_dist_km
FROM test_orbits a, test_orbits b
WHERE a.id < b.id
ORDER BY a.name, b.name;
-- GiST index scan: find all LEO sats (overlap with ISS)
-- GiST index scan: find all sats overlapping ISS (altitude AND inclination)
-- Should return only ISS itself (Equatorial-LEO is at same altitude but wrong inclination)
SET enable_seqscan = off;
SELECT name FROM test_orbits
WHERE tle && (SELECT tle FROM test_orbits WHERE name = 'ISS')