C++ Module: lambertSolver

Auxiliary Files

• Unit tests

Example usage

• scenarioLambertSolver

Executive Summary

This module solves Lambert’s problem using either the algorithm by R.H. Gooding (longer, extensive report available here ) or the algorithm by D. Izzo. Given two position vectors \(\mathbf{r}_{1}\) and \(\mathbf{r}_{2}\) as well as a requested time-of-flight \(t\), the solution to Lambert’s problem provides the corresponding velocity vectors v1_N and v2_N of the transfer orbit. The information about the type of algorithm (Gooding or Izzo), the position vectors, the time-of-flight as well as the number of complete revolutions around the central body, is provided by the LambertProblemMsgPayload input message. The velocity vectors are computed in the same frame as the input position vectors.

The input message’s transferDirection selects PROGRADE or RETROGRADE relative to the positive z-axis of that frame. The zero/default value is PROGRADE, preserving the behavior of existing callers. The position endpoints and flight time alone do not specify a preferred direction. Both algorithms honor this selection for zero- and multi-revolution transfers.

The algorithms by Gooding and Izzo provide solutions for elliptic, parabolic and hyperbolic transfer orbits. In the zero-revolution case, exactly one solution exists. In the multi-revolution case (meaning that the transfer orbit completes at least one complete revolution about the central body), either two or zero solutions exist, depending on whether the requested transfer time is greater or less than the minimum time-of-flight for the given number of revolutions. At the minimum time, within floating-point roundoff, the two solutions coalesce. Both output slots then contain the same orbit with both validity flags set to 1, and the root-finder iteration counts and step errors are zero. Above the minimum time, solution 1 has \(x < x_M\) and solution 2 has \(x > x_M\), where \(x_M\) is the location of the minimum of the time-of-flight curve. The LambertSolutionMsgPayload output message includes two solutions. If a solution does not exist, the corresponding velocity vectors in the message are equal to the zero vector, and the “validity” flag is equal to zero. A validity flag is set to 1 only after the root converges and both endpoint velocities are finite. If the minimum-time search cannot converge to a finite minimum, the module logs a warning and returns zero velocities with both validity flags set to zero.

Izzo Zero-Revolution Starter

For non-dimensional flight times between \(T_1=\frac{2}{3}(1-\lambda^3)\) and \(T_{00}=\arccos\lambda+\lambda\sqrt{1-\lambda^2}\), the Izzo starter interpolates linearly in \(\log(1+x)\) versus \(\log T\):

\[x_0 = \left(\frac{T_{00}}{T}\right)^{\log 2/\log(T_{00}/T_1)}-1, \qquad T_1 < T < T_{00}.\]

The positive exponent keeps the starter between zero and one and makes it continuous with the adjacent branches: \(x_0=0\) at \(T=T_{00}\) and \(x_0=1\) at \(T=T_1\). This avoids unnecessary safeguarded iterations from a negative starter for a positive root. The root convergence checks below still determine whether the resulting solution is valid.

Root Convergence

The Gooding Halley steps and Izzo Householder steps are safeguarded by a bracket for each monotonic branch of the flight-time curve. A nonfinite step, a step outside that bracket, or a high-order step that does not shrink sufficiently relative to the preceding high-order step is replaced by bisection. Bisection clears this step-size history, allowing an admissible high-order step to resume. The first high-order step also starts with no step-size history. Multi-revolution roots remain on their respective sides of the minimum-time point. For a zero-revolution hyperbolic transfer, the upper bracket is expanded at most 64 times. A finite starter at that upper endpoint is retained. Each root then receives at most 80 refinement steps.

Convergence requires a finite root in its admissible domain, a final step or computed correction no larger than \(10^{-8}\), and a relative non-dimensional flight-time residual no larger than \(10^{-12}\). A small correction with an acceptable residual is recognized before the bracket safeguard, even when the correction rounds to zero. A small correction alone does not establish convergence. An exact finite flight-time match has a zero correction.

The analytic parabolic solution uses \(x=1\) with zero iterations and zero step error when it satisfies the same relative flight-time tolerance. Nearby requests use the smooth Battin flight-time formula rather than a constant parabolic approximation, so the residual check can distinguish them.

If numerical limits prevent convergence or velocity reconstruction produces a nonfinite value, the module logs a warning. The affected solution slot has a zero validity flag, zero velocities, and zero performance fields. The other branch can still be valid. Every update starts with fresh root results, so a later numerically resolvable request can recover after a root-search failure.

