Harden DE reader: layout validation, Chebyshev clamping, FD safety

Three independent reviews (manual, failure-mode analysis, JPL spec
cross-reference) confirmed the mathematical core is correct. This
commit addresses defensive coding and operational behavior:

- Fix header byte-offset comments (12 groups = 144 bytes, not 156)
- Add layout validation before Chebyshev interpolation (prevent
  buffer underread for bodies absent from a DE edition)
- Clamp Chebyshev argument to [-1,+1] with debug assertion for
  values beyond 1e-10 tolerance (catches structural normalization
  errors vs normal FP boundary rounding)
- Add O_CLOEXEC to prevent FD leaks to child processes
- Change GUC from PGC_SIGHUP to PGC_BACKEND to match actual
  one-shot initialization behavior
- Fix provider consistency: planet_velocity_de() now accepts a
  use_de flag to match the provider used for positions (rule 7)
- Optimize eph_de_moon() to use raw geocentric Moon (center=-1)
  instead of computing Earth just to subtract it back out
- Pre-compute obliquity trig constants (verified to full precision)
- Tighten canary check from 0.9-1.1 AU to 0.97-1.04 AU
- Return NAN for missing constants (0.0 was ambiguous)
- Add _Static_assert for sizeof(double) == 8
- Remove unused HDR_* macros
- Zero Datum values before setting null flags in ephemeris_info
- Replace magic numbers with DE_MOON/DE_SUN constants

All 13 regression tests pass. Zero compiler warnings.

src/astro_math.h

@@ -18,6 +18,11 @@
 #define RAD_TO_DEG  (180.0 / M_PI)
 #define ARCSEC_TO_RAD (M_PI / (180.0 * 3600.0))
 
+/* Pre-computed obliquity trig (IAU 1976, OBLIQUITY_J2000 = 0.40909280422232897 rad).
+ * Avoids recomputing cos/sin on every frame rotation call. */
+#define COS_OBLIQUITY_J2000  0.91748206206918181
+#define SIN_OBLIQUITY_J2000  0.39777715593191371
+
 /*
  * Greenwich Mean Sidereal Time from Julian date.
  * Vallado, "Fundamentals of Astrodynamics", Eq. 3-47.
@@ -128,12 +133,9 @@ equatorial_to_horizontal(double ha, double dec, double lat,
 static inline void
 ecliptic_to_equatorial(const double *ecl, double *equ)
 {
-    double cos_eps = cos(OBLIQUITY_J2000);
-    double sin_eps = sin(OBLIQUITY_J2000);
-
     equ[0] = ecl[0];
-    equ[1] = ecl[1] * cos_eps - ecl[2] * sin_eps;
+    equ[1] = ecl[1] * COS_OBLIQUITY_J2000 - ecl[2] * SIN_OBLIQUITY_J2000;
-    equ[2] = ecl[1] * sin_eps + ecl[2] * cos_eps;
+    equ[2] = ecl[1] * SIN_OBLIQUITY_J2000 + ecl[2] * COS_OBLIQUITY_J2000;
 }
 
 
@@ -144,12 +146,9 @@ ecliptic_to_equatorial(const double *ecl, double *equ)
 static inline void
 equatorial_to_ecliptic(const double *equ, double *ecl)
 {
-    double cos_eps = cos(OBLIQUITY_J2000);
-    double sin_eps = sin(OBLIQUITY_J2000);
-
     ecl[0] = equ[0];
-    ecl[1] = equ[1] * cos_eps + equ[2] * sin_eps;
+    ecl[1] = equ[1] * COS_OBLIQUITY_J2000 + equ[2] * SIN_OBLIQUITY_J2000;
-    ecl[2] = -equ[1] * sin_eps + equ[2] * cos_eps;
+    ecl[2] = -equ[1] * SIN_OBLIQUITY_J2000 + equ[2] * COS_OBLIQUITY_J2000;
 }
 
 

src/de_funcs.c

@@ -251,21 +251,32 @@ moon_observe_de(PG_FUNCTION_ARGS)
 
 /*
  * Compute planet heliocentric velocity via numerical differentiation.
- * Uses DE positions if available, VSOP87 otherwise.
+ *
+ * use_de: if true, use DE positions; if false, use VSOP87.
+ * Must match the provider used for the corresponding position query
+ * (rule 7: same provider for position and velocity).
  */
 static void
