Add v0.12.0: equatorial GiST operator class + DE moon equatorial functions

Feature A: GiST index for equatorial type with KNN ordering (<-> strategy 15). 24-byte float-precision spherical bounding box, RA-wrapping aware merge/split, Vincenty lower-bound distance for correct KNN pruning. Apollo-hardened with epsilon-widened bounds, circular-aware picksplit, compile-time size assertions. Feature B: 4 new DE moon equatorial functions (galilean_equatorial_de, saturn_moon_equatorial_de, uranus_moon_equatorial_de, mars_moon_equatorial_de). Same-provider rule enforced, transparent VSOP87 fallback. 120 -> 132 SQL objects. 22 regression suites passing.
2026-02-24 13:15:34 -07:00 · 2026-02-24 13:15:34 -07:00 · 84ce1f1b8d
commit 84ce1f1b8d
parent 608370c746
10 changed files with 3074 additions and 4 deletions
--- a/9
+++ b/9
@ -9,7 +9,8 @@ DATA = sql/pg_orrery--0.1.0.sql sql/pg_orrery--0.2.0.sql sql/pg_orrery--0.1.0--0
       sql/pg_orrery--0.8.0.sql sql/pg_orrery--0.7.0--0.8.0.sql \
       sql/pg_orrery--0.9.0.sql sql/pg_orrery--0.8.0--0.9.0.sql \
       sql/pg_orrery--0.10.0.sql sql/pg_orrery--0.9.0--0.10.0.sql \
-       sql/pg_orrery--0.11.0.sql sql/pg_orrery--0.10.0--0.11.0.sql
+       sql/pg_orrery--0.11.0.sql sql/pg_orrery--0.10.0--0.11.0.sql \
+       sql/pg_orrery--0.12.0.sql sql/pg_orrery--0.11.0--0.12.0.sql

 # Our extension C sources
 OBJS = src/pg_orrery.o src/tle_type.o src/eci_type.o src/observer_type.o \
@ -25,7 +26,8 @@ OBJS = src/pg_orrery.o src/tle_type.o src/eci_type.o src/observer_type.o \
       src/spgist_tle.o \
       src/orbital_elements_type.o \
       src/equatorial_funcs.o \
-       src/refraction_funcs.o
+       src/refraction_funcs.o \
+       src/gist_equatorial.o

 # Vendored SGP4/SDP4 sources (pure C, from Bill Gray's sat_code, MIT license)
 SGP4_DIR = src/sgp4
@ -41,7 +43,8 @@ OBJS += $(SGP4_OBJS)
 REGRESS = tle_parse sgp4_propagate coord_transforms pass_prediction gist_index convenience \
         star_observe kepler_comet planet_observe moon_observe lambert_transfer \
         de_ephemeris od_fit spgist_tle orbital_elements equatorial refraction \
-         aberration v011_features vallado_518
+         aberration v011_features vallado_518 \
+         gist_equatorial v012_features
 REGRESS_OPTS = --inputdir=test

 # Pure C — no C++ runtime needed. LAPACK for OD solver (dgelss_).
--- a/pg_orrery.control
+++ b/pg_orrery.control
@ -1,4 +1,4 @@
 comment = 'A database orrery — celestial mechanics types and functions for PostgreSQL'
-default_version = '0.11.0'
+default_version = '0.12.0'
 module_pathname = '$libdir/pg_orrery'
 relocatable = true
--- a/sql/pg_orrery--0.11.0--0.12.0.sql
+++ b/sql/pg_orrery--0.11.0--0.12.0.sql
@ -0,0 +1,72 @@
+-- pg_orrery 0.11.0 -> 0.12.0 migration
+--
+-- Adds equatorial GiST operator class for KNN sky queries
+-- and DE moon equatorial functions for all 4 planetary moon families.
+
+-- ============================================================
+-- GiST support functions for equatorial type
+-- ============================================================
+
+CREATE FUNCTION gist_eq_consistent(internal, equatorial, smallint, oid, internal) RETURNS bool
+    AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION gist_eq_union(internal, internal) RETURNS internal
+    AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION gist_eq_compress(internal) RETURNS internal
+    AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION gist_eq_decompress(internal) RETURNS internal
+    AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION gist_eq_penalty(internal, internal, internal) RETURNS internal
+    AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION gist_eq_picksplit(internal, internal) RETURNS internal
+    AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION gist_eq_same(internal, internal, internal) RETURNS internal
+    AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION gist_eq_distance(internal, equatorial, smallint, oid, internal) RETURNS float8
+    AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+-- ============================================================
+-- Equatorial GiST operator class (KNN ordering only)
+-- ============================================================
+
+CREATE OPERATOR CLASS eq_gist_ops
+    DEFAULT FOR TYPE equatorial USING gist AS
+        OPERATOR  15  <-> (equatorial, equatorial) FOR ORDER BY pg_catalog.float_ops,
+        FUNCTION  1   gist_eq_consistent(internal, equatorial, smallint, oid, internal),
+        FUNCTION  2   gist_eq_union(internal, internal),
+        FUNCTION  3   gist_eq_compress(internal),
+        FUNCTION  4   gist_eq_decompress(internal),
+        FUNCTION  5   gist_eq_penalty(internal, internal, internal),
+        FUNCTION  6   gist_eq_picksplit(internal, internal),
+        FUNCTION  7   gist_eq_same(internal, internal, internal),
+        FUNCTION  8   gist_eq_distance(internal, equatorial, smallint, oid, internal);
+
+-- ============================================================
+-- DE moon equatorial functions (STABLE, fall back to VSOP87)
+-- ============================================================
+
+CREATE FUNCTION galilean_equatorial_de(int4, timestamptz) RETURNS equatorial
+    AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION galilean_equatorial_de(int4, timestamptz) IS
+    'Geocentric RA/Dec of a Galilean moon via DE parent position (falls back to VSOP87). 0=Io, 1=Europa, 2=Ganymede, 3=Callisto.';
+
+CREATE FUNCTION saturn_moon_equatorial_de(int4, timestamptz) RETURNS equatorial
+    AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION saturn_moon_equatorial_de(int4, timestamptz) IS
+    'Geocentric RA/Dec of a Saturn moon via DE parent position (falls back to VSOP87). 0=Mimas..7=Hyperion.';
+
+CREATE FUNCTION uranus_moon_equatorial_de(int4, timestamptz) RETURNS equatorial
+    AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION uranus_moon_equatorial_de(int4, timestamptz) IS
+    'Geocentric RA/Dec of a Uranus moon via DE parent position (falls back to VSOP87). 0=Miranda..4=Oberon.';
+
+CREATE FUNCTION mars_moon_equatorial_de(int4, timestamptz) RETURNS equatorial
+    AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION mars_moon_equatorial_de(int4, timestamptz) IS
+    'Geocentric RA/Dec of a Mars moon via DE parent position (falls back to VSOP87). 0=Phobos, 1=Deimos.';
--- a/sql/pg_orrery--0.12.0.sql
+++ b/sql/pg_orrery--0.12.0.sql
--- a/src/de_funcs.c
+++ b/src/de_funcs.c
@ -56,6 +56,10 @@ PG_FUNCTION_INFO_V1(moon_observe_apparent_de);
 PG_FUNCTION_INFO_V1(planet_equatorial_apparent_de);
 PG_FUNCTION_INFO_V1(moon_equatorial_apparent_de);
 PG_FUNCTION_INFO_V1(small_body_observe_apparent_de);
+PG_FUNCTION_INFO_V1(galilean_equatorial_de);
+PG_FUNCTION_INFO_V1(saturn_moon_equatorial_de);
+PG_FUNCTION_INFO_V1(uranus_moon_equatorial_de);
+PG_FUNCTION_INFO_V1(mars_moon_equatorial_de);
 PG_FUNCTION_INFO_V1(pg_orrery_ephemeris_info);


@ -619,6 +623,147 @@ mars_moon_observe_de(PG_FUNCTION_ARGS)
 }