The performance message reports the number of root refinements, including bisection steps, and the final step or converged correction magnitude. A valid solution can require more iterations than older solver versions allowed. Downstream acceptance policies still apply: C++ Module: lambertValidator currently requires fewer than six root iterations when its constraint checks are enabled.

Message Connection Descriptions

The following diagram and table list all the module input and output messages. The module message connection is set by the user from Python. The message type contains a link to the message structure definition, while the description provides information on what this message is used for.

Module I/O Messages

Msg Variable Name	Msg Type	Description
lambertProblemInMsg	LambertProblemMsgPayload	lambert problem setup input message
lambertSolutionOutMsg	LambertSolutionMsgPayload	lambert problem solution output message
lambertPerformanceOutMsg	LambertPerformanceMsgPayload	lambert problem performance message (additional information about the solution process)

Module Assumptions and Limitations

The input message must provide finite, positive mu and transferTime, nonnegative numRevolutions, and finite position components in both r1_N and r2_N. The endpoint norms computed in double precision must also be finite and nonzero; zero vectors and vectors whose norm calculation underflows to zero or overflows are rejected before normalization. Invalid numeric inputs raise BasiliskError before solving or publishing a new output message. This includes NaN and infinite inputs, which previously could reach geometry or root-search calculations. Position components are classified individually so quiet NaNs are rejected without raising the FE_INVALID floating-point flag through a vectorized finiteness check. Finite, valid inputs that exceed the numerical resolution of the solver still follow the convergence-failure behavior above.

solverMethod must be IZZO (0, the default for a zero-initialized message) or GOODING (1). Unknown values raise BasiliskError before solving or publishing a new output. A corrected request can recover without resetting the module. The message field stores a raw integer, which is validated before conversion to the internal SolverMethod enum. C++ callers reading the field into an enum variable must validate the integer before explicitly casting it; assigning IZZO or GOODING to the message remains supported.

setAlignmentThreshold requires a finite value strictly between 0 and 90 degrees. The default is 1 degree. An invalid value raises BasiliskError and leaves the previous setting unchanged.

The transfer angle is measured along the requested direction of motion. PROGRADE selects an orbit with positive z-component of angular momentum, while RETROGRADE selects a negative z-component. This convention is relative to the input frame, so the direction must be selected consistently when changing frames. Neither direction implies a short transfer angle: either can require an angle greater than 180 degrees.

For polar geometry, the z-axis cannot distinguish the two directions. The deterministic convention is to use \(\mathbf{r}_{1}\times\mathbf{r}_{2}\) as the orbit normal for PROGRADE (the shorter transfer angle), and its negative for RETROGRADE (the longer transfer angle). This tie includes roundoff from normalizing the position vectors: with radial unit vectors \(\mathbf{u}_1\) and \(\mathbf{u}_2\), the solver treats the geometry as polar when

\[|(\mathbf{u}_1\times\mathbf{u}_2)_z| \le 8\epsilon\left(|u_{1x}u_{2y}| + |u_{1y}u_{2x}|\right),\]

where \(\epsilon\) is double-precision machine epsilon. Outside this bound, the sign selects the direction normally, including for near-polar transfers. The message stores transferDirection as an integer so unknown values raise a Basilisk error before conversion to the internal direction enum.

Position vectors near parallel or antiparallel alignment do not reliably define a transfer plane. The module rejects inputs when the magnitude of the cross product of the radial unit vectors is below the sine of alignmentThreshold (in degrees). The module also rejects a cross product when every component is within its cancellation-roundoff bound, using the same componentwise scaling as the polar bound above. This prevents normalization roundoff from defining a spurious plane for parallel or antiparallel endpoints, even with a tiny positive threshold. Small cross products with a resolvable component remain subject to the configured angular threshold. A computed zero cross-product norm is always rejected before normalization. These checks cover angles near both 0 and 180 degrees. A rejected update returns zero velocities, cleared validity flags, and zero performance fields in both solution slots. It also skips normalization of the rejected orbit normal and evaluation of lambda and the scaled flight time, avoiding invalid square roots or division by zero for rejected collinear geometry. The alignment check is evaluated on every update: a later acceptable geometry or a relaxed threshold can produce a solution on the same module instance without calling Reset or recreating the module.

