warehack.ing/pg_orrery
2
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
pg_orrery/src/lagrange.h
Ryan Malloy 024c0c1e0c Harden Newton-Raphson gamma bounds, improve v0.20.0 test coverage 
Add positive-gamma clamp in L1/L2/L3 Newton-Raphson iterations to
prevent divergence on extreme mass ratios. Add missing CREATE EXTENSION,
tighter L1/L2 precision checks (4 decimal places), lagrange_distance_oe
test with Ceres, L1-Earth-L2 ordering verification, and DE fallback
tests for planetary moon Lagrange functions.
2026-02-28 19:09:08 -07:00

493 lines
15 KiB
C
Raw Blame History

/*
* 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 (gamma_new <= 0.0)
gamma_new = gamma * 0.5; /* keep gamma positive */
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 (gamma_new <= 0.0)
gamma_new = gamma * 0.5; /* keep gamma positive */
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 (gamma_new <= 0.0)
gamma_new = gamma * 0.5; /* keep gamma positive */
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 */
Reference in New Issue View Git Blame Copy Permalink