1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
// SPDX-License-Identifier: AGPL-3.0-only
//! Externally validate the kshana CR3BP single-shooting differential corrector
//! against the **NASA/JPL Three-Body Periodic Orbit Database** (Solar System
//! Dynamics / Dynamical Systems group, `periodic_orbits.api` v1.0) for the
//! Earth-Moon **L2 Southern Halo (NRHO)** family — specifically the **9:2
//! lunar-synodic NRHO**, the NASA Gateway reference orbit, plus 4 catalog
//! neighbours spanning the 9:2 regime (6.53 .. 6.59 days).
//!
//! The CR3BP equations of motion are a one-parameter (mass ratio `mu`) family,
//! so the non-dimensional **period T**, **Jacobi constant C**, the perpendicular-
//! crossing **initial state** `{x0, z0, vy0}` and the non-dimensional **perilune
//! radius** are frame-independent invariants of the periodic orbit. kshana's
//! `mu = 0.012150585609624` matches the JPL system `mass_ratio` to 12 sig figs,
//! so the catalog row is a genuine external authority for those invariants.
//!
//! Each fixture row carries the JPL catalog state. The test feeds kshana's
//! `differential_correct_halo` the JPL `x0` (held fixed, as the corrector
//! requires) and a *perturbed* `{z0, vy0}` seed, lets the single-shooting STM
//! corrector converge, and compares the converged orbit's invariants against the
//! catalog. The perilune oracle in the fixture is computed by an **independent
//! integrator** (scipy DOP853, rtol=atol=1e-13) from the JPL state — so the
//! perilune is a two-integrator cross-check, not a self-check.
//!
//! UNITS (load-bearing for honesty): the perilune is validated in the JPL length
//! unit `lunit = 389703.265 km`. kshana::cr3bp hard-codes `EARTH_MOON_DIST_KM =
//! 384400 km`, a 1.4% *labelling* choice; the test converts kshana's perilune
//! back to non-dim (÷384400) and into the JPL unit before comparing, so the gap
//! measured is dynamics, not the unit convention.
//!
//! TOLERANCES (the residuals are the honest, reported fixed-step-RK4 /
//! single-shooting / finite-grid-perilune gap, NOT loosened to force a pass):
//! - Jacobi C: |dC| <= 5e-4 (measured worst ~1.5e-5)
//! - Period T: |dT|/T <= 1e-3 (measured worst ~9.5e-5)
//! - IC components: |dz0|,|dvy0|<= 1e-3 (measured worst ~1.4e-5)
//! - Perilune (JPL km): |dperi| <= 150 km (measured worst ~1.2 km)
//!
//! HONEST SCOPE: this validates the differential-correction CORE of
//! `kshana::cr3bp` against the published JPL catalog. It does NOT validate an
//! ephemeris (DE) cislunar model, the de-normalised MCI/MCMF transforms, or
//! station-keeping — those are separate follow-ons.
//!
//! Reference data, provenance, the JPL query URL and the committed generator
//! live in `tests/fixtures/cislunar_mission_analysis/`.
use kshana::cr3bp::{differential_correct_halo, Cr3bpState, EARTH_MOON_DIST_KM, EARTH_MOON_MU};
const REF: &str =
include_str!("fixtures/cislunar_mission_analysis/cislunar_mission_analysis_reference.txt");
// Tolerances on the non-dimensional CR3BP invariants (see header).
const TOL_C: f64 = 5e-4; // Jacobi constant, absolute
const TOL_T_REL: f64 = 1e-3; // period, relative
const TOL_IC: f64 = 1e-3; // {z0, vy0} initial-state components, absolute non-dim
const TOL_PERI_KM: f64 = 150.0; // perilune radius in the JPL length unit (km)
struct Member {
name: String,
x0: f64,
z0: f64,
vy0: f64,
c: f64,
t: f64,
peri_km_jpl: f64,
}
fn parse_system_lunit(text: &str) -> f64 {
for line in text.lines() {
if let Some(rest) = line.strip_prefix("# SYSTEM ") {
for tok in rest.split_whitespace() {
if let Some(v) = tok.strip_prefix("lunit_km=") {
return v.parse().expect("lunit_km parse");
}
}
}
}
panic!("no '# SYSTEM ... lunit_km=' line in fixture");
}
fn parse_members(text: &str) -> Vec<Member> {
let mut out = Vec::new();
for line in text.lines() {
if !line.starts_with("NRHO ") {
continue;
}
// NRHO name | x0 | z0 | vy0 | jacobi_C | period_T | period_days | peri_nondim | peri_km_jpl
let p: Vec<&str> = line.splitn(9, '|').collect();
assert_eq!(p.len(), 9, "NRHO row needs 9 |-fields: {line}");
let f = |s: &str| -> f64 { s.trim().parse().expect("float parse") };
out.push(Member {
name: p[0].trim_start_matches("NRHO").trim().to_string(),
x0: f(p[1]),
z0: f(p[2]),
vy0: f(p[3]),
c: f(p[4]),
t: f(p[5]),
peri_km_jpl: f(p[8]),
});
}
out
}
#[test]
fn cr3bp_corrector_reproduces_jpl_l2_southern_nrho_9to2() {
let mu = EARTH_MOON_MU;
let lunit_jpl = parse_system_lunit(REF);
let members = parse_members(REF);
assert!(
members.len() >= 3,
"expected >=3 NRHO reference members (9:2 + neighbours), got {}",
members.len()
);
let mut saw_92 = false;
let mut worst_c = 0.0_f64;
let mut worst_t_rel = 0.0_f64;
let mut worst_ic = 0.0_f64;
let mut worst_peri = 0.0_f64;
for m in &members {
if m.name == "L2S_NRHO_9to2" {
saw_92 = true;
}
// Seed with the JPL x0 held fixed (the corrector parametrises by x0) and a
// 2%-perturbed {z0, vy0} so the corrector genuinely has to converge, not
// merely echo the catalog state back.
let guess = Cr3bpState {
r: [m.x0, 0.0, m.z0 * 0.98],
v: [0.0, m.vy0 * 0.98, 0.0],
};
let orbit = differential_correct_halo(&guess, mu, 1e-11, 200)
.unwrap_or_else(|| panic!("{}: kshana corrector did not converge", m.name));
// x0 is held fixed by construction — assert it stayed put.
assert!(
(orbit.ic.r[0] - m.x0).abs() < 1e-12,
"{}: x0 drifted ({} vs {})",
m.name,
orbit.ic.r[0],
m.x0
);
// (1) Jacobi constant.
let dc = (orbit.jacobi - m.c).abs();
worst_c = worst_c.max(dc);
assert!(
dc <= TOL_C,
"{}: Jacobi C {:.10} vs JPL {:.10} (|dC|={:.2e} > {:.0e})",
m.name,
orbit.jacobi,
m.c,
dc,
TOL_C
);
// (2) Period (relative).
let dt_rel = (orbit.period - m.t).abs() / m.t;
worst_t_rel = worst_t_rel.max(dt_rel);
assert!(
dt_rel <= TOL_T_REL,
"{}: period T {:.10} vs JPL {:.10} (|dT|/T={:.2e} > {:.0e})",
m.name,
orbit.period,
m.t,
dt_rel,
TOL_T_REL
);
// (3) Initial-state components {z0, vy0} (x0 already checked).
for (lbl, got, want) in [("z0", orbit.ic.r[2], m.z0), ("vy0", orbit.ic.v[1], m.vy0)] {
let d = (got - want).abs();
worst_ic = worst_ic.max(d);
assert!(
d <= TOL_IC,
"{}: {} {:.10} vs JPL {:.10} (|d|={:.2e} > {:.0e})",
m.name,
lbl,
got,
want,
d,
TOL_IC
);
}
// (4) Perilune radius. kshana returns km in its 384400-km convention;
// convert to non-dim then into the JPL length unit before comparing to
// the independent (scipy DOP853) oracle value.
let peri_kshana_384400 = orbit.perilune_radius_km(mu, 1200);
let peri_nondim = peri_kshana_384400 / EARTH_MOON_DIST_KM;
let peri_km_jpl = peri_nondim * lunit_jpl;
let dperi = (peri_km_jpl - m.peri_km_jpl).abs();
worst_peri = worst_peri.max(dperi);
assert!(
dperi <= TOL_PERI_KM,
"{}: perilune {:.2} km vs JPL/scipy {:.2} km (|d|={:.2} km > {:.0} km)",
m.name,
peri_km_jpl,
m.peri_km_jpl,
dperi,
TOL_PERI_KM
);
}
assert!(
saw_92,
"the 9:2 NRHO primary-gate member (L2S_NRHO_9to2) must be present and validated"
);
eprintln!(
"cislunar_mission_analysis_reference: {} JPL L2-S NRHO members vs JPL catalog \
(9:2 primary gate) — worst |dC|={:.2e}, worst |dT|/T={:.2e}, worst |dIC|={:.2e}, \
worst |dperilune|={:.2} km (JPL lunit)",
members.len(),
worst_c,
worst_t_rel,
worst_ic,
worst_peri
);
}