From dfd085f17666c7f5b85c18e3ad0341311158b6a7 Mon Sep 17 00:00:00 2001 From: Ryan Malloy Date: Sat, 28 Feb 2026 14:01:55 -0700 Subject: [PATCH 1/2] Add CR3BP Lagrange point solver (pure math, no PG dependency) Quintic Newton-Raphson for L1/L2/L3, analytic L4/L5. Includes Sun-planet, Earth-Moon, and planet-moon mass ratio constants from IAU 2012 / JPL DE441. Co-rotating to ecliptic J2000 frame transform. Hill sphere and libration zone radius. 210/210 standalone tests pass. --- Makefile | 8 + src/lagrange.h | 486 +++++++++++++++++++++++++++++++++++++++++++ test/test_lagrange.c | 459 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 953 insertions(+) create mode 100644 src/lagrange.h create mode 100644 test/test_lagrange.c diff --git a/Makefile b/Makefile index b3d52cc..e600459 100644 --- a/Makefile +++ b/Makefile @@ -85,6 +85,14 @@ test-de-reader: test/test_de_reader.c src/de_reader.c src/de_reader.h .PHONY: test-de-reader +# ── Standalone Lagrange solver unit test (no PostgreSQL dependency) ── +# CR3BP quintic solver, co-rotating transform, Hill radius. +test-lagrange: test/test_lagrange.c src/lagrange.h + $(CC) -Wall -Werror -Isrc -o test/test_lagrange $< -lm + ./test/test_lagrange + +.PHONY: test-lagrange + # ── Standalone OD math unit test (no PostgreSQL dependency) ── # Element converters, inverse coordinate transforms, Brouwer/Kozai inverse. test-od-math: test/test_od_math.c src/od_math.c src/od_math.h diff --git a/src/lagrange.h b/src/lagrange.h new file mode 100644 index 0000000..efc1d3c --- /dev/null +++ b/src/lagrange.h @@ -0,0 +1,486 @@ +/* + * lagrange.h -- Circular restricted three-body problem (CR3BP) solver + * + * Computes the five Lagrange equilibrium points for any gravitational + * two-body system. The solver is pure C with no PostgreSQL dependency, + * no global state, and no memory allocation. + * + * The CR3BP uses the mass parameter mu = M_secondary / (M_primary + M_secondary). + * In the co-rotating frame normalized to unit separation, L1/L2/L3 lie + * on the x-axis and L4/L5 form equilateral triangles. + * + * L1/L2/L3 positions come from Newton-Raphson on the quintic + * equilibrium polynomial. L4/L5 are exact analytic. + * + * References: + * Szebehely V., "Theory of Orbits" (1967), Academic Press + * Murray & Dermott, "Solar System Dynamics" (1999), Cambridge + */ + +#ifndef PG_ORRERY_LAGRANGE_H +#define PG_ORRERY_LAGRANGE_H + +#include + +/* ── Lagrange point identifiers ────────────────────────── */ + +#define LAGRANGE_L1 1 +#define LAGRANGE_L2 2 +#define LAGRANGE_L3 3 +#define LAGRANGE_L4 4 +#define LAGRANGE_L5 5 + +/* ── Sun-planet mass ratios ────────────────────────────── */ + +/* + * GM_sun / GM_planet ratios. Convert to CR3BP mu via: + * mu = 1.0 / (1.0 + ratio) + * + * Sources: IAU 2012 nominal masses, JPL DE441 constants. + * The Earth ratio includes the Moon (Earth+Moon system barycenter). + */ +#define SUN_MERCURY_RATIO 6023682.155 +#define SUN_VENUS_RATIO 408523.7187 +#define SUN_EARTH_RATIO 332946.0487 /* Earth+Moon system */ +#define SUN_MARS_RATIO 3098703.59 +#define SUN_JUPITER_RATIO 1047.348644 +#define SUN_SATURN_RATIO 3497.9018 +#define SUN_URANUS_RATIO 22902.98 +#define SUN_NEPTUNE_RATIO 19412.26 + +/* ── Earth-Moon mass ratio ─────────────────────────────── */ + +/* + * M_earth / M_moon. From DE441 EMRAT constant. + * mu = 1.0 / (1.0 + EARTH_MOON_EMRAT) + */ +#define EARTH_MOON_EMRAT 81.300568 + +/* ── Planet-moon GM ratios ─────────────────────────────── */ + +/* + * GM_planet / GM_moon from spacecraft-derived values. + * mu = 1.0 / (1.0 + ratio) + * + * Galilean moons (Schubert et al. 2004, Anderson et al. 1996-2001): + */ +#define JUPITER_IO_RATIO 22423.9 /* GM_Jup / GM_Io */ +#define JUPITER_EUROPA_RATIO 39478.0 /* GM_Jup / GM_Europa */ +#define JUPITER_GANYMEDE_RATIO 12716.0 /* GM_Jup / GM_Ganymede */ +#define JUPITER_CALLISTO_RATIO 17350.0 /* GM_Jup / GM_Callisto */ + +/* + * Saturn moons (Jacobson et al. 2006): + */ +#define SATURN_MIMAS_RATIO 15108611.0 +#define SATURN_ENCELADUS_RATIO 4955938.0 +#define SATURN_TETHYS_RATIO 6137851.0 +#define SATURN_DIONE_RATIO 3430825.0 +#define SATURN_RHEA_RATIO 1629997.0 +#define SATURN_TITAN_RATIO 4226.5 /* Titan is massive */ +#define SATURN_IAPETUS_RATIO 3148296.0 +#define SATURN_HYPERION_RATIO 6.821e9 /* tiny */ + +/* + * Uranus moons (Jacobson et al. 1992): + */ +#define URANUS_MIRANDA_RATIO 1311870.0 +#define URANUS_ARIEL_RATIO 65229.0 +#define URANUS_UMBRIEL_RATIO 72449.0 +#define URANUS_TITANIA_RATIO 24399.0 +#define URANUS_OBERON_RATIO 25399.0 + +/* + * Mars moons (Jacobson 2014): + */ +#define MARS_PHOBOS_RATIO 5.8775e7 +#define MARS_DEIMOS_RATIO 3.919e8 + +/* ── Maximum Newton-Raphson iterations ─────────────────── */ + +#define LAGRANGE_MAX_ITER 50 + +/* ── Core API ──────────────────────────────────────────── */ + +/* + * Solve for a Lagrange point in the normalized co-rotating frame. + * + * mu: mass parameter = M2 / (M1 + M2), must be in (0, 0.5] + * point_id: LAGRANGE_L1 through LAGRANGE_L5 + * x, y: output co-rotating coordinates (normalized to unit separation) + * Origin at barycenter. Primary at (-mu, 0), secondary at (1-mu, 0). + * + * Returns 0 on success, -1 on invalid input or convergence failure. + */ +static inline int +lagrange_corotating(double mu, int point_id, double *x, double *y) +{ + double gamma, f, fp, gamma_new; + int i; + + if (mu <= 0.0 || mu > 0.5 || point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + return -1; + + switch (point_id) + { + case LAGRANGE_L1: + /* + * L1: between primary and secondary. + * Solve: gamma^5 - (3-mu)*gamma^4 + (3-2*mu)*gamma^3 + * - mu*gamma^2 + 2*mu*gamma - mu = 0 + * where gamma = distance from secondary toward primary. + * Initial guess: Hill sphere approximation. + */ + gamma = cbrt(mu / 3.0); + for (i = 0; i < LAGRANGE_MAX_ITER; i++) + { + double g2 = gamma * gamma; + double g3 = g2 * gamma; + double g4 = g3 * gamma; + double g5 = g4 * gamma; + + f = g5 - (3.0 - mu) * g4 + (3.0 - 2.0 * mu) * g3 + - mu * g2 + 2.0 * mu * gamma - mu; + fp = 5.0 * g4 - 4.0 * (3.0 - mu) * g3 + + 3.0 * (3.0 - 2.0 * mu) * g2 + - 2.0 * mu * gamma + 2.0 * mu; + + if (fabs(fp) < 1e-30) + return -1; + + gamma_new = gamma - f / fp; + if (fabs(gamma_new - gamma) < 1e-15) + break; + gamma = gamma_new; + } + if (i == LAGRANGE_MAX_ITER) + return -1; + + *x = 1.0 - mu - gamma; + *y = 0.0; + break; + + case LAGRANGE_L2: + /* + * L2: beyond secondary, away from primary. + * Solve: gamma^5 + (3-mu)*gamma^4 + (3-2*mu)*gamma^3 + * - mu*gamma^2 - 2*mu*gamma - mu = 0 + */ + gamma = cbrt(mu / 3.0); + for (i = 0; i < LAGRANGE_MAX_ITER; i++) + { + double g2 = gamma * gamma; + double g3 = g2 * gamma; + double g4 = g3 * gamma; + double g5 = g4 * gamma; + + f = g5 + (3.0 - mu) * g4 + (3.0 - 2.0 * mu) * g3 + - mu * g2 - 2.0 * mu * gamma - mu; + fp = 5.0 * g4 + 4.0 * (3.0 - mu) * g3 + + 3.0 * (3.0 - 2.0 * mu) * g2 + - 2.0 * mu * gamma - 2.0 * mu; + + if (fabs(fp) < 1e-30) + return -1; + + gamma_new = gamma - f / fp; + if (fabs(gamma_new - gamma) < 1e-15) + break; + gamma = gamma_new; + } + if (i == LAGRANGE_MAX_ITER) + return -1; + + *x = 1.0 - mu + gamma; + *y = 0.0; + break; + + case LAGRANGE_L3: + /* + * L3: opposite side from secondary, beyond primary. + * Solve: gamma^5 + (2+mu)*gamma^4 + (1+2*mu)*gamma^3 + * - (1-mu)*gamma^2 - 2*(1-mu)*gamma - (1-mu) = 0 + * where gamma = distance from primary. + */ + gamma = 1.0 - 7.0 * mu / 12.0; /* Szebehely approximation */ + for (i = 0; i < LAGRANGE_MAX_ITER; i++) + { + double g2 = gamma * gamma; + double g3 = g2 * gamma; + double g4 = g3 * gamma; + double g5 = g4 * gamma; + double one_minus_mu = 1.0 - mu; + + f = g5 + (2.0 + mu) * g4 + (1.0 + 2.0 * mu) * g3 + - one_minus_mu * g2 - 2.0 * one_minus_mu * gamma + - one_minus_mu; + fp = 5.0 * g4 + 4.0 * (2.0 + mu) * g3 + + 3.0 * (1.0 + 2.0 * mu) * g2 + - 2.0 * one_minus_mu * gamma + - 2.0 * one_minus_mu; + + if (fabs(fp) < 1e-30) + return -1; + + gamma_new = gamma - f / fp; + if (fabs(gamma_new - gamma) < 1e-15) + break; + gamma = gamma_new; + } + if (i == LAGRANGE_MAX_ITER) + return -1; + + *x = -mu - gamma; + *y = 0.0; + break; + + case LAGRANGE_L4: + /* Equilateral triangle, leading */ + *x = 0.5 - mu; + *y = sqrt(3.0) / 2.0; + break; + + case LAGRANGE_L5: + /* Equilateral triangle, trailing */ + *x = 0.5 - mu; + *y = -sqrt(3.0) / 2.0; + break; + + default: + return -1; + } + + return 0; +} + + +/* + * Transform a co-rotating Lagrange point to physical ecliptic J2000. + * + * The co-rotating frame has origin at the barycenter, x-axis along + * the primary→secondary direction, z-axis along the orbital angular + * momentum. We construct this frame from the instantaneous positions + * and velocity of the secondary relative to the primary. + * + * primary[3]: heliocentric position of primary (AU, ecliptic J2000) + * secondary[3]: heliocentric position of secondary (AU, ecliptic J2000) + * sec_vel[3]: velocity of secondary relative to primary (AU/day) + * mu: mass parameter M2/(M1+M2) + * point_id: LAGRANGE_L1..L5 + * result[3]: output heliocentric position (AU, ecliptic J2000) + * + * Returns 0 on success, -1 on failure. + */ +static inline int +lagrange_position(const double primary[3], const double secondary[3], + const double sec_vel[3], double mu, int point_id, + double result[3]) +{ + double d[3], sep, e_x[3], e_z[3], e_y[3]; + double hx, hy, hz, hmag; + double x_co, y_co; + int rc; + + /* Displacement: primary → secondary */ + d[0] = secondary[0] - primary[0]; + d[1] = secondary[1] - primary[1]; + d[2] = secondary[2] - primary[2]; + + sep = sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]); + if (sep < 1e-30) + return -1; + + /* Unit vector along primary→secondary */ + e_x[0] = d[0] / sep; + e_x[1] = d[1] / sep; + e_x[2] = d[2] / sep; + + /* Angular momentum direction: h = d x v */ + hx = d[1] * sec_vel[2] - d[2] * sec_vel[1]; + hy = d[2] * sec_vel[0] - d[0] * sec_vel[2]; + hz = d[0] * sec_vel[1] - d[1] * sec_vel[0]; + + hmag = sqrt(hx*hx + hy*hy + hz*hz); + if (hmag < 1e-30) + return -1; + + e_z[0] = hx / hmag; + e_z[1] = hy / hmag; + e_z[2] = hz / hmag; + + /* e_y = e_z x e_x (completes right-handed frame) */ + e_y[0] = e_z[1] * e_x[2] - e_z[2] * e_x[1]; + e_y[1] = e_z[2] * e_x[0] - e_z[0] * e_x[2]; + e_y[2] = e_z[0] * e_x[1] - e_z[1] * e_x[0]; + + /* Solve for co-rotating coordinates */ + rc = lagrange_corotating(mu, point_id, &x_co, &y_co); + if (rc != 0) + return -1; + + /* + * Physical position relative to barycenter: + * P_bary = primary + mu * d (barycenter location) + * L_phys = P_bary + sep * (x_co * e_x + y_co * e_y) + * + * But x_co is already relative to barycenter (origin in co-rotating + * frame), so: + * L_phys = primary + mu * d + sep * (x_co * e_x + y_co * e_y) + */ + result[0] = primary[0] + mu * d[0] + + sep * (x_co * e_x[0] + y_co * e_y[0]); + result[1] = primary[1] + mu * d[1] + + sep * (x_co * e_x[1] + y_co * e_y[1]); + result[2] = primary[2] + mu * d[2] + + sep * (x_co * e_x[2] + y_co * e_y[2]); + + return 0; +} + + +/* + * Hill sphere radius. + * + * separation_au: distance between primary and secondary (AU) + * mu: mass parameter M2/(M1+M2) + * + * Returns Hill radius in AU. + */ +static inline double +lagrange_hill_radius(double separation_au, double mu) +{ + return separation_au * cbrt(mu / 3.0); +} + + +/* + * Libration zone radius (approximate). + * + * For L1/L2: same as Hill radius (zone extends ~r_Hill from L-point). + * For L4/L5: horseshoe/tadpole width ~ separation * sqrt(mu) (Dermott 1981). + * For L3: ~ separation * (7/12) * mu (very narrow). + * + * separation_au: distance between primary and secondary (AU) + * mu: mass parameter + * point_id: LAGRANGE_L1..L5 + * + * Returns approximate zone radius in AU, or -1.0 on error. + */ +static inline double +lagrange_zone_radius(double separation_au, double mu, int point_id) +{ + switch (point_id) + { + case LAGRANGE_L1: + case LAGRANGE_L2: + return lagrange_hill_radius(separation_au, mu); + + case LAGRANGE_L3: + return separation_au * (7.0 / 12.0) * mu; + + case LAGRANGE_L4: + case LAGRANGE_L5: + return separation_au * sqrt(mu); + + default: + return -1.0; + } +} + + +/* + * Look up the Sun-planet mass ratio for a pg_orrery body_id. + * + * body_id: 1=Mercury..8=Neptune (pg_orrery convention) + * Returns the GM_sun/GM_planet ratio, or -1.0 for invalid body_id. + */ +static inline double +sun_planet_ratio(int body_id) +{ + switch (body_id) + { + case 1: return SUN_MERCURY_RATIO; + case 2: return SUN_VENUS_RATIO; + case 3: return SUN_EARTH_RATIO; + case 4: return SUN_MARS_RATIO; + case 5: return SUN_JUPITER_RATIO; + case 6: return SUN_SATURN_RATIO; + case 7: return SUN_URANUS_RATIO; + case 8: return SUN_NEPTUNE_RATIO; + default: return -1.0; + } +} + + +/* + * Compute mu from a Sun/planet GM ratio. + * mu = 1 / (1 + ratio) + */ +static inline double +mu_from_ratio(double ratio) +{ + return 1.0 / (1.0 + ratio); +} + + +/* + * Look up planet-moon GM ratio for a specific moon. + * + * family: 'g' (Galilean), 's' (Saturn), 'u' (Uranus), 'm' (Mars) + * moon_id: 0-based index within family + * Returns ratio, or -1.0 for invalid. + */ +static inline double +planet_moon_ratio(char family, int moon_id) +{ + switch (family) + { + case 'g': /* Galilean */ + switch (moon_id) + { + case 0: return JUPITER_IO_RATIO; + case 1: return JUPITER_EUROPA_RATIO; + case 2: return JUPITER_GANYMEDE_RATIO; + case 3: return JUPITER_CALLISTO_RATIO; + default: return -1.0; + } + + case 's': /* Saturn */ + switch (moon_id) + { + case 0: return SATURN_MIMAS_RATIO; + case 1: return SATURN_ENCELADUS_RATIO; + case 2: return SATURN_TETHYS_RATIO; + case 3: return SATURN_DIONE_RATIO; + case 4: return SATURN_RHEA_RATIO; + case 5: return SATURN_TITAN_RATIO; + case 6: return SATURN_IAPETUS_RATIO; + case 7: return SATURN_HYPERION_RATIO; + default: return -1.0; + } + + case 'u': /* Uranus */ + switch (moon_id) + { + case 0: return URANUS_MIRANDA_RATIO; + case 1: return URANUS_ARIEL_RATIO; + case 2: return URANUS_UMBRIEL_RATIO; + case 3: return URANUS_TITANIA_RATIO; + case 4: return URANUS_OBERON_RATIO; + default: return -1.0; + } + + case 'm': /* Mars */ + switch (moon_id) + { + case 0: return MARS_PHOBOS_RATIO; + case 1: return MARS_DEIMOS_RATIO; + default: return -1.0; + } + + default: + return -1.0; + } +} + +#endif /* PG_ORRERY_LAGRANGE_H */ diff --git a/test/test_lagrange.c b/test/test_lagrange.c new file mode 100644 index 0000000..8998aa8 --- /dev/null +++ b/test/test_lagrange.c @@ -0,0 +1,459 @@ +/* + * test_lagrange.c -- Standalone unit test for the Lagrange solver + * + * Verifies quintic solutions, L4/L5 geometry, Hill radius, + * zone radius, and co-rotating to physical frame transform. + * + * No PostgreSQL dependency. + * + * Build: cc -Wall -Werror -Isrc -o test/test_lagrange \ + * test/test_lagrange.c -lm + * Run: ./test/test_lagrange + */ + +#include "lagrange.h" + +#include +#include +#include + +/* ── 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) + +/* ── Tests ─────────────────────────────────────────────── */ + +/* + * Verify equilibrium: at a Lagrange point, the net force in the + * co-rotating frame should vanish. We check the effective potential + * gradient by evaluating the quintic polynomial. + */ +static void +test_equilibrium_check(double mu, int point_id, const char *label) +{ + double x, y; + int rc; + char buf[128]; + + rc = lagrange_corotating(mu, point_id, &x, &y); + snprintf(buf, sizeof(buf), "%s: convergence", label); + RUN(rc == 0, buf); + + if (rc != 0) + return; + + if (point_id <= LAGRANGE_L3) + { + /* + * For collinear points, verify equilibrium directly. + * At equilibrium on the x-axis: + * x - (1-mu)*(x+mu)/|x+mu|^3 - mu*(x-1+mu)/|x-1+mu|^3 = 0 + */ + double dx1 = x + mu; /* distance from primary (at -mu) */ + double dx2 = x - 1.0 + mu; /* distance from secondary (at 1-mu) */ + double r1 = fabs(dx1); + double r2 = fabs(dx2); + double residual; + + residual = x - (1.0 - mu) * dx1 / (r1 * r1 * r1) + - mu * dx2 / (r2 * r2 * r2); + + snprintf(buf, sizeof(buf), "%s: equilibrium residual", label); + CLOSE(residual, 0.0, 1e-12, buf); + } + else + { + /* L4/L5: equidistant from both primaries at unit distance */ + double r1 = sqrt((x + mu) * (x + mu) + y * y); + double r2 = sqrt((x - 1.0 + mu) * (x - 1.0 + mu) + y * y); + + snprintf(buf, sizeof(buf), "%s: distance to primary", label); + CLOSE(r1, 1.0, 1e-14, buf); + + snprintf(buf, sizeof(buf), "%s: distance to secondary", label); + CLOSE(r2, 1.0, 1e-14, buf); + } +} + + +static void +test_sun_earth(void) +{ + double mu = mu_from_ratio(SUN_EARTH_RATIO); + double x, y; + int rc; + + fprintf(stderr, "\n── Sun-Earth system (mu = %.6e) ──\n", mu); + + /* L1: between Sun and Earth, ~0.01 AU from Earth */ + rc = lagrange_corotating(mu, LAGRANGE_L1, &x, &y); + RUN(rc == 0, "L1 converges"); + /* L1 should be between barycenter and secondary */ + RUN(x > -mu && x < 1.0 - mu, "L1 between primaries"); + + /* Distance from secondary (Earth at 1-mu) */ + { + double d_from_earth = (1.0 - mu) - x; + CLOSE(d_from_earth, 0.01, 0.002, "L1 ~0.01 AU from Earth"); + } + + /* L2: beyond Earth, also ~0.01 AU */ + rc = lagrange_corotating(mu, LAGRANGE_L2, &x, &y); + RUN(rc == 0, "L2 converges"); + { + double d_from_earth = x - (1.0 - mu); + CLOSE(d_from_earth, 0.01, 0.002, "L2 ~0.01 AU from Earth"); + } + + /* L3: opposite side from Earth */ + rc = lagrange_corotating(mu, LAGRANGE_L3, &x, &y); + RUN(rc == 0, "L3 converges"); + RUN(x < -mu, "L3 beyond primary (opposite side)"); + + test_equilibrium_check(mu, LAGRANGE_L1, "Sun-Earth L1"); + test_equilibrium_check(mu, LAGRANGE_L2, "Sun-Earth L2"); + test_equilibrium_check(mu, LAGRANGE_L3, "Sun-Earth L3"); + test_equilibrium_check(mu, LAGRANGE_L4, "Sun-Earth L4"); + test_equilibrium_check(mu, LAGRANGE_L5, "Sun-Earth L5"); +} + + +static void +test_sun_jupiter(void) +{ + double mu = mu_from_ratio(SUN_JUPITER_RATIO); + double x, y; + int rc; + + fprintf(stderr, "\n── Sun-Jupiter system (mu = %.6e) ──\n", mu); + + /* L4/L5: should be at 60 degrees from Jupiter */ + rc = lagrange_corotating(mu, LAGRANGE_L4, &x, &y); + RUN(rc == 0, "L4 converges"); + { + /* Angle from secondary: atan2(y, x - (1-mu)) */ + double angle = atan2(y, x - (1.0 - mu)); + double angle_deg = angle * 180.0 / M_PI; + /* L4 leads secondary by ~60 degrees (but angle from barycenter) */ + /* Actually check equilateral property */ + double d_prim = sqrt((x + mu) * (x + mu) + y * y); + double d_sec = sqrt((x - 1.0 + mu) * (x - 1.0 + mu) + y * y); + CLOSE(d_prim, 1.0, 1e-14, "L4 unit distance from primary"); + CLOSE(d_sec, 1.0, 1e-14, "L4 unit distance from secondary"); + RUN(y > 0.0, "L4 above x-axis (leading)"); + (void)angle_deg; /* used implicitly via assertions */ + } + + test_equilibrium_check(mu, LAGRANGE_L1, "Sun-Jupiter L1"); + test_equilibrium_check(mu, LAGRANGE_L2, "Sun-Jupiter L2"); + test_equilibrium_check(mu, LAGRANGE_L3, "Sun-Jupiter L3"); + test_equilibrium_check(mu, LAGRANGE_L4, "Sun-Jupiter L4"); + test_equilibrium_check(mu, LAGRANGE_L5, "Sun-Jupiter L5"); +} + + +static void +test_earth_moon(void) +{ + double mu = mu_from_ratio(EARTH_MOON_EMRAT); + + fprintf(stderr, "\n── Earth-Moon system (mu = %.6e) ──\n", mu); + + test_equilibrium_check(mu, LAGRANGE_L1, "Earth-Moon L1"); + test_equilibrium_check(mu, LAGRANGE_L2, "Earth-Moon L2"); + test_equilibrium_check(mu, LAGRANGE_L3, "Earth-Moon L3"); + test_equilibrium_check(mu, LAGRANGE_L4, "Earth-Moon L4"); + test_equilibrium_check(mu, LAGRANGE_L5, "Earth-Moon L5"); + + /* Earth-Moon L1 should be ~326,000 km from Earth (~84.7% of separation) */ + { + double x, y; + int rc; + double earth_moon_km = 384400.0; /* mean distance */ + + rc = lagrange_corotating(mu, LAGRANGE_L1, &x, &y); + RUN(rc == 0, "Earth-Moon L1 converges"); + + /* In co-rotating frame, Earth is at -mu, Moon at 1-mu. + * L1 is between them. Distance from Earth = x + mu. */ + { + double frac = (x + mu); /* fraction of separation from Earth */ + double km_from_earth = frac * earth_moon_km; + CLOSE(km_from_earth, 326000.0, 5000.0, + "E-M L1 ~326,000 km from Earth"); + } + } +} + + +static void +test_l4_l5_symmetry(void) +{ + double mu = mu_from_ratio(SUN_JUPITER_RATIO); + double x4, y4, x5, y5; + int rc; + + fprintf(stderr, "\n── L4/L5 symmetry ──\n"); + + rc = lagrange_corotating(mu, LAGRANGE_L4, &x4, &y4); + RUN(rc == 0, "L4 converges"); + rc = lagrange_corotating(mu, LAGRANGE_L5, &x5, &y5); + RUN(rc == 0, "L5 converges"); + + CLOSE(x4, x5, 1e-15, "L4 and L5 same x-coordinate"); + CLOSE(y4, -y5, 1e-15, "L4 and L5 mirror in y"); +} + + +static void +test_l1_l2_ordering(void) +{ + double mu = mu_from_ratio(SUN_EARTH_RATIO); + double x1, y1, x2, y2, x3, y3; + int rc; + + fprintf(stderr, "\n── L1/L2/L3 ordering ──\n"); + + rc = lagrange_corotating(mu, LAGRANGE_L1, &x1, &y1); + RUN(rc == 0, "L1 converges"); + rc = lagrange_corotating(mu, LAGRANGE_L2, &x2, &y2); + RUN(rc == 0, "L2 converges"); + rc = lagrange_corotating(mu, LAGRANGE_L3, &x3, &y3); + RUN(rc == 0, "L3 converges"); + + /* Ordering: L3 < primary < L1 < secondary < L2 */ + RUN(x3 < -mu, "L3 < primary"); + RUN(x1 > -mu && x1 < 1.0 - mu, "L1 between primaries"); + RUN(x2 > 1.0 - mu, "L2 beyond secondary"); +} + + +static void +test_hill_radius(void) +{ + double mu_jup, mu_earth; + double hill_jup, hill_earth; + + fprintf(stderr, "\n── Hill radius ──\n"); + + mu_jup = mu_from_ratio(SUN_JUPITER_RATIO); + mu_earth = mu_from_ratio(SUN_EARTH_RATIO); + + /* Jupiter at ~5.2 AU */ + hill_jup = lagrange_hill_radius(5.2, mu_jup); + CLOSE(hill_jup, 0.355, 0.02, "Jupiter Hill radius ~0.35 AU"); + + /* Earth at ~1.0 AU */ + hill_earth = lagrange_hill_radius(1.0, mu_earth); + CLOSE(hill_earth, 0.01, 0.002, "Earth Hill radius ~0.01 AU"); +} + + +static void +test_zone_radius(void) +{ + double mu = mu_from_ratio(SUN_JUPITER_RATIO); + double zr; + + fprintf(stderr, "\n── Zone radius ──\n"); + + zr = lagrange_zone_radius(5.2, mu, LAGRANGE_L1); + RUN(zr > 0.0, "L1 zone radius positive"); + + zr = lagrange_zone_radius(5.2, mu, LAGRANGE_L4); + RUN(zr > 0.0, "L4 zone radius positive"); + + zr = lagrange_zone_radius(5.2, mu, 99); + RUN(zr < 0.0, "invalid point_id returns -1"); +} + + +static void +test_physical_transform(void) +{ + double primary[3] = {0.0, 0.0, 0.0}; /* Sun at origin */ + double secondary[3] = {1.0, 0.0, 0.0}; /* "planet" at 1 AU on x-axis */ + double sec_vel[3] = {0.0, 0.01720209895, 0.0}; /* ~Gauss constant, circular */ + double mu = 0.001; /* ~Jupiter-like */ + double result[3]; + int rc; + + fprintf(stderr, "\n── Physical frame transform ──\n"); + + /* L1: should be between Sun and planet, on x-axis */ + rc = lagrange_position(primary, secondary, sec_vel, mu, LAGRANGE_L1, result); + RUN(rc == 0, "L1 transform succeeds"); + RUN(result[0] > 0.0 && result[0] < 1.0, "L1 between Sun and planet on x-axis"); + CLOSE(result[1], 0.0, 1e-10, "L1 y-component ~0"); + CLOSE(result[2], 0.0, 1e-10, "L1 z-component ~0"); + + /* L2: beyond planet */ + rc = lagrange_position(primary, secondary, sec_vel, mu, LAGRANGE_L2, result); + RUN(rc == 0, "L2 transform succeeds"); + RUN(result[0] > 1.0, "L2 beyond planet"); + CLOSE(result[1], 0.0, 1e-10, "L2 y-component ~0"); + + /* L4: 60 degrees ahead, above x-axis in ecliptic plane */ + rc = lagrange_position(primary, secondary, sec_vel, mu, LAGRANGE_L4, result); + RUN(rc == 0, "L4 transform succeeds"); + { + double dist = sqrt(result[0]*result[0] + result[1]*result[1] + result[2]*result[2]); + /* L4 should be ~1 AU from Sun (equilateral triangle) */ + CLOSE(dist, 1.0, 0.01, "L4 ~1 AU from Sun"); + RUN(result[1] > 0.0, "L4 positive y (leading)"); + } + + /* L5: symmetric with L4 */ + rc = lagrange_position(primary, secondary, sec_vel, mu, LAGRANGE_L5, result); + RUN(rc == 0, "L5 transform succeeds"); + RUN(result[1] < 0.0, "L5 negative y (trailing)"); +} + + +static void +test_extreme_mass_ratios(void) +{ + double x, y; + int rc; + + fprintf(stderr, "\n── Extreme mass ratios ──\n"); + + /* Very small mu (like Mercury around the Sun) */ + { + double mu = mu_from_ratio(SUN_MERCURY_RATIO); /* ~1.66e-7 */ + rc = lagrange_corotating(mu, LAGRANGE_L1, &x, &y); + RUN(rc == 0, "tiny mu L1 converges"); + test_equilibrium_check(mu, LAGRANGE_L1, "Mercury L1"); + } + + /* Moderately large mu */ + { + double mu = 0.1; + rc = lagrange_corotating(mu, LAGRANGE_L1, &x, &y); + RUN(rc == 0, "mu=0.1 L1 converges"); + test_equilibrium_check(mu, LAGRANGE_L1, "mu=0.1 L1"); + test_equilibrium_check(mu, LAGRANGE_L2, "mu=0.1 L2"); + test_equilibrium_check(mu, LAGRANGE_L3, "mu=0.1 L3"); + } + + /* Equal mass (mu = 0.5, maximum) */ + { + double mu = 0.5; + rc = lagrange_corotating(mu, LAGRANGE_L1, &x, &y); + RUN(rc == 0, "mu=0.5 L1 converges"); + test_equilibrium_check(mu, LAGRANGE_L1, "mu=0.5 L1"); + test_equilibrium_check(mu, LAGRANGE_L2, "mu=0.5 L2"); + test_equilibrium_check(mu, LAGRANGE_L3, "mu=0.5 L3"); + /* L4/L5 at (0, +-sqrt(3)/2) for equal mass */ + rc = lagrange_corotating(mu, LAGRANGE_L4, &x, &y); + RUN(rc == 0, "mu=0.5 L4 converges"); + CLOSE(x, 0.0, 1e-15, "mu=0.5 L4 x=0"); + CLOSE(y, sqrt(3.0)/2.0, 1e-15, "mu=0.5 L4 y=sqrt(3)/2"); + } +} + + +static void +test_error_cases(void) +{ + double x, y; + int rc; + + fprintf(stderr, "\n── Error cases ──\n"); + + rc = lagrange_corotating(0.0, LAGRANGE_L1, &x, &y); + RUN(rc != 0, "mu=0 rejected"); + + rc = lagrange_corotating(-0.1, LAGRANGE_L1, &x, &y); + RUN(rc != 0, "negative mu rejected"); + + rc = lagrange_corotating(0.6, LAGRANGE_L1, &x, &y); + RUN(rc != 0, "mu>0.5 rejected"); + + rc = lagrange_corotating(0.01, 0, &x, &y); + RUN(rc != 0, "point_id=0 rejected"); + + rc = lagrange_corotating(0.01, 6, &x, &y); + RUN(rc != 0, "point_id=6 rejected"); + + /* Mass ratio lookups */ + RUN(sun_planet_ratio(1) > 0.0, "Mercury ratio valid"); + RUN(sun_planet_ratio(8) > 0.0, "Neptune ratio valid"); + RUN(sun_planet_ratio(0) < 0.0, "Sun ratio invalid"); + RUN(sun_planet_ratio(9) < 0.0, "body 9 invalid"); + + RUN(planet_moon_ratio('g', 0) > 0.0, "Io ratio valid"); + RUN(planet_moon_ratio('g', 4) < 0.0, "Galilean moon 4 invalid"); + RUN(planet_moon_ratio('s', 7) > 0.0, "Hyperion ratio valid"); + RUN(planet_moon_ratio('s', 8) < 0.0, "Saturn moon 8 invalid"); + RUN(planet_moon_ratio('u', 4) > 0.0, "Oberon ratio valid"); + RUN(planet_moon_ratio('u', 5) < 0.0, "Uranus moon 5 invalid"); + RUN(planet_moon_ratio('m', 1) > 0.0, "Deimos ratio valid"); + RUN(planet_moon_ratio('m', 2) < 0.0, "Mars moon 2 invalid"); + RUN(planet_moon_ratio('x', 0) < 0.0, "unknown family invalid"); +} + + +static void +test_all_planets(void) +{ + int body; + + fprintf(stderr, "\n── All planets equilibrium ──\n"); + + for (body = 1; body <= 8; body++) + { + double ratio = sun_planet_ratio(body); + double mu = mu_from_ratio(ratio); + char label[64]; + int pt; + + for (pt = LAGRANGE_L1; pt <= LAGRANGE_L5; pt++) + { + snprintf(label, sizeof(label), "body %d L%d", body, pt); + test_equilibrium_check(mu, pt, label); + } + } +} + + +/* ── Main ──────────────────────────────────────────────── */ + +int +main(void) +{ + fprintf(stderr, "Lagrange solver unit test\n"); + fprintf(stderr, "========================\n"); + + test_sun_earth(); + test_sun_jupiter(); + test_earth_moon(); + test_l4_l5_symmetry(); + test_l1_l2_ordering(); + test_hill_radius(); + test_zone_radius(); + test_physical_transform(); + test_extreme_mass_ratios(); + test_error_cases(); + test_all_planets(); + + fprintf(stderr, "\n%d/%d tests passed\n", n_pass, n_run); + + return (n_pass == n_run) ? 0 : 1; +} From dbc1f20a4607dac667a8c26b6e33f7f0822c2602 Mon Sep 17 00:00:00 2001 From: Ryan Malloy Date: Sat, 28 Feb 2026 14:21:28 -0700 Subject: [PATCH 2/2] Add v0.20.0: Lagrange point SQL functions, DE variants, regression tests MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 37 new SQL objects (188 → 225 total): - Sun-planet L1-L5: heliocentric, observe, equatorial, distance (5 IMMUTABLE) - Earth-Moon L1-L5: observe, equatorial via ELP2000-82B (2 IMMUTABLE) - Planetary moon L1-L5: Galilean/Saturn/Uranus/Mars families (8 IMMUTABLE) - Hill radius, zone radius, mass ratio, point name (5 IMMUTABLE) - DE variants with VSOP87/ELP2000-82B fallback (17 STABLE) All 31 regression tests pass. 210/210 standalone math tests pass. --- .gitignore | 1 + CLAUDE.md | 24 +- Makefile | 9 +- pg_orrery.control | 2 +- sql/pg_orrery--0.19.0--0.20.0.sql | 244 ++++ sql/pg_orrery--0.20.0.sql | 2202 +++++++++++++++++++++++++++++ src/de_funcs.c | 968 +++++++++++++ src/lagrange_funcs.c | 916 ++++++++++++ test/expected/v020_features.out | 323 +++++ test/sql/v020_features.sql | 209 +++ 10 files changed, 4885 insertions(+), 13 deletions(-) create mode 100644 sql/pg_orrery--0.19.0--0.20.0.sql create mode 100644 sql/pg_orrery--0.20.0.sql create mode 100644 src/lagrange_funcs.c create mode 100644 test/expected/v020_features.out create mode 100644 test/sql/v020_features.sql diff --git a/.gitignore b/.gitignore index e83e025..7a0bc18 100644 --- a/.gitignore +++ b/.gitignore @@ -24,6 +24,7 @@ test/test_de_reader test/test_od_math test/test_od_iod test/test_od_gauss +test/test_lagrange # Bench — downloaded TLE catalogs (large, ephemeral) # Already-tracked files (active.tle, spacetrack_full*.tle) are unaffected. diff --git a/CLAUDE.md b/CLAUDE.md index 29f42ca..bf607f5 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -1,9 +1,9 @@ # pg_orrery — A Database Orrery for PostgreSQL ## What This Is -A database orrery — celestial mechanics types and functions for PostgreSQL. Native C extension using PGXS, 188 SQL objects (172 user-visible functions + 16 GiST support), 9 custom types, covering satellites (SGP4/SDP4), planets (VSOP87 + optional JPL DE441), Moon (ELP2000-82B), 19 planetary moons (L1.2/TASS17/GUST86/MarsSat), stars (with proper motion and annual parallax), comets, asteroids (MPC catalog), Jupiter radio bursts, interplanetary Lambert transfers, equatorial RA/Dec coordinates with GiST-indexed angular separation, atmospheric refraction, annual stellar aberration, light-time correction, rise/set prediction (geometric + refracted + event windows + sun almanac) with status diagnostics, IAU constellation identification with full name lookup (Roman 1987), twilight dawn/dusk (civil/nautical/astronomical), lunar phase (angle, illumination, name, age), planet apparent magnitude with Saturn ring correction (Mallama & Hilton 2018), solar elongation, planet phase fraction, conjunction detection, satellite eclipse prediction (conical shadow with penumbral fraction), observing night quality assessment, lunar libration (optical + physical, Meeus Ch. 53 + p. 373), and angular separation rate. +A database orrery — celestial mechanics types and functions for PostgreSQL. Native C extension using PGXS, 225 SQL objects (209 user-visible functions + 16 GiST support), 9 custom types, covering satellites (SGP4/SDP4), planets (VSOP87 + optional JPL DE441), Moon (ELP2000-82B), 19 planetary moons (L1.2/TASS17/GUST86/MarsSat), stars (with proper motion and annual parallax), comets, asteroids (MPC catalog), Jupiter radio bursts, interplanetary Lambert transfers, equatorial RA/Dec coordinates with GiST-indexed angular separation, atmospheric refraction, annual stellar aberration, light-time correction, rise/set prediction (geometric + refracted + event windows + sun almanac) with status diagnostics, IAU constellation identification with full name lookup (Roman 1987), twilight dawn/dusk (civil/nautical/astronomical), lunar phase (angle, illumination, name, age), planet apparent magnitude with Saturn ring correction (Mallama & Hilton 2018), solar elongation, planet phase fraction, conjunction detection, satellite eclipse prediction (conical shadow with penumbral fraction), observing night quality assessment, lunar libration (optical + physical, Meeus Ch. 53 + p. 373), angular separation rate, and Lagrange point equilibrium positions (CR3BP L1-L5 for Sun-planet, Earth-Moon, and 19 planetary moon systems). -**Current version:** 0.19.0 +**Current version:** 0.20.0 **Repository:** https://git.supported.systems/warehack.ing/pg_orrery **Documentation:** https://pg-orrery.warehack.ing @@ -11,7 +11,7 @@ A database orrery — celestial mechanics types and functions for PostgreSQL. Na ```bash make PG_CONFIG=/usr/bin/pg_config # Compile with PGXS sudo make install PG_CONFIG=/usr/bin/pg_config # Install extension -make installcheck PG_CONFIG=/usr/bin/pg_config # Run 30 regression test suites +make installcheck PG_CONFIG=/usr/bin/pg_config # Run 31 regression test suites ``` Requires: PostgreSQL 17 development headers, GCC, Make. @@ -27,7 +27,7 @@ Image: `git.supported.systems/warehack.ing/pg_orrery:pg17` ## Project Layout ``` -pg_orrery.control # Extension metadata (version 0.19.0) +pg_orrery.control # Extension metadata (version 0.20.0) Makefile # PGXS build + Docker targets sql/ pg_orrery--0.1.0.sql # v0.1.0: satellite types/functions/operators @@ -49,6 +49,7 @@ sql/ pg_orrery--0.17.0.sql # v0.17.0: elongation, phase, eclipse, night quality, libration (174 objects) pg_orrery--0.18.0.sql # v0.18.0: ring tilt, penumbral eclipse, rise/set windows, angular rate (184 objects) pg_orrery--0.19.0.sql # v0.19.0: sun almanac, conjunctions, penumbral fraction, physical libration (188 objects) + pg_orrery--0.20.0.sql # v0.20.0: Lagrange point equilibrium positions (225 objects) pg_orrery--0.1.0--0.2.0.sql # Migration: v0.1.0 → v0.2.0 (adds solar system) pg_orrery--0.2.0--0.3.0.sql # Migration: v0.2.0 → v0.3.0 (adds DE ephemeris) pg_orrery--0.3.0--0.4.0.sql # Migration: v0.3.0 → v0.4.0 @@ -67,6 +68,7 @@ sql/ pg_orrery--0.16.0--0.17.0.sql # Migration: v0.16.0 → v0.17.0 (elongation, phase, eclipse, night quality, libration) pg_orrery--0.17.0--0.18.0.sql # Migration: v0.17.0 → v0.18.0 (ring tilt, penumbral eclipse, rise/set windows, angular rate) pg_orrery--0.18.0--0.19.0.sql # Migration: v0.18.0 → v0.19.0 (sun almanac, conjunctions, penumbral fraction, physical libration) + pg_orrery--0.19.0--0.20.0.sql # Migration: v0.19.0 → v0.20.0 (Lagrange points) src/ pg_orrery.c # PG_MODULE_MAGIC + _PG_init() (GUC registration) types.h # All struct definitions + constants + DE body ID mapping @@ -100,6 +102,8 @@ src/ magnitude_funcs.c # planet_magnitude() (with Saturn ring correction), solar_elongation(), planet_phase(), saturn_ring_tilt() eclipse_funcs.c # satellite eclipse prediction (conical shadow with penumbral fraction, Vallado §5.3) libration.h / libration_funcs.c # lunar libration (optical Meeus Ch. 53 + physical p. 373) + lagrange.h # CR3BP solver (header-only): quintic solver, co-rotating frame, Hill radius + lagrange_funcs.c # Lagrange point SQL functions (Sun-planet, Earth-Moon, planetary moons) l12.c / l12.h # L1.2 Galilean moon theory (Lieske 1998) tass17.c / tass17.h # TASS 1.7 Saturn moon theory (Vienne & Duriez 1995) gust86.c / gust86.h # GUST86 Uranus moon theory (Laskar & Jacobson 1987) @@ -124,7 +128,7 @@ src/ PROVENANCE.md # Vendoring decision, modifications, verification LICENSE # MIT license (Bill Gray / Project Pluto) test/ - sql/ # 30 regression test suites + sql/ # 31 regression test suites expected/ # Expected output data/vallado_518.json # 518 Vallado test vectors (AIAA 2006-6753-Rev1) docs/ @@ -151,7 +155,7 @@ All types are fixed-size, `STORAGE = plain`, `ALIGNMENT = double`. No TOAST over | `orbital_elements` | 72 | Classical Keplerian elements for comets/asteroids (epoch, q, e, inc, omega, Omega, tp, H, G) | | `equatorial` | 24 | Apparent RA (hours), Dec (degrees), distance (km) — of date | -## Function Domains (188 SQL objects) +## Function Domains (225 SQL objects) | Domain | Theory | Key Functions | Count | |--------|--------|---------------|-------| @@ -179,6 +183,7 @@ All types are fixed-size, `STORAGE = plain`, `ALIGNMENT = double`. No TOAST over | Observing quality | Composite (twilight+Moon) | `observing_night_quality()` | 1 | | Lunar libration | Meeus (1998) Ch. 53 + p. 373 | `moon_libration_longitude()`, `moon_libration()`, `moon_subsolar_longitude()`, `moon_physical_libration()` | 6 | | Constellation | Roman (1987) CDS VI/42 | `constellation()`, `constellation_full_name()` | 3 | +| Lagrange points | CR3BP quintic + VSOP87 | `lagrange_heliocentric()`, `lagrange_observe()`, `hill_radius()` | 37 | | Diagnostics | -- | `pg_orrery_ephemeris_info()` | 1 | All functions are `PARALLEL SAFE`. VSOP87/ELP82B functions are `IMMUTABLE` (compiled-in coefficients). DE functions are `STABLE` (external file dependency). @@ -312,7 +317,7 @@ All numerical logic is byte-identical to upstream. Verified against 518 Vallado ## Testing -30 regression test suites via `make installcheck`: +31 regression test suites via `make installcheck`: | Suite | What it tests | |-------|--------------| @@ -346,10 +351,11 @@ All numerical logic is byte-identical to upstream. Verified against 518 Vallado | v017_features | Solar elongation ranges/errors, planet phase ranges, satellite eclipse, observing night quality, lunar libration ranges, subsolar longitude | | v018_features | Saturn ring tilt range/variation, penumbral eclipse (shadow state, penumbra precedes umbra), rise/set event windows (Sun/Moon/planet, refracted vs geometric), angular separation rate (generic + planet convenience) | | v019_features | Sun almanac events (count/order/types/polar/refraction/window guard), conjunction detection (Jupiter-Saturn 2020, Moon-Venus, same-body error, threshold filter), penumbral fraction (range/bounds/eclipse consistency), physical libration (small corrections, time variation, total libration range) | +| v020_features | Lagrange L1-L5 heliocentric/observe/equatorial, Hill radius, zone radius, mass ratio, DE fallback, all planet + moon families, input validation | ### PG Version Matrix -Test all 30 regression suites + DE reader unit test across PostgreSQL 14-18 using Docker: +Test all 31 regression suites + DE reader unit test across PostgreSQL 14-18 using Docker: ```bash make test-matrix # Full matrix (PG 14-18) @@ -375,7 +381,7 @@ Logs saved to `test/matrix-logs/pg${ver}.log`. The script reuses the Dockerfile Starlight docs at `docs/` — 44+ MDX pages covering all domains. -Sections: Getting Started, Guides (9 domain walkthroughs incl. DE ephemeris), Workflow Translation (Skyfield/Horizons/GMAT/Radio Jupiter Pro comparisons), Reference (all 188 SQL objects incl. DE variants, equatorial GiST, refraction, rise/set, sun almanac, conjunction detection, constellation, twilight, lunar phase, planet magnitude, Saturn ring tilt, solar elongation, planet phase, satellite eclipse with penumbral fraction, observing quality, lunar libration with physical corrections, angular separation rate), Architecture (Hamilton's principles, constant custody, observation pipeline), Performance (benchmarks). +Sections: Getting Started, Guides (9 domain walkthroughs incl. DE ephemeris), Workflow Translation (Skyfield/Horizons/GMAT/Radio Jupiter Pro comparisons), Reference (all 225 SQL objects incl. DE variants, equatorial GiST, refraction, rise/set, sun almanac, conjunction detection, constellation, Lagrange points, twilight, lunar phase, planet magnitude, Saturn ring tilt, solar elongation, planet phase, satellite eclipse with penumbral fraction, observing quality, lunar libration with physical corrections, angular separation rate), Architecture (Hamilton's principles, constant custody, observation pipeline), Performance (benchmarks). ### Local Development ```bash diff --git a/Makefile b/Makefile index e600459..35206de 100644 --- a/Makefile +++ b/Makefile @@ -17,7 +17,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.16.0.sql sql/pg_orrery--0.15.0--0.16.0.sql \ sql/pg_orrery--0.17.0.sql sql/pg_orrery--0.16.0--0.17.0.sql \ sql/pg_orrery--0.18.0.sql sql/pg_orrery--0.17.0--0.18.0.sql \ - sql/pg_orrery--0.19.0.sql sql/pg_orrery--0.18.0--0.19.0.sql + sql/pg_orrery--0.19.0.sql sql/pg_orrery--0.18.0--0.19.0.sql \ + sql/pg_orrery--0.20.0.sql sql/pg_orrery--0.19.0--0.20.0.sql # Our extension C sources OBJS = src/pg_orrery.o src/tle_type.o src/eci_type.o src/observer_type.o \ @@ -38,7 +39,8 @@ OBJS = src/pg_orrery.o src/tle_type.o src/eci_type.o src/observer_type.o \ src/rise_set_funcs.o \ src/constellation_data.o src/constellation_funcs.o \ src/lunar_phase_funcs.o src/magnitude_funcs.o \ - src/eclipse_funcs.o src/libration_funcs.o + src/eclipse_funcs.o src/libration_funcs.o \ + src/lagrange_funcs.o # Vendored SGP4/SDP4 sources (pure C, from Bill Gray's sat_code, MIT license) SGP4_DIR = src/sgp4 @@ -62,7 +64,8 @@ REGRESS = tle_parse sgp4_propagate coord_transforms pass_prediction gist_index c v016_features \ v017_features \ v018_features \ - v019_features + v019_features \ + v020_features REGRESS_OPTS = --inputdir=test # Pure C — no C++ runtime needed. LAPACK for OD solver (dgelss_). diff --git a/pg_orrery.control b/pg_orrery.control index d6aed22..4060903 100644 --- 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.19.0' +default_version = '0.20.0' module_pathname = '$libdir/pg_orrery' relocatable = true diff --git a/sql/pg_orrery--0.19.0--0.20.0.sql b/sql/pg_orrery--0.19.0--0.20.0.sql new file mode 100644 index 0000000..1b69487 --- /dev/null +++ b/sql/pg_orrery--0.19.0--0.20.0.sql @@ -0,0 +1,244 @@ +-- pg_orrery 0.19.0 -> 0.20.0: Lagrange point support +-- CR3BP equilibrium positions for Sun-planet, Earth-Moon, and planetary moon systems. + +-- ============================================================ +-- Sun-planet Lagrange functions (5) +-- ============================================================ + +CREATE FUNCTION lagrange_heliocentric(int4, int4, timestamptz) RETURNS heliocentric + AS 'MODULE_PATHNAME', 'lagrange_heliocentric' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_heliocentric(int4, int4, timestamptz) IS + 'Heliocentric ecliptic J2000 position of a Sun-planet Lagrange point. body_id: 1-8 (Mercury-Neptune), point_id: 1-5 (L1-L5).'; + +CREATE FUNCTION lagrange_observe(int4, int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'lagrange_observe' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_observe(int4, int4, observer, timestamptz) IS + 'Observe a Sun-planet Lagrange point from a ground station. body_id: 1-8, point_id: 1-5.'; + +CREATE FUNCTION lagrange_equatorial(int4, int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'lagrange_equatorial' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_equatorial(int4, int4, timestamptz) IS + 'Geocentric RA/Dec of a Sun-planet Lagrange point. body_id: 1-8, point_id: 1-5.'; + +CREATE FUNCTION lagrange_distance(int4, int4, heliocentric, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'lagrange_distance' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_distance(int4, int4, heliocentric, timestamptz) IS + 'Distance (AU) from a heliocentric position to a Sun-planet Lagrange point.'; + +CREATE FUNCTION lagrange_distance_oe(int4, int4, orbital_elements, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'lagrange_distance_oe' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_distance_oe(int4, int4, orbital_elements, timestamptz) IS + 'Distance (AU) from an asteroid/comet (orbital_elements) to a Sun-planet Lagrange point.'; + +-- ============================================================ +-- Earth-Moon Lagrange functions (2) +-- ============================================================ + +CREATE FUNCTION lunar_lagrange_observe(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'lunar_lagrange_observe' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lunar_lagrange_observe(int4, observer, timestamptz) IS + 'Observe an Earth-Moon Lagrange point. point_id: 1-5 (L1-L5).'; + +CREATE FUNCTION lunar_lagrange_equatorial(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'lunar_lagrange_equatorial' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lunar_lagrange_equatorial(int4, timestamptz) IS + 'Geocentric RA/Dec of an Earth-Moon Lagrange point. point_id: 1-5 (L1-L5).'; + +-- ============================================================ +-- Planetary moon Lagrange functions (8) +-- ============================================================ + +CREATE FUNCTION galilean_lagrange_observe(int4, int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'galilean_lagrange_observe' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION galilean_lagrange_observe(int4, int4, observer, timestamptz) IS + 'Observe a Jupiter-Galilean moon Lagrange point. body_id: 0-3 (Io-Callisto), point_id: 1-5.'; + +CREATE FUNCTION galilean_lagrange_equatorial(int4, int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'galilean_lagrange_equatorial' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION galilean_lagrange_equatorial(int4, int4, timestamptz) IS + 'Geocentric RA/Dec of a Jupiter-Galilean moon Lagrange point. body_id: 0-3, point_id: 1-5.'; + +CREATE FUNCTION saturn_moon_lagrange_observe(int4, int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'saturn_moon_lagrange_observe' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION saturn_moon_lagrange_observe(int4, int4, observer, timestamptz) IS + 'Observe a Saturn moon Lagrange point. body_id: 0-7 (Mimas-Hyperion), point_id: 1-5.'; + +CREATE FUNCTION saturn_moon_lagrange_equatorial(int4, int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'saturn_moon_lagrange_equatorial' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION saturn_moon_lagrange_equatorial(int4, int4, timestamptz) IS + 'Geocentric RA/Dec of a Saturn moon Lagrange point. body_id: 0-7, point_id: 1-5.'; + +CREATE FUNCTION uranus_moon_lagrange_observe(int4, int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'uranus_moon_lagrange_observe' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION uranus_moon_lagrange_observe(int4, int4, observer, timestamptz) IS + 'Observe a Uranus moon Lagrange point. body_id: 0-4 (Miranda-Oberon), point_id: 1-5.'; + +CREATE FUNCTION uranus_moon_lagrange_equatorial(int4, int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'uranus_moon_lagrange_equatorial' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION uranus_moon_lagrange_equatorial(int4, int4, timestamptz) IS + 'Geocentric RA/Dec of a Uranus moon Lagrange point. body_id: 0-4, point_id: 1-5.'; + +CREATE FUNCTION mars_moon_lagrange_observe(int4, int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'mars_moon_lagrange_observe' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION mars_moon_lagrange_observe(int4, int4, observer, timestamptz) IS + 'Observe a Mars moon Lagrange point. body_id: 0-1 (Phobos-Deimos), point_id: 1-5.'; + +CREATE FUNCTION mars_moon_lagrange_equatorial(int4, int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'mars_moon_lagrange_equatorial' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION mars_moon_lagrange_equatorial(int4, int4, timestamptz) IS + 'Geocentric RA/Dec of a Mars moon Lagrange point. body_id: 0-1, point_id: 1-5.'; + +-- ============================================================ +-- Hill radius / zone / convenience (5) +-- ============================================================ + +CREATE FUNCTION hill_radius(int4, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'hill_radius_func' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION hill_radius(int4, timestamptz) IS + 'Hill sphere radius (AU) for a Sun-planet system. body_id: 1-8.'; + +CREATE FUNCTION hill_radius_lunar(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'hill_radius_lunar' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION hill_radius_lunar(timestamptz) IS + 'Hill sphere radius (AU) for the Earth-Moon system.'; + +CREATE FUNCTION lagrange_zone_radius(int4, int4, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'lagrange_zone_radius_func' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_zone_radius(int4, int4, timestamptz) IS + 'Approximate libration zone radius (AU) for a Sun-planet Lagrange point.'; + +CREATE FUNCTION lagrange_mass_ratio(int4) RETURNS float8 + AS 'MODULE_PATHNAME', 'lagrange_mass_ratio' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_mass_ratio(int4) IS + 'CR3BP mass parameter mu = M_planet / (M_sun + M_planet) for debugging. body_id: 1-8.'; + +CREATE FUNCTION lagrange_point_name(int4) RETURNS text + AS 'MODULE_PATHNAME', 'lagrange_point_name' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_point_name(int4) IS + 'Human-readable name for a Lagrange point ID (1->''L1'', ..., 5->''L5'').'; + +-- ============================================================ +-- DE variant functions (17) -- STABLE +-- ============================================================ + +CREATE FUNCTION lagrange_heliocentric_de(int4, int4, timestamptz) RETURNS heliocentric + AS 'MODULE_PATHNAME', 'lagrange_heliocentric_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_heliocentric_de(int4, int4, timestamptz) IS + 'DE variant of lagrange_heliocentric(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION lagrange_observe_de(int4, int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'lagrange_observe_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_observe_de(int4, int4, observer, timestamptz) IS + 'DE variant of lagrange_observe(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION lagrange_equatorial_de(int4, int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'lagrange_equatorial_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_equatorial_de(int4, int4, timestamptz) IS + 'DE variant of lagrange_equatorial(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION lagrange_distance_de(int4, int4, heliocentric, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'lagrange_distance_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_distance_de(int4, int4, heliocentric, timestamptz) IS + 'DE variant of lagrange_distance(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION lagrange_distance_oe_de(int4, int4, orbital_elements, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'lagrange_distance_oe_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_distance_oe_de(int4, int4, orbital_elements, timestamptz) IS + 'DE variant of lagrange_distance_oe(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION lunar_lagrange_observe_de(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'lunar_lagrange_observe_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lunar_lagrange_observe_de(int4, observer, timestamptz) IS + 'DE variant of lunar_lagrange_observe(). Falls back to ELP2000-82B if DE unavailable.'; + +CREATE FUNCTION lunar_lagrange_equatorial_de(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'lunar_lagrange_equatorial_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lunar_lagrange_equatorial_de(int4, timestamptz) IS + 'DE variant of lunar_lagrange_equatorial(). Falls back to ELP2000-82B if DE unavailable.'; + +CREATE FUNCTION galilean_lagrange_observe_de(int4, int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'galilean_lagrange_observe_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION galilean_lagrange_observe_de(int4, int4, observer, timestamptz) IS + 'DE variant of galilean_lagrange_observe(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION galilean_lagrange_equatorial_de(int4, int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'galilean_lagrange_equatorial_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION galilean_lagrange_equatorial_de(int4, int4, timestamptz) IS + 'DE variant of galilean_lagrange_equatorial(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION saturn_moon_lagrange_observe_de(int4, int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'saturn_moon_lagrange_observe_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION saturn_moon_lagrange_observe_de(int4, int4, observer, timestamptz) IS + 'DE variant of saturn_moon_lagrange_observe(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION saturn_moon_lagrange_equatorial_de(int4, int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'saturn_moon_lagrange_equatorial_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION saturn_moon_lagrange_equatorial_de(int4, int4, timestamptz) IS + 'DE variant of saturn_moon_lagrange_equatorial(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION uranus_moon_lagrange_observe_de(int4, int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'uranus_moon_lagrange_observe_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION uranus_moon_lagrange_observe_de(int4, int4, observer, timestamptz) IS + 'DE variant of uranus_moon_lagrange_observe(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION uranus_moon_lagrange_equatorial_de(int4, int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'uranus_moon_lagrange_equatorial_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION uranus_moon_lagrange_equatorial_de(int4, int4, timestamptz) IS + 'DE variant of uranus_moon_lagrange_equatorial(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION mars_moon_lagrange_observe_de(int4, int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'mars_moon_lagrange_observe_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION mars_moon_lagrange_observe_de(int4, int4, observer, timestamptz) IS + 'DE variant of mars_moon_lagrange_observe(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION mars_moon_lagrange_equatorial_de(int4, int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'mars_moon_lagrange_equatorial_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION mars_moon_lagrange_equatorial_de(int4, int4, timestamptz) IS + 'DE variant of mars_moon_lagrange_equatorial(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION hill_radius_de(int4, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'hill_radius_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION hill_radius_de(int4, timestamptz) IS + 'DE variant of hill_radius(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION lagrange_zone_radius_de(int4, int4, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'lagrange_zone_radius_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_zone_radius_de(int4, int4, timestamptz) IS + 'DE variant of lagrange_zone_radius(). Falls back to VSOP87 if DE unavailable.'; diff --git a/sql/pg_orrery--0.20.0.sql b/sql/pg_orrery--0.20.0.sql new file mode 100644 index 0000000..54204ba --- /dev/null +++ b/sql/pg_orrery--0.20.0.sql @@ -0,0 +1,2202 @@ +-- pg_orrery -- Orbital mechanics types and functions for PostgreSQL +-- +-- Types: tle, eci_position, geodetic, topocentric, observer, pass_event +-- Provides SGP4/SDP4 propagation, coordinate transforms, pass prediction, +-- and GiST indexing on altitude bands for conjunction screening. +-- +-- All propagation uses WGS-72 constants (matching TLE mean element fitting). +-- Coordinate output uses WGS-84 (matching modern geodetic standards). + +-- ============================================================ +-- Shell types (forward declarations) +-- ============================================================ + +CREATE TYPE tle; +CREATE TYPE eci_position; +CREATE TYPE geodetic; +CREATE TYPE topocentric; +CREATE TYPE observer; +CREATE TYPE pass_event; + + +-- ============================================================ +-- TLE type: Two-Line Element set +-- ============================================================ + +CREATE FUNCTION tle_in(cstring) RETURNS tle + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION tle_out(tle) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION tle_recv(internal) RETURNS tle + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION tle_send(tle) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE tle ( + INPUT = tle_in, + OUTPUT = tle_out, + RECEIVE = tle_recv, + SEND = tle_send, + INTERNALLENGTH = 112, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE tle IS 'Two-Line Element set — parsed mean orbital elements for SGP4/SDP4 propagation'; + +-- TLE accessor functions + +CREATE FUNCTION tle_epoch(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_epoch(tle) IS 'TLE epoch as Julian date (UTC)'; + +CREATE FUNCTION tle_norad_id(tle) RETURNS int4 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_norad_id(tle) IS 'NORAD catalog number'; + +CREATE FUNCTION tle_inclination(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_inclination(tle) IS 'Orbital inclination in degrees'; + +CREATE FUNCTION tle_eccentricity(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_eccentricity(tle) IS 'Orbital eccentricity (dimensionless)'; + +CREATE FUNCTION tle_raan(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_raan(tle) IS 'Right ascension of ascending node in degrees'; + +CREATE FUNCTION tle_arg_perigee(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_arg_perigee(tle) IS 'Argument of perigee in degrees'; + +CREATE FUNCTION tle_mean_anomaly(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_mean_anomaly(tle) IS 'Mean anomaly in degrees'; + +CREATE FUNCTION tle_mean_motion(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_mean_motion(tle) IS 'Mean motion in revolutions per day'; + +CREATE FUNCTION tle_bstar(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_bstar(tle) IS 'B* drag coefficient (1/earth-radii)'; + +CREATE FUNCTION tle_period(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_period(tle) IS 'Orbital period in minutes'; + +CREATE FUNCTION tle_age(tle, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_age(tle, timestamptz) IS 'TLE age in days (positive = past epoch, negative = before epoch)'; + +CREATE FUNCTION tle_perigee(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_perigee(tle) IS 'Perigee altitude in km above WGS-72 ellipsoid'; + +CREATE FUNCTION tle_apogee(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_apogee(tle) IS 'Apogee altitude in km above WGS-72 ellipsoid'; + +CREATE FUNCTION tle_intl_desig(tle) RETURNS text + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_intl_desig(tle) IS 'International designator (COSPAR ID)'; + +CREATE FUNCTION tle_from_lines(text, text) RETURNS tle + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_from_lines(text, text) IS + 'Construct TLE from separate line1/line2 text columns'; + + +-- ============================================================ +-- ECI position type: True Equator Mean Equinox (TEME) frame +-- ============================================================ + +CREATE FUNCTION eci_in(cstring) RETURNS eci_position + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_out(eci_position) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_recv(internal) RETURNS eci_position + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_send(eci_position) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE eci_position ( + INPUT = eci_in, + OUTPUT = eci_out, + RECEIVE = eci_recv, + SEND = eci_send, + INTERNALLENGTH = 48, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE eci_position IS 'Earth-Centered Inertial position and velocity in TEME frame (km, km/s)'; + +-- ECI accessor functions + +CREATE FUNCTION eci_x(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_y(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_z(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_vx(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_vy(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_vz(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_speed(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eci_speed(eci_position) IS 'Velocity magnitude in km/s'; + +CREATE FUNCTION eci_altitude(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eci_altitude(eci_position) IS 'Approximate geocentric altitude in km (radius - WGS72_AE)'; + + +-- ============================================================ +-- Geodetic type: WGS-84 latitude/longitude/altitude +-- ============================================================ + +CREATE FUNCTION geodetic_in(cstring) RETURNS geodetic + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION geodetic_out(geodetic) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION geodetic_recv(internal) RETURNS geodetic + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION geodetic_send(geodetic) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE geodetic ( + INPUT = geodetic_in, + OUTPUT = geodetic_out, + RECEIVE = geodetic_recv, + SEND = geodetic_send, + INTERNALLENGTH = 24, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE geodetic IS 'Geodetic coordinates on WGS-84 ellipsoid (lat/lon in degrees, altitude in km)'; + +CREATE FUNCTION geodetic_lat(geodetic) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION geodetic_lon(geodetic) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION geodetic_alt(geodetic) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + + +-- ============================================================ +-- Topocentric type: observer-relative az/el/range +-- ============================================================ + +CREATE FUNCTION topocentric_in(cstring) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION topocentric_out(topocentric) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION topocentric_recv(internal) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION topocentric_send(topocentric) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE topocentric ( + INPUT = topocentric_in, + OUTPUT = topocentric_out, + RECEIVE = topocentric_recv, + SEND = topocentric_send, + INTERNALLENGTH = 32, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE topocentric IS 'Topocentric coordinates relative to observer (azimuth, elevation, range, range rate)'; + +CREATE FUNCTION topo_azimuth(topocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION topo_azimuth(topocentric) IS 'Azimuth in degrees (0=N, 90=E, 180=S, 270=W)'; + +CREATE FUNCTION topo_elevation(topocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION topo_elevation(topocentric) IS 'Elevation in degrees (0=horizon, 90=zenith)'; + +CREATE FUNCTION topo_range(topocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION topo_range(topocentric) IS 'Slant range in km'; + +CREATE FUNCTION topo_range_rate(topocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION topo_range_rate(topocentric) IS 'Range rate in km/s (positive = receding)'; + + +-- ============================================================ +-- Observer type: ground station location +-- ============================================================ + +CREATE FUNCTION observer_in(cstring) RETURNS observer + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION observer_out(observer) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION observer_recv(internal) RETURNS observer + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION observer_send(observer) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE observer ( + INPUT = observer_in, + OUTPUT = observer_out, + RECEIVE = observer_recv, + SEND = observer_send, + INTERNALLENGTH = 24, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE observer IS 'Ground observer location (accepts: 40.0N 105.3W 1655m or decimal degrees)'; + +CREATE FUNCTION observer_lat(observer) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION observer_lat(observer) IS 'Latitude in degrees (positive = North)'; + +CREATE FUNCTION observer_lon(observer) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION observer_lon(observer) IS 'Longitude in degrees (positive = East)'; + +CREATE FUNCTION observer_alt(observer) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION observer_alt(observer) IS 'Altitude in meters above WGS-84 ellipsoid'; + +CREATE FUNCTION observer_from_geodetic(float8, float8, float8 DEFAULT 0.0) RETURNS observer + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION observer_from_geodetic(float8, float8, float8) IS + 'Construct observer from lat (deg), lon (deg), altitude (meters). Avoids text formatting round-trips.'; + + +-- ============================================================ +-- Pass event type: satellite visibility window +-- ============================================================ + +CREATE FUNCTION pass_event_in(cstring) RETURNS pass_event + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION pass_event_out(pass_event) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION pass_event_recv(internal) RETURNS pass_event + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION pass_event_send(pass_event) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE pass_event ( + INPUT = pass_event_in, + OUTPUT = pass_event_out, + RECEIVE = pass_event_recv, + SEND = pass_event_send, + INTERNALLENGTH = 48, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE pass_event IS 'Satellite pass event (AOS/MAX/LOS times, max elevation, AOS/LOS azimuths)'; + +CREATE FUNCTION pass_aos_time(pass_event) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_aos_time(pass_event) IS 'Acquisition of signal time'; + +CREATE FUNCTION pass_max_el_time(pass_event) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_max_el_time(pass_event) IS 'Maximum elevation time'; + +CREATE FUNCTION pass_los_time(pass_event) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_los_time(pass_event) IS 'Loss of signal time'; + +CREATE FUNCTION pass_max_elevation(pass_event) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_max_elevation(pass_event) IS 'Maximum elevation in degrees'; + +CREATE FUNCTION pass_aos_azimuth(pass_event) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_aos_azimuth(pass_event) IS 'AOS azimuth in degrees (0=N)'; + +CREATE FUNCTION pass_los_azimuth(pass_event) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_los_azimuth(pass_event) IS 'LOS azimuth in degrees (0=N)'; + +CREATE FUNCTION pass_duration(pass_event) RETURNS interval + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_duration(pass_event) IS 'Pass duration (LOS - AOS)'; + + +-- ============================================================ +-- SGP4/SDP4 propagation functions +-- ============================================================ + +CREATE FUNCTION sgp4_propagate(tle, timestamptz) RETURNS eci_position + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sgp4_propagate(tle, timestamptz) IS + 'Propagate TLE to a point in time using SGP4 (near-earth) or SDP4 (deep-space). Returns TEME ECI position/velocity.'; + +CREATE FUNCTION sgp4_propagate_safe(tle, timestamptz) RETURNS eci_position + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE PARALLEL SAFE; +COMMENT ON FUNCTION sgp4_propagate_safe(tle, timestamptz) IS + 'Like sgp4_propagate but returns NULL on error instead of raising an exception. For batch queries with potentially invalid TLEs.'; + +CREATE FUNCTION sgp4_propagate_series(tle, timestamptz, timestamptz, interval) + RETURNS TABLE(t timestamptz, x float8, y float8, z float8, vx float8, vy float8, vz float8) + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE + ROWS 100; +COMMENT ON FUNCTION sgp4_propagate_series(tle, timestamptz, timestamptz, interval) IS + 'Propagate TLE over a time range at regular intervals. Returns time series of TEME positions.'; + +CREATE FUNCTION tle_distance(tle, tle, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_distance(tle, tle, timestamptz) IS + 'Euclidean distance in km between two TLEs at a reference time'; + + +-- ============================================================ +-- Coordinate transform functions +-- ============================================================ + +CREATE FUNCTION eci_to_geodetic(eci_position, timestamptz) RETURNS geodetic + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eci_to_geodetic(eci_position, timestamptz) IS + 'Convert TEME ECI position to WGS-84 geodetic coordinates at given time'; + +CREATE FUNCTION eci_to_topocentric(eci_position, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eci_to_topocentric(eci_position, observer, timestamptz) IS + 'Convert TEME ECI position to topocentric (az/el/range) relative to observer'; + +CREATE FUNCTION subsatellite_point(tle, timestamptz) RETURNS geodetic + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION subsatellite_point(tle, timestamptz) IS + 'Subsatellite (nadir) point on WGS-84 ellipsoid at given time'; + +CREATE FUNCTION ground_track(tle, timestamptz, timestamptz, interval) + RETURNS TABLE(t timestamptz, lat float8, lon float8, alt float8) + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE + ROWS 100; +COMMENT ON FUNCTION ground_track(tle, timestamptz, timestamptz, interval) IS + 'Ground track as time series of subsatellite points (lat/lon in degrees, alt in km)'; + +CREATE FUNCTION observe(tle, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION observe(tle, observer, timestamptz) IS + 'Propagate TLE and compute observer-relative look angles in one call. Returns topocentric (az/el/range/range_rate).'; + +CREATE FUNCTION observe_safe(tle, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE PARALLEL SAFE; +COMMENT ON FUNCTION observe_safe(tle, observer, timestamptz) IS + 'Like observe() but returns NULL on propagation error. For batch queries with potentially invalid/decayed TLEs.'; + + +-- ============================================================ +-- Pass prediction functions +-- ============================================================ + +CREATE FUNCTION next_pass(tle, observer, timestamptz) RETURNS pass_event + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION next_pass(tle, observer, timestamptz) IS + 'Find the next satellite pass over observer (searches up to 7 days ahead)'; + +CREATE FUNCTION predict_passes(tle, observer, timestamptz, timestamptz, float8 DEFAULT 0.0) + RETURNS SETOF pass_event + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE + ROWS 10; +COMMENT ON FUNCTION predict_passes(tle, observer, timestamptz, timestamptz, float8) IS + 'Predict all satellite passes over observer in time window. Optional min_elevation in degrees.'; + +CREATE FUNCTION pass_visible(tle, observer, timestamptz, timestamptz) RETURNS boolean + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_visible(tle, observer, timestamptz, timestamptz) IS + 'True if any pass occurs over observer in the time window'; + + +-- ============================================================ +-- GiST operator support functions +-- ============================================================ + +-- Overlap operator: do orbital keys overlap in altitude AND inclination? +CREATE FUNCTION tle_overlap(tle, tle) RETURNS boolean + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +-- Altitude distance operator (altitude-only, for KNN ordering) +CREATE FUNCTION tle_alt_distance(tle, tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OPERATOR && ( + LEFTARG = tle, + RIGHTARG = tle, + FUNCTION = tle_overlap, + COMMUTATOR = &&, + RESTRICT = areasel, + JOIN = areajoinsel +); + +COMMENT ON OPERATOR && (tle, tle) IS 'Orbital key overlap (altitude band AND inclination range) — necessary condition for conjunction'; + +CREATE OPERATOR <-> ( + LEFTARG = tle, + RIGHTARG = tle, + FUNCTION = tle_alt_distance, + COMMUTATOR = <-> +); + +COMMENT ON OPERATOR <-> (tle, tle) IS '2-D orbital distance in km: L2 norm of altitude-band gap and inclination gap (radians × Earth radius). Returns 0 when both dimensions overlap.'; + + +-- ============================================================ +-- GiST operator class for 2-D orbital indexing (altitude + inclination) +-- ============================================================ + +-- GiST internal support functions +CREATE FUNCTION gist_tle_compress(internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_decompress(internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_consistent(internal, tle, smallint, oid, internal) RETURNS boolean + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_union(internal, internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_penalty(internal, internal, internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_picksplit(internal, internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_same(internal, internal, internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_distance(internal, tle, smallint, oid, internal) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OPERATOR CLASS tle_ops + DEFAULT FOR TYPE tle USING gist AS + OPERATOR 3 && , + OPERATOR 15 <-> (tle, tle) FOR ORDER BY float_ops, + FUNCTION 1 gist_tle_consistent(internal, tle, smallint, oid, internal), + FUNCTION 2 gist_tle_union(internal, internal), + FUNCTION 3 gist_tle_compress(internal), + FUNCTION 4 gist_tle_decompress(internal), + FUNCTION 5 gist_tle_penalty(internal, internal, internal), + FUNCTION 6 gist_tle_picksplit(internal, internal), + FUNCTION 7 gist_tle_same(internal, internal, internal), + FUNCTION 8 gist_tle_distance(internal, tle, smallint, oid, internal); + + +-- ============================================================ +-- Heliocentric type: ecliptic J2000 position in AU +-- ============================================================ + +CREATE TYPE heliocentric; + +CREATE FUNCTION heliocentric_in(cstring) RETURNS heliocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION heliocentric_out(heliocentric) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION heliocentric_recv(internal) RETURNS heliocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION heliocentric_send(heliocentric) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE heliocentric ( + INPUT = heliocentric_in, + OUTPUT = heliocentric_out, + RECEIVE = heliocentric_recv, + SEND = heliocentric_send, + INTERNALLENGTH = 24, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE heliocentric IS 'Heliocentric position in ecliptic J2000 frame (x, y, z in AU)'; + +CREATE FUNCTION helio_x(heliocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION helio_x(heliocentric) IS 'X component in AU (ecliptic J2000, vernal equinox direction)'; + +CREATE FUNCTION helio_y(heliocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION helio_y(heliocentric) IS 'Y component in AU (ecliptic J2000)'; + +CREATE FUNCTION helio_z(heliocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION helio_z(heliocentric) IS 'Z component in AU (ecliptic J2000, north ecliptic pole)'; + +CREATE FUNCTION helio_distance(heliocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION helio_distance(heliocentric) IS 'Heliocentric distance in AU'; + + +-- ============================================================ +-- Star observation functions +-- ============================================================ + +CREATE FUNCTION star_observe(float8, float8, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION star_observe(float8, float8, observer, timestamptz) IS + 'Observe a star from (ra_hours J2000, dec_degrees J2000, observer, time). Returns topocentric az/el. Range is 0 (infinite distance).'; + +CREATE FUNCTION star_observe_safe(float8, float8, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE PARALLEL SAFE; +COMMENT ON FUNCTION star_observe_safe(float8, float8, observer, timestamptz) IS + 'Like star_observe but returns NULL on invalid inputs. For batch queries over star catalogs.'; + + +-- ============================================================ +-- Keplerian propagation functions +-- ============================================================ + +CREATE FUNCTION kepler_propagate( + q_au float8, eccentricity float8, + inclination_deg float8, arg_perihelion_deg float8, + long_asc_node_deg float8, perihelion_jd float8, + t timestamptz +) RETURNS heliocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION kepler_propagate(float8, float8, float8, float8, float8, float8, timestamptz) IS + 'Two-body Keplerian propagation from classical orbital elements. Returns heliocentric ecliptic J2000 position in AU. Handles elliptic, parabolic, and hyperbolic orbits.'; + + +-- ============================================================ +-- Comet observation +-- ============================================================ + +CREATE FUNCTION comet_observe( + q_au float8, eccentricity float8, + inclination_deg float8, arg_perihelion_deg float8, + long_asc_node_deg float8, perihelion_jd float8, + earth_x_au float8, earth_y_au float8, earth_z_au float8, + obs observer, t timestamptz +) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION comet_observe(float8, float8, float8, float8, float8, float8, float8, float8, float8, observer, timestamptz) IS + 'Observe a comet/asteroid from orbital elements. Requires Earth heliocentric ecliptic J2000 position (AU). Returns topocentric az/el with geocentric range in km.'; + + +-- ============================================================ +-- VSOP87 planets, ELP82B Moon, Sun observation +-- ============================================================ + +CREATE FUNCTION planet_heliocentric(int4, timestamptz) RETURNS heliocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_heliocentric(int4, timestamptz) IS + 'VSOP87 heliocentric ecliptic J2000 position (AU). Body IDs: 0=Sun, 1=Mercury, ..., 8=Neptune.'; + +CREATE FUNCTION planet_observe(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_observe(int4, observer, timestamptz) IS + 'Observe a planet from (body_id 1-8, observer, time). Returns topocentric az/el with geocentric range in km.'; + +CREATE FUNCTION sun_observe(observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_observe(observer, timestamptz) IS + 'Observe the Sun from (observer, time). Returns topocentric az/el with geocentric range in km.'; + +CREATE FUNCTION moon_observe(observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_observe(observer, timestamptz) IS + 'Observe the Moon via ELP2000-82B from (observer, time). Returns topocentric az/el with geocentric range in km.'; + + +-- ============================================================ +-- Planetary moon observation +-- ============================================================ + +CREATE FUNCTION galilean_observe(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION galilean_observe(int4, observer, timestamptz) IS + 'Observe a Galilean moon of Jupiter via L1.2 theory. Body IDs: 0=Io, 1=Europa, 2=Ganymede, 3=Callisto.'; + +CREATE FUNCTION saturn_moon_observe(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION saturn_moon_observe(int4, observer, timestamptz) IS + 'Observe a Saturn moon via TASS 1.7. Body IDs: 0=Mimas, 1=Enceladus, 2=Tethys, 3=Dione, 4=Rhea, 5=Titan, 6=Iapetus, 7=Hyperion.'; + +CREATE FUNCTION uranus_moon_observe(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION uranus_moon_observe(int4, observer, timestamptz) IS + 'Observe a Uranus moon via GUST86. Body IDs: 0=Miranda, 1=Ariel, 2=Umbriel, 3=Titania, 4=Oberon.'; + +CREATE FUNCTION mars_moon_observe(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION mars_moon_observe(int4, observer, timestamptz) IS + 'Observe a Mars moon via MarsSat. Body IDs: 0=Phobos, 1=Deimos.'; + + +-- ============================================================ +-- Jupiter decametric radio burst prediction +-- ============================================================ + +CREATE FUNCTION io_phase_angle(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION io_phase_angle(timestamptz) IS + 'Io orbital phase angle in degrees [0,360). 0=superior conjunction (behind Jupiter). Standard Radio JOVE convention.'; + +CREATE FUNCTION jupiter_cml(observer, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION jupiter_cml(observer, timestamptz) IS + 'Jupiter Central Meridian Longitude, System III (1965.0), in degrees [0,360). Light-time corrected.'; + +CREATE FUNCTION jupiter_burst_probability(float8, float8) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION jupiter_burst_probability(float8, float8) IS + 'Estimated Jupiter decametric burst probability (0-1) from (io_phase_deg, cml_deg). Based on Carr et al. (1983) source regions A, B, C, D.'; + + +-- ============================================================ +-- Interplanetary transfer orbits (Lambert solver) +-- ============================================================ + +CREATE FUNCTION lambert_transfer( + dep_body_id int4, arr_body_id int4, + dep_time timestamptz, arr_time timestamptz, + OUT c3_departure float8, OUT c3_arrival float8, + OUT v_inf_departure float8, OUT v_inf_arrival float8, + OUT tof_days float8, OUT transfer_sma float8 +) RETURNS RECORD + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lambert_transfer(int4, int4, timestamptz, timestamptz) IS + 'Solve Lambert transfer between two planets. Returns C3 (km^2/s^2), v_infinity (km/s), TOF (days), SMA (AU). Body IDs 1-8.'; + +CREATE FUNCTION lambert_c3(int4, int4, timestamptz, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lambert_c3(int4, int4, timestamptz, timestamptz) IS + 'Departure C3 (km^2/s^2) for a Lambert transfer. Returns NULL if solver fails. For pork chop plots.'; + + +-- ============================================================ +-- DE ephemeris functions (optional high-precision) +-- ============================================================ + +CREATE FUNCTION planet_heliocentric_de(int4, timestamptz) RETURNS heliocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_heliocentric_de(int4, timestamptz) IS + 'Heliocentric ecliptic J2000 position via JPL DE (sub-arcsecond). Falls back to VSOP87 if DE unavailable. Body IDs: 0=Sun, 1-8=Mercury-Neptune.'; + +CREATE FUNCTION planet_observe_de(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_observe_de(int4, observer, timestamptz) IS + 'Observe planet via JPL DE. Falls back to VSOP87. Body IDs: 1-8 (Mercury-Neptune).'; + +CREATE FUNCTION sun_observe_de(observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_observe_de(observer, timestamptz) IS + 'Observe Sun via JPL DE. Falls back to VSOP87.'; + +CREATE FUNCTION moon_observe_de(observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_observe_de(observer, timestamptz) IS + 'Observe Moon via JPL DE. Falls back to ELP2000-82B.'; + +CREATE FUNCTION lambert_transfer_de( + dep_body_id int4, arr_body_id int4, + dep_time timestamptz, arr_time timestamptz, + OUT c3_departure float8, OUT c3_arrival float8, + OUT v_inf_departure float8, OUT v_inf_arrival float8, + OUT tof_days float8, OUT transfer_sma float8 +) RETURNS RECORD + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lambert_transfer_de(int4, int4, timestamptz, timestamptz) IS + 'Lambert transfer via JPL DE positions. Falls back to VSOP87. Body IDs 1-8.'; + +CREATE FUNCTION lambert_c3_de(int4, int4, timestamptz, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lambert_c3_de(int4, int4, timestamptz, timestamptz) IS + 'Departure C3 via JPL DE. Falls back to VSOP87. For pork chop plots.'; + +CREATE FUNCTION galilean_observe_de(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION galilean_observe_de(int4, observer, timestamptz) IS + 'Observe Galilean moon with JPL DE parent position. L1.2 moon offsets. Body IDs: 0-3 (Io-Callisto).'; + +CREATE FUNCTION saturn_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION saturn_moon_observe_de(int4, observer, timestamptz) IS + 'Observe Saturn moon with JPL DE parent position. TASS 1.7 moon offsets. Body IDs: 0-7 (Mimas-Hyperion).'; + +CREATE FUNCTION uranus_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION uranus_moon_observe_de(int4, observer, timestamptz) IS + 'Observe Uranus moon with JPL DE parent position. GUST86 moon offsets. Body IDs: 0-4 (Miranda-Oberon).'; + +CREATE FUNCTION mars_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION mars_moon_observe_de(int4, observer, timestamptz) IS + 'Observe Mars moon with JPL DE parent position. MarsSat moon offsets. Body IDs: 0-1 (Phobos-Deimos).'; + + +-- Diagnostic function + +CREATE FUNCTION pg_orrery_ephemeris_info( + OUT provider text, OUT file_path text, + OUT start_jd float8, OUT end_jd float8, + OUT version int4, OUT au_km float8 +) RETURNS RECORD + AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE; +COMMENT ON FUNCTION pg_orrery_ephemeris_info() IS + 'Returns current ephemeris provider status: VSOP87 or JPL_DE with file path, JD range, version, and AU value.'; + + +-- ============================================================ +-- Orbit determination (TLE fitting from observations) +-- ============================================================ + +-- Fit TLE from ECI position/velocity ephemeris +-- weights: per-observation weighting (NULL = uniform) + +CREATE FUNCTION tle_from_eci( + positions eci_position[], times timestamptz[], + seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false, + max_iter int4 DEFAULT 15, + weights float8[] DEFAULT NULL, + OUT fitted_tle tle, OUT iterations int4, + OUT rms_final float8, OUT rms_initial float8, OUT status text, + OUT condition_number float8, OUT covariance float8[], OUT nstate int4 +) RETURNS RECORD + AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE; +COMMENT ON FUNCTION tle_from_eci(eci_position[], timestamptz[], tle, boolean, int4, float8[]) IS + 'Fit a TLE from ECI position/velocity observations via differential correction. Optional per-observation weights for heterogeneous sensor fusion. Returns fitted TLE, RMS residuals, convergence status, condition number, and formal covariance matrix.'; + +-- Fit TLE from topocentric observations (az/el/range) — single observer +-- fit_range_rate: include range_rate as 4th residual component +-- weights: per-observation weighting (NULL = uniform) + +CREATE FUNCTION tle_from_topocentric( + observations topocentric[], times timestamptz[], + obs observer, + seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false, + max_iter int4 DEFAULT 15, + fit_range_rate boolean DEFAULT false, + weights float8[] DEFAULT NULL, + OUT fitted_tle tle, OUT iterations int4, + OUT rms_final float8, OUT rms_initial float8, OUT status text, + OUT condition_number float8, OUT covariance float8[], OUT nstate int4 +) RETURNS RECORD + AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE; +COMMENT ON FUNCTION tle_from_topocentric(topocentric[], timestamptz[], observer, tle, boolean, int4, boolean, float8[]) IS + 'Fit a TLE from topocentric (az/el/range) observations via differential correction. Optional range_rate fitting and per-observation weights. Returns fitted TLE, RMS residuals, convergence status, condition number, and formal covariance matrix.'; + +-- Fit TLE from topocentric observations — multiple observers + +CREATE FUNCTION tle_from_topocentric( + observations topocentric[], times timestamptz[], + observers observer[], observer_ids int4[], + seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false, + max_iter int4 DEFAULT 15, + fit_range_rate boolean DEFAULT false, + weights float8[] DEFAULT NULL, + OUT fitted_tle tle, OUT iterations int4, + OUT rms_final float8, OUT rms_initial float8, OUT status text, + OUT condition_number float8, OUT covariance float8[], OUT nstate int4 +) RETURNS RECORD + AS 'MODULE_PATHNAME', 'tle_from_topocentric_multi' + LANGUAGE C STABLE PARALLEL SAFE; +COMMENT ON FUNCTION tle_from_topocentric(topocentric[], timestamptz[], observer[], int4[], tle, boolean, int4, boolean, float8[]) IS + 'Fit a TLE from topocentric observations collected by multiple ground stations. observer_ids[i] indexes into observers[]. Optional range_rate fitting and per-observation weights.'; + +-- Per-observation residuals diagnostic + +CREATE FUNCTION tle_fit_residuals( + fitted tle, + positions eci_position[], + times timestamptz[] +) RETURNS TABLE ( + t timestamptz, + dx_km float8, + dy_km float8, + dz_km float8, + pos_err_km float8 +) + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_fit_residuals(tle, eci_position[], timestamptz[]) IS + 'Compute per-observation position residuals (km) between a TLE and ECI observations. Useful for fit quality assessment.'; + +-- Fit TLE from RA/Dec observations — single observer +-- Uses Gauss method for initial orbit determination when no seed is provided. +-- RA in hours [0,24), Dec in degrees [-90,90] (matches star_observe convention). +-- RMS output is in radians for angles-only mode. + +CREATE FUNCTION tle_from_angles( + ra_hours float8[], dec_degrees float8[], + times timestamptz[], + obs observer, + seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false, + max_iter int4 DEFAULT 15, + weights float8[] DEFAULT NULL, + OUT fitted_tle tle, OUT iterations int4, + OUT rms_final float8, OUT rms_initial float8, OUT status text, + OUT condition_number float8, OUT covariance float8[], OUT nstate int4 +) RETURNS RECORD AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE; +COMMENT ON FUNCTION tle_from_angles(float8[], float8[], timestamptz[], observer, tle, boolean, int4, float8[]) IS + 'Fit a TLE from angles-only (RA/Dec) observations via Gauss IOD + differential correction. RA in hours [0,24), Dec in degrees [-90,90]. RMS output in radians. Uses Gauss method for seed-free initial guess.'; + +-- Fit TLE from RA/Dec observations — multiple observers + +CREATE FUNCTION tle_from_angles( + ra_hours float8[], dec_degrees float8[], + times timestamptz[], + observers observer[], observer_ids int4[], + seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false, + max_iter int4 DEFAULT 15, + weights float8[] DEFAULT NULL, + OUT fitted_tle tle, OUT iterations int4, + OUT rms_final float8, OUT rms_initial float8, OUT status text, + OUT condition_number float8, OUT covariance float8[], OUT nstate int4 +) RETURNS RECORD + AS 'MODULE_PATHNAME', 'tle_from_angles_multi' + LANGUAGE C STABLE PARALLEL SAFE; +COMMENT ON FUNCTION tle_from_angles(float8[], float8[], timestamptz[], observer[], int4[], tle, boolean, int4, float8[]) IS + 'Fit a TLE from angles-only (RA/Dec) observations from multiple ground stations via Gauss IOD + differential correction.'; +-- pg_orrery 0.6.0 -> 0.7.0 migration +-- +-- Adds SP-GiST orbital trie index for satellite pass prediction. +-- 2-level trie: SMA (L0) + inclination (L1) with query-time RAAN filter. +-- The &? operator answers "might this satellite be visible?" + +-- ============================================================ +-- observer_window composite type (query parameter bundle) +-- ============================================================ + +CREATE TYPE observer_window AS ( + obs observer, + t_start timestamptz, + t_end timestamptz, + min_el float8 +); + +COMMENT ON TYPE observer_window IS + 'Observation query parameters: observer location, time window, and minimum elevation angle (degrees). Used with the &? visibility cone operator.'; + +-- ============================================================ +-- Visibility cone operator function +-- ============================================================ + +CREATE FUNCTION tle_visibility_possible(tle, observer_window) RETURNS boolean + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; + +COMMENT ON FUNCTION tle_visibility_possible(tle, observer_window) IS + 'Could this satellite be visible from the observer during the time window? Combines altitude, inclination, and RAAN checks. Conservative superset — survivors need SGP4 propagation for ground truth.'; + +-- ============================================================ +-- &? operator (visibility cone check) +-- ============================================================ +-- The indexed column (tle) MUST be the left argument so PostgreSQL +-- can form a ScanKey and pass it to inner_consistent for pruning. + +CREATE OPERATOR &? ( + LEFTARG = tle, + RIGHTARG = observer_window, + FUNCTION = tle_visibility_possible, + RESTRICT = contsel, + JOIN = contjoinsel +); + +COMMENT ON OPERATOR &? (tle, observer_window) IS + 'Visibility cone check: could this satellite be visible from the observer during the time window? Index-accelerated via SP-GiST orbital trie.'; + +-- ============================================================ +-- SP-GiST support functions +-- ============================================================ + +CREATE FUNCTION spgist_tle_config(internal, internal) RETURNS void + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION spgist_tle_choose(internal, internal) RETURNS void + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION spgist_tle_picksplit(internal, internal) RETURNS void + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION spgist_tle_inner_consistent(internal, internal) RETURNS void + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION spgist_tle_leaf_consistent(internal, internal) RETURNS void + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +-- ============================================================ +-- SP-GiST operator class (opt-in, not DEFAULT) +-- ============================================================ + +CREATE OPERATOR CLASS tle_spgist_ops + FOR TYPE tle USING spgist AS + OPERATOR 1 &? (tle, observer_window), + FUNCTION 1 spgist_tle_config(internal, internal), + FUNCTION 2 spgist_tle_choose(internal, internal), + FUNCTION 3 spgist_tle_picksplit(internal, internal), + FUNCTION 4 spgist_tle_inner_consistent(internal, internal), + FUNCTION 5 spgist_tle_leaf_consistent(internal, internal); +-- pg_orrery 0.7.0 -> 0.8.0 migration +-- +-- Adds orbital_elements type for comets/asteroids, MPC MPCORB.DAT parser, +-- and small_body_observe()/small_body_heliocentric() observation functions. + +-- ============================================================ +-- orbital_elements type +-- ============================================================ + +CREATE TYPE orbital_elements; + +CREATE FUNCTION orbital_elements_in(cstring) RETURNS orbital_elements + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION orbital_elements_out(orbital_elements) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION orbital_elements_recv(internal) RETURNS orbital_elements + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION orbital_elements_send(orbital_elements) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE orbital_elements ( + INPUT = orbital_elements_in, + OUTPUT = orbital_elements_out, + RECEIVE = orbital_elements_recv, + SEND = orbital_elements_send, + INTERNALLENGTH = 72, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE orbital_elements IS + 'Classical Keplerian orbital elements for comets and asteroids (epoch, q, e, inc, omega, Omega, tp, H, G). 72 bytes, fixed-size.'; + + +-- ============================================================ +-- Accessor functions +-- ============================================================ + +CREATE FUNCTION oe_epoch(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_epoch(orbital_elements) IS 'Osculation epoch (Julian date)'; + +CREATE FUNCTION oe_perihelion(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_perihelion(orbital_elements) IS 'Perihelion distance q (AU)'; + +CREATE FUNCTION oe_eccentricity(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_eccentricity(orbital_elements) IS 'Eccentricity'; + +CREATE FUNCTION oe_inclination(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_inclination(orbital_elements) IS 'Inclination (degrees)'; + +CREATE FUNCTION oe_arg_perihelion(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_arg_perihelion(orbital_elements) IS 'Argument of perihelion (degrees)'; + +CREATE FUNCTION oe_raan(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_raan(orbital_elements) IS 'Longitude of ascending node (degrees)'; + +CREATE FUNCTION oe_tp(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_tp(orbital_elements) IS 'Time of perihelion passage (Julian date)'; + +CREATE FUNCTION oe_h_mag(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_h_mag(orbital_elements) IS 'Absolute magnitude H (NaN if unknown)'; + +CREATE FUNCTION oe_g_slope(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_g_slope(orbital_elements) IS 'Slope parameter G (NaN if unknown)'; + +CREATE FUNCTION oe_semi_major_axis(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_semi_major_axis(orbital_elements) IS 'Semi-major axis a = q/(1-e) in AU. NULL for parabolic/hyperbolic orbits (e >= 1).'; + +CREATE FUNCTION oe_period_years(orbital_elements) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_period_years(orbital_elements) IS 'Orbital period in years = a^1.5 (Kepler third law). NULL for parabolic/hyperbolic orbits (e >= 1).'; + + +-- ============================================================ +-- MPC MPCORB.DAT parser +-- ============================================================ + +CREATE FUNCTION oe_from_mpc(text) RETURNS orbital_elements + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION oe_from_mpc(text) IS + 'Parse one MPCORB.DAT fixed-width line into orbital_elements. Converts MPC packed epoch, computes perihelion distance and tp from (a, e, M).'; + + +-- ============================================================ +-- Observation functions +-- ============================================================ + +CREATE FUNCTION small_body_heliocentric(orbital_elements, timestamptz) RETURNS heliocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION small_body_heliocentric(orbital_elements, timestamptz) IS + 'Heliocentric ecliptic J2000 position of a comet/asteroid from its orbital elements at a given time.'; + +CREATE FUNCTION small_body_observe(orbital_elements, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION small_body_observe(orbital_elements, observer, timestamptz) IS + 'Observe a comet/asteroid from orbital elements. Auto-fetches Earth via VSOP87. Returns topocentric az/el with geocentric range in km.'; +-- pg_orrery 0.8.0 -> 0.9.0 migration +-- +-- Adds equatorial type (apparent RA/Dec of date), atmospheric refraction, +-- stellar proper motion, and light-time corrected _apparent() functions. + +-- ============================================================ +-- equatorial type — apparent RA/Dec of date +-- ============================================================ + +CREATE TYPE equatorial; + +CREATE FUNCTION equatorial_in(cstring) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION equatorial_out(equatorial) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION equatorial_recv(internal) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION equatorial_send(equatorial) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE equatorial ( + INPUT = equatorial_in, + OUTPUT = equatorial_out, + RECEIVE = equatorial_recv, + SEND = equatorial_send, + INTERNALLENGTH = 24, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE equatorial IS + 'Apparent equatorial coordinates of date: RA (hours), Dec (degrees), distance (km). Solar system: J2000 precessed via IAU 1976. Satellites: TEME frame (~of-date to ~arcsecond). 24 bytes, fixed-size.'; + + +-- ============================================================ +-- Equatorial accessor functions +-- ============================================================ + +CREATE FUNCTION eq_ra(equatorial) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eq_ra(equatorial) IS 'Right ascension in hours [0, 24)'; + +CREATE FUNCTION eq_dec(equatorial) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eq_dec(equatorial) IS 'Declination in degrees [-90, 90]'; + +CREATE FUNCTION eq_distance(equatorial) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eq_distance(equatorial) IS 'Distance in km (0 for stars without parallax)'; + + +-- ============================================================ +-- Satellite RA/Dec functions +-- ============================================================ + +CREATE FUNCTION eci_to_equatorial(eci_position, observer, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eci_to_equatorial(eci_position, observer, timestamptz) IS + 'Topocentric apparent RA/Dec from ECI position. Observer parallax-corrected — LEO parallax is ~1 degree.'; + +CREATE FUNCTION eci_to_equatorial_geo(eci_position, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eci_to_equatorial_geo(eci_position, timestamptz) IS + 'Geocentric apparent RA/Dec from ECI position. Observer-independent — the direction of the TEME position vector.'; + + +-- ============================================================ +-- Solar system equatorial functions (VSOP87) +-- ============================================================ + +CREATE FUNCTION planet_equatorial(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_equatorial(int4, timestamptz) IS + 'Geocentric apparent RA/Dec of a planet via VSOP87. Body IDs: 1=Mercury through 8=Neptune.'; + +CREATE FUNCTION sun_equatorial(timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_equatorial(timestamptz) IS + 'Geocentric apparent RA/Dec of the Sun via VSOP87.'; + +CREATE FUNCTION moon_equatorial(timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_equatorial(timestamptz) IS + 'Geocentric apparent RA/Dec of the Moon via ELP2000-82B.'; + +CREATE FUNCTION small_body_equatorial(orbital_elements, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION small_body_equatorial(orbital_elements, timestamptz) IS + 'Geocentric apparent RA/Dec of a comet/asteroid from orbital elements. Earth via VSOP87.'; + +CREATE FUNCTION star_equatorial(float8, float8, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION star_equatorial(float8, float8, timestamptz) IS + 'Apparent RA/Dec of a star at a given time. Precesses J2000 catalog coordinates (RA hours, Dec degrees) to date via IAU 1976.'; + + +-- ============================================================ +-- Atmospheric refraction (Bennett 1982) +-- ============================================================ + +CREATE FUNCTION atmospheric_refraction(float8) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION atmospheric_refraction(float8) IS + 'Atmospheric refraction correction in degrees for a given geometric elevation (degrees). Standard atmosphere: P=1010 mbar, T=10C. Bennett (1982) formula with domain guard at -1 degree.'; + +CREATE FUNCTION atmospheric_refraction_ext(float8, float8, float8) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION atmospheric_refraction_ext(float8, float8, float8) IS + 'Atmospheric refraction with pressure/temperature correction. Args: elevation_deg, pressure_mbar, temperature_celsius. Meeus P/T factor applied to Bennett formula.'; + +CREATE FUNCTION topo_elevation_apparent(topocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION topo_elevation_apparent(topocentric) IS + 'Apparent elevation in degrees — geometric elevation plus atmospheric refraction correction.'; + + +-- ============================================================ +-- Refracted pass prediction +-- ============================================================ + +CREATE FUNCTION predict_passes_refracted( + tle, observer, timestamptz, timestamptz, float8 DEFAULT 0.0 +) RETURNS SETOF pass_event + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE + ROWS 20; +COMMENT ON FUNCTION predict_passes_refracted(tle, observer, timestamptz, timestamptz, float8) IS + 'Predict satellite passes using a refracted horizon threshold (-0.569 deg geometric). Atmospheric refraction makes satellites visible ~35 seconds earlier at AOS and later at LOS.'; + + +-- ============================================================ +-- Stellar proper motion +-- ============================================================ + +CREATE FUNCTION star_observe_pm( + float8, float8, float8, float8, float8, float8, observer, timestamptz +) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION star_observe_pm(float8, float8, float8, float8, float8, float8, observer, timestamptz) IS + 'Observe a star with proper motion. Args: ra_hours, dec_deg, pm_ra_masyr (mu_alpha*cos(delta)), pm_dec_masyr, parallax_mas, rv_kms, observer, time. Hipparcos/Gaia convention for pm_ra.'; + +CREATE FUNCTION star_equatorial_pm( + float8, float8, float8, float8, float8, float8, timestamptz +) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION star_equatorial_pm(float8, float8, float8, float8, float8, float8, timestamptz) IS + 'Apparent RA/Dec of a star with proper motion. Args: ra_hours, dec_deg, pm_ra_masyr, pm_dec_masyr, parallax_mas, rv_kms, time. Distance from parallax if > 0.'; + + +-- ============================================================ +-- Light-time corrected observation functions +-- ============================================================ + +CREATE FUNCTION planet_observe_apparent(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_observe_apparent(int4, observer, timestamptz) IS + 'Observe a planet with single-iteration light-time correction. Body at retarded time, Earth at observation time. VSOP87.'; + +CREATE FUNCTION sun_observe_apparent(observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_observe_apparent(observer, timestamptz) IS + 'Observe the Sun with light-time correction (~8.3 min). VSOP87.'; + +CREATE FUNCTION small_body_observe_apparent(orbital_elements, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION small_body_observe_apparent(orbital_elements, observer, timestamptz) IS + 'Observe a comet/asteroid with single-iteration light-time correction. Kepler propagation at retarded time, Earth via VSOP87 at observation time.'; + + +-- ============================================================ +-- Light-time corrected equatorial functions +-- ============================================================ + +CREATE FUNCTION planet_equatorial_apparent(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_equatorial_apparent(int4, timestamptz) IS + 'Geocentric apparent RA/Dec of a planet with light-time correction. VSOP87.'; + +CREATE FUNCTION moon_equatorial_apparent(timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_equatorial_apparent(timestamptz) IS + 'Geocentric apparent RA/Dec of the Moon with light-time correction (~1.3 sec). ELP2000-82B.'; + +CREATE FUNCTION small_body_equatorial_apparent(orbital_elements, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION small_body_equatorial_apparent(orbital_elements, timestamptz) IS + 'Geocentric apparent RA/Dec of a comet/asteroid with light-time correction.'; + + +-- ============================================================ +-- DE ephemeris equatorial variants (STABLE) +-- ============================================================ + +CREATE FUNCTION planet_equatorial_de(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_equatorial_de(int4, timestamptz) IS + 'Geocentric apparent RA/Dec of a planet via JPL DE ephemeris (falls back to VSOP87 + equatorial).'; + +CREATE FUNCTION moon_equatorial_de(timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_equatorial_de(timestamptz) IS + 'Geocentric apparent RA/Dec of the Moon via JPL DE ephemeris (falls back to ELP2000-82B + equatorial).'; +-- pg_orrery 0.9.0 -> 0.10.0 migration +-- +-- Adds annual aberration to existing _apparent() functions, +-- 6 new _apparent_de() variants, equatorial angular separation +-- operator and cone predicate, and stellar annual parallax. + +-- ============================================================ +-- Equatorial angular distance and cone search +-- ============================================================ + +CREATE FUNCTION eq_angular_distance(equatorial, equatorial) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eq_angular_distance(equatorial, equatorial) IS + 'Angular separation in degrees between two equatorial positions. Vincenty formula (stable at 0 and 180 degrees).'; + +CREATE FUNCTION eq_within_cone(equatorial, equatorial, float8) RETURNS bool + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eq_within_cone(equatorial, equatorial, float8) IS + 'True if first position is within radius_deg of second position. Cosine shortcut for fast rejection.'; + +CREATE OPERATOR <-> ( + LEFTARG = equatorial, + RIGHTARG = equatorial, + FUNCTION = eq_angular_distance, + COMMUTATOR = <-> +); +COMMENT ON OPERATOR <-> (equatorial, equatorial) IS + 'Angular separation in degrees between two equatorial positions.'; + + +-- ============================================================ +-- DE apparent observation functions (STABLE, light-time + aberration) +-- ============================================================ + +CREATE FUNCTION planet_observe_apparent_de(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_observe_apparent_de(int4, observer, timestamptz) IS + 'Observe a planet with light-time correction and annual aberration via JPL DE (falls back to VSOP87).'; + +CREATE FUNCTION sun_observe_apparent_de(observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_observe_apparent_de(observer, timestamptz) IS + 'Observe the Sun with aberration via JPL DE (falls back to VSOP87).'; + +CREATE FUNCTION moon_observe_apparent_de(observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_observe_apparent_de(observer, timestamptz) IS + 'Observe the Moon with light-time correction and annual aberration via JPL DE (falls back to ELP2000-82B).'; + +CREATE FUNCTION planet_equatorial_apparent_de(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_equatorial_apparent_de(int4, timestamptz) IS + 'Geocentric apparent RA/Dec of a planet with light-time correction and annual aberration via JPL DE (falls back to VSOP87).'; + +CREATE FUNCTION moon_equatorial_apparent_de(timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_equatorial_apparent_de(timestamptz) IS + 'Geocentric apparent RA/Dec of the Moon with light-time correction and annual aberration via JPL DE (falls back to ELP2000-82B).'; + +CREATE FUNCTION small_body_observe_apparent_de(orbital_elements, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION small_body_observe_apparent_de(orbital_elements, observer, timestamptz) IS + 'Observe a comet/asteroid with light-time correction and annual aberration. Earth position via JPL DE (falls back to VSOP87).'; +-- pg_orrery 0.10.0 -> 0.11.0 migration +-- +-- Adds make_orbital_elements() constructors and +-- geocentric equatorial functions for planetary moons. + +-- ============================================================ +-- orbital_elements constructors +-- ============================================================ + +CREATE FUNCTION make_orbital_elements( + epoch_jd float8, q_au float8, e float8, + inc_rad float8, omega_rad float8, node_rad float8, + tp_jd float8, h_mag float8, g_slope float8 +) RETURNS orbital_elements + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION make_orbital_elements(float8,float8,float8,float8,float8,float8,float8,float8,float8) IS + 'Construct orbital_elements from 9 floats (angular elements in radians).'; + +CREATE FUNCTION make_orbital_elements_deg( + epoch_jd float8, q_au float8, e float8, + inc_deg float8, omega_deg float8, node_deg float8, + tp_jd float8, h_mag float8, g_slope float8 +) RETURNS orbital_elements + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION make_orbital_elements_deg(float8,float8,float8,float8,float8,float8,float8,float8,float8) IS + 'Construct orbital_elements from 9 floats (angular elements in degrees). Matches text I/O and most catalog column layouts.'; + + +-- ============================================================ +-- Planetary moon equatorial functions +-- ============================================================ + +CREATE FUNCTION galilean_equatorial(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION galilean_equatorial(int4, timestamptz) IS + 'Geometric geocentric RA/Dec of a Galilean moon (0=Io, 1=Europa, 2=Ganymede, 3=Callisto). L1.2 theory + VSOP87. No light-time or aberration correction.'; + +CREATE FUNCTION saturn_moon_equatorial(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION saturn_moon_equatorial(int4, timestamptz) IS + 'Geometric geocentric RA/Dec of a Saturn moon (0=Mimas..7=Hyperion). TASS17 theory + VSOP87. No light-time or aberration correction.'; + +CREATE FUNCTION uranus_moon_equatorial(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION uranus_moon_equatorial(int4, timestamptz) IS + 'Geometric geocentric RA/Dec of a Uranus moon (0=Miranda..4=Oberon). GUST86 theory + VSOP87. No light-time or aberration correction.'; + +CREATE FUNCTION mars_moon_equatorial(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION mars_moon_equatorial(int4, timestamptz) IS + 'Geometric geocentric RA/Dec of a Mars moon (0=Phobos, 1=Deimos). MarsSat theory + VSOP87. No light-time or aberration correction.'; +-- 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.'; + + +-- ============================================================ +-- v0.13.0: make_equatorial() constructor +-- ============================================================ + +CREATE FUNCTION make_equatorial(ra_hours float8, dec_deg float8, distance_km float8) + RETURNS equatorial + AS 'MODULE_PATHNAME', 'make_equatorial' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +COMMENT ON FUNCTION make_equatorial(float8, float8, float8) IS + 'Construct equatorial from RA (hours [0,24)), Dec (degrees [-90,90]), distance (km).'; + + +-- ============================================================ +-- v0.13.0: Rise/set prediction functions +-- ============================================================ + +CREATE FUNCTION planet_next_rise(body_id int4, obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_next_rise(int4, observer, timestamptz) IS + 'Next geometric rise time for a planet. Returns NULL if no rise within 7 days. body_id: 1=Mercury..8=Neptune.'; + +CREATE FUNCTION planet_next_set(body_id int4, obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_next_set(int4, observer, timestamptz) IS + 'Next geometric set time for a planet. Returns NULL if no set within 7 days. body_id: 1=Mercury..8=Neptune.'; + +CREATE FUNCTION sun_next_rise(obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_next_rise(observer, timestamptz) IS + 'Next geometric sunrise. Returns NULL if Sun does not rise within 7 days (polar night).'; + +CREATE FUNCTION sun_next_set(obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_next_set(observer, timestamptz) IS + 'Next geometric sunset. Returns NULL if Sun does not set within 7 days (midnight sun).'; + +CREATE FUNCTION moon_next_rise(obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_next_rise(observer, timestamptz) IS + 'Next geometric moonrise. Returns NULL if Moon does not rise within 7 days.'; + +CREATE FUNCTION moon_next_set(obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_next_set(observer, timestamptz) IS + 'Next geometric moonset. Returns NULL if Moon does not set within 7 days.'; + +CREATE FUNCTION sun_next_rise_refracted(obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_next_rise_refracted(observer, timestamptz) IS + 'Next refracted sunrise (-0.833 deg threshold: refraction + semidiameter). Earlier than geometric by ~4 min.'; + +CREATE FUNCTION sun_next_set_refracted(obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_next_set_refracted(observer, timestamptz) IS + 'Next refracted sunset (-0.833 deg threshold: refraction + semidiameter). Later than geometric by ~4 min.'; + + +-- ============================================================ +-- v0.14.0: Refracted planet/moon rise/set +-- ============================================================ + +CREATE FUNCTION planet_next_rise_refracted(body_id int4, obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_next_rise_refracted(int4, observer, timestamptz) IS + 'Next refracted rise time for a planet (-0.569 deg threshold: atmospheric refraction only). Earlier than geometric.'; + +CREATE FUNCTION planet_next_set_refracted(body_id int4, obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_next_set_refracted(int4, observer, timestamptz) IS + 'Next refracted set time for a planet (-0.569 deg threshold: atmospheric refraction only). Later than geometric.'; + +CREATE FUNCTION moon_next_rise_refracted(obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_next_rise_refracted(observer, timestamptz) IS + 'Next refracted moonrise (-0.833 deg threshold: refraction + semidiameter). Earlier than geometric.'; + +CREATE FUNCTION moon_next_set_refracted(obs observer, t timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_next_set_refracted(observer, timestamptz) IS + 'Next refracted moonset (-0.833 deg threshold: refraction + semidiameter). Later than geometric.'; + + +-- ============================================================ +-- v0.14.0: Constellation identification (Roman 1987, CDS VI/42) +-- ============================================================ + +CREATE FUNCTION constellation(eq equatorial) RETURNS text + AS 'MODULE_PATHNAME', 'constellation_from_equatorial' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION constellation(equatorial) IS + 'IAU constellation abbreviation (3 letters) from equatorial coordinates (Roman 1987).'; + +CREATE FUNCTION constellation(ra_hours float8, dec_deg float8) RETURNS text + AS 'MODULE_PATHNAME', 'constellation_from_radec' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION constellation(float8, float8) IS + 'IAU constellation from J2000 RA (hours [0,24)) and Dec (degrees [-90,90]).'; +-- pg_orrery 0.14.0 -> 0.15.0 migration +-- +-- Adds: constellation_full_name (1 function), +-- rise/set status diagnostics (3 functions). + +-- ============================================================ +-- Constellation full name lookup +-- ============================================================ + +CREATE FUNCTION constellation_full_name(abbr text) RETURNS text + AS 'MODULE_PATHNAME', 'constellation_full_name_from_abbr' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION constellation_full_name(text) IS + 'Full IAU constellation name from 3-letter abbreviation. Returns NULL for invalid abbreviation.'; + +-- ============================================================ +-- Rise/set status diagnostics +-- ============================================================ + +CREATE FUNCTION sun_rise_set_status(obs observer, t timestamptz) RETURNS text + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_rise_set_status(observer, timestamptz) IS + 'Classify Sun visibility: rises_and_sets, circumpolar, or never_rises.'; + +CREATE FUNCTION moon_rise_set_status(obs observer, t timestamptz) RETURNS text + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_rise_set_status(observer, timestamptz) IS + 'Classify Moon visibility: rises_and_sets, circumpolar, or never_rises.'; + +CREATE FUNCTION planet_rise_set_status(body_id int4, obs observer, t timestamptz) RETURNS text + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_rise_set_status(int4, observer, timestamptz) IS + 'Classify planet visibility: rises_and_sets, circumpolar, or never_rises. Body IDs 1-8 (Mercury-Neptune).'; +-- pg_orrery 0.15.0 -> 0.16.0: twilight, lunar phase, planet magnitude + +-- ============================================================ +-- Twilight functions (6) +-- ============================================================ + +CREATE FUNCTION sun_civil_dawn(observer, timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME', 'sun_civil_dawn' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_civil_dawn(observer, timestamptz) IS + 'Next civil dawn (Sun crosses -6 deg rising). Outdoor activities without artificial light.'; + +CREATE FUNCTION sun_civil_dusk(observer, timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME', 'sun_civil_dusk' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_civil_dusk(observer, timestamptz) IS + 'Next civil dusk (Sun crosses -6 deg setting). Artificial light needed.'; + +CREATE FUNCTION sun_nautical_dawn(observer, timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME', 'sun_nautical_dawn' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_nautical_dawn(observer, timestamptz) IS + 'Next nautical dawn (Sun crosses -12 deg rising). Horizon visible at sea.'; + +CREATE FUNCTION sun_nautical_dusk(observer, timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME', 'sun_nautical_dusk' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_nautical_dusk(observer, timestamptz) IS + 'Next nautical dusk (Sun crosses -12 deg setting). Horizon no longer visible at sea.'; + +CREATE FUNCTION sun_astronomical_dawn(observer, timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME', 'sun_astronomical_dawn' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_astronomical_dawn(observer, timestamptz) IS + 'Next astronomical dawn (Sun crosses -18 deg rising). Sky was fully dark.'; + +CREATE FUNCTION sun_astronomical_dusk(observer, timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME', 'sun_astronomical_dusk' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_astronomical_dusk(observer, timestamptz) IS + 'Next astronomical dusk (Sun crosses -18 deg setting). Sky becomes fully dark.'; + +-- ============================================================ +-- Lunar phase functions (4) +-- ============================================================ + +CREATE FUNCTION moon_phase_angle(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'moon_phase_angle' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_phase_angle(timestamptz) IS + 'Sun-Earth-Moon phase angle in degrees [0,360). 0=new, 90=first quarter, 180=full, 270=last quarter.'; + +CREATE FUNCTION moon_illumination(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'moon_illumination' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_illumination(timestamptz) IS + 'Illuminated fraction of the Moon disk [0.0, 1.0].'; + +CREATE FUNCTION moon_phase_name(timestamptz) RETURNS text + AS 'MODULE_PATHNAME', 'moon_phase_name' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_phase_name(timestamptz) IS + 'Moon phase name: new_moon, waxing_crescent, first_quarter, waxing_gibbous, full_moon, waning_gibbous, last_quarter, waning_crescent.'; + +CREATE FUNCTION moon_age(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'moon_age' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_age(timestamptz) IS + 'Days since last new moon [0, ~29.53), approximated from phase angle.'; + +-- ============================================================ +-- Planet magnitude (1) +-- ============================================================ + +CREATE FUNCTION planet_magnitude(int4, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'planet_magnitude' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_magnitude(int4, timestamptz) IS + 'Apparent visual magnitude of a planet (Mallama & Hilton 2018). Body IDs 1-8. Saturn ring tilt not modeled.'; +-- pg_orrery 0.16.0 -> 0.17.0: solar elongation, planet phase, satellite eclipse, +-- observing night quality, lunar libration + +-- ============================================================ +-- Solar elongation (1) +-- ============================================================ + +CREATE FUNCTION solar_elongation(int4, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'solar_elongation' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION solar_elongation(int4, timestamptz) IS + 'Sun-Earth-Planet angle in degrees [0, 180]. How far a planet appears from the Sun. Body IDs 1-8.'; + +-- ============================================================ +-- Planet phase fraction (1) +-- ============================================================ + +CREATE FUNCTION planet_phase(int4, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'planet_phase' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_phase(int4, timestamptz) IS + 'Illuminated fraction of a planet disk as seen from Earth [0.0, 1.0]. Body IDs 1-8.'; + +-- ============================================================ +-- Satellite eclipse prediction (4) +-- ============================================================ + +CREATE FUNCTION satellite_is_eclipsed(tle, timestamptz) RETURNS bool + AS 'MODULE_PATHNAME', 'satellite_is_eclipsed' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION satellite_is_eclipsed(tle, timestamptz) IS + 'True if the satellite is in Earth cylindrical shadow at the given time.'; + +CREATE FUNCTION satellite_next_eclipse_entry(tle, timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME', 'satellite_next_eclipse_entry' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION satellite_next_eclipse_entry(tle, timestamptz) IS + 'Next time the satellite enters Earth shadow (up to 7-day search). NULL if none found.'; + +CREATE FUNCTION satellite_next_eclipse_exit(tle, timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME', 'satellite_next_eclipse_exit' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION satellite_next_eclipse_exit(tle, timestamptz) IS + 'Next time the satellite exits Earth shadow (up to 7-day search). NULL if none found.'; + +CREATE FUNCTION satellite_eclipse_fraction(tle, timestamptz, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'satellite_eclipse_fraction' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION satellite_eclipse_fraction(tle, timestamptz, timestamptz) IS + 'Fraction of the given time window the satellite spends in eclipse [0.0, 1.0].'; + +-- ============================================================ +-- Observing night quality (1) +-- ============================================================ + +CREATE FUNCTION observing_night_quality(observer, timestamptz DEFAULT NOW()) +RETURNS text AS $$ +DECLARE + astro_dusk timestamptz; + astro_dawn timestamptz; + dark_hours float8; + illum float8; + moon_up bool; + score int := 100; +BEGIN + -- Astronomical darkness window + astro_dusk := sun_astronomical_dusk($1, $2); + IF astro_dusk IS NULL THEN + RETURN 'poor'; -- No astronomical darkness (polar summer) + END IF; + astro_dawn := sun_astronomical_dawn($1, astro_dusk); + IF astro_dawn IS NULL THEN + RETURN 'poor'; + END IF; + + dark_hours := extract(epoch FROM astro_dawn - astro_dusk) / 3600.0; + + -- Short dark window penalty + IF dark_hours < 2.0 THEN score := score - 40; + ELSIF dark_hours < 4.0 THEN score := score - 20; + ELSIF dark_hours < 6.0 THEN score := score - 10; + END IF; + + -- Moon illumination penalty + illum := moon_illumination(astro_dusk); + IF illum > 0.75 THEN + -- Check if Moon is above horizon during darkness + moon_up := topo_elevation(moon_observe($1, astro_dusk)) > 0 + OR topo_elevation(moon_observe($1, astro_dusk + (astro_dawn - astro_dusk) / 2)) > 0; + IF moon_up THEN + score := score - (illum * 30)::int; -- Up to -30 for full moon + END IF; + END IF; + + -- Classify + IF score >= 80 THEN RETURN 'excellent'; + ELSIF score >= 60 THEN RETURN 'good'; + ELSIF score >= 40 THEN RETURN 'fair'; + ELSE RETURN 'poor'; + END IF; +END; +$$ LANGUAGE plpgsql STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION observing_night_quality(observer, timestamptz) IS + 'Composite observing quality assessment: excellent/good/fair/poor based on darkness duration and Moon interference.'; + +-- ============================================================ +-- Lunar libration (5) +-- ============================================================ + +CREATE FUNCTION moon_libration_longitude(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'moon_libration_longitude' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_libration_longitude(timestamptz) IS + 'Optical libration in longitude (degrees, typically [-8, +8]). Meeus Ch. 53.'; + +CREATE FUNCTION moon_libration_latitude(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'moon_libration_latitude' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_libration_latitude(timestamptz) IS + 'Optical libration in latitude (degrees, typically [-7, +7]). Meeus Ch. 53.'; + +CREATE FUNCTION moon_libration_position_angle(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'moon_libration_position_angle' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_libration_position_angle(timestamptz) IS + 'Position angle of the Moon axis (degrees, [0, 360)). Meeus Ch. 53.'; + +CREATE FUNCTION moon_libration(timestamptz, + OUT l float8, OUT b float8, OUT p float8) RETURNS record + AS 'MODULE_PATHNAME', 'moon_libration' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_libration(timestamptz) IS + 'All three libration values: longitude (l), latitude (b), position angle (p) in degrees.'; + +CREATE FUNCTION moon_subsolar_longitude(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'moon_subsolar_longitude' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_subsolar_longitude(timestamptz) IS + 'Selenographic longitude of the sub-solar point (degrees, [0, 360)). Determines the lunar terminator position.'; +-- pg_orrery 0.17.0 -> 0.18.0: Saturn ring tilt, penumbral eclipse, +-- rise/set event windows, angular separation rate + +-- ============================================================ +-- Saturn ring tilt (1) +-- ============================================================ + +CREATE FUNCTION saturn_ring_tilt(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'saturn_ring_tilt' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION saturn_ring_tilt(timestamptz) IS + 'Sub-observer latitude B'' of Earth relative to Saturn ring plane (degrees, [-27, +27]). Uses IAU 2000 pole direction.'; + +-- ============================================================ +-- Penumbral eclipse prediction (4) +-- ============================================================ + +CREATE FUNCTION satellite_in_penumbra(tle, timestamptz) RETURNS bool + AS 'MODULE_PATHNAME', 'satellite_in_penumbra' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION satellite_in_penumbra(tle, timestamptz) IS + 'True if the satellite is in Earth penumbral shadow (partial sunlight) at the given time.'; + +CREATE FUNCTION satellite_shadow_state(tle, timestamptz) RETURNS text + AS 'MODULE_PATHNAME', 'satellite_shadow_state' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION satellite_shadow_state(tle, timestamptz) IS + 'Shadow state of satellite: ''sunlit'', ''penumbra'', or ''umbra''. Uses conical shadow model.'; + +CREATE FUNCTION satellite_next_penumbra_entry(tle, timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME', 'satellite_next_penumbra_entry' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION satellite_next_penumbra_entry(tle, timestamptz) IS + 'Next time the satellite enters Earth penumbral shadow (up to 7-day search). NULL if none found.'; + +CREATE FUNCTION satellite_next_penumbra_exit(tle, timestamptz) RETURNS timestamptz + AS 'MODULE_PATHNAME', 'satellite_next_penumbra_exit' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION satellite_next_penumbra_exit(tle, timestamptz) IS + 'Next time the satellite exits Earth penumbral shadow (up to 7-day search). NULL if none found.'; + +-- ============================================================ +-- Rise/set event windows (3 SRFs) +-- ============================================================ + +CREATE FUNCTION planet_rise_set_events( + body_id int4, observer, start timestamptz, stop timestamptz, + refracted bool DEFAULT false +) RETURNS TABLE(event_time timestamptz, event_type text) + AS 'MODULE_PATHNAME', 'planet_rise_set_events' + LANGUAGE C STABLE STRICT PARALLEL SAFE + ROWS 10; +COMMENT ON FUNCTION planet_rise_set_events(int4, observer, timestamptz, timestamptz, bool) IS + 'All rise and set events for a planet within a time window. Returns TABLE(event_time, event_type). Max 366-day window.'; + +CREATE FUNCTION sun_rise_set_events( + observer, start timestamptz, stop timestamptz, + refracted bool DEFAULT false +) RETURNS TABLE(event_time timestamptz, event_type text) + AS 'MODULE_PATHNAME', 'sun_rise_set_events' + LANGUAGE C STABLE STRICT PARALLEL SAFE + ROWS 10; +COMMENT ON FUNCTION sun_rise_set_events(observer, timestamptz, timestamptz, bool) IS + 'All rise and set events for the Sun within a time window. Returns TABLE(event_time, event_type). Max 366-day window.'; + +CREATE FUNCTION moon_rise_set_events( + observer, start timestamptz, stop timestamptz, + refracted bool DEFAULT false +) RETURNS TABLE(event_time timestamptz, event_type text) + AS 'MODULE_PATHNAME', 'moon_rise_set_events' + LANGUAGE C STABLE STRICT PARALLEL SAFE + ROWS 10; +COMMENT ON FUNCTION moon_rise_set_events(observer, timestamptz, timestamptz, bool) IS + 'All rise and set events for the Moon within a time window. Returns TABLE(event_time, event_type). Max 366-day window.'; + +-- ============================================================ +-- Angular separation rate (2) +-- ============================================================ + +CREATE FUNCTION eq_angular_rate( + equatorial, equatorial, equatorial, equatorial, float8 +) RETURNS float8 + AS 'MODULE_PATHNAME', 'eq_angular_rate' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eq_angular_rate(equatorial, equatorial, equatorial, equatorial, float8) IS + 'Rate of change of angular separation (deg/hr). Args: pos1_t0, pos2_t0, pos1_t1, pos2_t1, dt_seconds. Positive = separating, negative = approaching.'; + +CREATE FUNCTION planet_angular_rate(int4, int4, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'planet_angular_rate' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_angular_rate(int4, int4, timestamptz) IS + 'Rate of angular separation change between two bodies (deg/hr). Body IDs: 0=Sun, 1-8=planets, 10=Moon. Uses 1-minute finite difference.'; +-- pg_orrery 0.18.0 -> 0.19.0: Sun almanac SRF, conjunction detection, +-- penumbral fraction, physical libration + +-- ============================================================ +-- Sun almanac events SRF (1) +-- ============================================================ + +CREATE FUNCTION sun_almanac_events( + observer, start timestamptz, stop timestamptz, + refracted bool DEFAULT false +) RETURNS TABLE(event_time timestamptz, event_type text) + AS 'MODULE_PATHNAME', 'sun_almanac_events' + LANGUAGE C STABLE STRICT PARALLEL SAFE + ROWS 50; +COMMENT ON FUNCTION sun_almanac_events(observer, timestamptz, timestamptz, bool) IS + 'All Sun events (rise, set, civil/nautical/astronomical dawn and dusk) within a time window, sorted chronologically. Replaces chained individual twilight queries. Max 366-day window.'; + +-- ============================================================ +-- Conjunction detection SRF (1) +-- ============================================================ + +CREATE FUNCTION planet_conjunctions( + int4, int4, timestamptz, timestamptz, + max_separation float8 DEFAULT 10.0 +) RETURNS TABLE(conjunction_time timestamptz, separation_deg float8) + AS 'MODULE_PATHNAME', 'planet_conjunctions' + LANGUAGE C STABLE STRICT PARALLEL SAFE + ROWS 10; +COMMENT ON FUNCTION planet_conjunctions(int4, int4, timestamptz, timestamptz, float8) IS + 'Finds conjunctions (angular separation minima) between two solar system bodies. Body IDs: 0=Sun, 1-8=planets, 10=Moon. max_separation filters results (degrees, default 10). Max 3660-day (10-year) window.'; + +-- ============================================================ +-- Penumbral fraction (1) +-- ============================================================ + +CREATE FUNCTION satellite_penumbral_fraction(tle, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'satellite_penumbral_fraction' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION satellite_penumbral_fraction(tle, timestamptz) IS + 'Continuous shadow depth: 0.0 = full sunlight, 1.0 = full umbra. Linear interpolation in penumbral zone.'; + +-- ============================================================ +-- Physical libration (1) +-- ============================================================ + +CREATE FUNCTION moon_physical_libration( + timestamptz, + OUT tau float8, OUT rho float8 +) RETURNS record + AS 'MODULE_PATHNAME', 'moon_physical_libration' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_physical_libration(timestamptz) IS + 'Physical libration corrections (Meeus p. 373): tau = longitude correction, rho = latitude correction (both in degrees, typically |value| < 0.1).'; +-- pg_orrery 0.19.0 -> 0.20.0: Lagrange point support +-- CR3BP equilibrium positions for Sun-planet, Earth-Moon, and planetary moon systems. + +-- ============================================================ +-- Sun-planet Lagrange functions (5) +-- ============================================================ + +CREATE FUNCTION lagrange_heliocentric(int4, int4, timestamptz) RETURNS heliocentric + AS 'MODULE_PATHNAME', 'lagrange_heliocentric' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_heliocentric(int4, int4, timestamptz) IS + 'Heliocentric ecliptic J2000 position of a Sun-planet Lagrange point. body_id: 1-8 (Mercury-Neptune), point_id: 1-5 (L1-L5).'; + +CREATE FUNCTION lagrange_observe(int4, int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'lagrange_observe' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_observe(int4, int4, observer, timestamptz) IS + 'Observe a Sun-planet Lagrange point from a ground station. body_id: 1-8, point_id: 1-5.'; + +CREATE FUNCTION lagrange_equatorial(int4, int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'lagrange_equatorial' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_equatorial(int4, int4, timestamptz) IS + 'Geocentric RA/Dec of a Sun-planet Lagrange point. body_id: 1-8, point_id: 1-5.'; + +CREATE FUNCTION lagrange_distance(int4, int4, heliocentric, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'lagrange_distance' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_distance(int4, int4, heliocentric, timestamptz) IS + 'Distance (AU) from a heliocentric position to a Sun-planet Lagrange point.'; + +CREATE FUNCTION lagrange_distance_oe(int4, int4, orbital_elements, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'lagrange_distance_oe' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_distance_oe(int4, int4, orbital_elements, timestamptz) IS + 'Distance (AU) from an asteroid/comet (orbital_elements) to a Sun-planet Lagrange point.'; + +-- ============================================================ +-- Earth-Moon Lagrange functions (2) +-- ============================================================ + +CREATE FUNCTION lunar_lagrange_observe(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'lunar_lagrange_observe' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lunar_lagrange_observe(int4, observer, timestamptz) IS + 'Observe an Earth-Moon Lagrange point. point_id: 1-5 (L1-L5).'; + +CREATE FUNCTION lunar_lagrange_equatorial(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'lunar_lagrange_equatorial' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lunar_lagrange_equatorial(int4, timestamptz) IS + 'Geocentric RA/Dec of an Earth-Moon Lagrange point. point_id: 1-5 (L1-L5).'; + +-- ============================================================ +-- Planetary moon Lagrange functions (8) +-- ============================================================ + +CREATE FUNCTION galilean_lagrange_observe(int4, int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'galilean_lagrange_observe' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION galilean_lagrange_observe(int4, int4, observer, timestamptz) IS + 'Observe a Jupiter-Galilean moon Lagrange point. body_id: 0-3 (Io-Callisto), point_id: 1-5.'; + +CREATE FUNCTION galilean_lagrange_equatorial(int4, int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'galilean_lagrange_equatorial' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION galilean_lagrange_equatorial(int4, int4, timestamptz) IS + 'Geocentric RA/Dec of a Jupiter-Galilean moon Lagrange point. body_id: 0-3, point_id: 1-5.'; + +CREATE FUNCTION saturn_moon_lagrange_observe(int4, int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'saturn_moon_lagrange_observe' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION saturn_moon_lagrange_observe(int4, int4, observer, timestamptz) IS + 'Observe a Saturn moon Lagrange point. body_id: 0-7 (Mimas-Hyperion), point_id: 1-5.'; + +CREATE FUNCTION saturn_moon_lagrange_equatorial(int4, int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'saturn_moon_lagrange_equatorial' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION saturn_moon_lagrange_equatorial(int4, int4, timestamptz) IS + 'Geocentric RA/Dec of a Saturn moon Lagrange point. body_id: 0-7, point_id: 1-5.'; + +CREATE FUNCTION uranus_moon_lagrange_observe(int4, int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'uranus_moon_lagrange_observe' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION uranus_moon_lagrange_observe(int4, int4, observer, timestamptz) IS + 'Observe a Uranus moon Lagrange point. body_id: 0-4 (Miranda-Oberon), point_id: 1-5.'; + +CREATE FUNCTION uranus_moon_lagrange_equatorial(int4, int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'uranus_moon_lagrange_equatorial' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION uranus_moon_lagrange_equatorial(int4, int4, timestamptz) IS + 'Geocentric RA/Dec of a Uranus moon Lagrange point. body_id: 0-4, point_id: 1-5.'; + +CREATE FUNCTION mars_moon_lagrange_observe(int4, int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'mars_moon_lagrange_observe' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION mars_moon_lagrange_observe(int4, int4, observer, timestamptz) IS + 'Observe a Mars moon Lagrange point. body_id: 0-1 (Phobos-Deimos), point_id: 1-5.'; + +CREATE FUNCTION mars_moon_lagrange_equatorial(int4, int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'mars_moon_lagrange_equatorial' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION mars_moon_lagrange_equatorial(int4, int4, timestamptz) IS + 'Geocentric RA/Dec of a Mars moon Lagrange point. body_id: 0-1, point_id: 1-5.'; + +-- ============================================================ +-- Hill radius / zone / convenience (5) +-- ============================================================ + +CREATE FUNCTION hill_radius(int4, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'hill_radius_func' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION hill_radius(int4, timestamptz) IS + 'Hill sphere radius (AU) for a Sun-planet system. body_id: 1-8.'; + +CREATE FUNCTION hill_radius_lunar(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'hill_radius_lunar' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION hill_radius_lunar(timestamptz) IS + 'Hill sphere radius (AU) for the Earth-Moon system.'; + +CREATE FUNCTION lagrange_zone_radius(int4, int4, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'lagrange_zone_radius_func' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_zone_radius(int4, int4, timestamptz) IS + 'Approximate libration zone radius (AU) for a Sun-planet Lagrange point.'; + +CREATE FUNCTION lagrange_mass_ratio(int4) RETURNS float8 + AS 'MODULE_PATHNAME', 'lagrange_mass_ratio' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_mass_ratio(int4) IS + 'CR3BP mass parameter mu = M_planet / (M_sun + M_planet) for debugging. body_id: 1-8.'; + +CREATE FUNCTION lagrange_point_name(int4) RETURNS text + AS 'MODULE_PATHNAME', 'lagrange_point_name' + LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_point_name(int4) IS + 'Human-readable name for a Lagrange point ID (1->''L1'', ..., 5->''L5'').'; + +-- ============================================================ +-- DE variant functions (17) -- STABLE +-- ============================================================ + +CREATE FUNCTION lagrange_heliocentric_de(int4, int4, timestamptz) RETURNS heliocentric + AS 'MODULE_PATHNAME', 'lagrange_heliocentric_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_heliocentric_de(int4, int4, timestamptz) IS + 'DE variant of lagrange_heliocentric(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION lagrange_observe_de(int4, int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'lagrange_observe_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_observe_de(int4, int4, observer, timestamptz) IS + 'DE variant of lagrange_observe(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION lagrange_equatorial_de(int4, int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'lagrange_equatorial_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_equatorial_de(int4, int4, timestamptz) IS + 'DE variant of lagrange_equatorial(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION lagrange_distance_de(int4, int4, heliocentric, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'lagrange_distance_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_distance_de(int4, int4, heliocentric, timestamptz) IS + 'DE variant of lagrange_distance(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION lagrange_distance_oe_de(int4, int4, orbital_elements, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'lagrange_distance_oe_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_distance_oe_de(int4, int4, orbital_elements, timestamptz) IS + 'DE variant of lagrange_distance_oe(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION lunar_lagrange_observe_de(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'lunar_lagrange_observe_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lunar_lagrange_observe_de(int4, observer, timestamptz) IS + 'DE variant of lunar_lagrange_observe(). Falls back to ELP2000-82B if DE unavailable.'; + +CREATE FUNCTION lunar_lagrange_equatorial_de(int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'lunar_lagrange_equatorial_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lunar_lagrange_equatorial_de(int4, timestamptz) IS + 'DE variant of lunar_lagrange_equatorial(). Falls back to ELP2000-82B if DE unavailable.'; + +CREATE FUNCTION galilean_lagrange_observe_de(int4, int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'galilean_lagrange_observe_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION galilean_lagrange_observe_de(int4, int4, observer, timestamptz) IS + 'DE variant of galilean_lagrange_observe(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION galilean_lagrange_equatorial_de(int4, int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'galilean_lagrange_equatorial_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION galilean_lagrange_equatorial_de(int4, int4, timestamptz) IS + 'DE variant of galilean_lagrange_equatorial(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION saturn_moon_lagrange_observe_de(int4, int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'saturn_moon_lagrange_observe_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION saturn_moon_lagrange_observe_de(int4, int4, observer, timestamptz) IS + 'DE variant of saturn_moon_lagrange_observe(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION saturn_moon_lagrange_equatorial_de(int4, int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'saturn_moon_lagrange_equatorial_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION saturn_moon_lagrange_equatorial_de(int4, int4, timestamptz) IS + 'DE variant of saturn_moon_lagrange_equatorial(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION uranus_moon_lagrange_observe_de(int4, int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'uranus_moon_lagrange_observe_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION uranus_moon_lagrange_observe_de(int4, int4, observer, timestamptz) IS + 'DE variant of uranus_moon_lagrange_observe(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION uranus_moon_lagrange_equatorial_de(int4, int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'uranus_moon_lagrange_equatorial_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION uranus_moon_lagrange_equatorial_de(int4, int4, timestamptz) IS + 'DE variant of uranus_moon_lagrange_equatorial(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION mars_moon_lagrange_observe_de(int4, int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME', 'mars_moon_lagrange_observe_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION mars_moon_lagrange_observe_de(int4, int4, observer, timestamptz) IS + 'DE variant of mars_moon_lagrange_observe(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION mars_moon_lagrange_equatorial_de(int4, int4, timestamptz) RETURNS equatorial + AS 'MODULE_PATHNAME', 'mars_moon_lagrange_equatorial_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION mars_moon_lagrange_equatorial_de(int4, int4, timestamptz) IS + 'DE variant of mars_moon_lagrange_equatorial(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION hill_radius_de(int4, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'hill_radius_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION hill_radius_de(int4, timestamptz) IS + 'DE variant of hill_radius(). Falls back to VSOP87 if DE unavailable.'; + +CREATE FUNCTION lagrange_zone_radius_de(int4, int4, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME', 'lagrange_zone_radius_de' + LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lagrange_zone_radius_de(int4, int4, timestamptz) IS + 'DE variant of lagrange_zone_radius(). Falls back to VSOP87 if DE unavailable.'; diff --git a/src/de_funcs.c b/src/de_funcs.c index d588cf8..d289833 100644 --- a/src/de_funcs.c +++ b/src/de_funcs.c @@ -33,6 +33,7 @@ #include "tass17.h" #include "gust86.h" #include "marssat.h" +#include "lagrange.h" #include #include @@ -62,6 +63,25 @@ PG_FUNCTION_INFO_V1(uranus_moon_equatorial_de); PG_FUNCTION_INFO_V1(mars_moon_equatorial_de); PG_FUNCTION_INFO_V1(pg_orrery_ephemeris_info); +/* Lagrange DE variants */ +PG_FUNCTION_INFO_V1(lagrange_heliocentric_de); +PG_FUNCTION_INFO_V1(lagrange_observe_de); +PG_FUNCTION_INFO_V1(lagrange_equatorial_de); +PG_FUNCTION_INFO_V1(lagrange_distance_de); +PG_FUNCTION_INFO_V1(lagrange_distance_oe_de); +PG_FUNCTION_INFO_V1(lunar_lagrange_observe_de); +PG_FUNCTION_INFO_V1(lunar_lagrange_equatorial_de); +PG_FUNCTION_INFO_V1(galilean_lagrange_observe_de); +PG_FUNCTION_INFO_V1(galilean_lagrange_equatorial_de); +PG_FUNCTION_INFO_V1(saturn_moon_lagrange_observe_de); +PG_FUNCTION_INFO_V1(saturn_moon_lagrange_equatorial_de); +PG_FUNCTION_INFO_V1(uranus_moon_lagrange_observe_de); +PG_FUNCTION_INFO_V1(uranus_moon_lagrange_equatorial_de); +PG_FUNCTION_INFO_V1(mars_moon_lagrange_observe_de); +PG_FUNCTION_INFO_V1(mars_moon_lagrange_equatorial_de); +PG_FUNCTION_INFO_V1(hill_radius_de); +PG_FUNCTION_INFO_V1(lagrange_zone_radius_de); + /* ================================================================ * planet_heliocentric_de(body_id int, timestamptz) -> heliocentric @@ -1306,3 +1326,951 @@ pg_orrery_ephemeris_info(PG_FUNCTION_ARGS) tuple = heap_form_tuple(tupdesc, values, nulls); PG_RETURN_DATUM(HeapTupleGetDatum(tuple)); } + + +/* ================================================================ + * Lagrange point functions with DE ephemeris + * + * DE variants of the Lagrange functions in lagrange_funcs.c. + * Each uses DE for planet/Earth positions where possible, + * falling back to VSOP87/ELP2000-82B on any DE failure. + * Rule 7 always holds: both target and Earth from the same provider. + * ================================================================ + */ + + +/* + * Validate body_id and point_id for Sun-planet Lagrange functions. + * Duplicated from lagrange_funcs.c (no cross-TU symbols). + */ +static void +validate_lagrange_args_de(int body_id, int point_id, const char *func_name) +{ + if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("%s: body_id %d must be 1-8 (Mercury-Neptune)", + func_name, body_id))); + + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("%s: point_id %d must be 1-5 (L1-L5)", + func_name, point_id))); +} + + +/* + * Compute Sun-planet Lagrange point using DE for planet position+velocity. + * Falls back to VSOP87 if DE unavailable (rule 7: consistent provider). + * + * Returns 0 on success, -1 on solver failure. + * Sets *used_de to indicate which provider was used. + */ +static int +compute_sun_planet_lagrange_de(int body_id, int point_id, double jd, + double result[3], bool *used_de) +{ + double planet_xyz[6]; + double sun[3] = {0.0, 0.0, 0.0}; + double planet_vel[3]; + double ratio, mu; + + /* Try DE for planet position */ + if (eph_de_planet(body_id, jd, planet_xyz)) + { + /* Velocity via central difference (DE provides position only) */ + double pos_fwd[6], pos_bwd[6]; + double dt = 0.01; /* days */ + + 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) + { + planet_vel[0] = (pos_fwd[0] - pos_bwd[0]) / (2.0 * dt); + planet_vel[1] = (pos_fwd[1] - pos_bwd[1]) / (2.0 * dt); + planet_vel[2] = (pos_fwd[2] - pos_bwd[2]) / (2.0 * dt); + *used_de = true; + } + else + { + /* DE boundary straddled — fall through to VSOP87 */ + goto vsop87_fallback; + } + } + else + { +vsop87_fallback: + { + int vsop_body = body_id - 1; + + if (eph_de_available()) + ereport(NOTICE, + (errmsg("DE ephemeris unavailable for this query, falling back to VSOP87"))); + + GetVsop87Coor(jd, vsop_body, planet_xyz); + /* VSOP87 provides analytic velocity in xyz[3..5] */ + planet_vel[0] = planet_xyz[3]; + planet_vel[1] = planet_xyz[4]; + planet_vel[2] = planet_xyz[5]; + *used_de = false; + } + } + + ratio = sun_planet_ratio(body_id); + mu = mu_from_ratio(ratio); + + return lagrange_position(sun, planet_xyz, planet_vel, mu, point_id, result); +} + + +/* + * Compute Earth-Moon Lagrange point using DE for Moon position. + * Falls back to ELP2000-82B if DE unavailable. + * + * Moon velocity always via central difference (neither DE nor ELP + * provide direct velocity). Result is geocentric ecliptic J2000. + */ +static int +compute_lunar_lagrange_de(int point_id, double jd, double result_geo[3], + bool *used_de) +{ + double moon_pos[3], moon_fwd[3], moon_bwd[3]; + double moon_vel[3]; + double earth[3] = {0.0, 0.0, 0.0}; + double mu; + double dt = 0.001; /* days */ + + if (eph_de_moon(jd, moon_pos)) + { + bool got_fwd = eph_de_moon(jd + dt, moon_fwd); + bool got_bwd = eph_de_moon(jd - dt, moon_bwd); + + if (got_fwd && got_bwd) + { + *used_de = true; + } + else + { + /* DE boundary straddled */ + goto elp_fallback; + } + } + else + { +elp_fallback: + if (eph_de_available()) + ereport(NOTICE, + (errmsg("DE ephemeris unavailable, falling back to ELP2000-82B"))); + + GetElp82bCoor(jd, moon_pos); + GetElp82bCoor(jd + dt, moon_fwd); + GetElp82bCoor(jd - dt, moon_bwd); + *used_de = false; + } + + moon_vel[0] = (moon_fwd[0] - moon_bwd[0]) / (2.0 * dt); + moon_vel[1] = (moon_fwd[1] - moon_bwd[1]) / (2.0 * dt); + moon_vel[2] = (moon_fwd[2] - moon_bwd[2]) / (2.0 * dt); + + mu = mu_from_ratio(EARTH_MOON_EMRAT); + + return lagrange_position(earth, moon_pos, moon_vel, mu, point_id, + result_geo); +} + + +/* + * Observe a planetary moon Lagrange point using DE for parent planet + * and Earth positions. Falls back to VSOP87 if DE unavailable. + * + * lp_rel[3]: L-point offset relative to parent (ecliptic J2000, AU) + * parent_body_id: pg_orrery body ID (5=Jupiter, 6=Saturn, etc.) + */ +static void +observe_pmoon_lagrange_de(const double lp_rel[3], int parent_body_id, + double jd, const pg_observer *obs, + pg_topocentric *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); + } + + geo_ecl[0] = (parent_xyz[0] + lp_rel[0]) - earth_xyz[0]; + geo_ecl[1] = (parent_xyz[1] + lp_rel[1]) - earth_xyz[1]; + geo_ecl[2] = (parent_xyz[2] + lp_rel[2]) - earth_xyz[2]; + + observe_from_geocentric(geo_ecl, jd, obs, result); +} + + +static void +equatorial_pmoon_lagrange_de(const double lp_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; + + 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); + } + + geo_ecl[0] = (parent_xyz[0] + lp_rel[0]) - earth_xyz[0]; + geo_ecl[1] = (parent_xyz[1] + lp_rel[1]) - earth_xyz[1]; + geo_ecl[2] = (parent_xyz[2] + lp_rel[2]) - earth_xyz[2]; + + geocentric_to_equatorial(geo_ecl, jd, result); +} + + +/* + * Compute planetary moon Lagrange offset (no cross-TU call). + * Duplicated from lagrange_funcs.c. + */ +static int +compute_planetary_moon_lagrange_de(const double moon_rel[3], + const double moon_vel[3], + char family, int moon_id, + int point_id, + double lp_rel[3]) +{ + double planet_origin[3] = {0.0, 0.0, 0.0}; + double ratio, mu; + + ratio = planet_moon_ratio(family, moon_id); + if (ratio < 0.0) + return -1; + + mu = mu_from_ratio(ratio); + + return lagrange_position(planet_origin, moon_rel, moon_vel, mu, + point_id, lp_rel); +} + + +/* ================================================================ + * 1. lagrange_heliocentric_de(body_id, point_id, timestamptz) + * -> heliocentric + * + * DE variant of lagrange_heliocentric(). STABLE. + * ================================================================ + */ +Datum +lagrange_heliocentric_de(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + int64 ts = PG_GETARG_INT64(2); + double jd; + double lp[3]; + bool used_de; + pg_heliocentric *result; + + validate_lagrange_args_de(body_id, point_id, "lagrange_heliocentric_de"); + + jd = timestamptz_to_jd(ts); + + if (compute_sun_planet_lagrange_de(body_id, point_id, jd, lp, &used_de) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("lagrange_heliocentric_de: solver failed for body %d L%d", + body_id, point_id))); + + result = (pg_heliocentric *) palloc(sizeof(pg_heliocentric)); + result->x = lp[0]; + result->y = lp[1]; + result->z = lp[2]; + + PG_RETURN_POINTER(result); +} + + +/* ================================================================ + * 2. lagrange_observe_de(body_id, point_id, observer, timestamptz) + * -> topocentric + * + * DE variant of lagrange_observe(). STABLE. + * Rule 7: L-point + Earth from same provider. + * ================================================================ + */ +Datum +lagrange_observe_de(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2); + int64 ts = PG_GETARG_INT64(3); + double jd; + double lp[3], earth_xyz[6], geo_ecl[3]; + bool used_de; + pg_topocentric *result; + + validate_lagrange_args_de(body_id, point_id, "lagrange_observe_de"); + + jd = timestamptz_to_jd(ts); + + if (compute_sun_planet_lagrange_de(body_id, point_id, jd, lp, &used_de) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("lagrange_observe_de: solver failed"))); + + /* Earth from same provider as the L-point (rule 7) */ + if (used_de && eph_de_earth(jd, earth_xyz)) + { + /* DE succeeded for both */ + } + else + { + GetVsop87Coor(jd, 2, earth_xyz); + } + + geo_ecl[0] = lp[0] - earth_xyz[0]; + geo_ecl[1] = lp[1] - earth_xyz[1]; + geo_ecl[2] = lp[2] - earth_xyz[2]; + + result = (pg_topocentric *) palloc(sizeof(pg_topocentric)); + observe_from_geocentric(geo_ecl, jd, obs, result); + + PG_RETURN_POINTER(result); +} + + +/* ================================================================ + * 3. lagrange_equatorial_de(body_id, point_id, timestamptz) + * -> equatorial + * + * DE variant of lagrange_equatorial(). STABLE. + * ================================================================ + */ +Datum +lagrange_equatorial_de(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + int64 ts = PG_GETARG_INT64(2); + double jd; + double lp[3], earth_xyz[6], geo_ecl[3]; + bool used_de; + pg_equatorial *result; + + validate_lagrange_args_de(body_id, point_id, "lagrange_equatorial_de"); + + jd = timestamptz_to_jd(ts); + + if (compute_sun_planet_lagrange_de(body_id, point_id, jd, lp, &used_de) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("lagrange_equatorial_de: solver failed"))); + + if (used_de && eph_de_earth(jd, earth_xyz)) + { + /* DE succeeded */ + } + else + { + GetVsop87Coor(jd, 2, earth_xyz); + } + + geo_ecl[0] = lp[0] - earth_xyz[0]; + geo_ecl[1] = lp[1] - earth_xyz[1]; + geo_ecl[2] = lp[2] - earth_xyz[2]; + + result = (pg_equatorial *) palloc(sizeof(pg_equatorial)); + geocentric_to_equatorial(geo_ecl, jd, result); + + PG_RETURN_POINTER(result); +} + + +/* ================================================================ + * 4. lagrange_distance_de(body_id, point_id, heliocentric, timestamptz) + * -> float8 + * + * Distance from a heliocentric position to a Sun-planet L-point (AU). + * Uses DE for planet position. STABLE. + * ================================================================ + */ +Datum +lagrange_distance_de(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + pg_heliocentric *pos = (pg_heliocentric *) PG_GETARG_POINTER(2); + int64 ts = PG_GETARG_INT64(3); + double jd; + double lp[3]; + double dx, dy, dz, dist; + bool used_de; + + validate_lagrange_args_de(body_id, point_id, "lagrange_distance_de"); + + jd = timestamptz_to_jd(ts); + + if (compute_sun_planet_lagrange_de(body_id, point_id, jd, lp, &used_de) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("lagrange_distance_de: solver failed"))); + + dx = pos->x - lp[0]; + dy = pos->y - lp[1]; + dz = pos->z - lp[2]; + dist = sqrt(dx*dx + dy*dy + dz*dz); + + PG_RETURN_FLOAT8(dist); +} + + +/* ================================================================ + * 5. lagrange_distance_oe_de(body_id, point_id, orbital_elements, timestamptz) + * -> float8 + * + * Distance from an asteroid/comet to a Sun-planet L-point (AU). + * Uses DE for L-point computation. STABLE. + * ================================================================ + */ +Datum +lagrange_distance_oe_de(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + pg_orbital_elements *oe = (pg_orbital_elements *) PG_GETARG_POINTER(2); + int64 ts = PG_GETARG_INT64(3); + double jd; + double lp[3], body_pos[3]; + double dx, dy, dz, dist; + bool used_de; + + validate_lagrange_args_de(body_id, point_id, "lagrange_distance_oe_de"); + + jd = timestamptz_to_jd(ts); + + if (compute_sun_planet_lagrange_de(body_id, point_id, jd, lp, &used_de) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("lagrange_distance_oe_de: solver failed"))); + + kepler_position(oe->q, oe->e, oe->inc, oe->arg_peri, oe->raan, + oe->tp, jd, body_pos); + + dx = body_pos[0] - lp[0]; + dy = body_pos[1] - lp[1]; + dz = body_pos[2] - lp[2]; + dist = sqrt(dx*dx + dy*dy + dz*dz); + + PG_RETURN_FLOAT8(dist); +} + + +/* ================================================================ + * 6. lunar_lagrange_observe_de(point_id, observer, timestamptz) + * -> topocentric + * + * Earth-Moon L-point using DE Moon position. STABLE. + * ================================================================ + */ +Datum +lunar_lagrange_observe_de(PG_FUNCTION_ARGS) +{ + int32 point_id = PG_GETARG_INT32(0); + pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(1); + int64 ts = PG_GETARG_INT64(2); + double jd; + double lp_geo[3]; + bool used_de; + pg_topocentric *result; + + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("lunar_lagrange_observe_de: point_id %d must be 1-5", + point_id))); + + jd = timestamptz_to_jd(ts); + + if (compute_lunar_lagrange_de(point_id, jd, lp_geo, &used_de) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("lunar_lagrange_observe_de: solver failed"))); + + /* lp_geo is already geocentric ecliptic J2000 (AU) */ + result = (pg_topocentric *) palloc(sizeof(pg_topocentric)); + observe_from_geocentric(lp_geo, jd, obs, result); + + PG_RETURN_POINTER(result); +} + + +/* ================================================================ + * 7. lunar_lagrange_equatorial_de(point_id, timestamptz) -> equatorial + * + * Earth-Moon L-point using DE Moon position. STABLE. + * ================================================================ + */ +Datum +lunar_lagrange_equatorial_de(PG_FUNCTION_ARGS) +{ + int32 point_id = PG_GETARG_INT32(0); + int64 ts = PG_GETARG_INT64(1); + double jd; + double lp_geo[3]; + bool used_de; + pg_equatorial *result; + + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("lunar_lagrange_equatorial_de: point_id %d must be 1-5", + point_id))); + + jd = timestamptz_to_jd(ts); + + if (compute_lunar_lagrange_de(point_id, jd, lp_geo, &used_de) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("lunar_lagrange_equatorial_de: solver failed"))); + + result = (pg_equatorial *) palloc(sizeof(pg_equatorial)); + geocentric_to_equatorial(lp_geo, jd, result); + + PG_RETURN_POINTER(result); +} + + +/* ================================================================ + * Planetary moon Lagrange points with DE parent positions + * + * Moon-theory offset (L1.2, TASS17, GUST86, MarsSat) computed + * relative to parent, then Lagrange solver applied. Parent planet + * and Earth positions from DE (with VSOP87 fallback). + * ================================================================ + */ + + +/* ── 8. Galilean Lagrange observe DE ─────────────────────── */ + +Datum +galilean_lagrange_observe_de(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2); + int64 ts = PG_GETARG_INT64(3); + double jd; + double moon_xyz[3], moon_vel[3], lp_rel[3]; + pg_topocentric *result; + + if (body_id < L12_IO || body_id > L12_CALLISTO) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("galilean_lagrange_observe_de: body_id %d must be 0-3", + body_id))); + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("galilean_lagrange_observe_de: point_id %d must be 1-5", + point_id))); + + jd = timestamptz_to_jd(ts); + GetL12Coor(jd, body_id, moon_xyz, moon_vel); + + if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 'g', body_id, + point_id, lp_rel) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("galilean_lagrange_observe_de: solver failed"))); + + result = (pg_topocentric *) palloc(sizeof(pg_topocentric)); + observe_pmoon_lagrange_de(lp_rel, BODY_JUPITER, jd, obs, result); + + PG_RETURN_POINTER(result); +} + + +/* ── 9. Galilean Lagrange equatorial DE ──────────────────── */ + +Datum +galilean_lagrange_equatorial_de(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + int64 ts = PG_GETARG_INT64(2); + double jd; + double moon_xyz[3], moon_vel[3], lp_rel[3]; + pg_equatorial *result; + + if (body_id < L12_IO || body_id > L12_CALLISTO) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("galilean_lagrange_equatorial_de: body_id %d must be 0-3", + body_id))); + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("galilean_lagrange_equatorial_de: point_id %d must be 1-5", + point_id))); + + jd = timestamptz_to_jd(ts); + GetL12Coor(jd, body_id, moon_xyz, moon_vel); + + if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 'g', body_id, + point_id, lp_rel) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("galilean_lagrange_equatorial_de: solver failed"))); + + result = (pg_equatorial *) palloc(sizeof(pg_equatorial)); + equatorial_pmoon_lagrange_de(lp_rel, BODY_JUPITER, jd, result); + + PG_RETURN_POINTER(result); +} + + +/* ── 10. Saturn moon Lagrange observe DE ─────────────────── */ + +Datum +saturn_moon_lagrange_observe_de(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2); + int64 ts = PG_GETARG_INT64(3); + double jd; + double moon_xyz[3], moon_vel[3], lp_rel[3]; + pg_topocentric *result; + + if (body_id < TASS17_MIMAS || body_id > TASS17_HYPERION) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("saturn_moon_lagrange_observe_de: body_id %d must be 0-7", + body_id))); + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("saturn_moon_lagrange_observe_de: point_id %d must be 1-5", + point_id))); + + jd = timestamptz_to_jd(ts); + GetTass17Coor(jd, body_id, moon_xyz, moon_vel); + + if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 's', body_id, + point_id, lp_rel) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("saturn_moon_lagrange_observe_de: solver failed"))); + + result = (pg_topocentric *) palloc(sizeof(pg_topocentric)); + observe_pmoon_lagrange_de(lp_rel, BODY_SATURN, jd, obs, result); + + PG_RETURN_POINTER(result); +} + + +/* ── 11. Saturn moon Lagrange equatorial DE ──────────────── */ + +Datum +saturn_moon_lagrange_equatorial_de(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + int64 ts = PG_GETARG_INT64(2); + double jd; + double moon_xyz[3], moon_vel[3], lp_rel[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_lagrange_equatorial_de: body_id %d must be 0-7", + body_id))); + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("saturn_moon_lagrange_equatorial_de: point_id %d must be 1-5", + point_id))); + + jd = timestamptz_to_jd(ts); + GetTass17Coor(jd, body_id, moon_xyz, moon_vel); + + if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 's', body_id, + point_id, lp_rel) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("saturn_moon_lagrange_equatorial_de: solver failed"))); + + result = (pg_equatorial *) palloc(sizeof(pg_equatorial)); + equatorial_pmoon_lagrange_de(lp_rel, BODY_SATURN, jd, result); + + PG_RETURN_POINTER(result); +} + + +/* ── 12. Uranus moon Lagrange observe DE ─────────────────── */ + +Datum +uranus_moon_lagrange_observe_de(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2); + int64 ts = PG_GETARG_INT64(3); + double jd; + double moon_xyz[3], moon_vel[3], lp_rel[3]; + pg_topocentric *result; + + if (body_id < GUST86_MIRANDA || body_id > GUST86_OBERON) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("uranus_moon_lagrange_observe_de: body_id %d must be 0-4", + body_id))); + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("uranus_moon_lagrange_observe_de: point_id %d must be 1-5", + point_id))); + + jd = timestamptz_to_jd(ts); + GetGust86Coor(jd, body_id, moon_xyz, moon_vel); + + if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 'u', body_id, + point_id, lp_rel) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("uranus_moon_lagrange_observe_de: solver failed"))); + + result = (pg_topocentric *) palloc(sizeof(pg_topocentric)); + observe_pmoon_lagrange_de(lp_rel, BODY_URANUS, jd, obs, result); + + PG_RETURN_POINTER(result); +} + + +/* ── 13. Uranus moon Lagrange equatorial DE ──────────────── */ + +Datum +uranus_moon_lagrange_equatorial_de(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + int64 ts = PG_GETARG_INT64(2); + double jd; + double moon_xyz[3], moon_vel[3], lp_rel[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_lagrange_equatorial_de: body_id %d must be 0-4", + body_id))); + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("uranus_moon_lagrange_equatorial_de: point_id %d must be 1-5", + point_id))); + + jd = timestamptz_to_jd(ts); + GetGust86Coor(jd, body_id, moon_xyz, moon_vel); + + if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 'u', body_id, + point_id, lp_rel) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("uranus_moon_lagrange_equatorial_de: solver failed"))); + + result = (pg_equatorial *) palloc(sizeof(pg_equatorial)); + equatorial_pmoon_lagrange_de(lp_rel, BODY_URANUS, jd, result); + + PG_RETURN_POINTER(result); +} + + +/* ── 14. Mars moon Lagrange observe DE ───────────────────── */ + +Datum +mars_moon_lagrange_observe_de(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2); + int64 ts = PG_GETARG_INT64(3); + double jd; + double moon_xyz[3], moon_vel[3], lp_rel[3]; + pg_topocentric *result; + + if (body_id < MARS_SAT_PHOBOS || body_id > MARS_SAT_DEIMOS) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("mars_moon_lagrange_observe_de: body_id %d must be 0-1", + body_id))); + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("mars_moon_lagrange_observe_de: point_id %d must be 1-5", + point_id))); + + jd = timestamptz_to_jd(ts); + GetMarsSatCoor(jd, body_id, moon_xyz, moon_vel); + + if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 'm', body_id, + point_id, lp_rel) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("mars_moon_lagrange_observe_de: solver failed"))); + + result = (pg_topocentric *) palloc(sizeof(pg_topocentric)); + observe_pmoon_lagrange_de(lp_rel, BODY_MARS, jd, obs, result); + + PG_RETURN_POINTER(result); +} + + +/* ── 15. Mars moon Lagrange equatorial DE ────────────────── */ + +Datum +mars_moon_lagrange_equatorial_de(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + int64 ts = PG_GETARG_INT64(2); + double jd; + double moon_xyz[3], moon_vel[3], lp_rel[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_lagrange_equatorial_de: body_id %d must be 0-1", + body_id))); + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("mars_moon_lagrange_equatorial_de: point_id %d must be 1-5", + point_id))); + + jd = timestamptz_to_jd(ts); + GetMarsSatCoor(jd, body_id, moon_xyz, moon_vel); + + if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 'm', body_id, + point_id, lp_rel) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("mars_moon_lagrange_equatorial_de: solver failed"))); + + result = (pg_equatorial *) palloc(sizeof(pg_equatorial)); + equatorial_pmoon_lagrange_de(lp_rel, BODY_MARS, jd, result); + + PG_RETURN_POINTER(result); +} + + +/* ================================================================ + * 16. hill_radius_de(body_id, timestamptz) -> float8 (AU) + * + * Hill radius using DE for planet heliocentric distance. STABLE. + * ================================================================ + */ +Datum +hill_radius_de(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int64 ts = PG_GETARG_INT64(1); + double jd; + double planet_xyz[6], sep, ratio, mu; + + if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("hill_radius_de: body_id %d must be 1-8", body_id))); + + jd = timestamptz_to_jd(ts); + + if (!eph_de_planet(body_id, jd, planet_xyz)) + { + if (eph_de_available()) + ereport(NOTICE, + (errmsg("DE ephemeris unavailable, falling back to VSOP87"))); + + GetVsop87Coor(jd, body_id - 1, planet_xyz); + } + + sep = sqrt(planet_xyz[0]*planet_xyz[0] + + planet_xyz[1]*planet_xyz[1] + + planet_xyz[2]*planet_xyz[2]); + + ratio = sun_planet_ratio(body_id); + mu = mu_from_ratio(ratio); + + PG_RETURN_FLOAT8(lagrange_hill_radius(sep, mu)); +} + + +/* ================================================================ + * 17. lagrange_zone_radius_de(body_id, point_id, timestamptz) + * -> float8 (AU) + * + * Libration zone radius using DE for planet distance. STABLE. + * ================================================================ + */ +Datum +lagrange_zone_radius_de(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + int64 ts = PG_GETARG_INT64(2); + double jd; + double planet_xyz[6], sep, ratio, mu, zr; + + validate_lagrange_args_de(body_id, point_id, "lagrange_zone_radius_de"); + + jd = timestamptz_to_jd(ts); + + if (!eph_de_planet(body_id, jd, planet_xyz)) + { + if (eph_de_available()) + ereport(NOTICE, + (errmsg("DE ephemeris unavailable, falling back to VSOP87"))); + + GetVsop87Coor(jd, body_id - 1, planet_xyz); + } + + sep = sqrt(planet_xyz[0]*planet_xyz[0] + + planet_xyz[1]*planet_xyz[1] + + planet_xyz[2]*planet_xyz[2]); + + ratio = sun_planet_ratio(body_id); + mu = mu_from_ratio(ratio); + + zr = lagrange_zone_radius(sep, mu, point_id); + + PG_RETURN_FLOAT8(zr); +} diff --git a/src/lagrange_funcs.c b/src/lagrange_funcs.c new file mode 100644 index 0000000..d6bd28d --- /dev/null +++ b/src/lagrange_funcs.c @@ -0,0 +1,916 @@ +/* + * lagrange_funcs.c -- Lagrange point observation functions + * + * SQL functions for computing positions and observations of the five + * Lagrange equilibrium points in the circular restricted three-body + * problem (CR3BP). Covers Sun-planet, Earth-Moon, and planetary moon + * systems. + * + * The pipeline mirrors planet_funcs.c / moon_funcs.c: + * 1. Primary + secondary heliocentric positions (VSOP87) + * 2. Secondary velocity via analytic (VSOP87) or central difference + * 3. lagrange_position() -> L-point heliocentric ecliptic J2000 + * 4. Geocentric ecliptic (subtract Earth) + * 5. observe_from_geocentric() or geocentric_to_equatorial() + * + * All functions are IMMUTABLE STRICT PARALLEL SAFE (compiled-in + * mass ratios, VSOP87/ELP82B coefficients). + */ + +#include "postgres.h" +#include "fmgr.h" +#include "funcapi.h" +#include "utils/timestamp.h" +#include "utils/builtins.h" +#include "types.h" +#include "astro_math.h" +#include "lagrange.h" +#include "vsop87.h" +#include "kepler.h" +#include "elp82b.h" +#include "l12.h" +#include "tass17.h" +#include "gust86.h" +#include "marssat.h" +#include + +/* ── Forward declarations ──────────────────────────────── */ + +/* Sun-planet */ +PG_FUNCTION_INFO_V1(lagrange_heliocentric); +PG_FUNCTION_INFO_V1(lagrange_observe); +PG_FUNCTION_INFO_V1(lagrange_equatorial); +PG_FUNCTION_INFO_V1(lagrange_distance); +PG_FUNCTION_INFO_V1(lagrange_distance_oe); + +/* Earth-Moon */ +PG_FUNCTION_INFO_V1(lunar_lagrange_observe); +PG_FUNCTION_INFO_V1(lunar_lagrange_equatorial); + +/* Planetary moons */ +PG_FUNCTION_INFO_V1(galilean_lagrange_observe); +PG_FUNCTION_INFO_V1(galilean_lagrange_equatorial); +PG_FUNCTION_INFO_V1(saturn_moon_lagrange_observe); +PG_FUNCTION_INFO_V1(saturn_moon_lagrange_equatorial); +PG_FUNCTION_INFO_V1(uranus_moon_lagrange_observe); +PG_FUNCTION_INFO_V1(uranus_moon_lagrange_equatorial); +PG_FUNCTION_INFO_V1(mars_moon_lagrange_observe); +PG_FUNCTION_INFO_V1(mars_moon_lagrange_equatorial); + +/* Hill/zone/convenience */ +PG_FUNCTION_INFO_V1(hill_radius_func); +PG_FUNCTION_INFO_V1(hill_radius_lunar); +PG_FUNCTION_INFO_V1(lagrange_zone_radius_func); +PG_FUNCTION_INFO_V1(lagrange_mass_ratio); +PG_FUNCTION_INFO_V1(lagrange_point_name); + + +/* ── Internal helpers ──────────────────────────────────── */ + +/* + * Validate body_id and point_id common to Sun-planet functions. + */ +static void +validate_lagrange_args(int body_id, int point_id, const char *func_name) +{ + if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("%s: body_id %d must be 1-8 (Mercury-Neptune)", + func_name, body_id))); + + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("%s: point_id %d must be 1-5 (L1-L5)", + func_name, point_id))); +} + + +/* + * Compute Sun-planet Lagrange point in heliocentric ecliptic J2000 (AU). + * + * Uses VSOP87 for planet position and velocity. The Sun is at (0,0,0). + */ +static int +compute_sun_planet_lagrange(int body_id, int point_id, double jd, + double result[3]) +{ + double planet_xyz[6]; + double sun[3] = {0.0, 0.0, 0.0}; + double ratio, mu; + int vsop_body = body_id - 1; + + GetVsop87Coor(jd, vsop_body, planet_xyz); + + ratio = sun_planet_ratio(body_id); + mu = mu_from_ratio(ratio); + + /* planet_xyz[3..5] is velocity in AU/day from VSOP87 */ + return lagrange_position(sun, planet_xyz, &planet_xyz[3], mu, point_id, + result); +} + + +/* + * Compute Earth-Moon Lagrange point in geocentric ecliptic J2000 (AU). + * + * ELP2000-82B provides position only. Velocity via central difference + * (dt=0.001 day ~ 86 seconds). Error ~40 km at lunar distance — negligible + * vs libration zone scale (~60,000 km for L4/L5). + * + * Returns position relative to Earth (geocentric), not heliocentric. + */ +static int +compute_lunar_lagrange(int point_id, double jd, double result_geo[3]) +{ + double moon_pos[3], moon_fwd[3], moon_bwd[3]; + double moon_vel[3]; + double earth[3] = {0.0, 0.0, 0.0}; /* geocentric frame: Earth at origin */ + double mu; + double dt = 0.001; /* days */ + int rc; + + /* Moon geocentric position */ + GetElp82bCoor(jd, moon_pos); + + /* Velocity via central difference */ + GetElp82bCoor(jd + dt, moon_fwd); + GetElp82bCoor(jd - dt, moon_bwd); + moon_vel[0] = (moon_fwd[0] - moon_bwd[0]) / (2.0 * dt); + moon_vel[1] = (moon_fwd[1] - moon_bwd[1]) / (2.0 * dt); + moon_vel[2] = (moon_fwd[2] - moon_bwd[2]) / (2.0 * dt); + + mu = mu_from_ratio(EARTH_MOON_EMRAT); + + rc = lagrange_position(earth, moon_pos, moon_vel, mu, point_id, + result_geo); + return rc; +} + + +/* + * Compute planet-moon Lagrange point heliocentric ecliptic J2000 (AU). + * + * The moon theory provides both position and velocity relative to the + * parent planet. We work in a frame centered on the planet. + */ +static int +compute_planetary_moon_lagrange(const double moon_rel[3], + const double moon_vel[3], + char family, int moon_id, + int point_id, + double lp_rel[3]) +{ + double planet_origin[3] = {0.0, 0.0, 0.0}; + double ratio, mu; + + ratio = planet_moon_ratio(family, moon_id); + if (ratio < 0.0) + return -1; + + mu = mu_from_ratio(ratio); + + return lagrange_position(planet_origin, moon_rel, moon_vel, mu, + point_id, lp_rel); +} + + +/* ================================================================ + * lagrange_heliocentric(body_id int4, point_id int4, timestamptz) + * -> heliocentric + * + * Sun-planet Lagrange point in heliocentric ecliptic J2000. + * ================================================================ + */ +Datum +lagrange_heliocentric(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + int64 ts = PG_GETARG_INT64(2); + double jd; + double lp[3]; + pg_heliocentric *result; + + validate_lagrange_args(body_id, point_id, "lagrange_heliocentric"); + + jd = timestamptz_to_jd(ts); + + if (compute_sun_planet_lagrange(body_id, point_id, jd, lp) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("lagrange_heliocentric: solver failed for body %d L%d", + body_id, point_id))); + + result = (pg_heliocentric *) palloc(sizeof(pg_heliocentric)); + result->x = lp[0]; + result->y = lp[1]; + result->z = lp[2]; + + PG_RETURN_POINTER(result); +} + + +/* ================================================================ + * lagrange_observe(body_id, point_id, observer, timestamptz) + * -> topocentric + * + * Full pipeline: L-point heliocentric -> geocentric -> az/el. + * ================================================================ + */ +Datum +lagrange_observe(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2); + int64 ts = PG_GETARG_INT64(3); + double jd; + double lp[3], earth_xyz[6], geo_ecl[3]; + pg_topocentric *result; + + validate_lagrange_args(body_id, point_id, "lagrange_observe"); + + jd = timestamptz_to_jd(ts); + + if (compute_sun_planet_lagrange(body_id, point_id, jd, lp) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("lagrange_observe: solver failed"))); + + /* Geocentric */ + GetVsop87Coor(jd, 2, earth_xyz); + geo_ecl[0] = lp[0] - earth_xyz[0]; + geo_ecl[1] = lp[1] - earth_xyz[1]; + geo_ecl[2] = lp[2] - earth_xyz[2]; + + result = (pg_topocentric *) palloc(sizeof(pg_topocentric)); + observe_from_geocentric(geo_ecl, jd, obs, result); + + PG_RETURN_POINTER(result); +} + + +/* ================================================================ + * lagrange_equatorial(body_id, point_id, timestamptz) -> equatorial + * ================================================================ + */ +Datum +lagrange_equatorial(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + int64 ts = PG_GETARG_INT64(2); + double jd; + double lp[3], earth_xyz[6], geo_ecl[3]; + pg_equatorial *result; + + validate_lagrange_args(body_id, point_id, "lagrange_equatorial"); + + jd = timestamptz_to_jd(ts); + + if (compute_sun_planet_lagrange(body_id, point_id, jd, lp) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("lagrange_equatorial: solver failed"))); + + GetVsop87Coor(jd, 2, earth_xyz); + geo_ecl[0] = lp[0] - earth_xyz[0]; + geo_ecl[1] = lp[1] - earth_xyz[1]; + geo_ecl[2] = lp[2] - earth_xyz[2]; + + result = (pg_equatorial *) palloc(sizeof(pg_equatorial)); + geocentric_to_equatorial(geo_ecl, jd, result); + + PG_RETURN_POINTER(result); +} + + +/* ================================================================ + * lagrange_distance(body_id, point_id, heliocentric, timestamptz) + * -> float8 + * + * Distance from a heliocentric position to a Sun-planet L-point (AU). + * ================================================================ + */ +Datum +lagrange_distance(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + pg_heliocentric *pos = (pg_heliocentric *) PG_GETARG_POINTER(2); + int64 ts = PG_GETARG_INT64(3); + double jd; + double lp[3]; + double dx, dy, dz, dist; + + validate_lagrange_args(body_id, point_id, "lagrange_distance"); + + jd = timestamptz_to_jd(ts); + + if (compute_sun_planet_lagrange(body_id, point_id, jd, lp) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("lagrange_distance: solver failed"))); + + dx = pos->x - lp[0]; + dy = pos->y - lp[1]; + dz = pos->z - lp[2]; + dist = sqrt(dx*dx + dy*dy + dz*dz); + + PG_RETURN_FLOAT8(dist); +} + + +/* ================================================================ + * lagrange_distance_oe(body_id, point_id, orbital_elements, timestamptz) + * -> float8 + * + * Distance from an asteroid/comet to a Sun-planet L-point (AU). + * ================================================================ + */ +Datum +lagrange_distance_oe(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + pg_orbital_elements *oe = (pg_orbital_elements *) PG_GETARG_POINTER(2); + int64 ts = PG_GETARG_INT64(3); + double jd; + double lp[3], body_pos[3]; + double dx, dy, dz, dist; + + validate_lagrange_args(body_id, point_id, "lagrange_distance_oe"); + + jd = timestamptz_to_jd(ts); + + if (compute_sun_planet_lagrange(body_id, point_id, jd, lp) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("lagrange_distance_oe: solver failed"))); + + kepler_position(oe->q, oe->e, oe->inc, oe->arg_peri, oe->raan, + oe->tp, jd, body_pos); + + dx = body_pos[0] - lp[0]; + dy = body_pos[1] - lp[1]; + dz = body_pos[2] - lp[2]; + dist = sqrt(dx*dx + dy*dy + dz*dz); + + PG_RETURN_FLOAT8(dist); +} + + +/* ================================================================ + * lunar_lagrange_observe(point_id, observer, timestamptz) -> topocentric + * + * Earth-Moon Lagrange point observation. + * ================================================================ + */ +Datum +lunar_lagrange_observe(PG_FUNCTION_ARGS) +{ + int32 point_id = PG_GETARG_INT32(0); + pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(1); + int64 ts = PG_GETARG_INT64(2); + double jd; + double lp_geo[3]; + pg_topocentric *result; + + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("lunar_lagrange_observe: point_id %d must be 1-5", + point_id))); + + jd = timestamptz_to_jd(ts); + + if (compute_lunar_lagrange(point_id, jd, lp_geo) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("lunar_lagrange_observe: solver failed"))); + + /* lp_geo is already geocentric ecliptic J2000 (AU) */ + result = (pg_topocentric *) palloc(sizeof(pg_topocentric)); + observe_from_geocentric(lp_geo, jd, obs, result); + + PG_RETURN_POINTER(result); +} + + +/* ================================================================ + * lunar_lagrange_equatorial(point_id, timestamptz) -> equatorial + * ================================================================ + */ +Datum +lunar_lagrange_equatorial(PG_FUNCTION_ARGS) +{ + int32 point_id = PG_GETARG_INT32(0); + int64 ts = PG_GETARG_INT64(1); + double jd; + double lp_geo[3]; + pg_equatorial *result; + + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("lunar_lagrange_equatorial: point_id %d must be 1-5", + point_id))); + + jd = timestamptz_to_jd(ts); + + if (compute_lunar_lagrange(point_id, jd, lp_geo) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("lunar_lagrange_equatorial: solver failed"))); + + result = (pg_equatorial *) palloc(sizeof(pg_equatorial)); + geocentric_to_equatorial(lp_geo, jd, result); + + PG_RETURN_POINTER(result); +} + + +/* ================================================================ + * Planetary moon Lagrange points + * + * Each family duplicates the observe/equatorial static helpers from + * moon_funcs.c (per no-cross-TU-symbol convention), but routes through + * the Lagrange solver instead of returning the moon position directly. + * ================================================================ + */ + +/* + * Internal: observe a planetary moon Lagrange point. + * L-point offset is relative to planet, same as moon offset. + */ +static void +observe_pmoon_lagrange(const double lp_rel[3], int vsop_parent, + double jd, const pg_observer *obs, + pg_topocentric *result) +{ + double parent_xyz[6]; + double earth_xyz[6]; + double geo_ecl[3]; + + GetVsop87Coor(jd, vsop_parent, parent_xyz); + GetVsop87Coor(jd, 2, earth_xyz); + + geo_ecl[0] = (parent_xyz[0] + lp_rel[0]) - earth_xyz[0]; + geo_ecl[1] = (parent_xyz[1] + lp_rel[1]) - earth_xyz[1]; + geo_ecl[2] = (parent_xyz[2] + lp_rel[2]) - earth_xyz[2]; + + observe_from_geocentric(geo_ecl, jd, obs, result); +} + +static void +equatorial_pmoon_lagrange(const double lp_rel[3], int vsop_parent, + double jd, pg_equatorial *result) +{ + double parent_xyz[6]; + double earth_xyz[6]; + double geo_ecl[3]; + + GetVsop87Coor(jd, vsop_parent, parent_xyz); + GetVsop87Coor(jd, 2, earth_xyz); + + geo_ecl[0] = (parent_xyz[0] + lp_rel[0]) - earth_xyz[0]; + geo_ecl[1] = (parent_xyz[1] + lp_rel[1]) - earth_xyz[1]; + geo_ecl[2] = (parent_xyz[2] + lp_rel[2]) - earth_xyz[2]; + + geocentric_to_equatorial(geo_ecl, jd, result); +} + + +/* ── Galilean moons ────────────────────────────────────── */ + +Datum +galilean_lagrange_observe(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2); + int64 ts = PG_GETARG_INT64(3); + double jd; + double moon_xyz[3], moon_vel[3], lp_rel[3]; + pg_topocentric *result; + + if (body_id < L12_IO || body_id > L12_CALLISTO) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("galilean_lagrange_observe: body_id %d must be 0-3", + body_id))); + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("galilean_lagrange_observe: point_id %d must be 1-5", + point_id))); + + jd = timestamptz_to_jd(ts); + GetL12Coor(jd, body_id, moon_xyz, moon_vel); + + if (compute_planetary_moon_lagrange(moon_xyz, moon_vel, 'g', body_id, + point_id, lp_rel) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("galilean_lagrange_observe: solver failed"))); + + result = (pg_topocentric *) palloc(sizeof(pg_topocentric)); + observe_pmoon_lagrange(lp_rel, 4, jd, obs, result); + + PG_RETURN_POINTER(result); +} + + +Datum +galilean_lagrange_equatorial(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + int64 ts = PG_GETARG_INT64(2); + double jd; + double moon_xyz[3], moon_vel[3], lp_rel[3]; + pg_equatorial *result; + + if (body_id < L12_IO || body_id > L12_CALLISTO) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("galilean_lagrange_equatorial: body_id %d must be 0-3", + body_id))); + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("galilean_lagrange_equatorial: point_id %d must be 1-5", + point_id))); + + jd = timestamptz_to_jd(ts); + GetL12Coor(jd, body_id, moon_xyz, moon_vel); + + if (compute_planetary_moon_lagrange(moon_xyz, moon_vel, 'g', body_id, + point_id, lp_rel) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("galilean_lagrange_equatorial: solver failed"))); + + result = (pg_equatorial *) palloc(sizeof(pg_equatorial)); + equatorial_pmoon_lagrange(lp_rel, 4, jd, result); + + PG_RETURN_POINTER(result); +} + + +/* ── Saturn moons ──────────────────────────────────────── */ + +Datum +saturn_moon_lagrange_observe(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2); + int64 ts = PG_GETARG_INT64(3); + double jd; + double moon_xyz[3], moon_vel[3], lp_rel[3]; + pg_topocentric *result; + + if (body_id < TASS17_MIMAS || body_id > TASS17_HYPERION) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("saturn_moon_lagrange_observe: body_id %d must be 0-7", + body_id))); + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("saturn_moon_lagrange_observe: point_id %d must be 1-5", + point_id))); + + jd = timestamptz_to_jd(ts); + GetTass17Coor(jd, body_id, moon_xyz, moon_vel); + + if (compute_planetary_moon_lagrange(moon_xyz, moon_vel, 's', body_id, + point_id, lp_rel) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("saturn_moon_lagrange_observe: solver failed"))); + + result = (pg_topocentric *) palloc(sizeof(pg_topocentric)); + observe_pmoon_lagrange(lp_rel, 5, jd, obs, result); + + PG_RETURN_POINTER(result); +} + + +Datum +saturn_moon_lagrange_equatorial(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + int64 ts = PG_GETARG_INT64(2); + double jd; + double moon_xyz[3], moon_vel[3], lp_rel[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_lagrange_equatorial: body_id %d must be 0-7", + body_id))); + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("saturn_moon_lagrange_equatorial: point_id %d must be 1-5", + point_id))); + + jd = timestamptz_to_jd(ts); + GetTass17Coor(jd, body_id, moon_xyz, moon_vel); + + if (compute_planetary_moon_lagrange(moon_xyz, moon_vel, 's', body_id, + point_id, lp_rel) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("saturn_moon_lagrange_equatorial: solver failed"))); + + result = (pg_equatorial *) palloc(sizeof(pg_equatorial)); + equatorial_pmoon_lagrange(lp_rel, 5, jd, result); + + PG_RETURN_POINTER(result); +} + + +/* ── Uranus moons ──────────────────────────────────────── */ + +Datum +uranus_moon_lagrange_observe(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2); + int64 ts = PG_GETARG_INT64(3); + double jd; + double moon_xyz[3], moon_vel[3], lp_rel[3]; + pg_topocentric *result; + + if (body_id < GUST86_MIRANDA || body_id > GUST86_OBERON) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("uranus_moon_lagrange_observe: body_id %d must be 0-4", + body_id))); + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("uranus_moon_lagrange_observe: point_id %d must be 1-5", + point_id))); + + jd = timestamptz_to_jd(ts); + GetGust86Coor(jd, body_id, moon_xyz, moon_vel); + + if (compute_planetary_moon_lagrange(moon_xyz, moon_vel, 'u', body_id, + point_id, lp_rel) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("uranus_moon_lagrange_observe: solver failed"))); + + result = (pg_topocentric *) palloc(sizeof(pg_topocentric)); + observe_pmoon_lagrange(lp_rel, 6, jd, obs, result); + + PG_RETURN_POINTER(result); +} + + +Datum +uranus_moon_lagrange_equatorial(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + int64 ts = PG_GETARG_INT64(2); + double jd; + double moon_xyz[3], moon_vel[3], lp_rel[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_lagrange_equatorial: body_id %d must be 0-4", + body_id))); + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("uranus_moon_lagrange_equatorial: point_id %d must be 1-5", + point_id))); + + jd = timestamptz_to_jd(ts); + GetGust86Coor(jd, body_id, moon_xyz, moon_vel); + + if (compute_planetary_moon_lagrange(moon_xyz, moon_vel, 'u', body_id, + point_id, lp_rel) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("uranus_moon_lagrange_equatorial: solver failed"))); + + result = (pg_equatorial *) palloc(sizeof(pg_equatorial)); + equatorial_pmoon_lagrange(lp_rel, 6, jd, result); + + PG_RETURN_POINTER(result); +} + + +/* ── Mars moons ────────────────────────────────────────── */ + +Datum +mars_moon_lagrange_observe(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2); + int64 ts = PG_GETARG_INT64(3); + double jd; + double moon_xyz[3], moon_vel[3], lp_rel[3]; + pg_topocentric *result; + + if (body_id < MARS_SAT_PHOBOS || body_id > MARS_SAT_DEIMOS) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("mars_moon_lagrange_observe: body_id %d must be 0-1", + body_id))); + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("mars_moon_lagrange_observe: point_id %d must be 1-5", + point_id))); + + jd = timestamptz_to_jd(ts); + GetMarsSatCoor(jd, body_id, moon_xyz, moon_vel); + + if (compute_planetary_moon_lagrange(moon_xyz, moon_vel, 'm', body_id, + point_id, lp_rel) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("mars_moon_lagrange_observe: solver failed"))); + + result = (pg_topocentric *) palloc(sizeof(pg_topocentric)); + observe_pmoon_lagrange(lp_rel, 3, jd, obs, result); + + PG_RETURN_POINTER(result); +} + + +Datum +mars_moon_lagrange_equatorial(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + int64 ts = PG_GETARG_INT64(2); + double jd; + double moon_xyz[3], moon_vel[3], lp_rel[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_lagrange_equatorial: body_id %d must be 0-1", + body_id))); + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("mars_moon_lagrange_equatorial: point_id %d must be 1-5", + point_id))); + + jd = timestamptz_to_jd(ts); + GetMarsSatCoor(jd, body_id, moon_xyz, moon_vel); + + if (compute_planetary_moon_lagrange(moon_xyz, moon_vel, 'm', body_id, + point_id, lp_rel) != 0) + ereport(ERROR, + (errcode(ERRCODE_INTERNAL_ERROR), + errmsg("mars_moon_lagrange_equatorial: solver failed"))); + + result = (pg_equatorial *) palloc(sizeof(pg_equatorial)); + equatorial_pmoon_lagrange(lp_rel, 3, jd, result); + + PG_RETURN_POINTER(result); +} + + +/* ================================================================ + * Hill radius / zone / convenience functions + * ================================================================ + */ + +/* hill_radius(body_id int4, timestamptz) -> float8 (AU) */ +Datum +hill_radius_func(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int64 ts = PG_GETARG_INT64(1); + double jd; + double planet_xyz[6], sep, ratio, mu; + + if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("hill_radius: body_id %d must be 1-8", body_id))); + + jd = timestamptz_to_jd(ts); + GetVsop87Coor(jd, body_id - 1, planet_xyz); + + sep = sqrt(planet_xyz[0]*planet_xyz[0] + + planet_xyz[1]*planet_xyz[1] + + planet_xyz[2]*planet_xyz[2]); + + ratio = sun_planet_ratio(body_id); + mu = mu_from_ratio(ratio); + + PG_RETURN_FLOAT8(lagrange_hill_radius(sep, mu)); +} + + +/* hill_radius_lunar(timestamptz) -> float8 (AU) */ +Datum +hill_radius_lunar(PG_FUNCTION_ARGS) +{ + int64 ts = PG_GETARG_INT64(0); + double jd; + double moon_pos[3], sep, mu; + + jd = timestamptz_to_jd(ts); + GetElp82bCoor(jd, moon_pos); + + sep = sqrt(moon_pos[0]*moon_pos[0] + + moon_pos[1]*moon_pos[1] + + moon_pos[2]*moon_pos[2]); + + mu = mu_from_ratio(EARTH_MOON_EMRAT); + + PG_RETURN_FLOAT8(lagrange_hill_radius(sep, mu)); +} + + +/* lagrange_zone_radius(body_id, point_id, timestamptz) -> float8 (AU) */ +Datum +lagrange_zone_radius_func(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + int32 point_id = PG_GETARG_INT32(1); + int64 ts = PG_GETARG_INT64(2); + double jd; + double planet_xyz[6], sep, ratio, mu, zr; + + validate_lagrange_args(body_id, point_id, "lagrange_zone_radius"); + + jd = timestamptz_to_jd(ts); + GetVsop87Coor(jd, body_id - 1, planet_xyz); + + sep = sqrt(planet_xyz[0]*planet_xyz[0] + + planet_xyz[1]*planet_xyz[1] + + planet_xyz[2]*planet_xyz[2]); + + ratio = sun_planet_ratio(body_id); + mu = mu_from_ratio(ratio); + + zr = lagrange_zone_radius(sep, mu, point_id); + + PG_RETURN_FLOAT8(zr); +} + + +/* lagrange_mass_ratio(body_id int4) -> float8 */ +Datum +lagrange_mass_ratio(PG_FUNCTION_ARGS) +{ + int32 body_id = PG_GETARG_INT32(0); + double ratio; + + if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("lagrange_mass_ratio: body_id %d must be 1-8", + body_id))); + + ratio = sun_planet_ratio(body_id); + + PG_RETURN_FLOAT8(mu_from_ratio(ratio)); +} + + +/* lagrange_point_name(point_id int4) -> text */ +Datum +lagrange_point_name(PG_FUNCTION_ARGS) +{ + int32 point_id = PG_GETARG_INT32(0); + + if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("lagrange_point_name: point_id %d must be 1-5", + point_id))); + + switch (point_id) + { + case LAGRANGE_L1: PG_RETURN_TEXT_P(cstring_to_text("L1")); + case LAGRANGE_L2: PG_RETURN_TEXT_P(cstring_to_text("L2")); + case LAGRANGE_L3: PG_RETURN_TEXT_P(cstring_to_text("L3")); + case LAGRANGE_L4: PG_RETURN_TEXT_P(cstring_to_text("L4")); + case LAGRANGE_L5: PG_RETURN_TEXT_P(cstring_to_text("L5")); + default: PG_RETURN_NULL(); + } +} diff --git a/test/expected/v020_features.out b/test/expected/v020_features.out new file mode 100644 index 0000000..a2a934a --- /dev/null +++ b/test/expected/v020_features.out @@ -0,0 +1,323 @@ +-- v020_features: Lagrange point support +-- Tests Sun-planet, Earth-Moon, planetary moon Lagrange points, +-- Hill radius, zone radius, DE fallback, and input validation. +-- Reference observer: Greenwich, UK +\set obs '''(51.4769,-0.0005,0)''' +-- Reference time: J2000 epoch (2000-01-01 12:00:00 UTC) +\set t '''2000-01-01 12:00:00+00''' +-- ============================================================ +-- Sun-Earth L1/L2: should be ~0.01 AU from Earth (~1.5 million km) +-- SOHO is at L1, JWST at L2. +-- ============================================================ +-- L1 heliocentric: should be close to Earth's heliocentric (~1 AU from Sun) +SELECT + round(helio_distance(lagrange_heliocentric(3, 1, :t ::timestamptz))::numeric, 2) AS sun_dist_au; + sun_dist_au +------------- + 0.97 +(1 row) + +-- L2 heliocentric: also ~1 AU from Sun, slightly further than L1 +SELECT + round(helio_distance(lagrange_heliocentric(3, 2, :t ::timestamptz))::numeric, 2) AS sun_dist_au; + sun_dist_au +------------- + 0.99 +(1 row) + +-- L1 between Sun and Earth (closer to Sun than L2) +SELECT + helio_distance(lagrange_heliocentric(3, 1, :t ::timestamptz)) + < + helio_distance(lagrange_heliocentric(3, 2, :t ::timestamptz)) + AS l1_closer_than_l2; + l1_closer_than_l2 +------------------- + t +(1 row) + +-- ============================================================ +-- Sun-Jupiter L4/L5: ~60 degrees from Jupiter, ~5.2 AU from Sun +-- These are the Trojan asteroid zones. +-- ============================================================ +-- L4/L5 should be ~5.2 AU from Sun +SELECT + round(helio_distance(lagrange_heliocentric(5, 4, :t ::timestamptz))::numeric, 1) AS l4_sun_dist; + l4_sun_dist +------------- + 5.0 +(1 row) + +SELECT + round(helio_distance(lagrange_heliocentric(5, 5, :t ::timestamptz))::numeric, 1) AS l5_sun_dist; + l5_sun_dist +------------- + 5.0 +(1 row) + +-- L4 and L5 equidistant from Sun (within 0.001 AU) +SELECT + abs( + helio_distance(lagrange_heliocentric(5, 4, :t ::timestamptz)) + - + helio_distance(lagrange_heliocentric(5, 5, :t ::timestamptz)) + ) < 0.001 AS l4_l5_equidistant; + l4_l5_equidistant +------------------- + t +(1 row) + +-- ============================================================ +-- Earth-Moon L1: ~326,000 km from Earth +-- ============================================================ +-- lunar_lagrange_equatorial returns distance in km +SELECT + round(eq_distance(lunar_lagrange_equatorial(1, :t ::timestamptz))::numeric, -3) + BETWEEN 300000 AND 360000 AS em_l1_in_range; + em_l1_in_range +---------------- + t +(1 row) + +-- ============================================================ +-- lagrange_observe returns valid az/el +-- ============================================================ +SELECT + topo_elevation(lagrange_observe(3, 2, :obs ::observer, :t ::timestamptz)) + BETWEEN -90 AND 90 AS valid_elevation; + valid_elevation +----------------- + t +(1 row) + +-- lagrange_equatorial returns valid RA/Dec +SELECT + eq_ra(lagrange_equatorial(3, 1, :t ::timestamptz)) BETWEEN 0 AND 24 AS valid_ra, + eq_dec(lagrange_equatorial(3, 1, :t ::timestamptz)) BETWEEN -90 AND 90 AS valid_dec; + valid_ra | valid_dec +----------+----------- + t | t +(1 row) + +-- ============================================================ +-- lagrange_distance self-test: L-point distance to itself ≈ 0 +-- ============================================================ +SELECT + round(lagrange_distance( + 5, 4, + lagrange_heliocentric(5, 4, :t ::timestamptz), + :t ::timestamptz + )::numeric, 10) AS self_distance; + self_distance +--------------- + 0.0000000000 +(1 row) + +-- ============================================================ +-- Hill radius +-- ============================================================ +-- Jupiter Hill radius ~0.35 AU +SELECT + round(hill_radius(5, :t ::timestamptz)::numeric, 2) + BETWEEN 0.30 AND 0.40 AS jupiter_hill_ok; + jupiter_hill_ok +----------------- + t +(1 row) + +-- Earth Hill radius ~0.01 AU +SELECT + round(hill_radius(3, :t ::timestamptz)::numeric, 3) + BETWEEN 0.008 AND 0.012 AS earth_hill_ok; + earth_hill_ok +--------------- + t +(1 row) + +-- Lunar Hill radius (much smaller, AU) +SELECT + hill_radius_lunar(:t ::timestamptz) > 0 AS lunar_hill_positive; + lunar_hill_positive +--------------------- + t +(1 row) + +-- ============================================================ +-- Zone radius +-- ============================================================ +SELECT + lagrange_zone_radius(5, 4, :t ::timestamptz) > 0 AS jup_l4_zone_positive; + jup_l4_zone_positive +---------------------- + t +(1 row) + +SELECT + lagrange_zone_radius(5, 1, :t ::timestamptz) > 0 AS jup_l1_zone_positive; + jup_l1_zone_positive +---------------------- + t +(1 row) + +-- ============================================================ +-- Convenience functions +-- ============================================================ +-- lagrange_mass_ratio returns small positive number +SELECT + lagrange_mass_ratio(5) > 0 AND lagrange_mass_ratio(5) < 0.01 AS jupiter_mu_ok; + jupiter_mu_ok +--------------- + t +(1 row) + +SELECT + lagrange_mass_ratio(3) > 0 AND lagrange_mass_ratio(3) < 0.001 AS earth_mu_ok; + earth_mu_ok +------------- + t +(1 row) + +-- lagrange_point_name +SELECT lagrange_point_name(1) AS l1_name; + l1_name +--------- + L1 +(1 row) + +SELECT lagrange_point_name(5) AS l5_name; + l5_name +--------- + L5 +(1 row) + +-- ============================================================ +-- All planets produce valid results +-- ============================================================ +SELECT body_id, + round(helio_distance(lagrange_heliocentric(body_id, 1, :t ::timestamptz))::numeric, 2) AS sun_dist_au +FROM generate_series(1, 8) AS body_id +ORDER BY body_id; + body_id | sun_dist_au +---------+------------- + 1 | 0.46 + 2 | 0.71 + 3 | 0.97 + 4 | 1.38 + 5 | 4.63 + 6 | 8.77 + 7 | 19.44 + 8 | 29.35 +(8 rows) + +-- ============================================================ +-- Planetary moon Lagrange points +-- ============================================================ +-- Galilean: Io L4 (body=0, point=4) +SELECT + eq_ra(galilean_lagrange_equatorial(0, 4, :t ::timestamptz)) BETWEEN 0 AND 24 + AS io_l4_valid_ra; + io_l4_valid_ra +---------------- + t +(1 row) + +-- Saturn: Titan L1 (body=5, point=1) +SELECT + eq_ra(saturn_moon_lagrange_equatorial(5, 1, :t ::timestamptz)) BETWEEN 0 AND 24 + AS titan_l1_valid_ra; + titan_l1_valid_ra +------------------- + t +(1 row) + +-- Uranus: Titania L2 (body=3, point=2) +SELECT + eq_ra(uranus_moon_lagrange_equatorial(3, 2, :t ::timestamptz)) BETWEEN 0 AND 24 + AS titania_l2_valid_ra; + titania_l2_valid_ra +--------------------- + t +(1 row) + +-- Mars: Phobos L5 (body=0, point=5) +SELECT + eq_ra(mars_moon_lagrange_equatorial(0, 5, :t ::timestamptz)) BETWEEN 0 AND 24 + AS phobos_l5_valid_ra; + phobos_l5_valid_ra +-------------------- + t +(1 row) + +-- Galilean observe returns valid topocentric +SELECT + topo_elevation(galilean_lagrange_observe(2, 3, :obs ::observer, :t ::timestamptz)) + BETWEEN -90 AND 90 AS ganymede_l3_valid_el; + ganymede_l3_valid_el +---------------------- + t +(1 row) + +-- ============================================================ +-- DE fallback (no DE loaded, should produce same results as VSOP87) +-- ============================================================ +SELECT + round(helio_distance(lagrange_heliocentric_de(3, 1, :t ::timestamptz))::numeric, 4) AS de_l1_dist; + de_l1_dist +------------ + 0.9735 +(1 row) + +SELECT + round(helio_distance(lagrange_heliocentric(3, 1, :t ::timestamptz))::numeric, 4) AS vsop_l1_dist; + vsop_l1_dist +-------------- + 0.9735 +(1 row) + +-- DE hill_radius fallback +SELECT + round(hill_radius_de(5, :t ::timestamptz)::numeric, 4) = + round(hill_radius(5, :t ::timestamptz)::numeric, 4) + AS hill_de_matches_vsop; + hill_de_matches_vsop +---------------------- + t +(1 row) + +-- ============================================================ +-- Input validation +-- ============================================================ +-- Bad body_id +SELECT lagrange_heliocentric(0, 1, :t ::timestamptz); -- Sun not valid +ERROR: lagrange_heliocentric: body_id 0 must be 1-8 (Mercury-Neptune) +SELECT lagrange_heliocentric(9, 1, :t ::timestamptz); -- body 9 invalid +ERROR: lagrange_heliocentric: body_id 9 must be 1-8 (Mercury-Neptune) +-- Bad point_id +SELECT lagrange_heliocentric(3, 0, :t ::timestamptz); -- point 0 invalid +ERROR: lagrange_heliocentric: point_id 0 must be 1-5 (L1-L5) +SELECT lagrange_heliocentric(3, 6, :t ::timestamptz); -- point 6 invalid +ERROR: lagrange_heliocentric: point_id 6 must be 1-5 (L1-L5) +-- Bad lunar point_id +SELECT lunar_lagrange_equatorial(0, :t ::timestamptz); -- point 0 invalid +ERROR: lunar_lagrange_equatorial: point_id 0 must be 1-5 +SELECT lunar_lagrange_equatorial(6, :t ::timestamptz); -- point 6 invalid +ERROR: lunar_lagrange_equatorial: point_id 6 must be 1-5 +-- Bad planetary moon body_id +SELECT galilean_lagrange_equatorial(4, 1, :t ::timestamptz); -- Galilean 4 invalid +ERROR: galilean_lagrange_equatorial: body_id 4 must be 0-3 +SELECT saturn_moon_lagrange_equatorial(8, 1, :t ::timestamptz); -- Saturn 8 invalid +ERROR: saturn_moon_lagrange_equatorial: body_id 8 must be 0-7 +SELECT uranus_moon_lagrange_equatorial(5, 1, :t ::timestamptz); -- Uranus 5 invalid +ERROR: uranus_moon_lagrange_equatorial: body_id 5 must be 0-4 +SELECT mars_moon_lagrange_equatorial(2, 1, :t ::timestamptz); -- Mars 2 invalid +ERROR: mars_moon_lagrange_equatorial: body_id 2 must be 0-1 +-- lagrange_mass_ratio bad body +SELECT lagrange_mass_ratio(0); +ERROR: lagrange_mass_ratio: body_id 0 must be 1-8 +SELECT lagrange_mass_ratio(9); +ERROR: lagrange_mass_ratio: body_id 9 must be 1-8 +-- lagrange_point_name bad id +SELECT lagrange_point_name(0); +ERROR: lagrange_point_name: point_id 0 must be 1-5 +SELECT lagrange_point_name(6); +ERROR: lagrange_point_name: point_id 6 must be 1-5 diff --git a/test/sql/v020_features.sql b/test/sql/v020_features.sql new file mode 100644 index 0000000..a578ba3 --- /dev/null +++ b/test/sql/v020_features.sql @@ -0,0 +1,209 @@ +-- v020_features: Lagrange point support +-- Tests Sun-planet, Earth-Moon, planetary moon Lagrange points, +-- Hill radius, zone radius, DE fallback, and input validation. + +-- Reference observer: Greenwich, UK +\set obs '''(51.4769,-0.0005,0)''' + +-- Reference time: J2000 epoch (2000-01-01 12:00:00 UTC) +\set t '''2000-01-01 12:00:00+00''' + +-- ============================================================ +-- Sun-Earth L1/L2: should be ~0.01 AU from Earth (~1.5 million km) +-- SOHO is at L1, JWST at L2. +-- ============================================================ + +-- L1 heliocentric: should be close to Earth's heliocentric (~1 AU from Sun) +SELECT + round(helio_distance(lagrange_heliocentric(3, 1, :t ::timestamptz))::numeric, 2) AS sun_dist_au; + +-- L2 heliocentric: also ~1 AU from Sun, slightly further than L1 +SELECT + round(helio_distance(lagrange_heliocentric(3, 2, :t ::timestamptz))::numeric, 2) AS sun_dist_au; + +-- L1 between Sun and Earth (closer to Sun than L2) +SELECT + helio_distance(lagrange_heliocentric(3, 1, :t ::timestamptz)) + < + helio_distance(lagrange_heliocentric(3, 2, :t ::timestamptz)) + AS l1_closer_than_l2; + +-- ============================================================ +-- Sun-Jupiter L4/L5: ~60 degrees from Jupiter, ~5.2 AU from Sun +-- These are the Trojan asteroid zones. +-- ============================================================ + +-- L4/L5 should be ~5.2 AU from Sun +SELECT + round(helio_distance(lagrange_heliocentric(5, 4, :t ::timestamptz))::numeric, 1) AS l4_sun_dist; + +SELECT + round(helio_distance(lagrange_heliocentric(5, 5, :t ::timestamptz))::numeric, 1) AS l5_sun_dist; + +-- L4 and L5 equidistant from Sun (within 0.001 AU) +SELECT + abs( + helio_distance(lagrange_heliocentric(5, 4, :t ::timestamptz)) + - + helio_distance(lagrange_heliocentric(5, 5, :t ::timestamptz)) + ) < 0.001 AS l4_l5_equidistant; + +-- ============================================================ +-- Earth-Moon L1: ~326,000 km from Earth +-- ============================================================ + +-- lunar_lagrange_equatorial returns distance in km +SELECT + round(eq_distance(lunar_lagrange_equatorial(1, :t ::timestamptz))::numeric, -3) + BETWEEN 300000 AND 360000 AS em_l1_in_range; + +-- ============================================================ +-- lagrange_observe returns valid az/el +-- ============================================================ + +SELECT + topo_elevation(lagrange_observe(3, 2, :obs ::observer, :t ::timestamptz)) + BETWEEN -90 AND 90 AS valid_elevation; + +-- lagrange_equatorial returns valid RA/Dec +SELECT + eq_ra(lagrange_equatorial(3, 1, :t ::timestamptz)) BETWEEN 0 AND 24 AS valid_ra, + eq_dec(lagrange_equatorial(3, 1, :t ::timestamptz)) BETWEEN -90 AND 90 AS valid_dec; + +-- ============================================================ +-- lagrange_distance self-test: L-point distance to itself ≈ 0 +-- ============================================================ + +SELECT + round(lagrange_distance( + 5, 4, + lagrange_heliocentric(5, 4, :t ::timestamptz), + :t ::timestamptz + )::numeric, 10) AS self_distance; + +-- ============================================================ +-- Hill radius +-- ============================================================ + +-- Jupiter Hill radius ~0.35 AU +SELECT + round(hill_radius(5, :t ::timestamptz)::numeric, 2) + BETWEEN 0.30 AND 0.40 AS jupiter_hill_ok; + +-- Earth Hill radius ~0.01 AU +SELECT + round(hill_radius(3, :t ::timestamptz)::numeric, 3) + BETWEEN 0.008 AND 0.012 AS earth_hill_ok; + +-- Lunar Hill radius (much smaller, AU) +SELECT + hill_radius_lunar(:t ::timestamptz) > 0 AS lunar_hill_positive; + +-- ============================================================ +-- Zone radius +-- ============================================================ + +SELECT + lagrange_zone_radius(5, 4, :t ::timestamptz) > 0 AS jup_l4_zone_positive; + +SELECT + lagrange_zone_radius(5, 1, :t ::timestamptz) > 0 AS jup_l1_zone_positive; + +-- ============================================================ +-- Convenience functions +-- ============================================================ + +-- lagrange_mass_ratio returns small positive number +SELECT + lagrange_mass_ratio(5) > 0 AND lagrange_mass_ratio(5) < 0.01 AS jupiter_mu_ok; + +SELECT + lagrange_mass_ratio(3) > 0 AND lagrange_mass_ratio(3) < 0.001 AS earth_mu_ok; + +-- lagrange_point_name +SELECT lagrange_point_name(1) AS l1_name; +SELECT lagrange_point_name(5) AS l5_name; + +-- ============================================================ +-- All planets produce valid results +-- ============================================================ + +SELECT body_id, + round(helio_distance(lagrange_heliocentric(body_id, 1, :t ::timestamptz))::numeric, 2) AS sun_dist_au +FROM generate_series(1, 8) AS body_id +ORDER BY body_id; + +-- ============================================================ +-- Planetary moon Lagrange points +-- ============================================================ + +-- Galilean: Io L4 (body=0, point=4) +SELECT + eq_ra(galilean_lagrange_equatorial(0, 4, :t ::timestamptz)) BETWEEN 0 AND 24 + AS io_l4_valid_ra; + +-- Saturn: Titan L1 (body=5, point=1) +SELECT + eq_ra(saturn_moon_lagrange_equatorial(5, 1, :t ::timestamptz)) BETWEEN 0 AND 24 + AS titan_l1_valid_ra; + +-- Uranus: Titania L2 (body=3, point=2) +SELECT + eq_ra(uranus_moon_lagrange_equatorial(3, 2, :t ::timestamptz)) BETWEEN 0 AND 24 + AS titania_l2_valid_ra; + +-- Mars: Phobos L5 (body=0, point=5) +SELECT + eq_ra(mars_moon_lagrange_equatorial(0, 5, :t ::timestamptz)) BETWEEN 0 AND 24 + AS phobos_l5_valid_ra; + +-- Galilean observe returns valid topocentric +SELECT + topo_elevation(galilean_lagrange_observe(2, 3, :obs ::observer, :t ::timestamptz)) + BETWEEN -90 AND 90 AS ganymede_l3_valid_el; + +-- ============================================================ +-- DE fallback (no DE loaded, should produce same results as VSOP87) +-- ============================================================ + +SELECT + round(helio_distance(lagrange_heliocentric_de(3, 1, :t ::timestamptz))::numeric, 4) AS de_l1_dist; + +SELECT + round(helio_distance(lagrange_heliocentric(3, 1, :t ::timestamptz))::numeric, 4) AS vsop_l1_dist; + +-- DE hill_radius fallback +SELECT + round(hill_radius_de(5, :t ::timestamptz)::numeric, 4) = + round(hill_radius(5, :t ::timestamptz)::numeric, 4) + AS hill_de_matches_vsop; + +-- ============================================================ +-- Input validation +-- ============================================================ + +-- Bad body_id +SELECT lagrange_heliocentric(0, 1, :t ::timestamptz); -- Sun not valid +SELECT lagrange_heliocentric(9, 1, :t ::timestamptz); -- body 9 invalid + +-- Bad point_id +SELECT lagrange_heliocentric(3, 0, :t ::timestamptz); -- point 0 invalid +SELECT lagrange_heliocentric(3, 6, :t ::timestamptz); -- point 6 invalid + +-- Bad lunar point_id +SELECT lunar_lagrange_equatorial(0, :t ::timestamptz); -- point 0 invalid +SELECT lunar_lagrange_equatorial(6, :t ::timestamptz); -- point 6 invalid + +-- Bad planetary moon body_id +SELECT galilean_lagrange_equatorial(4, 1, :t ::timestamptz); -- Galilean 4 invalid +SELECT saturn_moon_lagrange_equatorial(8, 1, :t ::timestamptz); -- Saturn 8 invalid +SELECT uranus_moon_lagrange_equatorial(5, 1, :t ::timestamptz); -- Uranus 5 invalid +SELECT mars_moon_lagrange_equatorial(2, 1, :t ::timestamptz); -- Mars 2 invalid + +-- lagrange_mass_ratio bad body +SELECT lagrange_mass_ratio(0); +SELECT lagrange_mass_ratio(9); + +-- lagrange_point_name bad id +SELECT lagrange_point_name(0); +SELECT lagrange_point_name(6);