A code map of lambert_izzo: how a single call to lambert(&input)
flows through the crate, and where each stage lives in the source.
The crate solves Lambert's two-point boundary-value problem under
two-body gravity using D. Izzo's revisited algorithm
(arXiv:1403.2705; PDF at
docs/izzo.pdf). For an intro to what the problem
is and why three TOF regimes exist, read
docs/concepts.md first — this document picks up
where that one leaves off and points at code.
The public surface is plain f64 and [f64; 3] — no math-library
dependency. The math is dimensionally homogeneous and frame-invariant:
pass r1 / r2 in any consistent inertial frame and consistent units,
and the returned velocities come back in that same frame and units.
A call to lambert(&input) runs three top-level stages — input
validation + geometry pre-compute (from_inputs), the find_xy core
(which itself nests multi-rev enumeration → Householder iteration → TOF
evaluation), and velocity reconstruction (build_solutions).
LambertInput
│
▼
1. Input validation → geometry.rs (Geometry::from_inputs)
│
▼
2. Geometry pre-compute → geometry.rs (lambda, big_t, gamma, ρ, σ, …)
│
▼
┌──── find_xy: branch loop ────────────────────────────────────┐
│ 3. Multi-rev enumeration → root_finding.rs (find_xy loop) │
│ 4. Root finding (per M) → root_finding.rs (Householder + │
│ Halley T_min) │
│ 5. TOF evaluation (in 4) → tof.rs (Battin / Lancaster / │
│ Lagrange) │
└──────────────────────────────────────────────────────────────┘
│
▼
6. Velocity reconstruction → lib.rs (reconstruct → build_solutions)
│
▼
LambertSolutions
The driver is lambert in
crates/lambert_izzo/src/lib.rs.
Three lines:
let geom = Geometry::from_inputs(input.r1, input.r2, input.tof, input.mu, input.way)?;
let roots = find_xy(&geom, input.revolutions)?;
Ok(build_solutions(&geom, &roots))Public errors are defined in error.rs; the validation logic itself
runs at the top of Geometry::from_inputs. Each check returns a typed
LambertError variant with structured field data — pattern-match on
the value, never parse a string.
LambertErrorenum:crates/lambert_izzo/src/error.rs. Variants:NonFiniteInput,NonPositiveTimeOfFlight,NonPositiveMu,DegeneratePositionVector,CollinearGeometry,NoConvergence,SingularDenominator.- Endpoint and parameter identity (so callers can pattern-match
rather than read messages):
Position::{R1, R2}andNonFiniteParameter, both inerror.rs. - Construction-time bound on revolution counts (a separate error type,
raised before the solver is even called):
RevsOutOfRangeinbounded_revs.rs. - Validation site: the top of
Geometry::from_inputsingeometry.rs— finiteness, sign, and norm-floor checks; collinearity is checked once the cross product is available.
Geometry::from_inputs is the only constructor; the resulting
Geometry struct is then threaded read-only through stages 3–6 so the
chord, semi-perimeter, and λ are never re-derived in two places.
- Struct definition:
Geometryingeometry.rs. Fields are paper symbols (lambda,big_t,gamma,rho,sigma, plus the unit-vector basisir1,ir2,it1,it2and normsr1n,r2n). - Build site:
Geometry::from_inputsingeometry.rs. - λ sign convention (Izzo Eq. 7) and short/long-way handling: inside
from_inputs. - Round-off guards — both apply
.max(0.0)before asqrtto absorb the few-ulp negatives that appear whencrounds slightly aboves, orρslightly above1.
The Izzo formulation admits the always-present single-revolution
solution plus up to ⌊T/π⌋ multi-revolution branches, each of which
yields a long-period and a short-period root. The crate caps
enumeration at BoundedRevs::MAX = 32, type-enforced at construction
time so the bounded return collection always fits without truncation.
RevolutionBudgetenum (SingleOnly/UpTo(BoundedRevs)):lib.rs.BoundedRevsnewtype +MAX = 32:bounded_revs.rs.- Iterator over
1..=Mas validatedBoundedRevsvalues, used to drive the loop in stage 4:RevolutionBudget::iter_revsinlib.rs. - Branch loop site: inside
find_xyinroot_finding.rs. - Silent-skip semantics (the load-bearing detail): higher-
Mbranches are dropped from the returnedmultiset whenT_min(M) > tof, and the loop stops at the first infeasibleM(sinceT_minis monotone inM). Two checks gate this:- Quick reject
tof < M·π: every branch needs at leastM·πnon-dimensional time. - Boundary-zone Halley
T_minconfirmation whentoflies below the analytic minimum atx = 0.
- Quick reject
- Caller-side observability — the highest
Mactually solved:LambertSolutions::max_feasible_revsinlib.rs.
Householder's third-order method on T(x) − T_target = 0, started
from analytic guesses (Izzo Eq. 30 single-rev, Eq. 31 multi-rev). The
iteration consumes the analytic derivatives of T(x) (Eq. 22)
produced by stage 5 — see tof_derivatives_with_y below.
Eq. 30 typo caveat. The paper's Eq. 30 displays the middle-branch initial guess as
(T0/T)^(log2(T1/T0)) − 1, which fails the boundaryx = 1atT = T1. The textual derivation just before Eq. 30 (page 12) gives the correct form(T0/T)^(ln 2 / ln(T0/T1)) − 1, and the reference PyKEP implementation uses the same corrected form.lambert_izzofollows that corrected form — see the in-source comment ininitial_guess_single_rev.
- Driver:
find_xyinroot_finding.rs. Returns the always-present single-revRootplus anArrayVec<RootPair, MAX_MULTI_REV_PAIRS>of feasible multi-rev pairs. - Single-rev initial guess (Eq. 30 plus the derivative-matched
hyperbolic starter for the
T ≤ T1branch):initial_guess_single_revinroot_finding.rs. - Multi-rev initial guesses (Eq. 31, returns the long-period and
short-period asymptote starters as
(x0l, x0r)):initial_guess_multi_revinroot_finding.rs. - Householder iteration (Eq. 22):
root_finding.rs. ReturnsLambertError::SingularDenominatorif the denominator vanishes algebraically, orLambertError::NoConvergenceif theHOUSEHOLDER_MAX_ITERS = 15cap is reached. - Halley iteration on
dT/dx = 0for the multi-revT_minfeasibility check used by stage 3:halley_t_mininroot_finding.rs. Cubic convergence reachesf64precision well insideHALLEY_MAX_ITERS = 12. - All convergence tolerances and iteration caps:
constants.rs.
T(x, λ, M) is evaluated in three regimes that blend at fixed
thresholds in |x − 1|. The dispatcher picks one; the Householder
loop never sees the choice. For the why behind the regimes (which
formula loses precision where, and what the thresholds are tuned
against), read the matching section in
docs/concepts.md.
- Dispatcher:
x_to_tof_with_yincrates/lambert_izzo/src/tof.rs. - Battin hypergeometric series (Eq. 20),
|x − 1| ≤ BATTIN_THRESHOLD, single-rev only. - Lancaster–Blanchard (Eq. 18),
BATTIN_THRESHOLD < |x − 1| ≤ LAGRANGE_THRESHOLD. - Lagrange (Eq. 9), elsewhere.
- Threshold constants:
BATTIN_THRESHOLD,LAGRANGE_THRESHOLDinconstants.rs. - Analytic derivatives
(dT/dx, d²T/dx², d³T/dx³)(Eq. 22) consumed by Householder in stage 4:tof_derivatives_with_yintof.rs.
Each converged (x, y) root becomes a (v1, v2) pair via Izzo
Algorithm 1 — radial and tangential components computed from λ, ρ, σ,
γ, then projected onto the in-plane basis built in stage 2.
- Per-root reconstruction:
reconstructincrates/lambert_izzo/src/lib.rs. - Driver wrapping
reconstructover the single-rev root and every multi-rev pair, plus the parallelLambertDiagnosticsassembly:build_solutionsinlib.rs. - The output type and its diagnostics counterpart:
LambertSolutions(withmax_feasible_revs) andLambertDiagnostics, both inlib.rs.
docs/concepts.md— what Lambert's problem is, the role of the dimensionlessx, and why three TOF regimes exist.docs/izzo.pdf— D. Izzo (2014), the algorithm reference. Inline comments cite this paper asEq. N/Algorithm N.