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.
This commit is contained in:
parent
0cf55f28ac
commit
df9863dcc2
8
Makefile
8
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
|
||||
|
||||
486
src/lagrange.h
Normal file
486
src/lagrange.h
Normal file
@ -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 <math.h>
|
||||
|
||||
/* ── 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 */
|
||||
459
test/test_lagrange.c
Normal file
459
test/test_lagrange.c
Normal file
@ -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 <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;
|
||||
}
|
||||
Loading…
x
Reference in New Issue
Block a user