-planet_velocity_de(int body_id, double jd, double vel[3])
+planet_velocity_de(int body_id, double jd, bool use_de, double vel[3])
 {
     double pos_fwd[6], pos_bwd[6];
     double dt = 0.01;  /* days */
-    bool got_fwd, got_bwd;
 
-    got_fwd = eph_de_planet(body_id, jd + dt, pos_fwd);
+    if (use_de)
-    got_bwd = eph_de_planet(body_id, jd - dt, pos_bwd);
+    {
-
+        bool got_fwd = eph_de_planet(body_id, jd + dt, pos_fwd);
-    if (!got_fwd || !got_bwd)
+        bool got_bwd = eph_de_planet(body_id, jd - dt, pos_bwd);
+
+        if (!got_fwd || !got_bwd)
+        {
+            /* DE boundary straddled — use VSOP87 for both to stay consistent */
+            int vsop_body = body_id - 1;
+            GetVsop87Coor(jd + dt, vsop_body, pos_fwd);
+            GetVsop87Coor(jd - dt, vsop_body, pos_bwd);
+        }
+    }
+    else
     {
-        /* Fall back to VSOP87 */
         int vsop_body = body_id - 1;
         GetVsop87Coor(jd + dt, vsop_body, pos_fwd);
         GetVsop87Coor(jd - dt, vsop_body, pos_bwd);
@@ -341,9 +352,9 @@ lambert_transfer_de(PG_FUNCTION_ARGS)
     if (!lambert_solve_uv(r1, r2, tof_days, GAUSS_K2, 1, &lr))
         PG_RETURN_NULL();
 
-    /* Planet velocities */
+    /* Planet velocities (same provider as positions — rule 7) */
-    planet_velocity_de(dep_body, dep_jd, v_planet_dep);
+    planet_velocity_de(dep_body, dep_jd, have_de, v_planet_dep);
-    planet_velocity_de(arr_body, arr_jd, v_planet_arr);
+    planet_velocity_de(arr_body, arr_jd, have_de, v_planet_arr);
 
     au_per_day_to_km_per_s = AU_KM / 86400.0;
 
@@ -427,7 +438,7 @@ lambert_c3_de(PG_FUNCTION_ARGS)
     if (!lambert_solve_uv(r1, r2, tof_days, GAUSS_K2, 1, &lr))
         PG_RETURN_NULL();
 
-    planet_velocity_de(dep_body, dep_jd, v_planet_dep);
+    planet_velocity_de(dep_body, dep_jd, have_de, v_planet_dep);
 
     au_per_day_to_km_per_s = AU_KM / 86400.0;
 
@@ -638,6 +649,10 @@ pg_orbit_ephemeris_info(PG_FUNCTION_ARGS)
     {
         values[0] = CStringGetTextDatum("VSOP87");
 
+        values[1] = (Datum) 0;
+        values[2] = (Datum) 0;
+        values[3] = (Datum) 0;
+        values[4] = (Datum) 0;
         nulls[1] = true;   /* no file path */
         nulls[2] = true;   /* no start_jd */
         nulls[3] = true;   /* no end_jd */

src/de_reader.c

@@ -22,19 +22,19 @@
 
 #include "de_reader.h"
 
+#include <assert.h>
 #include <stdlib.h>
 #include <string.h>
 #include <math.h>
 #include <unistd.h>
 #include <fcntl.h>
 
+/* swap_double() assumes 8-byte doubles (IEEE 754 on all platforms PG supports) */
+_Static_assert(sizeof(double) == 8, "DE reader requires 8-byte doubles");
+
 /* Known AU value for byte-order detection */
 #define DE_AU_KNOWN     149597870.700
 
-/* Header offsets (in doubles, 0-based within record 1) */
-#define HDR_TITLE_WORDS     84          /* 3 title lines x 84 chars = 252 chars = ~32 doubles */
-#define HDR_CNAME_WORDS     400         /* constant names */
-#define HDR_SS_OFFSET       (HDR_TITLE_WORDS / sizeof(double))  /* not used directly */
 
 /*
  * Byte-swap a double (for big-endian DE files on little-endian hosts).
@@ -145,6 +145,22 @@ interp_component(de_handle *h, int body, int comp, double jd)
     /* Normalize to [-1, +1] within the sub-interval */
     x = 2.0 * (t_sub - sub_idx * sub_length) / sub_length - 1.0;
 
+    /* Clamp to [-1, +1]: floating-point arithmetic at interval
+     * boundaries can produce |x| slightly > 1.0, which causes
+     * Chebyshev polynomials to diverge. Legitimate rounding should
+     * never exceed ~1e-14; anything larger indicates a structural
+     * arithmetic error in the normalization above. */
+    if (x > 1.0)
+    {
+        assert(x < 1.0 + 1e-10);
+        x = 1.0;
+    }
+    if (x < -1.0)
+    {
+        assert(x > -1.0 - 1e-10);
+        x = -1.0;
+    }
+
     /*
      * Coefficient offset in the record buffer.
      * Layout offset is 1-based (Fortran convention), convert to 0-based.
@@ -200,6 +216,8 @@ load_record(de_handle *h, double jd)
 
 /*
  * Parse the header and validate the ephemeris file.
+ * Supports DE405 and later (DE430, DE440, DE441).
+ * Earlier formats have different header sizes and are not supported.
  */
 static int
 parse_header(de_handle *h)
@@ -216,9 +234,9 @@ parse_header(de_handle *h)
      *   2676..2679: NCON (int32)
      *   2680..2687: AU (double)
      *   2688..2695: EMRAT (double)
-     *   2696..2851: IPT[12][3] = 12 body groups x 3 ints (offset, ncoeff, nsub)
+     *   2696..2839: IPT[12][3] = 12 body groups x 3 ints (offset, ncoeff, nsub)
-     *   2852..2855: DE version number (int32)
+     *   2840..2843: DE version number (int32)
-     *   2856..2867: IPT[12][3] for 13th body (librations)
+     *   2844..2855: IPT[12][3] for 13th body (librations)
      *
      * After IPT parsing, we know ncoeff and can size the record buffer.
      *
@@ -226,7 +244,7 @@ parse_header(de_handle *h)
      * The record size in bytes = ncoeff * 8.
      */
 
-    unsigned char buf[4096];    /* enough for the header area */
+    unsigned char buf[4096];    /* DE405+ header fits in 4096 bytes */
     double  ss[3];
     double  au_val;
     int32_t ncon;
@@ -388,10 +406,11 @@ canary_check(de_handle *h)
     if (err != DE_OK)
         return DE_ERR_CANARY;
 
-    /* Rough check: Earth should be ~1 AU from Sun */
+    /* Earth-Sun distance varies 0.983-1.017 AU over the year.
+     * Tight tolerance catches garbled files that 0.9-1.1 would miss. */
     {
         double dist = sqrt(pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2]);
-        if (dist < 0.9 || dist > 1.1)
+        if (dist < 0.97 || dist > 1.04)
             return DE_ERR_CANARY;
     }
 
@@ -416,8 +435,9 @@ de_reader_open(const char *path, int *errcode)
     h->cached_recno = -1;
     h->record_buf = NULL;
 
-    /* Open the ephemeris file (read-only, no O_CLOEXEC dependency) */
+    /* Open the ephemeris file (read-only, close-on-exec to prevent FD
-    h->fd = open(path, O_RDONLY);
+     * leaks to child processes such as COPY ... PROGRAM or archiving) */
+    h->fd = open(path, O_RDONLY | O_CLOEXEC);
     if (h->fd < 0)
     {
         *errcode = DE_ERR_OPEN;
@@ -477,6 +497,16 @@ get_raw_pos(de_handle *h, double jd, int body, double pos[3])
     if (body == DE_NUTATION || body == DE_LIBRATION)
         return DE_ERR_BODY;
 
+    /* Body exists in the DE body table but may have no coefficients
+     * in this particular DE edition (e.g., Pluto omitted from some files).
+     * Returns DE_ERR_BODY since callers handle both "invalid index" and
+     * "not present" by falling back to VSOP87. */
+    {
+        de_body_layout *lay = &h->layout[body];
+        if (lay->offset <= 0 || lay->ncoeff <= 0 || lay->nsub <= 0)
+            return DE_ERR_BODY;
+    }
+
     err = load_record(h, jd);
     if (err != DE_OK)
         return err;
@@ -651,12 +681,14 @@ de_reader_get_const(de_handle *h, const char *name)
      * that are already parsed from the header.
      */
     if (!h || !name)
-        return 0.0;
+        return NAN;
 
     if (strcmp(name, "AU") == 0)
         return h->au_km;
     if (strcmp(name, "EMRAT") == 0)
         return h->emrat;
 
-    return 0.0;
+    /* NAN sentinel: callers use isnan() to distinguish "not found"
+     * from a legitimate zero-valued constant. */
+    return NAN;
 }

src/de_reader.h

@@ -133,7 +133,7 @@ void de_reader_close(de_handle *h);
 
 /*
  * Look up a named constant from the DE header.
- * Returns the value, or 0.0 if not found.
+ * Returns the value, or NAN if not found (use isnan() to check).
  */
 double de_reader_get_const(de_handle *h, const char *name);
 

src/eph_provider.c

@@ -54,7 +54,10 @@ eph_register_gucs(void)
         "Relative paths are relative to $PGDATA.",
         &eph_path_guc,
         "",                     /* default: empty = disabled */
-        PGC_SIGHUP,             /* changeable by SIGHUP reload */
+        PGC_BACKEND,            /* Set at backend startup only.
+                                 * Live reload (PGC_SIGHUP) would require an
+                                 * assign_hook that resets de_init_attempted,
+                                 * closes any open handle, and re-initializes. */
         GUC_SUPERUSER_ONLY,     /* only superuser can set */
         NULL,                   /* check_hook */
         NULL,                   /* assign_hook */
@@ -179,7 +182,10 @@ de_get_helio_ecliptic(int de_target, double jd, double ecl[3])
 
 /*
  * Internal: get Earth's heliocentric position from DE.
- * Earth = EMB - Moon/(1+EMRAT), already handled by de_reader.
+ *
+ * Earth = EMB(helio) - Moon(geocentric) / (1 + EMRAT).
+ * The EMB-to-Earth correction is also in de_reader.c:get_earth_pos();
+ * both must use the same formula.
  */
 static bool
 de_get_earth_helio_ecliptic(double jd, double ecl[3])
@@ -190,21 +196,20 @@ de_get_earth_helio_ecliptic(double jd, double ecl[3])
     if (!de_handle_ptr)
         return false;
 
-    /* Earth relative to Sun: use EMB and subtract Moon component */
+    /* EMB heliocentric, then correct to Earth below */
-    err = de_reader_get_pos(de_handle_ptr, jd, DE_EMB, 10, icrs);
+    err = de_reader_get_pos(de_handle_ptr, jd, DE_EMB, DE_SUN, icrs);
     if (err != DE_OK)
         return false;
 
-    /* Correct EMB to Earth: subtract Moon/(1+EMRAT) */
+    /* Correct EMB-helio to Earth-helio: subtract Moon/(1+EMRAT) */
     {
         double moon_icrs[3];
         double factor;
 
-        err = de_reader_get_pos(de_handle_ptr, jd, 9, -1, moon_icrs);
+        err = de_reader_get_pos(de_handle_ptr, jd, DE_MOON, -1, moon_icrs);
         if (err != DE_OK)
             return false;
 
-        /* Moon in DE is geocentric in km->AU. Factor for EMB->Earth */
         factor = 1.0 / (1.0 + de_handle_ptr->emrat);
         icrs[0] -= moon_icrs[0] * factor;
         icrs[1] -= moon_icrs[1] * factor;
@@ -283,8 +288,11 @@ eph_de_moon(double jd, double xyz[3])
     if (!eph_de_available())
         return false;
 
-    /* Moon geocentric: target=Moon(9), center=Earth(99) */
+    /* Moon geocentric: use center=-1 for raw geocentric Moon directly.
-    err = de_reader_get_pos(de_handle_ptr, jd, 9, 99, icrs);
+     * The Moon is stored geocentric in the DE file, so center=-1 avoids
+     * the roundabout of computing Earth (EMB - Moon/(1+EMRAT)) just to
+     * subtract it back out (~4x fewer interpolations). */
+    err = de_reader_get_pos(de_handle_ptr, jd, DE_MOON, -1, icrs);
     if (err != DE_OK)
         return false;