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; +}