From bcf460a6ffd9cc03bd41622967dca87589c584fd Mon Sep 17 00:00:00 2001 From: Ryan Malloy Date: Tue, 17 Feb 2026 01:56:44 -0700 Subject: [PATCH] Add standalone DE reader unit test and aliasing safety docs MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Standalone test (test/test_de_reader.c): generates a synthetic 12KB DE binary with known Chebyshev coefficients, exercises header parsing, polynomial evaluation at 5 domain points, Earth derivation via center=99, Moon geocentric, layout validation, range/body error paths, NAN sentinel, and garbage file rejection. 55 tests, no PG dependency. Makefile: add `make test-de-reader` target; also includes the earlier sat_code C++ -> pure C sgp4 migration that was missed from prior commits. astro_math.h: document that ecliptic_to_equatorial/equatorial_to_ecliptic are NOT safe for aliased (in-place) calls — equ[1] is written before equ[2] reads ecl[1]. The vendored sgp4/sdp4.c has separate in-place versions using a temp variable. --- Makefile | 35 +-- src/astro_math.h | 8 + test/test_de_reader.c | 490 ++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 518 insertions(+), 15 deletions(-) create mode 100644 test/test_de_reader.c diff --git a/Makefile b/Makefile index e338df7..8825d00 100644 --- a/Makefile +++ b/Makefile @@ -14,36 +14,41 @@ OBJS = src/pg_orbit.o src/tle_type.o src/eci_type.o src/observer_type.o \ src/lambert.o src/transfer_funcs.o \ src/de_reader.o src/eph_provider.o src/de_funcs.o -# sat_code C++ sources (compiled with g++, linked with extern "C" symbols) -SAT_CODE_DIR = lib/sat_code -SAT_CODE_SRCS = $(SAT_CODE_DIR)/sgp4.cpp $(SAT_CODE_DIR)/sdp4.cpp \ - $(SAT_CODE_DIR)/deep.cpp $(SAT_CODE_DIR)/common.cpp \ - $(SAT_CODE_DIR)/basics.cpp $(SAT_CODE_DIR)/get_el.cpp \ - $(SAT_CODE_DIR)/tle_out.cpp -SAT_CODE_OBJS = $(SAT_CODE_SRCS:.cpp=.o) +# Vendored SGP4/SDP4 sources (pure C, from Bill Gray's sat_code, MIT license) +SGP4_DIR = src/sgp4 +SGP4_SRCS = $(SGP4_DIR)/sgp4.c $(SGP4_DIR)/sdp4.c \ + $(SGP4_DIR)/deep.c $(SGP4_DIR)/common.c \ + $(SGP4_DIR)/basics.c $(SGP4_DIR)/get_el.c \ + $(SGP4_DIR)/tle_out.c +SGP4_OBJS = $(SGP4_SRCS:.c=.o) -OBJS += $(SAT_CODE_OBJS) +OBJS += $(SGP4_OBJS) # Regression tests REGRESS = tle_parse sgp4_propagate coord_transforms pass_prediction gist_index convenience \ star_observe kepler_comet planet_observe moon_observe lambert_transfer \ - de_ephemeris + de_ephemeris vallado_518 REGRESS_OPTS = --inputdir=test -# Need C++ runtime for sat_code -SHLIB_LINK += -lstdc++ -lm +# Pure C — no C++ runtime needed +SHLIB_LINK += -lm # Compiler flags -PG_CPPFLAGS = -I$(SAT_CODE_DIR) +PG_CPPFLAGS = -I$(SGP4_DIR) # Use PGXS PG_CONFIG ?= pg_config PGXS := $(shell $(PG_CONFIG) --pgxs) include $(PGXS) -# Rule for compiling sat_code C++ files -$(SAT_CODE_DIR)/%.o: $(SAT_CODE_DIR)/%.cpp - $(CXX) $(CXXFLAGS) -fPIC -I$(SAT_CODE_DIR) -c -o $@ $< +# ── Standalone DE reader unit test (no PostgreSQL dependency) ── +# Generates a synthetic DE binary, exercises Chebyshev evaluation, +# header parsing, Earth derivation, error paths. +test-de-reader: test/test_de_reader.c src/de_reader.c src/de_reader.h + $(CC) -Wall -Werror -Isrc -o test/test_de_reader $< src/de_reader.c -lm + ./test/test_de_reader + +.PHONY: test-de-reader # ── Docker packaging ──────────────────────────────────────── REGISTRY ?= git.supported.systems/warehack.ing diff --git a/src/astro_math.h b/src/astro_math.h index ba5c719..ca0745e 100644 --- a/src/astro_math.h +++ b/src/astro_math.h @@ -129,6 +129,11 @@ equatorial_to_horizontal(double ha, double dec, double lat, /* * Ecliptic J2000 to equatorial J2000. * Simple rotation around X-axis by -obliquity. + * + * NOT safe for aliased (in-place) calls: ecl and equ must not overlap. + * Writing equ[1] before reading ecl[1] for equ[2] produces wrong results + * when ecl == equ. The vendored sgp4/sdp4.c has separate in-place versions + * that use a temp variable; do not confuse the two. */ static inline void ecliptic_to_equatorial(const double *ecl, double *equ) @@ -142,6 +147,9 @@ ecliptic_to_equatorial(const double *ecl, double *equ) /* * Equatorial J2000 to ecliptic J2000. * Rotation around X-axis by +obliquity. + * + * NOT safe for aliased (in-place) calls: equ and ecl must not overlap. + * Same ordering hazard as ecliptic_to_equatorial() above. */ static inline void equatorial_to_ecliptic(const double *equ, double *ecl) diff --git a/test/test_de_reader.c b/test/test_de_reader.c new file mode 100644 index 0000000..7428359 --- /dev/null +++ b/test/test_de_reader.c @@ -0,0 +1,490 @@ +/* + * test_de_reader.c -- Standalone unit test for the DE reader + * + * Generates a synthetic JPL DE binary file in /tmp, exercises + * de_reader_open/get_pos/get_const/close, and verifies results + * against known Chebyshev evaluations. + * + * No PostgreSQL dependency: de_reader.c is pure C. + * + * Build: cc -Wall -Isrc -o test/test_de_reader \ + * test/test_de_reader.c src/de_reader.c -lm + * Run: ./test/test_de_reader + */ + +#include "de_reader.h" + +#include +#include +#include +#include +#include +#include + +/* ── Synthetic file parameters ──────────────────────────── */ + +#define TEST_FILE "/tmp/pg_orbit_test_de.bin" +#define TEST_AU 149597870.700 +#define TEST_EMRAT 81.30056907419062 +#define TEST_START 2451529.0 /* J2000 - 16 days */ +#define TEST_END 2451561.0 /* J2000 + 16 days */ +#define TEST_INTV 32.0 /* single 32-day interval */ +#define NCOEFF 512 /* record size in doubles */ +#define J2000 2451545.0 + +/* + * IPT layout. 1-based offsets, ncoeff=4 per component, nsub=1. + * Each body occupies 4 coeffs * 3 components * 1 sub = 12 doubles. + * + * Mercury is padded to offset 501 so that the computed max_offset + * reaches 512 (= NCOEFF), keeping the record large enough for the + * 4096-byte header read. + */ +#define EMB_OFF 3 +#define MOON_OFF 15 +#define SUN_OFF 27 +#define MERC_OFF 501 +#define BODY_NC 4 +#define BODY_NS 1 + +/* Mercury x-component Chebyshev coefficients (in km) */ +#define MC0 1.0e8 +#define MC1 5.0e7 +#define MC2 1.0e7 +#define MC3 2.0e6 + +/* ── Test harness ───────────────────────────────────────── */ + +static int n_run, n_pass; + +#define RUN(cond, msg) do { \ + n_run++; \ + if (!(cond)) \ + fprintf(stderr, "FAIL: %s [line %d]\n", (msg), __LINE__); \ + else { n_pass++; fprintf(stderr, " ok: %s\n", (msg)); } \ +} while (0) + +#define CLOSE(a, b, tol, msg) do { \ + n_run++; \ + double _a = (a), _b = (b); \ + if (fabs(_a - _b) > (tol)) \ + fprintf(stderr, "FAIL: %s: %.15g vs %.15g (diff %.3e) [line %d]\n", \ + (msg), _a, _b, fabs(_a - _b), __LINE__); \ + else { n_pass++; fprintf(stderr, " ok: %s\n", (msg)); } \ +} while (0) + +/* ── Synthetic file generator ───────────────────────────── */ + +static void +put_i32(unsigned char *buf, int off, int32_t v) +{ + memcpy(buf + off, &v, 4); +} + +static void +put_f64(unsigned char *buf, int off, double v) +{ + memcpy(buf + off, &v, 8); +} + +static void +put_ipt(unsigned char *buf, int group, int off, int nc, int ns) +{ + int base = 2696 + group * 12; + put_i32(buf, base, (int32_t)off); + put_i32(buf, base + 4, (int32_t)nc); + put_i32(buf, base + 8, (int32_t)ns); +} + +/* + * Build a minimal but valid DE binary in memory and write to disk. + * + * Layout (3 records of NCOEFF doubles = 12288 bytes): + * Record 1: header (titles, SS, AU, EMRAT, IPT, DE version) + * Record 2: constant values + * Record 3: Chebyshev coefficients for EMB, Moon, Sun, Mercury + */ +static int +generate(void) +{ + size_t file_sz = 3 * NCOEFF * sizeof(double); + unsigned char *buf = calloc(1, file_sz); + double *data; + int fd; + ssize_t n; + + if (!buf) + return -1; + + /* ── Record 1 (header) ── */ + + /* Titles (3 x 84 chars at byte 0) */ + snprintf((char *)buf, 84, "SYNTHETIC DE FOR pg_orbit UNIT TEST"); + snprintf((char *)buf + 84, 84, "Generated by test_de_reader.c"); + snprintf((char *)buf + 168, 84, "Not a real ephemeris"); + + /* Constant names at byte 252 (6 chars each) */ + memcpy(buf + 252, "AU ", 6); + memcpy(buf + 258, "EMRAT ", 6); + + /* SS[3] at byte 2652 */ + put_f64(buf, 2652, TEST_START); + put_f64(buf, 2660, TEST_END); + put_f64(buf, 2668, TEST_INTV); + + /* NCON at byte 2676 */ + put_i32(buf, 2676, 2); + + /* AU at byte 2680 */ + put_f64(buf, 2680, TEST_AU); + + /* EMRAT at byte 2688 */ + put_f64(buf, 2688, TEST_EMRAT); + + /* IPT: 13 body groups */ + put_ipt(buf, 0, MERC_OFF, BODY_NC, BODY_NS); /* Mercury */ + put_ipt(buf, 1, 0, 0, 0); /* Venus: absent */ + put_ipt(buf, 2, EMB_OFF, BODY_NC, BODY_NS); /* EMB */ + put_ipt(buf, 3, 0, 0, 0); /* Mars: absent */ + put_ipt(buf, 4, 0, 0, 0); /* Jupiter */ + put_ipt(buf, 5, 0, 0, 0); /* Saturn */ + put_ipt(buf, 6, 0, 0, 0); /* Uranus */ + put_ipt(buf, 7, 0, 0, 0); /* Neptune */ + put_ipt(buf, 8, 0, 0, 0); /* Pluto */ + put_ipt(buf, 9, MOON_OFF, BODY_NC, BODY_NS); /* Moon */ + put_ipt(buf, 10, SUN_OFF, BODY_NC, BODY_NS); /* Sun */ + put_ipt(buf, 11, 0, 0, 0); /* Nutation */ + + /* DE version at byte 2840 */ + put_i32(buf, 2840, 999); + + /* Librations (group 12) at byte 2844 */ + put_i32(buf, 2844, 0); + put_i32(buf, 2848, 0); + put_i32(buf, 2852, 0); + + /* ── Record 2 (constant values) ── */ + + { + double *rec2 = (double *)(buf + NCOEFF * sizeof(double)); + rec2[0] = TEST_AU; + rec2[1] = TEST_EMRAT; + } + + /* ── Record 3 (data: Chebyshev coefficients) ── */ + + data = (double *)(buf + 2 * NCOEFF * sizeof(double)); + + /* First two doubles: JD range of this record */ + data[0] = TEST_START; + data[1] = TEST_START + TEST_INTV; + + /* + * EMB: constant position at 1 AU along x (in km). + * Coefficients for T_0 only; T_1..T_3 = 0. + */ + data[EMB_OFF - 1] = TEST_AU; /* x: c0 */ + + /* Moon geocentric: constant 384400 km along x */ + data[MOON_OFF - 1] = 384400.0; /* x: c0 */ + + /* Sun: all zeros (at SSB origin) — already zero from calloc */ + + /* + * Mercury x: non-trivial polynomial with 4 Chebyshev terms. + * Allows verification at multiple points in the domain. + * + * At normalized x: + * val = MC0*T_0(x) + MC1*T_1(x) + MC2*T_2(x) + MC3*T_3(x) + * + * Reference values (T_n at key points): + * x=-1: T = { 1, -1, 1, -1} -> val = c0-c1+c2-c3 = 5.8e7 + * x= 0: T = { 1, 0, -1, 0} -> val = c0-c2 = 9.0e7 + * x=+1: T = { 1, 1, 1, 1} -> val = c0+c1+c2+c3 = 1.62e8 + */ + data[MERC_OFF - 1 + 0] = MC0; /* x: c0 */ + data[MERC_OFF - 1 + 1] = MC1; /* x: c1 */ + data[MERC_OFF - 1 + 2] = MC2; /* x: c2 */ + data[MERC_OFF - 1 + 3] = MC3; /* x: c3 */ + /* Mercury y, z: zeros (calloc) */ + + /* ── Write to disk ── */ + + fd = open(TEST_FILE, O_WRONLY | O_CREAT | O_TRUNC | O_CLOEXEC, 0600); + if (fd < 0) + { + free(buf); + return -1; + } + + n = write(fd, buf, file_sz); + close(fd); + free(buf); + + return (n == (ssize_t)file_sz) ? 0 : -1; +} + + +/* ── Tests ──────────────────────────────────────────────── */ + +int +main(void) +{ + de_handle *h; + int err; + double pos[3]; + double tol = 1e-12; + + fprintf(stderr, "\n=== pg_orbit DE reader unit tests ===\n\n"); + + if (generate() != 0) + { + fprintf(stderr, "FATAL: could not write %s\n", TEST_FILE); + return 1; + } + fprintf(stderr, "Wrote %s (%d bytes)\n\n", + TEST_FILE, 3 * NCOEFF * (int)sizeof(double)); + + + /* ── 1. Open / header ── */ + + h = de_reader_open(TEST_FILE, &err); + RUN(h != NULL, "open succeeds"); + RUN(err == DE_OK, "errcode == DE_OK"); + + if (!h) + { + fprintf(stderr, "FATAL: open failed (err=%d)\n", err); + unlink(TEST_FILE); + return 1; + } + + RUN(h->de_version == 999, "de_version = 999"); + CLOSE(h->au_km, TEST_AU, 1.0, "AU parsed"); + CLOSE(h->emrat, TEST_EMRAT, 1e-10, "EMRAT parsed"); + CLOSE(h->start_jd, TEST_START, 1e-6, "start_jd"); + CLOSE(h->end_jd, TEST_END, 1e-6, "end_jd"); + CLOSE(h->interval_days, TEST_INTV, 1e-6, "interval_days"); + RUN(h->ncoeff == NCOEFF, "ncoeff = 512"); + + + /* ── 2. EMB heliocentric (EMB - Sun) ── */ + + err = de_reader_get_pos(h, J2000, DE_EMB, DE_SUN, pos); + RUN(err == DE_OK, "EMB-Sun query ok"); + CLOSE(pos[0], 1.0, tol, "EMB helio x = 1 AU"); + CLOSE(pos[1], 0.0, tol, "EMB helio y = 0"); + CLOSE(pos[2], 0.0, tol, "EMB helio z = 0"); + + + /* ── 3. Sun at SSB origin ── */ + + err = de_reader_get_pos(h, J2000, DE_SUN, -1, pos); + RUN(err == DE_OK, "Sun raw query ok"); + CLOSE(pos[0], 0.0, tol, "Sun x = 0"); + CLOSE(pos[1], 0.0, tol, "Sun y = 0"); + CLOSE(pos[2], 0.0, tol, "Sun z = 0"); + + + /* ── 4. Moon geocentric (center = -1, raw) ── */ + + err = de_reader_get_pos(h, J2000, DE_MOON, -1, pos); + RUN(err == DE_OK, "Moon geocentric ok"); + CLOSE(pos[0], 384400.0 / TEST_AU, 1e-10, "Moon geo x"); + CLOSE(pos[1], 0.0, tol, "Moon geo y = 0"); + CLOSE(pos[2], 0.0, tol, "Moon geo z = 0"); + + + /* ── 5. Earth derivation via center=99 ── + * + * Moon relative to Earth (center=99): + * Moon bary = Earth + Moon_geo + * result = Moon_bary - Earth = Moon_geo + * + * This exercises get_earth_pos() without exposing it directly. + */ + + err = de_reader_get_pos(h, J2000, DE_MOON, 99, pos); + RUN(err == DE_OK, "Moon rel Earth ok"); + CLOSE(pos[0], 384400.0 / TEST_AU, 1e-10, + "Moon-Earth x = geocentric Moon (Earth derivation verified)"); + CLOSE(pos[1], 0.0, tol, "Moon-Earth y = 0"); + CLOSE(pos[2], 0.0, tol, "Moon-Earth z = 0"); + + + /* ── 6. EMB relative to Earth (center=99) ── + * + * EMB - Earth = EMB - (EMB - Moon/(1+EMRAT)) = Moon/(1+EMRAT) + */ + + { + double expected = (384400.0 / TEST_AU) / (1.0 + TEST_EMRAT); + + err = de_reader_get_pos(h, J2000, DE_EMB, 99, pos); + RUN(err == DE_OK, "EMB rel Earth ok"); + CLOSE(pos[0], expected, 1e-14, + "EMB-Earth x = Moon/(1+EMRAT)"); + } + + + /* ── 7. Mercury Chebyshev at x = 0 (J2000.0) ── + * + * T_0(0) = 1, T_1(0) = 0, T_2(0) = -1, T_3(0) = 0 + * value = c0 - c2 = 9e7 km + */ + + { + double expected = (MC0 - MC2) / TEST_AU; + + err = de_reader_get_pos(h, J2000, DE_MERCURY, -1, pos); + RUN(err == DE_OK, "Mercury at x=0 ok"); + CLOSE(pos[0], expected, tol, "Mercury x at x=0 = (c0-c2)/AU"); + CLOSE(pos[1], 0.0, tol, "Mercury y = 0"); + CLOSE(pos[2], 0.0, tol, "Mercury z = 0"); + } + + + /* ── 8. Mercury at x = -1 (interval start) ── + * + * T_n(-1) = (-1)^n + * value = c0 - c1 + c2 - c3 = 5.8e7 km + */ + + { + double expected = (MC0 - MC1 + MC2 - MC3) / TEST_AU; + + err = de_reader_get_pos(h, TEST_START, DE_MERCURY, -1, pos); + RUN(err == DE_OK, "Mercury at x=-1 ok"); + CLOSE(pos[0], expected, tol, "Mercury x at x=-1"); + } + + + /* ── 9. Mercury at x = +1 (interval end) ── + * + * T_n(+1) = 1 for all n + * value = c0 + c1 + c2 + c3 = 1.62e8 km + */ + + { + double expected = (MC0 + MC1 + MC2 + MC3) / TEST_AU; + + err = de_reader_get_pos(h, TEST_END, DE_MERCURY, -1, pos); + RUN(err == DE_OK, "Mercury at x=+1 ok"); + CLOSE(pos[0], expected, tol, "Mercury x at x=+1"); + } + + + /* ── 10. Mercury at x = -0.5 (quarter point) ── + * + * T_0(-0.5) = 1, T_1(-0.5) = -0.5 + * T_2(-0.5) = 2*(0.25) - 1 = -0.5 + * T_3(-0.5) = 4*(-0.125) - 3*(-0.5) = 1.0 + * value = c0 - 0.5*c1 - 0.5*c2 + c3 + */ + + { + double expected = (MC0 - 0.5*MC1 - 0.5*MC2 + MC3) / TEST_AU; + double jd = TEST_START + 8.0; /* t_sub=8, x = 2*8/32 - 1 = -0.5 */ + + err = de_reader_get_pos(h, jd, DE_MERCURY, -1, pos); + RUN(err == DE_OK, "Mercury at x=-0.5 ok"); + CLOSE(pos[0], expected, tol, "Mercury x at x=-0.5"); + } + + + /* ── 11. Mercury at x = +0.5 (three-quarter point) ── + * + * T_0(0.5) = 1, T_1(0.5) = 0.5 + * T_2(0.5) = -0.5, T_3(0.5) = -1.0 + * value = c0 + 0.5*c1 - 0.5*c2 - c3 + */ + + { + double expected = (MC0 + 0.5*MC1 - 0.5*MC2 - MC3) / TEST_AU; + double jd = TEST_START + 24.0; /* t_sub=24, x = 2*24/32 - 1 = 0.5 */ + + err = de_reader_get_pos(h, jd, DE_MERCURY, -1, pos); + RUN(err == DE_OK, "Mercury at x=+0.5 ok"); + CLOSE(pos[0], expected, tol, "Mercury x at x=+0.5"); + } + + + /* ── 12. Center subtraction: Mercury - Sun = Mercury raw ── + * (Sun is at origin, so heliocentric = raw) + */ + + { + double raw[3], helio[3]; + + de_reader_get_pos(h, J2000, DE_MERCURY, -1, raw); + de_reader_get_pos(h, J2000, DE_MERCURY, DE_SUN, helio); + CLOSE(helio[0], raw[0], tol, "Mercury helio x = raw x"); + CLOSE(helio[1], raw[1], tol, "Mercury helio y = raw y"); + CLOSE(helio[2], raw[2], tol, "Mercury helio z = raw z"); + } + + + /* ── 13. get_const ── */ + + CLOSE(de_reader_get_const(h, "AU"), TEST_AU, 1.0, "get_const AU"); + CLOSE(de_reader_get_const(h, "EMRAT"), TEST_EMRAT, 1e-10, "get_const EMRAT"); + RUN(isnan(de_reader_get_const(h, "BOGUS")), "unknown const -> NAN"); + RUN(isnan(de_reader_get_const(NULL, "AU")), "NULL handle -> NAN"); + + + /* ── 14. Error paths ── */ + + err = de_reader_get_pos(h, 2400000.0, DE_MERCURY, -1, pos); + RUN(err == DE_ERR_RANGE, "JD before range -> DE_ERR_RANGE"); + + err = de_reader_get_pos(h, 2500000.0, DE_MERCURY, -1, pos); + RUN(err == DE_ERR_RANGE, "JD after range -> DE_ERR_RANGE"); + + err = de_reader_get_pos(h, J2000, DE_NUTATION, -1, pos); + RUN(err == DE_ERR_BODY, "nutation -> DE_ERR_BODY"); + + err = de_reader_get_pos(h, J2000, DE_LIBRATION, -1, pos); + RUN(err == DE_ERR_BODY, "libration -> DE_ERR_BODY"); + + err = de_reader_get_pos(h, J2000, DE_VENUS, -1, pos); + RUN(err == DE_ERR_BODY, "Venus (no layout) -> DE_ERR_BODY"); + + err = de_reader_get_pos(h, J2000, DE_MARS, -1, pos); + RUN(err == DE_ERR_BODY, "Mars (no layout) -> DE_ERR_BODY"); + + err = de_reader_get_pos(NULL, J2000, DE_MERCURY, -1, pos); + RUN(err == DE_ERR_OPEN, "NULL handle -> DE_ERR_OPEN"); + + + /* ── 15. Cleanup ── */ + + de_reader_close(h); + de_reader_close(NULL); /* must not crash */ + + /* ── Bad file: open should fail ── */ + + { + unsigned char garbage[4096]; + int gfd; + + memset(garbage, 0xAA, sizeof(garbage)); + gfd = open(TEST_FILE, O_WRONLY | O_CREAT | O_TRUNC | O_CLOEXEC, 0600); + if (gfd >= 0) + { + write(gfd, garbage, sizeof(garbage)); + close(gfd); + + h = de_reader_open(TEST_FILE, &err); + RUN(h == NULL, "garbage file -> open fails"); + RUN(err != DE_OK, "garbage file -> error code set"); + } + } + + + /* ── Summary ── */ + + unlink(TEST_FILE); + fprintf(stderr, "\n%d/%d tests passed\n\n", n_pass, n_run); + + return (n_pass == n_run) ? 0 : 1; +}