+/* ================================================================
+ * Planetary moon equatorial coordinates with DE parent positions
+ *
+ * Same DE/VSOP87 dispatch as observe_moon_de(), but returns
+ * geocentric RA/Dec instead of topocentric az/el.
+ * ================================================================
+ */
+static void
+equatorial_moon_de(const double moon_rel[3], int parent_body_id,
+                   double jd, pg_equatorial *result)
+{
+    double parent_xyz[6];
+    double earth_xyz[6];
+    double geo_ecl[3];
+    bool   have_de;
+
+    /* Rule 7: both parent and Earth from same provider */
+    have_de = eph_de_planet(parent_body_id, jd, parent_xyz) &&
+              eph_de_earth(jd, earth_xyz);
+
+    if (!have_de)
+    {
+        int vsop_parent = parent_body_id - 1;
+
+        if (eph_de_available())
+            ereport(NOTICE,
+                    (errmsg("DE ephemeris unavailable, falling back to VSOP87")));
+
+        GetVsop87Coor(jd, vsop_parent, parent_xyz);
+        GetVsop87Coor(jd, 2, earth_xyz);   /* VSOP87 body 2 = Earth */
+    }
+
+    /* Moon geocentric = (parent + moon_relative) - Earth */
+    geo_ecl[0] = (parent_xyz[0] + moon_rel[0]) - earth_xyz[0];
+    geo_ecl[1] = (parent_xyz[1] + moon_rel[1]) - earth_xyz[1];
+    geo_ecl[2] = (parent_xyz[2] + moon_rel[2]) - earth_xyz[2];
+
+    geocentric_to_equatorial(geo_ecl, jd, result);
+}
+
+
+Datum
+galilean_equatorial_de(PG_FUNCTION_ARGS)
+{
+    int32          body_id = PG_GETARG_INT32(0);
+    int64          ts      = PG_GETARG_INT64(1);
+    double         jd;
+    double         moon_xyz[3];
+    pg_equatorial *result;
+
+    if (body_id < L12_IO || body_id > L12_CALLISTO)
+        ereport(ERROR,
+                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+                 errmsg("galilean_equatorial_de: body_id %d must be 0-3 (Io-Callisto)",
+                        body_id)));
+
+    jd = timestamptz_to_jd(ts);
+    GetL12Coor(jd, body_id, moon_xyz, NULL);
+
+    result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
+    equatorial_moon_de(moon_xyz, BODY_JUPITER, jd, result);
+
+    PG_RETURN_POINTER(result);
+}
+
+
+Datum
+saturn_moon_equatorial_de(PG_FUNCTION_ARGS)
+{
+    int32          body_id = PG_GETARG_INT32(0);
+    int64          ts      = PG_GETARG_INT64(1);
+    double         jd;
+    double         moon_xyz[3];
+    pg_equatorial *result;
+
+    if (body_id < TASS17_MIMAS || body_id > TASS17_HYPERION)
+        ereport(ERROR,
+                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+                 errmsg("saturn_moon_equatorial_de: body_id %d must be 0-7 (Mimas-Hyperion)",
+                        body_id)));
+
+    jd = timestamptz_to_jd(ts);
+    GetTass17Coor(jd, body_id, moon_xyz, NULL);
+
+    result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
+    equatorial_moon_de(moon_xyz, BODY_SATURN, jd, result);
+
+    PG_RETURN_POINTER(result);
+}
+
+
+Datum
+uranus_moon_equatorial_de(PG_FUNCTION_ARGS)
+{
+    int32          body_id = PG_GETARG_INT32(0);
+    int64          ts      = PG_GETARG_INT64(1);
+    double         jd;
+    double         moon_xyz[3];
+    pg_equatorial *result;
+
+    if (body_id < GUST86_MIRANDA || body_id > GUST86_OBERON)
+        ereport(ERROR,
+                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+                 errmsg("uranus_moon_equatorial_de: body_id %d must be 0-4 (Miranda-Oberon)",
+                        body_id)));
+
+    jd = timestamptz_to_jd(ts);
+    GetGust86Coor(jd, body_id, moon_xyz, NULL);
+
+    result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
+    equatorial_moon_de(moon_xyz, BODY_URANUS, jd, result);
+
+    PG_RETURN_POINTER(result);
+}
+
+
+Datum
+mars_moon_equatorial_de(PG_FUNCTION_ARGS)
+{
+    int32          body_id = PG_GETARG_INT32(0);
+    int64          ts      = PG_GETARG_INT64(1);
+    double         jd;
+    double         moon_xyz[3];
+    pg_equatorial *result;
+
+    if (body_id < MARS_SAT_PHOBOS || body_id > MARS_SAT_DEIMOS)
+        ereport(ERROR,
+                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+                 errmsg("mars_moon_equatorial_de: body_id %d must be 0-1 (Phobos-Deimos)",
+                        body_id)));
+
+    jd = timestamptz_to_jd(ts);
+    GetMarsSatCoor(jd, body_id, moon_xyz, NULL);
+
+    result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
+    equatorial_moon_de(moon_xyz, BODY_MARS, jd, result);
+
+    PG_RETURN_POINTER(result);
+}
+
+
 /* ================================================================
 * planet_equatorial_de(body_id int, timestamptz) -> equatorial
 *
--- a/src/gist_equatorial.c
+++ b/src/gist_equatorial.c
@ -0,0 +1,782 @@
+/*
+ * gist_equatorial.c -- GiST operator class for equatorial RA/Dec indexing
+ *
+ * Enables KNN nearest-neighbor queries on equatorial coordinates using
+ * the existing <-> (angular distance) operator.  Internal nodes store
+ * spherical bounding boxes in float precision; leaf nodes extract a
+ * point from the pg_equatorial datum.
+ *
+ * Key design: 24-byte float-based spherical box that fits inside
+ * sizeof(pg_equatorial).  GiST's index_form_tuple() copies typlen
+ * bytes from the datum pointer, so the key must be <= 24 bytes.
+ *
+ * RA wrapping: when ra_low > ra_high, the box wraps across 0h,
+ * covering [ra_low, 2*pi) union [0, ra_high].
+ */
+
+#include "postgres.h"
+#include "fmgr.h"
+#include "access/gist.h"
+#include "access/stratnum.h"
+#include "utils/float.h"
+#include "types.h"
+#include <math.h>
+#include <float.h>
+
+PG_FUNCTION_INFO_V1(gist_eq_consistent);
+PG_FUNCTION_INFO_V1(gist_eq_union);
+PG_FUNCTION_INFO_V1(gist_eq_compress);
+PG_FUNCTION_INFO_V1(gist_eq_decompress);
+PG_FUNCTION_INFO_V1(gist_eq_penalty);
+PG_FUNCTION_INFO_V1(gist_eq_picksplit);
+PG_FUNCTION_INFO_V1(gist_eq_same);
+PG_FUNCTION_INFO_V1(gist_eq_distance);
+
+/* Float comparison tolerance */
+#define EQ_KEY_EPSILON  1.0e-7f
+
+#define TWO_PI_F  ((float)(2.0 * M_PI))
+#define PI_F      ((float)M_PI)
+#define HALF_PI_F ((float)(M_PI / 2.0))
+
+/*
+ * Spherical bounding box in float precision.
+ *
+ * RA in radians [0, 2*pi).  When ra_low > ra_high the box wraps
+ * across the 0h meridian: covers [ra_low, 2*pi) union [0, ra_high].
+ * Dec in radians [-pi/2, pi/2].
+ *
+ * 6 floats = 24 bytes == sizeof(pg_equatorial).
+ */
+typedef struct eq_spherical_key
+{
+	float	ra_low;
+	float	ra_high;
+	float	dec_low;
+	float	dec_high;
+	float	_pad[2];	/* zero-fill to 24 bytes */
+} eq_spherical_key;
+
+/* GiST copies typlen bytes via index_form_tuple -- sizes must match exactly */
+StaticAssertDecl(sizeof(eq_spherical_key) == sizeof(pg_equatorial),
+                 "eq_spherical_key must match pg_equatorial size for GiST");
+
+
+/* Does this key wrap across the 0h meridian? */
+static inline bool
+key_wraps_ra(const eq_spherical_key *k)
+{
+	return k->ra_low > k->ra_high;
+}
+
+
+/* Is an RA value (radians) inside the key's RA interval? */
+static inline bool
+ra_in_range(float ra, const eq_spherical_key *k)
+{
+	if (key_wraps_ra(k))
+		return ra >= k->ra_low || ra <= k->ra_high;
+	else
+		return ra >= k->ra_low && ra <= k->ra_high;
+}
+
+
+/* RA interval span in radians (wrapping-aware) */
+static inline float
+ra_span(const eq_spherical_key *k)
+{
+	if (key_wraps_ra(k))
+		return (TWO_PI_F - k->ra_low) + k->ra_high;
+	else
+		return k->ra_high - k->ra_low;
+}
+
+
+/* Dec interval span */
+static inline float
+dec_span(const eq_spherical_key *k)
+{
+	return k->dec_high - k->dec_low;
+}
+
+
+/* Do two RA intervals overlap? (wrapping-aware) */
+static inline bool
+ra_ranges_overlap(const eq_spherical_key *a, const eq_spherical_key *b)
+{
+	/* If either wraps, test whether any endpoint of one is inside the other */
+	if (key_wraps_ra(a) || key_wraps_ra(b))
+	{
+		return ra_in_range(a->ra_low, b)
+			|| ra_in_range(a->ra_high, b)
+			|| ra_in_range(b->ra_low, a)
+			|| ra_in_range(b->ra_high, a);
+	}
+	/* Neither wraps: simple interval overlap */
+	return a->ra_low <= b->ra_high && b->ra_low <= a->ra_high;
+}
+
+
+/* Do two keys overlap in both RA and Dec? */
+static inline bool
+keys_overlap(const eq_spherical_key *a, const eq_spherical_key *b)
+{
+	return ra_ranges_overlap(a, b)
+		&& (a->dec_low <= b->dec_high)
+		&& (b->dec_low <= a->dec_high);
+}
+
+
+/*
+ * Merge src into dst, producing the smallest bounding box that
+ * contains both.  For Dec: simple min/max.  For RA: choose the
+ * smaller of the two possible merged intervals (may or may not wrap).
+ */
+static void
+key_merge_eq(eq_spherical_key *dst, const eq_spherical_key *src)
+{
+	float new_low, new_high;
+	float span_nowrap;
+
+	/* Dec: straightforward */
+	if (src->dec_low < dst->dec_low)
+		dst->dec_low = src->dec_low;
+	if (src->dec_high > dst->dec_high)
+		dst->dec_high = src->dec_high;
+
+	/*
+	 * RA: compute both possible merged intervals and pick the smaller.
+	 *
+	 * Option A (no-wrap): [min(all lows), max(all highs)]
+	 * Option B (wrap):    [max of the "upper" edges, min of the "lower" edges]
+	 *
+	 * We expand both dst and src to their covered RA values, then find
+	 * the tightest enclosing interval.
+	 */
+
+	/* Collect the four RA boundary points */
+	{
+		float pts[4];
+		float min_pt, max_pt;
+		int   i;
+
+		pts[0] = dst->ra_low;
+		pts[1] = dst->ra_high;
+		pts[2] = src->ra_low;
+		pts[3] = src->ra_high;
+
+		/* If either wraps, we can't just take min/max naively */
+		/* Option A: no-wrap interval */
+		min_pt = pts[0];
+		max_pt = pts[0];
+		for (i = 1; i < 4; i++)
+		{
+			if (pts[i] < min_pt) min_pt = pts[i];
+			if (pts[i] > max_pt) max_pt = pts[i];
+		}
+		span_nowrap = max_pt - min_pt;
+
+		/* But we need to check all four points are actually covered */
+		new_low = min_pt;
+		new_high = max_pt;
+	}
+
+	/*
+	 * If the no-wrap span > pi, a wrapping interval might be tighter.
+	 * Otherwise the no-wrap interval is already the smallest enclosing.
+	 */
+	if (span_nowrap > PI_F)
+	{
+		/*
+		 * Wrapping might be tighter.  But we need the actual wrapped
+		 * interval that contains both source intervals.
+		 *
+		 * Strategy: step through each of the 4 possible "gap" positions
+		 * (between consecutive sorted points) and find which gap, when
+		 * placed at the wrap point, gives the smallest enclosing interval.
+		 */
+		eq_spherical_key test;
+		float best_span = TWO_PI_F;
+		float best_low = 0, best_high = TWO_PI_F;
+
+		/* Sort the 4 points */
+		float pts[4];
+		int   i, j;
+
+		pts[0] = dst->ra_low;
+		pts[1] = dst->ra_high;
+		pts[2] = src->ra_low;
+		pts[3] = src->ra_high;
+
+		/* Simple insertion sort for 4 elements */
+		for (i = 1; i < 4; i++)
+		{
+			float tmp = pts[i];
+			j = i - 1;
+			while (j >= 0 && pts[j] > tmp)
+			{
+				pts[j + 1] = pts[j];
+				j--;
+			}
+			pts[j + 1] = tmp;
+		}
+
+		/* Try each gap as the "empty" region */
+		for (i = 0; i < 4; i++)
+		{
+			float gap_start = pts[i];
+			float gap_end   = pts[(i + 1) % 4];
+			float try_low, try_high, try_span;
+
+			if (i < 3)
+			{
+				/* Gap between pts[i] and pts[i+1] */
+				try_low  = gap_end;		/* interval starts after gap */
+				try_high = gap_start;	/* interval ends at gap (wrapping) */
+			}
+			else
+			{
+				/* Gap between pts[3] and pts[0] (wrapping through 0) */
+				try_low  = pts[0];
+				try_high = pts[3];
+			}
+
+			/* Compute span */
+			if (try_low > try_high)
+				try_span = (TWO_PI_F - try_low) + try_high;
+			else
+				try_span = try_high - try_low;
+
+			/* Verify both original intervals are contained */
+			test.ra_low = try_low;
+			test.ra_high = try_high;
+
+			if (ra_in_range(dst->ra_low, &test) &&
+				ra_in_range(dst->ra_high, &test) &&
+				ra_in_range(src->ra_low, &test) &&
+				ra_in_range(src->ra_high, &test) &&
+				try_span < best_span)
+			{
+				best_span = try_span;
+				best_low  = try_low;
+				best_high = try_high;
+			}
+		}
+
+		/* Safety: if no tight interval found, use full circle */
+		if (best_span > TWO_PI_F - EQ_KEY_EPSILON)
+		{
+			dst->ra_low  = 0.0f;
+			dst->ra_high = TWO_PI_F - EQ_KEY_EPSILON;
+		}
+		else
+		{
+			dst->ra_low  = best_low;
+			dst->ra_high = best_high;
+		}
+	}
+	else
+	{
+		/* No-wrap is tight enough */
+		dst->ra_low  = new_low;
+		dst->ra_high = new_high;
+	}
+}
+
+
+/*
+ * Vincenty angular distance between two points on the sphere.
+ * All angles in radians, result in radians.
+ */
+static inline double
+vincenty_rad(double ra1, double dec1, double ra2, double dec2)
+{
+	double d_ra, cos_d_ra, sin_d_ra;
+	double sin_d1, cos_d1, sin_d2, cos_d2;
+	double num1, num2, num, den;
+
+	d_ra     = ra2 - ra1;
+	cos_d_ra = cos(d_ra);
+	sin_d_ra = sin(d_ra);
+
+	sin_d1 = sin(dec1);
+	cos_d1 = cos(dec1);
+	sin_d2 = sin(dec2);
+	cos_d2 = cos(dec2);
+
+	num1 = cos_d2 * sin_d_ra;
+	num2 = cos_d1 * sin_d2 - sin_d1 * cos_d2 * cos_d_ra;
+	num  = sqrt(num1 * num1 + num2 * num2);
+	den  = sin_d1 * sin_d2 + cos_d1 * cos_d2 * cos_d_ra;
+
+	return atan2(num, den);
+}
+
+
+/*
+ * Lower-bound angular distance from a query point to any point inside
+ * a bounding box.  For leaf nodes (point = point) this is exact.
+ * For internal nodes, clamp the query to the nearest box boundary
+ * and compute great-circle distance to the clamped point.
+ *
+ * Box boundaries are widened by EQ_KEY_EPSILON before comparison to
+ * guarantee the lower-bound property under float-to-double promotion.
+ * The key stores float precision but query coordinates are doubles;
+ * without widening, rounding could make the function think the query
+ * is outside the box when it's on the boundary, producing a distance
+ * that overestimates the true minimum.  Overestimates violate the
+ * GiST KNN contract (distance must be a lower bound).
+ */
+static double
+distance_point_to_box(double q_ra, double q_dec,
+                      const eq_spherical_key *box)
+{
+	double clamp_ra, clamp_dec;
+	double ra_lo, ra_hi, dec_lo, dec_hi;
+	bool   ra_inside, dec_inside;
+
+	/* Widen box by float epsilon to guarantee lower bound */
+	ra_lo  = (double)box->ra_low  - (double)EQ_KEY_EPSILON;
+	ra_hi  = (double)box->ra_high + (double)EQ_KEY_EPSILON;
+	dec_lo = (double)box->dec_low - (double)EQ_KEY_EPSILON;
+	dec_hi = (double)box->dec_high + (double)EQ_KEY_EPSILON;
+
+	/* Clamp Dec limits to valid range */
+	if (dec_lo < -M_PI / 2.0) dec_lo = -M_PI / 2.0;
+	if (dec_hi >  M_PI / 2.0) dec_hi =  M_PI / 2.0;
+
+	/* Check RA containment using widened bounds */
+	if (box->ra_low > box->ra_high)
+		ra_inside = (q_ra >= ra_lo || q_ra <= ra_hi);  /* wraps */
+	else
+		ra_inside = (q_ra >= ra_lo && q_ra <= ra_hi);
+
+	dec_inside = (q_dec >= dec_lo && q_dec <= dec_hi);
+
+	/* If query is inside the widened box, distance lower bound is 0 */
+	if (ra_inside && dec_inside)
+		return 0.0;
+
+	/* Clamp Dec to widened box range */
+	clamp_dec = q_dec;
+	if (clamp_dec < dec_lo) clamp_dec = dec_lo;
+	if (clamp_dec > dec_hi) clamp_dec = dec_hi;
+
+	/* Clamp RA to nearest widened box edge */
+	if (ra_inside)
+	{
+		clamp_ra = q_ra;
+	}
+	else
+	{
+		double d_low, d_high;
+
+		d_low = fabs(q_ra - ra_lo);
+		if (d_low > M_PI) d_low = 2.0 * M_PI - d_low;
+
+		d_high = fabs(q_ra - ra_hi);
+		if (d_high > M_PI) d_high = 2.0 * M_PI - d_high;
+
+		clamp_ra = (d_low <= d_high) ? ra_lo : ra_hi;
+	}
+
+	return vincenty_rad(q_ra, q_dec, clamp_ra, clamp_dec);
+}
+
+
+/* Circular midpoint: average of two RA values on the circle */
+static inline float
+ra_midpoint(float a, float b)
+{
+	float mid;
+	if (fabsf(b - a) > PI_F)
+	{
+		/* Wrap: average through 0h */
+		mid = (a + b) / 2.0f + PI_F;
+		if (mid >= TWO_PI_F) mid -= TWO_PI_F;
+	}
+	else
+	{
+		mid = (a + b) / 2.0f;
+	}
+	return mid;
+}
+
+
+/* Circular distance between two RA values */
+static inline float
+ra_circular_dist(float a, float b)
+{
+	float d = fabsf(a - b);
+	return (d > PI_F) ? TWO_PI_F - d : d;
+}
+
+
+/* ================================================================
+ * GiST support function 1: consistent
+ *
+ * For KNN (strategy 15 / <->): always return true, letting the
+ * distance function handle pruning.  GiST uses distance to order
+ * and prune; consistent just needs to not reject valid candidates.
+ * ================================================================
+ */
+Datum
+gist_eq_consistent(PG_FUNCTION_ARGS)
+{
+	StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
+	bool *recheck = (bool *) PG_GETARG_POINTER(4);
+
+	/* Strategy 15 = KNN ordering (<->) -- the only operator in this class */
+	if (strategy != 15)
+		elog(ERROR, "gist_eq_consistent: unsupported strategy %d", strategy);
+
+	/* KNN distance handles all pruning; consistent is permissive */
+	*recheck = true;
+	PG_RETURN_BOOL(true);
+}
+
+
+/* ================================================================
+ * GiST support function 2: union
+ *
+ * Merge all entries into a single bounding box.
+ * ================================================================
+ */
+Datum
+gist_eq_union(PG_FUNCTION_ARGS)
+{
+	GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0);
+	int             *sizep    = (int *) PG_GETARG_POINTER(1);
+	int              i;
+	eq_spherical_key *result;
+	eq_spherical_key *cur;
+
+	/* Allocate sizeof(pg_equatorial) = 24 bytes, matching typlen */
+	result = (eq_spherical_key *) palloc0(sizeof(pg_equatorial));
+	cur = (eq_spherical_key *) DatumGetPointer(entryvec->vector[0].key);
+	*result = *cur;
+
+	for (i = 1; i < entryvec->n; i++)
+	{
+		cur = (eq_spherical_key *) DatumGetPointer(entryvec->vector[i].key);
+		key_merge_eq(result, cur);
+	}
+
+	*sizep = sizeof(pg_equatorial);
+
+	PG_RETURN_POINTER(result);
+}
+
+
+/* ================================================================
+ * GiST support function 3: compress
+ *
+ * Leaf: extract RA/Dec from pg_equatorial into an eq_spherical_key
+ * point (ra_low == ra_high, dec_low == dec_high).
+ * Internal: identity (already an eq_spherical_key from union).
+ * ================================================================
+ */
+Datum
+gist_eq_compress(PG_FUNCTION_ARGS)
+{
+	GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
+	GISTENTRY *retval;
+
+	if (entry->leafkey)
+	{
+		pg_equatorial    *eq  = (pg_equatorial *) DatumGetPointer(entry->key);
+		eq_spherical_key *key = (eq_spherical_key *) palloc0(sizeof(pg_equatorial));
+
+		key->ra_low   = (float)eq->ra;
+		key->ra_high  = (float)eq->ra;
+		key->dec_low  = (float)eq->dec;
+		key->dec_high = (float)eq->dec;
+
+		retval = (GISTENTRY *) palloc(sizeof(GISTENTRY));
+		gistentryinit(*retval, PointerGetDatum(key),
+		              entry->rel, entry->page, entry->offset, false);
+	}
+	else
+	{
+		retval = entry;
+	}
+
+	PG_RETURN_POINTER(retval);
+}
+
+
+/* ================================================================
+ * GiST support function 4: decompress -- identity
+ * ================================================================
+ */
+Datum
+gist_eq_decompress(PG_FUNCTION_ARGS)
+{
+	PG_RETURN_POINTER(PG_GETARG_POINTER(0));
+}
+
+
+/* ================================================================
+ * GiST support function 5: penalty
+ *
+ * Cost of inserting a new entry into an existing subtree.
+ * Uses half-perimeter (margin) on the sphere: RA span normalized
+ * by 2*pi + Dec span normalized by pi.
+ * ================================================================
+ */
+Datum
+gist_eq_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);
+	eq_spherical_key *orig = (eq_spherical_key *) DatumGetPointer(origentry->key);
+	eq_spherical_key *add  = (eq_spherical_key *) DatumGetPointer(newentry->key);
+	eq_spherical_key  merged;
+	float             orig_margin, merged_margin;
+
+	orig_margin = ra_span(orig) / TWO_PI_F + dec_span(orig) / PI_F;
+
+	merged = *orig;
+	key_merge_eq(&merged, add);
+	merged_margin = ra_span(&merged) / TWO_PI_F + dec_span(&merged) / PI_F;
+
+	*penalty = merged_margin - orig_margin;
+
+	PG_RETURN_POINTER(penalty);
+}
+
+
+/* ================================================================
+ * GiST support function 6: picksplit
+ *
+ * Split an overfull page.  Compute spread in RA and Dec, split
+ * along the dimension with greater spread.  RA spread uses
+ * circular distance.
+ *
+ * GiST picksplit convention: entryvec->vector[] is 1-based.
+ * ================================================================
+ */
+
+typedef struct
+{
+	int    index;
+	float  sortval;
+} eq_picksplit_item;
+
+static int
+eq_picksplit_cmp(const void *a, const void *b)
+{
+	float ma = ((const eq_picksplit_item *) a)->sortval;
+	float mb = ((const eq_picksplit_item *) b)->sortval;
+
+	if (ma < mb) return -1;
+	if (ma > mb) return 1;
+	return 0;
+}
+
+Datum
+gist_eq_picksplit(PG_FUNCTION_ARGS)
+{
+	GistEntryVector  *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0);
+	GIST_SPLITVEC    *splitvec = (GIST_SPLITVEC *) PG_GETARG_POINTER(1);
+	OffsetNumber      maxoff   = entryvec->n - 1;
+	int               nentries = maxoff - FirstOffsetNumber + 1;
+	eq_picksplit_item *items;
+	eq_spherical_key  *left_union, *right_union;
+	eq_spherical_key  *cur;
+	int               split_at, i;
+	OffsetNumber      off;
+	float             ra_spread, dec_spread;
+	float             ra_min, ra_max, dec_min, dec_max;
+	bool              split_on_ra;
+
+	/* First pass: compute spread in both dimensions */
+	cur = (eq_spherical_key *) DatumGetPointer(
+			entryvec->vector[FirstOffsetNumber].key);
+
+	dec_min = (cur->dec_low + cur->dec_high) / 2.0f;
+	dec_max = dec_min;
+
+	/* For RA, find centroid spread using circular distance from first point */
+	ra_min = ra_midpoint(cur->ra_low, cur->ra_high);
+	ra_max = ra_min;
+
+	for (off = FirstOffsetNumber + 1; off <= maxoff; off++)
+	{
+		float ra_mid, dec_mid;
+
+		cur = (eq_spherical_key *) DatumGetPointer(entryvec->vector[off].key);
+		ra_mid  = ra_midpoint(cur->ra_low, cur->ra_high);
+		dec_mid = (cur->dec_low + cur->dec_high) / 2.0f;
+
+		if (dec_mid < dec_min) dec_min = dec_mid;
+		if (dec_mid > dec_max) dec_max = dec_mid;
+
+		/*
+		 * For RA spread, track min/max of the raw midpoint values.
+		 * This is an approximation -- the actual circular spread
+		 * could be larger if data clusters near 0h, but it works
+		 * well enough for the split heuristic.
+		 */
+		if (ra_mid < ra_min) ra_min = ra_mid;
+		if (ra_mid > ra_max) ra_max = ra_mid;
+	}
+
+	/* Normalize: RA by 2*pi, Dec by pi */
+	ra_spread  = (ra_max - ra_min) / TWO_PI_F;
+	dec_spread = (dec_max - dec_min) / PI_F;
+
+	/*
+	 * If RA points span > half the circle, the linear min/max
+	 * underestimates the true spread.  Use circular max distance.
+	 */
+	if (ra_spread > 0.5f)
+	{
+		/* Recompute with circular distances from the first point */
+		float ref_ra = ra_midpoint(
+			((eq_spherical_key *)DatumGetPointer(
+				entryvec->vector[FirstOffsetNumber].key))->ra_low,
+			((eq_spherical_key *)DatumGetPointer(
+				entryvec->vector[FirstOffsetNumber].key))->ra_high);
+		float max_dist = 0;
+
+		for (off = FirstOffsetNumber; off <= maxoff; off++)
+		{
+			float d;
+			cur = (eq_spherical_key *) DatumGetPointer(entryvec->vector[off].key);
+			d = ra_circular_dist(ref_ra, ra_midpoint(cur->ra_low, cur->ra_high));
+			if (d > max_dist) max_dist = d;
+		}
+		ra_spread = max_dist / TWO_PI_F;
+	}
+
+	split_on_ra = (ra_spread >= dec_spread);
+
+	/* Second pass: compute sort values.
+	 * For RA, rotate all midpoints relative to a reference so that
+	 * clusters straddling 0h are contiguous in sort order. */
+	items = (eq_picksplit_item *) palloc(sizeof(eq_picksplit_item) * nentries);
+	{
+		float ref_ra = 0.0f;
+
+		if (split_on_ra)
+		{
+			cur = (eq_spherical_key *) DatumGetPointer(
+					entryvec->vector[FirstOffsetNumber].key);
+			ref_ra = ra_midpoint(cur->ra_low, cur->ra_high);
+		}
+
+		for (i = 0, off = FirstOffsetNumber; off <= maxoff; i++, off++)
+		{
+			cur = (eq_spherical_key *) DatumGetPointer(entryvec->vector[off].key);
+			items[i].index = off;
+			if (split_on_ra)
+			{
+				float mid = ra_midpoint(cur->ra_low, cur->ra_high);
+				float rotated = mid - ref_ra;
+				if (rotated < 0.0f) rotated += TWO_PI_F;
+				items[i].sortval = rotated;
+			}
+			else
+			{
+				items[i].sortval = (cur->dec_low + cur->dec_high) / 2.0f;
+			}
+		}
+	}
+
+	qsort(items, nentries, sizeof(eq_picksplit_item), eq_picksplit_cmp);
+
+	split_at = nentries / 2;
+
+	splitvec->spl_left  = (OffsetNumber *) palloc(sizeof(OffsetNumber) * nentries);
+	splitvec->spl_right = (OffsetNumber *) palloc(sizeof(OffsetNumber) * nentries);
+	splitvec->spl_nleft  = 0;
+	splitvec->spl_nright = 0;
+
+	left_union  = (eq_spherical_key *) palloc0(sizeof(pg_equatorial));
+	right_union = (eq_spherical_key *) palloc0(sizeof(pg_equatorial));
+
+	cur = (eq_spherical_key *) DatumGetPointer(
+			entryvec->vector[items[0].index].key);
+	*left_union = *cur;
+
+	cur = (eq_spherical_key *) DatumGetPointer(
+			entryvec->vector[items[split_at].index].key);
+	*right_union = *cur;
+
+	for (i = 0; i < nentries; i++)
+	{
+		OffsetNumber idx = items[i].index;
+
+		cur = (eq_spherical_key *) DatumGetPointer(entryvec->vector[idx].key);
+
+		if (i < split_at)
+		{
+			splitvec->spl_left[splitvec->spl_nleft++] = idx;
+			key_merge_eq(left_union, cur);
+		}
+		else
+		{
+			splitvec->spl_right[splitvec->spl_nright++] = idx;
+			key_merge_eq(right_union, cur);
+		}
+	}
+
+	splitvec->spl_ldatum = PointerGetDatum(left_union);
+	splitvec->spl_rdatum = PointerGetDatum(right_union);
+
+	pfree(items);
+
+	PG_RETURN_POINTER(splitvec);
+}
+
+
+/* ================================================================
+ * GiST support function 7: same
+ *
+ * Equality test on compressed keys.
+ * ================================================================
+ */
+Datum
+gist_eq_same(PG_FUNCTION_ARGS)
+{
+	eq_spherical_key *a = (eq_spherical_key *) PG_GETARG_POINTER(0);
+	eq_spherical_key *b = (eq_spherical_key *) PG_GETARG_POINTER(1);
+	bool             *result = (bool *) PG_GETARG_POINTER(2);
+
+	*result = (fabsf(a->ra_low  - b->ra_low)  < EQ_KEY_EPSILON
+			&& fabsf(a->ra_high - b->ra_high) < EQ_KEY_EPSILON
+			&& fabsf(a->dec_low  - b->dec_low)  < EQ_KEY_EPSILON
+			&& fabsf(a->dec_high - b->dec_high) < EQ_KEY_EPSILON);
+
+	PG_RETURN_POINTER(result);
+}
+
+
+/* ================================================================
+ * GiST support function 8: distance
+ *
+ * Lower-bound angular distance (degrees) from a query equatorial
+ * position to any point inside the bounding box.
+ *
+ * Leaf nodes: exact Vincenty distance.
+ * Internal nodes: clamp query point to nearest box boundary,
+ * compute great-circle distance to clamped point.
+ *
+ * The lower-bound property satisfies GiST's KNN distance contract.
+ * ================================================================
+ */
+Datum
+gist_eq_distance(PG_FUNCTION_ARGS)
+{
+	GISTENTRY        *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
+	pg_equatorial    *query = (pg_equatorial *) PG_GETARG_POINTER(1);
+	eq_spherical_key *key   = (eq_spherical_key *) DatumGetPointer(entry->key);
+	double            dist_rad;
+
+	dist_rad = distance_point_to_box(query->ra, query->dec, key);
+
+	/* Return degrees, matching the <-> operator */
+	PG_RETURN_FLOAT8(dist_rad * (180.0 / M_PI));
+}
--- a/test/expected/gist_equatorial.out
+++ b/test/expected/gist_equatorial.out
@ -0,0 +1,206 @@
+-- Test equatorial GiST index: KNN ordering, RA wrapping, cone search
+CREATE EXTENSION IF NOT EXISTS pg_orrery;
+NOTICE:  extension "pg_orrery" already exists, skipping
+-- ============================================================
+-- Test table: known sky positions
+-- ============================================================
+CREATE TABLE sky_test (
+    id serial,
+    name text,
+    eq equatorial
+);
+-- Planets and Sun at a fixed epoch
+INSERT INTO sky_test (name, eq) VALUES
+    ('Jupiter', planet_equatorial_apparent(5, '2024-06-15 12:00:00+00')),
+    ('Saturn',  planet_equatorial_apparent(6, '2024-06-15 12:00:00+00')),
+    ('Mars',    planet_equatorial_apparent(4, '2024-06-15 12:00:00+00')),
+    ('Venus',   planet_equatorial_apparent(2, '2024-06-15 12:00:00+00')),
+    ('Mercury', planet_equatorial_apparent(1, '2024-06-15 12:00:00+00')),
+    ('Sun',     sun_equatorial('2024-06-15 12:00:00+00')),
+    ('Moon',    moon_equatorial('2024-06-15 12:00:00+00'));
+-- Bright stars at well-known positions
+INSERT INTO sky_test (name, eq) VALUES
+    ('Polaris',  star_equatorial(2.530, 89.264, '2024-06-15 12:00:00+00')),
+    ('Sirius',   star_equatorial(6.752, -16.716, '2024-06-15 12:00:00+00')),
+    ('Vega',     star_equatorial(18.616, 38.784, '2024-06-15 12:00:00+00')),
+    ('Canopus',  star_equatorial(6.399, -52.696, '2024-06-15 12:00:00+00')),
+    ('Arcturus', star_equatorial(14.261, 19.182, '2024-06-15 12:00:00+00'));
+-- RA-wrapping test: objects near 0h and 23.9h
+INSERT INTO sky_test (name, eq) VALUES
+    ('NearZeroH',  '(0.10000000,15.00000000,0.000)'::equatorial),
+    ('Near24H',    '(23.90000000,15.00000000,0.000)'::equatorial);
+-- ============================================================
+-- Test 1: Create GiST index
+-- ============================================================
+CREATE INDEX idx_sky_gist ON sky_test USING gist (eq);
+-- ============================================================
+-- Test 2: KNN correctness -- seqscan vs index scan
+-- Query: 5 nearest to Jupiter
+-- ============================================================
+-- First get seqscan ordering
+SET enable_indexscan = off;
+SET enable_bitmapscan = off;
+SELECT 'knn_seq' AS test, name,
+       round((eq <-> planet_equatorial_apparent(5, '2024-06-15 12:00:00+00'))::numeric, 4) AS dist
+FROM sky_test
+WHERE name != 'Jupiter'
+ORDER BY eq <-> planet_equatorial_apparent(5, '2024-06-15 12:00:00+00')
+LIMIT 5;
+  test   |  name   |  dist   
+---------+---------+---------
+ knn_seq | Sun     | 20.1241
+ knn_seq | Mercury | 21.1875
+ knn_seq | Venus   | 23.0885
+ knn_seq | Mars    | 30.0971
+ knn_seq | Sirius  | 53.0538
+(5 rows)
+
+RESET enable_indexscan;
+RESET enable_bitmapscan;
+-- Now force index scan
+SET enable_seqscan = off;
+SELECT 'knn_idx' AS test, name,
+       round((eq <-> planet_equatorial_apparent(5, '2024-06-15 12:00:00+00'))::numeric, 4) AS dist
+FROM sky_test
+WHERE name != 'Jupiter'
+ORDER BY eq <-> planet_equatorial_apparent(5, '2024-06-15 12:00:00+00')
+LIMIT 5;
+  test   |  name   |  dist   
+---------+---------+---------
+ knn_idx | Sun     | 20.1241
+ knn_idx | Mercury | 21.1875
+ knn_idx | Venus   | 23.0885
+ knn_idx | Mars    | 30.0971
+ knn_idx | Sirius  | 53.0538
+(5 rows)
+
+RESET enable_seqscan;
+-- ============================================================
+-- Test 3: KNN near Polaris (high declination)
+-- ============================================================
+SET enable_seqscan = off;
+SELECT 'knn_polaris' AS test, name,
+       round((eq <-> star_equatorial(2.530, 89.264, '2024-06-15 12:00:00+00'))::numeric, 2) AS dist
+FROM sky_test
+ORDER BY eq <-> star_equatorial(2.530, 89.264, '2024-06-15 12:00:00+00')
+LIMIT 3;
+    test     |  name   | dist  
+-------------+---------+-------
+ knn_polaris | Polaris |  0.00
+ knn_polaris | Vega    | 51.57
+ knn_polaris | Mercury | 65.08
+(3 rows)
+
+RESET enable_seqscan;
+-- ============================================================
+-- Test 4: RA wrapping -- NearZeroH and Near24H should be neighbors
+-- (They are only 0.2h * 15 deg/h * cos(15) ~ 2.9 deg apart)
+-- ============================================================
+SET enable_seqscan = off;
+SELECT 'ra_wrap' AS test, name,
+       round((eq <-> '(0.10000000,15.00000000,0.000)'::equatorial)::numeric, 2) AS dist
+FROM sky_test
+ORDER BY eq <-> '(0.10000000,15.00000000,0.000)'::equatorial
+LIMIT 3;
+  test   |   name    | dist  
+---------+-----------+-------
+ ra_wrap | NearZeroH |  0.00
+ ra_wrap | Near24H   |  2.90
+ ra_wrap | Saturn    | 23.50
+(3 rows)
+
+RESET enable_seqscan;
+-- ============================================================
+-- Test 5: Cone search -- everything within 15 degrees of Vega
+-- ============================================================
+SET enable_seqscan = off;
+SELECT 'cone_vega' AS test, name,
+       round((eq <-> star_equatorial(18.616, 38.784, '2024-06-15 12:00:00+00'))::numeric, 2) AS dist
+FROM sky_test
+WHERE eq_within_cone(eq, star_equatorial(18.616, 38.784, '2024-06-15 12:00:00+00'), 15.0)
+ORDER BY eq <-> star_equatorial(18.616, 38.784, '2024-06-15 12:00:00+00');
+   test    | name | dist 
+-----------+------+------
+ cone_vega | Vega | 0.00
+(1 row)
+
+RESET enable_seqscan;
+-- ============================================================
+-- Test 6: EXPLAIN shows Index Scan
+-- ============================================================
+SET enable_seqscan = off;
+EXPLAIN (COSTS OFF)
+SELECT name FROM sky_test
+ORDER BY eq <-> '(12.00000000,0.00000000,0.000)'::equatorial
+LIMIT 3;
+                               QUERY PLAN                                
+-------------------------------------------------------------------------
+ Limit
+   ->  Index Scan using idx_sky_gist on sky_test
+         Order By: (eq <-> '(12.00000000,0.00000000,0.000)'::equatorial)
+(3 rows)
+
+RESET enable_seqscan;
+-- ============================================================
+-- Test 7: Empty table doesn't crash
+-- ============================================================
+CREATE TABLE sky_empty (eq equatorial);
+CREATE INDEX idx_sky_empty ON sky_empty USING gist (eq);
+SELECT 'empty_knn' AS test, count(*) AS n
+FROM (
+    SELECT eq FROM sky_empty
+    ORDER BY eq <-> '(12.00000000,0.00000000,0.000)'::equatorial
+) sub;
+   test    | n 
+-----------+---
+ empty_knn | 0
+(1 row)
+
+DROP TABLE sky_empty;
+-- ============================================================
+-- Test 8: Single row
+-- ============================================================
+CREATE TABLE sky_single (eq equatorial);
+INSERT INTO sky_single VALUES ('(6.00000000,30.00000000,1000.000)'::equatorial);
+CREATE INDEX idx_sky_single ON sky_single USING gist (eq);
+SET enable_seqscan = off;
+SELECT 'single_knn' AS test,
+       round((eq <-> '(12.00000000,0.00000000,0.000)'::equatorial)::numeric, 2) AS dist
+FROM sky_single
+ORDER BY eq <-> '(12.00000000,0.00000000,0.000)'::equatorial
+LIMIT 1;
+    test    | dist  
+------------+-------
+ single_knn | 90.00
+(1 row)
+
+RESET enable_seqscan;
+DROP TABLE sky_single;
+-- ============================================================
+-- Test 9: Larger batch -- verify no crashes on tree rebalancing
+-- ============================================================
+CREATE TABLE sky_batch (eq equatorial);
+INSERT INTO sky_batch
+SELECT planet_equatorial_apparent(
+    (i % 7) + 1 + (CASE WHEN (i % 7) + 1 >= 3 THEN 1 ELSE 0 END),
+    '2024-01-01 00:00:00+00'::timestamptz + (i || ' hours')::interval
+)
+FROM generate_series(1, 100) AS i;
+CREATE INDEX idx_sky_batch ON sky_batch USING gist (eq);
+SET enable_seqscan = off;
+SELECT 'batch_knn' AS test, count(*) AS n
+FROM (
+    SELECT eq
+    FROM sky_batch
+    ORDER BY eq <-> '(12.00000000,0.00000000,0.000)'::equatorial
+    LIMIT 10
+) sub;
+   test    | n  
+-----------+----
+ batch_knn | 10
+(1 row)
+
+RESET enable_seqscan;
+DROP TABLE sky_batch;
+-- Cleanup
+DROP TABLE sky_test;
--- a/test/expected/v012_features.out
+++ b/test/expected/v012_features.out
@ -0,0 +1,133 @@
+-- v0.12.0 feature tests: DE moon equatorial functions
+-- ============================================================
+-- Test 1: galilean_equatorial_de fallback matches VSOP87 variant
+-- Without DE configured, DE variant should produce identical results
+-- ============================================================
+SELECT 'galilean_eq_de_fallback' AS test,
+       moon_id,
+       round(eq_ra(galilean_equatorial_de(moon_id, '2024-06-15 12:00:00+00'))::numeric, 4) AS de_ra,
+       round(eq_ra(galilean_equatorial(moon_id, '2024-06-15 12:00:00+00'))::numeric, 4) AS vsop_ra,
+       round(eq_ra(galilean_equatorial_de(moon_id, '2024-06-15 12:00:00+00'))::numeric, 4) =
+       round(eq_ra(galilean_equatorial(moon_id, '2024-06-15 12:00:00+00'))::numeric, 4) AS match
+FROM generate_series(0, 3) AS moon_id
+ORDER BY moon_id;
+          test           | moon_id | de_ra  | vsop_ra | match 
+-------------------------+---------+--------+---------+-------
+ galilean_eq_de_fallback |       0 | 4.1957 |  4.1957 | t
+ galilean_eq_de_fallback |       1 | 4.1950 |  4.1950 | t
+ galilean_eq_de_fallback |       2 | 4.1937 |  4.1937 | t
+ galilean_eq_de_fallback |       3 | 4.2057 |  4.2057 | t
+(4 rows)
+
+-- ============================================================
+-- Test 2: saturn_moon_equatorial_de fallback (Titan, id=5)
+-- ============================================================
+SELECT 'saturn_eq_de_fallback' AS test,
+       round(eq_ra(saturn_moon_equatorial_de(5, '2024-06-15 12:00:00+00'))::numeric, 4) AS de_ra,
+       round(eq_ra(saturn_moon_equatorial(5, '2024-06-15 12:00:00+00'))::numeric, 4) AS vsop_ra,
+       round(eq_ra(saturn_moon_equatorial_de(5, '2024-06-15 12:00:00+00'))::numeric, 4) =
+       round(eq_ra(saturn_moon_equatorial(5, '2024-06-15 12:00:00+00'))::numeric, 4) AS match;
+         test          |  de_ra  | vsop_ra | match 
+-----------------------+---------+---------+-------
+ saturn_eq_de_fallback | 23.3909 | 23.3909 | t
+(1 row)
+
+-- ============================================================
+-- Test 3: uranus_moon_equatorial_de fallback (Titania, id=3)
+-- ============================================================
+SELECT 'uranus_eq_de_fallback' AS test,
+       round(eq_ra(uranus_moon_equatorial_de(3, '2024-06-15 12:00:00+00'))::numeric, 4) AS de_ra,
+       round(eq_ra(uranus_moon_equatorial(3, '2024-06-15 12:00:00+00'))::numeric, 4) AS vsop_ra,
+       round(eq_ra(uranus_moon_equatorial_de(3, '2024-06-15 12:00:00+00'))::numeric, 4) =
+       round(eq_ra(uranus_moon_equatorial(3, '2024-06-15 12:00:00+00'))::numeric, 4) AS match;
+         test          | de_ra  | vsop_ra | match 
+-----------------------+--------+---------+-------
+ uranus_eq_de_fallback | 3.5124 |  3.5124 | t
+(1 row)
+
+-- ============================================================
+-- Test 4: mars_moon_equatorial_de fallback (Phobos + Deimos)
+-- ============================================================
+SELECT 'mars_eq_de_fallback' AS test,
+       moon_id,
+       round(eq_ra(mars_moon_equatorial_de(moon_id, '2024-06-15 12:00:00+00'))::numeric, 4) AS de_ra,
+       round(eq_ra(mars_moon_equatorial(moon_id, '2024-06-15 12:00:00+00'))::numeric, 4) AS vsop_ra,
+       round(eq_ra(mars_moon_equatorial_de(moon_id, '2024-06-15 12:00:00+00'))::numeric, 4) =
+       round(eq_ra(mars_moon_equatorial(moon_id, '2024-06-15 12:00:00+00'))::numeric, 4) AS match
+FROM generate_series(0, 1) AS moon_id
+ORDER BY moon_id;
+        test         | moon_id | de_ra  | vsop_ra | match 
+---------------------+---------+--------+---------+-------
+ mars_eq_de_fallback |       0 | 2.1851 |  2.1851 | t
+ mars_eq_de_fallback |       1 | 2.1851 |  2.1851 | t
+(2 rows)
+
+-- ============================================================
+-- Test 5: All DE moon equatorial return valid RA/Dec ranges
+-- ============================================================
+SELECT 'de_moon_eq_valid' AS test,
+       'galilean' AS family,
+       moon_id,
+       eq_ra(galilean_equatorial_de(moon_id, '2024-06-15 12:00:00+00')) BETWEEN 0 AND 24 AS ra_valid,
+       eq_dec(galilean_equatorial_de(moon_id, '2024-06-15 12:00:00+00')) BETWEEN -90 AND 90 AS dec_valid
+FROM generate_series(0, 3) AS moon_id
+ORDER BY moon_id;
+       test       |  family  | moon_id | ra_valid | dec_valid 
+------------------+----------+---------+----------+-----------
+ de_moon_eq_valid | galilean |       0 | t        | t
+ de_moon_eq_valid | galilean |       1 | t        | t
+ de_moon_eq_valid | galilean |       2 | t        | t
+ de_moon_eq_valid | galilean |       3 | t        | t
+(4 rows)
+
+-- ============================================================
+-- Test 6: Invalid body_id rejection for all 4 families
+-- ============================================================
+DO $$
+BEGIN
+    PERFORM galilean_equatorial_de(5, '2024-06-15 12:00:00+00');
+    RAISE EXCEPTION 'should have failed';
+EXCEPTION WHEN numeric_value_out_of_range THEN
+    RAISE NOTICE 'galilean_eq_de_invalid: correctly rejected';
+END;
+$$;
+NOTICE:  galilean_eq_de_invalid: correctly rejected
+DO $$
+BEGIN
+    PERFORM saturn_moon_equatorial_de(8, '2024-06-15 12:00:00+00');
+    RAISE EXCEPTION 'should have failed';
+EXCEPTION WHEN numeric_value_out_of_range THEN
+    RAISE NOTICE 'saturn_eq_de_invalid: correctly rejected';
+END;
+$$;
+NOTICE:  saturn_eq_de_invalid: correctly rejected
+DO $$
+BEGIN
+    PERFORM uranus_moon_equatorial_de(5, '2024-06-15 12:00:00+00');
+    RAISE EXCEPTION 'should have failed';
+EXCEPTION WHEN numeric_value_out_of_range THEN
+    RAISE NOTICE 'uranus_eq_de_invalid: correctly rejected';
+END;
+$$;
+NOTICE:  uranus_eq_de_invalid: correctly rejected
+DO $$
+BEGIN
+    PERFORM mars_moon_equatorial_de(2, '2024-06-15 12:00:00+00');
+    RAISE EXCEPTION 'should have failed';
+EXCEPTION WHEN numeric_value_out_of_range THEN
+    RAISE NOTICE 'mars_eq_de_invalid: correctly rejected';
+END;
+$$;
+NOTICE:  mars_eq_de_invalid: correctly rejected
+-- ============================================================
+-- Test 7: Negative body_id rejection
+-- ============================================================
+DO $$
+BEGIN
+    PERFORM galilean_equatorial_de(-1, '2024-06-15 12:00:00+00');
+    RAISE EXCEPTION 'should have failed';
+EXCEPTION WHEN numeric_value_out_of_range THEN
+    RAISE NOTICE 'galilean_eq_de_negative: correctly rejected';
+END;
+$$;
+NOTICE:  galilean_eq_de_negative: correctly rejected
--- a/test/sql/gist_equatorial.sql
+++ b/test/sql/gist_equatorial.sql
@ -0,0 +1,161 @@
+-- Test equatorial GiST index: KNN ordering, RA wrapping, cone search
+CREATE EXTENSION IF NOT EXISTS pg_orrery;
+
+-- ============================================================
+-- Test table: known sky positions
+-- ============================================================
+CREATE TABLE sky_test (
+    id serial,
+    name text,
+    eq equatorial
+);
+
+-- Planets and Sun at a fixed epoch
+INSERT INTO sky_test (name, eq) VALUES
+    ('Jupiter', planet_equatorial_apparent(5, '2024-06-15 12:00:00+00')),
+    ('Saturn',  planet_equatorial_apparent(6, '2024-06-15 12:00:00+00')),
+    ('Mars',    planet_equatorial_apparent(4, '2024-06-15 12:00:00+00')),
+    ('Venus',   planet_equatorial_apparent(2, '2024-06-15 12:00:00+00')),
+    ('Mercury', planet_equatorial_apparent(1, '2024-06-15 12:00:00+00')),
+    ('Sun',     sun_equatorial('2024-06-15 12:00:00+00')),
+    ('Moon',    moon_equatorial('2024-06-15 12:00:00+00'));
+
+-- Bright stars at well-known positions
+INSERT INTO sky_test (name, eq) VALUES
+    ('Polaris',  star_equatorial(2.530, 89.264, '2024-06-15 12:00:00+00')),
+    ('Sirius',   star_equatorial(6.752, -16.716, '2024-06-15 12:00:00+00')),
+    ('Vega',     star_equatorial(18.616, 38.784, '2024-06-15 12:00:00+00')),
+    ('Canopus',  star_equatorial(6.399, -52.696, '2024-06-15 12:00:00+00')),
+    ('Arcturus', star_equatorial(14.261, 19.182, '2024-06-15 12:00:00+00'));
+
+-- RA-wrapping test: objects near 0h and 23.9h
+INSERT INTO sky_test (name, eq) VALUES
+    ('NearZeroH',  '(0.10000000,15.00000000,0.000)'::equatorial),
+    ('Near24H',    '(23.90000000,15.00000000,0.000)'::equatorial);
+
+-- ============================================================
+-- Test 1: Create GiST index
+-- ============================================================
+CREATE INDEX idx_sky_gist ON sky_test USING gist (eq);
+
+-- ============================================================
+-- Test 2: KNN correctness -- seqscan vs index scan
+-- Query: 5 nearest to Jupiter
+-- ============================================================
+-- First get seqscan ordering
+SET enable_indexscan = off;
+SET enable_bitmapscan = off;
+SELECT 'knn_seq' AS test, name,
+       round((eq <-> planet_equatorial_apparent(5, '2024-06-15 12:00:00+00'))::numeric, 4) AS dist
+FROM sky_test
+WHERE name != 'Jupiter'
+ORDER BY eq <-> planet_equatorial_apparent(5, '2024-06-15 12:00:00+00')
+LIMIT 5;
+RESET enable_indexscan;
+RESET enable_bitmapscan;
+
+-- Now force index scan
+SET enable_seqscan = off;
+SELECT 'knn_idx' AS test, name,
+       round((eq <-> planet_equatorial_apparent(5, '2024-06-15 12:00:00+00'))::numeric, 4) AS dist
+FROM sky_test
+WHERE name != 'Jupiter'
+ORDER BY eq <-> planet_equatorial_apparent(5, '2024-06-15 12:00:00+00')
+LIMIT 5;
+RESET enable_seqscan;
+
+-- ============================================================
+-- Test 3: KNN near Polaris (high declination)
+-- ============================================================
+SET enable_seqscan = off;
+SELECT 'knn_polaris' AS test, name,
+       round((eq <-> star_equatorial(2.530, 89.264, '2024-06-15 12:00:00+00'))::numeric, 2) AS dist
+FROM sky_test
+ORDER BY eq <-> star_equatorial(2.530, 89.264, '2024-06-15 12:00:00+00')
+LIMIT 3;
+RESET enable_seqscan;
+
+-- ============================================================
+-- Test 4: RA wrapping -- NearZeroH and Near24H should be neighbors
+-- (They are only 0.2h * 15 deg/h * cos(15) ~ 2.9 deg apart)
+-- ============================================================
+SET enable_seqscan = off;
+SELECT 'ra_wrap' AS test, name,
+       round((eq <-> '(0.10000000,15.00000000,0.000)'::equatorial)::numeric, 2) AS dist
+FROM sky_test
+ORDER BY eq <-> '(0.10000000,15.00000000,0.000)'::equatorial
+LIMIT 3;
+RESET enable_seqscan;
+
+-- ============================================================
+-- Test 5: Cone search -- everything within 15 degrees of Vega
+-- ============================================================
+SET enable_seqscan = off;
+SELECT 'cone_vega' AS test, name,
+       round((eq <-> star_equatorial(18.616, 38.784, '2024-06-15 12:00:00+00'))::numeric, 2) AS dist
+FROM sky_test
+WHERE eq_within_cone(eq, star_equatorial(18.616, 38.784, '2024-06-15 12:00:00+00'), 15.0)
+ORDER BY eq <-> star_equatorial(18.616, 38.784, '2024-06-15 12:00:00+00');
+RESET enable_seqscan;
+
+-- ============================================================
+-- Test 6: EXPLAIN shows Index Scan
+-- ============================================================
+SET enable_seqscan = off;
+EXPLAIN (COSTS OFF)
+SELECT name FROM sky_test
+ORDER BY eq <-> '(12.00000000,0.00000000,0.000)'::equatorial
+LIMIT 3;
+RESET enable_seqscan;
+
+-- ============================================================
+-- Test 7: Empty table doesn't crash
+-- ============================================================
+CREATE TABLE sky_empty (eq equatorial);
+CREATE INDEX idx_sky_empty ON sky_empty USING gist (eq);
+SELECT 'empty_knn' AS test, count(*) AS n
+FROM (
+    SELECT eq FROM sky_empty
+    ORDER BY eq <-> '(12.00000000,0.00000000,0.000)'::equatorial
+) sub;
+DROP TABLE sky_empty;
+
+-- ============================================================
+-- Test 8: Single row
+-- ============================================================
+CREATE TABLE sky_single (eq equatorial);
+INSERT INTO sky_single VALUES ('(6.00000000,30.00000000,1000.000)'::equatorial);
+CREATE INDEX idx_sky_single ON sky_single USING gist (eq);
+SET enable_seqscan = off;
+SELECT 'single_knn' AS test,
+       round((eq <-> '(12.00000000,0.00000000,0.000)'::equatorial)::numeric, 2) AS dist
+FROM sky_single
+ORDER BY eq <-> '(12.00000000,0.00000000,0.000)'::equatorial
+LIMIT 1;
+RESET enable_seqscan;
+DROP TABLE sky_single;
+
+-- ============================================================
+-- Test 9: Larger batch -- verify no crashes on tree rebalancing
+-- ============================================================
+CREATE TABLE sky_batch (eq equatorial);
+INSERT INTO sky_batch
+SELECT planet_equatorial_apparent(
+    (i % 7) + 1 + (CASE WHEN (i % 7) + 1 >= 3 THEN 1 ELSE 0 END),
+    '2024-01-01 00:00:00+00'::timestamptz + (i || ' hours')::interval
+)
+FROM generate_series(1, 100) AS i;
+CREATE INDEX idx_sky_batch ON sky_batch USING gist (eq);
+SET enable_seqscan = off;
+SELECT 'batch_knn' AS test, count(*) AS n
+FROM (
+    SELECT eq
+    FROM sky_batch
+    ORDER BY eq <-> '(12.00000000,0.00000000,0.000)'::equatorial
+    LIMIT 10
+) sub;
+RESET enable_seqscan;
+DROP TABLE sky_batch;
+
+-- Cleanup
+DROP TABLE sky_test;
--- a/test/sql/v012_features.sql
+++ b/test/sql/v012_features.sql
@ -0,0 +1,106 @@
+-- v0.12.0 feature tests: DE moon equatorial functions
+
+-- ============================================================
+-- Test 1: galilean_equatorial_de fallback matches VSOP87 variant
+-- Without DE configured, DE variant should produce identical results
+-- ============================================================
+SELECT 'galilean_eq_de_fallback' AS test,
+       moon_id,
+       round(eq_ra(galilean_equatorial_de(moon_id, '2024-06-15 12:00:00+00'))::numeric, 4) AS de_ra,
+       round(eq_ra(galilean_equatorial(moon_id, '2024-06-15 12:00:00+00'))::numeric, 4) AS vsop_ra,
+       round(eq_ra(galilean_equatorial_de(moon_id, '2024-06-15 12:00:00+00'))::numeric, 4) =
+       round(eq_ra(galilean_equatorial(moon_id, '2024-06-15 12:00:00+00'))::numeric, 4) AS match
+FROM generate_series(0, 3) AS moon_id
+ORDER BY moon_id;
+
+-- ============================================================
+-- Test 2: saturn_moon_equatorial_de fallback (Titan, id=5)
+-- ============================================================
+SELECT 'saturn_eq_de_fallback' AS test,
+       round(eq_ra(saturn_moon_equatorial_de(5, '2024-06-15 12:00:00+00'))::numeric, 4) AS de_ra,
+       round(eq_ra(saturn_moon_equatorial(5, '2024-06-15 12:00:00+00'))::numeric, 4) AS vsop_ra,
+       round(eq_ra(saturn_moon_equatorial_de(5, '2024-06-15 12:00:00+00'))::numeric, 4) =
+       round(eq_ra(saturn_moon_equatorial(5, '2024-06-15 12:00:00+00'))::numeric, 4) AS match;
+
+-- ============================================================
+-- Test 3: uranus_moon_equatorial_de fallback (Titania, id=3)
+-- ============================================================
+SELECT 'uranus_eq_de_fallback' AS test,
+       round(eq_ra(uranus_moon_equatorial_de(3, '2024-06-15 12:00:00+00'))::numeric, 4) AS de_ra,
+       round(eq_ra(uranus_moon_equatorial(3, '2024-06-15 12:00:00+00'))::numeric, 4) AS vsop_ra,
+       round(eq_ra(uranus_moon_equatorial_de(3, '2024-06-15 12:00:00+00'))::numeric, 4) =
+       round(eq_ra(uranus_moon_equatorial(3, '2024-06-15 12:00:00+00'))::numeric, 4) AS match;
+
+-- ============================================================
+-- Test 4: mars_moon_equatorial_de fallback (Phobos + Deimos)
+-- ============================================================
+SELECT 'mars_eq_de_fallback' AS test,
+       moon_id,
+       round(eq_ra(mars_moon_equatorial_de(moon_id, '2024-06-15 12:00:00+00'))::numeric, 4) AS de_ra,
+       round(eq_ra(mars_moon_equatorial(moon_id, '2024-06-15 12:00:00+00'))::numeric, 4) AS vsop_ra,
+       round(eq_ra(mars_moon_equatorial_de(moon_id, '2024-06-15 12:00:00+00'))::numeric, 4) =
+       round(eq_ra(mars_moon_equatorial(moon_id, '2024-06-15 12:00:00+00'))::numeric, 4) AS match
+FROM generate_series(0, 1) AS moon_id
+ORDER BY moon_id;
+
+-- ============================================================
+-- Test 5: All DE moon equatorial return valid RA/Dec ranges
+-- ============================================================
+SELECT 'de_moon_eq_valid' AS test,
+       'galilean' AS family,
+       moon_id,
+       eq_ra(galilean_equatorial_de(moon_id, '2024-06-15 12:00:00+00')) BETWEEN 0 AND 24 AS ra_valid,
+       eq_dec(galilean_equatorial_de(moon_id, '2024-06-15 12:00:00+00')) BETWEEN -90 AND 90 AS dec_valid
+FROM generate_series(0, 3) AS moon_id
+ORDER BY moon_id;
+
+-- ============================================================
+-- Test 6: Invalid body_id rejection for all 4 families
+-- ============================================================
+DO $$
+BEGIN
+    PERFORM galilean_equatorial_de(5, '2024-06-15 12:00:00+00');
+    RAISE EXCEPTION 'should have failed';
+EXCEPTION WHEN numeric_value_out_of_range THEN
+    RAISE NOTICE 'galilean_eq_de_invalid: correctly rejected';
+END;
+$$;
+
+DO $$
+BEGIN
+    PERFORM saturn_moon_equatorial_de(8, '2024-06-15 12:00:00+00');
+    RAISE EXCEPTION 'should have failed';
+EXCEPTION WHEN numeric_value_out_of_range THEN
+    RAISE NOTICE 'saturn_eq_de_invalid: correctly rejected';
+END;
+$$;
+
+DO $$
+BEGIN
+    PERFORM uranus_moon_equatorial_de(5, '2024-06-15 12:00:00+00');
+    RAISE EXCEPTION 'should have failed';
+EXCEPTION WHEN numeric_value_out_of_range THEN
+    RAISE NOTICE 'uranus_eq_de_invalid: correctly rejected';
+END;
+$$;
+
+DO $$
+BEGIN
+    PERFORM mars_moon_equatorial_de(2, '2024-06-15 12:00:00+00');
+    RAISE EXCEPTION 'should have failed';
+EXCEPTION WHEN numeric_value_out_of_range THEN
+    RAISE NOTICE 'mars_eq_de_invalid: correctly rejected';
+END;
+$$;
+
+-- ============================================================
+-- Test 7: Negative body_id rejection
+-- ============================================================
+DO $$
+BEGIN
+    PERFORM galilean_equatorial_de(-1, '2024-06-15 12:00:00+00');
+    RAISE EXCEPTION 'should have failed';
+EXCEPTION WHEN numeric_value_out_of_range THEN
+    RAISE NOTICE 'galilean_eq_de_negative: correctly rejected';
+END;
+$$;