/* * 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_orrery_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_orrery 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_orrery 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; }