The convergence tolerances apply to the non-dimensional root and flight-time equation; they are not absolute bounds on position or velocity error. Geometry and floating-point resolution can limit accuracy. If the flight-time evaluator cannot resolve a request to the convergence tolerances, the result is invalid.

User Guide

The module is first initialized as follows:

lambertModule = lambertSolver.LambertSolver()
lambertModule.ModelTag = "lambertSolver"
lambertModule.setAlignmentThreshold(1.0) # module defaults this value to 1.0 degrees if not specified
unitTestSim.AddModelToTask(unitTaskName, lambertModule)

The lambert problem input message is either created as a standalone message in python

lambertProblemInMsgData = messaging.LambertProblemMsgPayload()
lambertProblemInMsgData.solverMethod = messaging.IZZO
lambertProblemInMsgData.r1_N = np.array([10000. * 1000, 0. ,0.])
lambertProblemInMsgData.r2_N = np.array([0., 8000. * 1000,0.])
lambertProblemInMsgData.transferTime = 10000.
lambertProblemInMsgData.mu = 3.986004418e14
lambertProblemInMsgData.numRevolutions = 0
lambertProblemInMsg = messaging.LambertProblemMsg().write(lambertProblemInMsgData)

or obtained from another FSW module. The lambert problem input message is then connected.

lambertModule.lambertProblemInMsg.subscribeTo(lambertProblemInMsg)

To request a retrograde transfer, set the direction and write the updated payload to the connected problem message:

lambertProblemInMsgData.transferDirection = messaging.RETROGRADE
lambertProblemInMsg.write(lambertProblemInMsgData)

Use messaging.PROGRADE to select the default explicitly. A zero-initialized LambertProblemMsgPayload already selects prograde motion. C and C++ producers must initialize the new transferDirection field along with the other inputs.

---

class LambertSolver : public SysModel

#include <lambertSolver.h>

This module solves Lambert’s problem using either the Gooding or the Izzo algorithm.

Public Functions

LambertSolver()

This is the constructor for the module class. It sets default variable values and initializes the various parts of the model

~LambertSolver()

Module Destructor

void Reset(uint64_t currentSimNanos) override

This method is used to reset the module and checks that required input messages are connected.

Parameters:

currentSimNanos – current simulation time in nano-seconds

void UpdateState(uint64_t currentSimNanos) override

This is the main method that gets called every time the module is updated. It computes the solution of Lambert’s problem.

Parameters:

currentSimNanos – current simulation time in nano-seconds

void setAlignmentThreshold(const double value)

Set the minimum angular separation from parallel or antiparallel alignment.

Note

A rejected value leaves the previous threshold unchanged.

Parameters:

value – Threshold in degrees, finite and strictly between 0 and 90.

Throws:

BasiliskError – If the threshold is outside its valid range.

inline double getlignmentThreshold() const

getter for alignmentThreshold

Public Members

ReadFunctor<LambertProblemMsgPayload> lambertProblemInMsg

lambert problem input message

Message<LambertSolutionMsgPayload> lambertSolutionOutMsg

lambert solution output message

Message<LambertPerformanceMsgPayload> lambertPerformanceOutMsg

lambert performance output message

BSKLogger bskLogger

BSK Logging.

Private Functions

void readMessages()

This method reads the input messages each call of updateState. It also checks if the message contents are valid for this module.

void writeMessages(uint64_t currentSimNanos)

This method writes the output messages each call of updateState

Parameters:

currentSimNanos – current simulation time in nano-seconds

void problemGeometry()

This method computes the problem geometry for the given parameters of Lambert’s problem. The orbit frame is also determined.

std::array<double, 2> goodingInitialGuess(double lambda, double T)

This method computes the initial guess for the free variable x using the Gooding procedure

Parameters:

• lam – lambda parameter that defines the problem geometry

• T – non-dimensional time-of-flight

Returns:

std::array<double, 2>

std::array<double, 2> izzoInitialGuess(double lambda, double T)

This method computes the initial guess for the free variable x using the Izzo procedure

Parameters:

• lam – lambda parameter that defines the problem geometry

• T – non-dimensional time-of-flight

Returns:

std::array<double, 2>

void findx()

Find the free variable on each applicable branch of the time-of-flight curve.

std::array<Eigen::Vector3d, 2> computeVelocities(double x)

This method computes the velocities at the initial and final position for a given free variable x

