diff --git a/src/astro_math.h b/src/astro_math.h index 5cc27e9..ba5c719 100644 --- a/src/astro_math.h +++ b/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[2] = ecl[1] * sin_eps + ecl[2] * cos_eps; + equ[1] = ecl[1] * COS_OBLIQUITY_J2000 - ecl[2] * SIN_OBLIQUITY_J2000; + 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[2] = -equ[1] * sin_eps + equ[2] * cos_eps; + ecl[1] = equ[1] * COS_OBLIQUITY_J2000 + equ[2] * SIN_OBLIQUITY_J2000; + ecl[2] = -equ[1] * SIN_OBLIQUITY_J2000 + equ[2] * COS_OBLIQUITY_J2000; } diff --git a/src/de_funcs.c b/src/de_funcs.c index 75e97b8..39b1682 100644 --- a/src/de_funcs.c +++ b/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); - got_bwd = eph_de_planet(body_id, jd - dt, pos_bwd); - - if (!got_fwd || !got_bwd) + if (use_de) + { + bool got_fwd = eph_de_planet(body_id, jd + dt, pos_fwd); + 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_velocity_de(dep_body, dep_jd, v_planet_dep); - planet_velocity_de(arr_body, arr_jd, v_planet_arr); + /* Planet velocities (same provider as positions — rule 7) */ + planet_velocity_de(dep_body, dep_jd, have_de, v_planet_dep); + 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 */ diff --git a/src/de_reader.c b/src/de_reader.c index a11fb1a..0695d5e 100644 --- a/src/de_reader.c +++ b/src/de_reader.c @@ -22,19 +22,19 @@ #include "de_reader.h" +#include #include #include #include #include #include +/* 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) - * 2852..2855: DE version number (int32) - * 2856..2867: IPT[12][3] for 13th body (librations) + * 2696..2839: IPT[12][3] = 12 body groups x 3 ints (offset, ncoeff, nsub) + * 2840..2843: DE version number (int32) + * 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) */ - h->fd = open(path, O_RDONLY); + /* Open the ephemeris file (read-only, close-on-exec to prevent FD + * 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; } diff --git a/src/de_reader.h b/src/de_reader.h index f0ffeeb..52db174 100644 --- a/src/de_reader.h +++ b/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); diff --git a/src/eph_provider.c b/src/eph_provider.c index a42a8e1..3f4366a 100644 --- a/src/eph_provider.c +++ b/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 */ - err = de_reader_get_pos(de_handle_ptr, jd, DE_EMB, 10, icrs); + /* EMB heliocentric, then correct to Earth below */ + 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) */ - err = de_reader_get_pos(de_handle_ptr, jd, 9, 99, icrs); + /* Moon geocentric: use center=-1 for raw geocentric Moon directly. + * 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;