Lambert's problem is a fundamental challenge in orbital mechanics that involves determining the transfer orbit between two position vectors given a specific travel time. This problem is essential for mission planning, intercept maneuvers, and interplanetary trajectory design. The most common numerical approach utilizes universal variables and Stumpff functions to solve the governing transcendental equations.

Mathematical Foundations

Given an initial position vector $\mathbf{r}_1$, a final position vector $\mathbf{r}_2$, and the time of flight $\Delta t$, we seek the velocity vectors $\mathbf{v}_1$ and $\mathbf{v}_2$. The solution typically involves finding the root of the time-of-flight equation:

$$\sqrt{\mu}\Delta t = \frac{x^3 S(z) + A\sqrt{y}}{C(z)}$$

Where $S(z)$ and $C(z)$ are Stumpff functions, and $A$ is a geometric parameter determined by the transfer angle.

MATLAB Implementation

The following implementation uses a Newton-Raphson iteration to solve for the universal variable $z$, which characterizes the energy of the orbit.

function [v_start, v_end] = lambert_solver(r1, r2, dt, mu, orbit_direction)
    % r1, r2: Position vectors (1x3)
    % dt: Time of flight (seconds)
    % mu: Gravitational parameter
    % orbit_direction: 'prograde' or 'retrograde'

    mag_r1 = norm(r1);
    mag_r2 = norm(r2);
    
    % Calculate the change in true anomaly
    cos_delta_theta = dot(r1, r2) / (mag_r1 * mag_r2);
    
    % Ensure the value is within acos range
    cos_delta_theta = max(-1, min(1, cos_delta_theta));
    d_theta = acos(cos_delta_theta);
    
    % Adjust for retrograde motion
    if strcmp(orbit_direction, 'retrograde')
        d_theta = 2*pi - d_theta;
    end
    
    % Geometric parameter A
    A = sin(d_theta) * sqrt(mag_r1 * mag_r2 / (1 - cos(d_theta)));
    
    % Initial guess for z
    z_val = 0;
    tolerance = 1e-10;
    max_steps = 100;
    
    for i = 1:max_steps
        [S, C] = calculate_stumpff(z_val);
        
        % Calculate y(z)
        y_val = mag_r1 + mag_r2 + A * (z_val * S - 1) / sqrt(C);
        
        % Ensure y is positive for valid square roots
        if A > 0 && y_val < 0
            % Adjust z if the guess leads to an invalid region
            z_val = z_val + 0.5;
            continue;
        end
        
        % Universal variable x
        x_val = sqrt(y_val / C);
        
        % Calculated time of flight
        t_calc = (x_val^3 * S + A * sqrt(y_val)) / sqrt(mu);
        
        % Convergence check
        error = t_calc - dt;
        if abs(error) < tolerance
            break;
        end
        
        % Derivative calculation (Newton-Raphson)
        if z_val == 0
            derivative = (sqrt(2)/40) * y_val^1.5 + (A/8) * (sqrt(y_val) + A*sqrt(1/(2*y_val)));
        else
            derivative = (x_val^3 * (1/(2*z_val) * (C - 3*S/(2*C))) + 3*S*x_val^2/(2*C)) / sqrt(mu);
        end
        
        % Update z
        z_val = z_val - error / derivative;
    end
    
    % Compute Lagrange coefficients
    f_coeff = 1 - y_val / mag_r1;
    g_coeff = A * sqrt(y_val / mu);
    g_dot_coeff = 1 - y_val / mag_r2;
    
    % Resulting velocities
    v_start = (r2 - f_coeff * r1) / g_coeff;
    v_end = (g_dot_coeff * r2 - r1) / g_coeff;
end

function [S, C] = calculate_stumpff(z)
    % Computes Stumpff functions S(z) and C(z)
    if z > 0
        sz = sqrt(z);
        S = (sz - sin(sz)) / (sz^3);
        C = (1 - cos(sz)) / z;
    elseif z < 0
        sz = sqrt(-z);
        S = (sinh(sz) - sz) / (sz^3);
        C = (cosh(sz) - 1) / (-z);
    else
        S = 1/6;
        C = 1/2;
    end
end

Core Components

1. Stumpff Functions

The Stumpff functions $S(z)$ and $C(z)$ are series expansions that generalize the behvaior of orbits. They allow the solver to handle elliptical ($z > 0$), parabolic ($z = 0$), and hyperbolic ($z < 0$) trajectories within a single mathematical framework, avoiding singularities at the transition points.

2. Iterative Root-Finding

The solver iterates on the variable $z$. While the Newton-Raphson method is used here for its quadratic convergence, more robust solvers might employ bisection or Brent's method to handle cases where the derivative becomes unstable, particularly near the parabolic limit.

3. Lagrange Coefficients

Once the universal variable $z$ (and consequently $y$) is found, the velocity vectors are derived using the $f$ and $g$ functions. These functions represent the mapping between the state vectors at two different time epochs under two-body dynamics.

Practical Applications

Hohmann Transfer Verification

Lambert solvers are frequently used to verify classical maneuvers. For an Earth-to-Mars transfer, the algorithm calculates the required injection velocity given the planetary ephemerides at departure and arrival.

Satellite Intercepts

For Rendezvous, Proximity Operations, and Docking (RPOD), Lambert's problem calculates the impulse required to shift an interceptor from its current orbit to a target location in a fixed amount of time.

Optimization Considerations

• Multi-Revolution Solutions: The standard algorithm solves for the "short way" or "long way" within a single orbit. Advanced versions can find multiple solutions for transfers spanning several complete revolutions.
• Coordinate Systems: Ensure position vectors are in an inertial frame (e.g., J2000 or ECI) before performing the calculation.
• Numerical Stabiltiy: For very small transfer angles, the calculation of $A$ can become sensitive. Using alternative formulations, such as those by Battin or Gooding, can improve robustness in edge cases.