Parameters:

x – free variable of Lambert’s problem that satisfies the given time of flight

Returns:

std::array<Eigen::Vector3d, 2>

double x2tof(double x, int N, double lam)

This method computes the non-dimensional time of flight (TOF) for a given x

Parameters:

• x – free variable of Lambert’s problem that satisfies the given time of flight

• N – number of revolutions

• lam – lambda parameter that defines the problem geometry

Returns:

double

std::array<double, 3> dTdx(double x, double T, double lam)

This method computes the derivatives of the time of flight curve T(x)

Parameters:

• x – free variable of Lambert’s problem that satisfies the given time of flight

• T – requested non-dimensional time-of-flight

• lam – lambda parameter that defines the problem geometry

Returns:

std::array<double, 3>

RootResult solveRoot(double T, double x0, int N, double xLower, double xUpper, bool increasing)

Solve one monotonic flight-time branch with safeguarded Halley or Householder steps.

Note

Bisection replaces nonfinite, out-of-bracket, or insufficiently shrinking high-order steps. The shrink check compares consecutive high-order steps; bisection clears that history. Success requires both a small step or correction and a small relative flight-time residual, or an exact finite flight-time match. Numerical failure never constitutes convergence.

Parameters:

• T – Requested positive non-dimensional flight time.

• x0 – Initial root estimate; an out-of-bracket estimate is replaced by the midpoint.

• N – Number of complete revolutions.

• xLower – Lower bound of the root’s branch, including the open elliptic endpoint -1.

• xUpper – Upper bound, including the open endpoint 1 for a multi-revolution branch.

• increasing – True for the right multi-revolution branch; false for decreasing branches.

Returns:

Last root estimate, iteration count, step or correction magnitude, and explicit convergence status.

std::optional<std::array<double, 2>> getTmin(double T0M, int N)

Locate the multi-revolution minimum time using safeguarded Halley iterations.

Note

Halley steps remain inside a derivative-sign bracket. Bisection is used when a step leaves the bracket or does not sufficiently decrease in size. Convergence requires a finite time and a derivative residual consistent with the tolerance on x_M.

Parameters:

• T0M – Non-dimensional time of flight at x = 0 for N revolutions.

• N – Number of revolutions, which must be positive.

Returns:

The converged location x_M and minimum time T(x_M), or std::nullopt if the search fails.

double hypergeometricF(double z)

This method computes the hypergeometric function 2F1(a,b,c,z)

Parameters:

z – argument of hypergeometric function

Returns:

double

Private Members

double alignmentThreshold = 1.0

[deg] minimum angle between position vectors so they are not too aligned.

SolverMethod solverMethod

lambert solver algorithm (GOODING or IZZO)

LambertTransferDirection transferDirection = PROGRADE

Requested direction about the input frame’s +z axis.

Eigen::Vector3d r1_N

position vector at t0

Eigen::Vector3d r2_N

position vector at t1

double transferTime = {}

time of flight between r1_N and r2_N (t1-t0)

double mu = {}

gravitational parameter

int numberOfRevolutions = {}

number of revolutions

double TOF = {}

non-dimensional time-of-flight constraint

double lambda = {}

parameter of Lambert”s problem that defines problem geometry

bool multiRevSolution = {}

boolean flag if multi-revolution solutions exist or not

double xMin = {}

[-] Minimum-time root location separating the multi-revolution branches.

bool noSolution = {}

boolean flag if no solution should be returned (in case of 180 deg transfer angle)

std::array<Eigen::Vector3d, 3> Oframe1

array containing the orbit frame unit vectors at t0

std::array<Eigen::Vector3d, 3> Oframe2

array containing the orbit frame unit vectors at t1

std::array<Eigen::Vector3d, 2> vvecs

array containing the velocity vector solutions at t0 and t1

std::array<Eigen::Vector3d, 2> vvecsSol2

array containing the velocity vector solutions at t0 and t1 (sol 2)

RootResult solution = {}

Numerical result for the first solution.

RootResult solutionSol2 = {}

Numerical result for the second solution.

struct RootResult

Numerical result for one root of the flight-time equation.

Public Members

double x = {}

[-] Last evaluated root candidate.

int iterations = {}

Number of root iterations.

double error = {}

[-] Magnitude of the final step or converged correction in x.

bool converged = {}

True only after finite-domain, step, and residual checks pass.