Ryan Malloy
df9863dcc2
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.
460 lines
14 KiB
C
460 lines
14 KiB
C
/*
|
|
* 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 <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <math.h>
|
|
|
|
/* ── 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;
|
|
}
|