lambert-rs 0.0.18

Rust + PyO3 implementation of the Izzo Lambert solver

lambert_rs

High-performance Rust + PyO3 implementation of the Izzo Lambert solver with advanced optimization capabilities.

Installation

From PyPI

pip install lambert-rs

Dependencies: No required dependencies. The package optionally uses pandas if available for optimize_lambert_nm_multi_array (returns DataFrame), otherwise returns a dictionary.

Features

• Fast Lambert Problem Solving: Vectorized batch processing with parallel execution
• Multiple Solution Methods: Single solutions, grid search, and Nelder-Mead optimization
• Delta-V Calculations: Automatic computation of transfer impulses
• Two-Body Propagation: Universal variable formulation for orbit propagation
• Transfer Optimization: Find optimal departure/arrival times minimizing delta-v
• Rendezvous Optimization: Optimize for full rendezvous (two burns) or transfer only (one burn)
• Coordinate Frame Conversions: ECI ↔ RIC conversions with direction cosine matrices and angular velocity calculations

Quick Start

Basic Lambert Problem

Solve a single Lambert problem:

import numpy as np
import lambert_rs

# Initial and final positions (km)
r1 = np.array([7000., 0., 0.])
r2 = np.array([0., -8000., 0.])

# Time of flight (seconds)
tof = 3600.0

# Solve Lambert problem
sol = lambert_rs.lambert_izzo_single(
    r1, r2, tof,
    max_rev=0,        # Maximum number of revolutions
    retrograde=False, # Prograde transfer
    mu=398600.4418,  # Earth's gravitational parameter (km^3/s^2)
    tol=1e-8,        # Tolerance
    maxiter=50       # Maximum iterations
)

# Get the first valid solution
v1 = sol.v1[0]
v2 = sol.v2[0]
print(f"v1: {v1} km/s")
print(f"v2: {v2} km/s")
# Access Keplerian elements of the transfer orbit
print(f"Transfer orbit SMA: {sol.sma[0]:.2f} km")
print(f"Eccentricity: {sol.ecc[0]:.4f}")
print(f"Inclination: {sol.inc[0]:.2f} deg")
print(f"True anomaly at departure: {sol.ta_depart[0]:.2f} deg")
print(f"True anomaly at arrival: {sol.ta_arrive[0]:.2f} deg")

Batch Processing

Process multiple Lambert problems in parallel:

import numpy as np
import lambert_rs

# Multiple position pairs: shape (N, 3)
r1_batch = np.array([
    [7000., 0., 0.],
    [8000., 0., 0.],
    [9000., 0., 0.],
])

r2_batch = np.array([
    [0., -8000., 0.],
    [0., -9000., 0.],
    [0., -10000., 0.],
])

# Multiple time-of-flight values
tof_array = np.array([3600., 7200., 10800.])

# Solve all combinations
sol = lambert_rs.lambert_izzo_vec(
    r1_batch, r2_batch, tof_array,
    max_rev=0,
    retrograde=False,
    mu=398600.4418,
    tol=1e-8,
    maxiter=50
)

# Results shape: (N, num_tof, num_solutions, 3) for v1/v2
# Keplerian elements have same shape as tof: (N, num_tof, num_solutions)
print(f"v1 shape: {sol.v1.shape}")
print(f"v2 shape: {sol.v2.shape}")
print(f"sma shape: {sol.sma.shape}")
print(f"ecc shape: {sol.ecc.shape}")
print(f"inc shape: {sol.inc.shape}")

Delta-V Calculations

Calculate transfer impulses for rendezvous:

import numpy as np
import lambert_rs

# Initial and target states: [x, y, z, vx, vy, vz]
s1 = np.array([7000., 0., 0., 0., 7.5, 0.])  # Initial state
s2 = np.array([0., -8000., 0., 8.0, 0., 0.])  # Target state

# Time of flight
tof = 3600.0

# Calculate delta-v for all solutions
dv1_mag, dv2_mag, dv_total, valid, dv1_vec, dv2_vec = lambert_rs.lambert_izzo_vec_dv(
    s1, s2, tof,
    max_rev=0,
    retrograde=False,
    mu=398600.4418,
    tol=1e-8,
    maxiter=50
)

# Find best solution (minimum total delta-v)
if np.any(valid):
    valid_indices = np.where(valid)[0]
    best_idx = valid_indices[np.nanargmin(dv_total[valid])]
    print(f"Best delta-v: {dv_total[best_idx]:.3f} km/s")
    print(f"First burn: {dv1_vec[best_idx]} km/s")
    print(f"Second burn: {dv2_vec[best_idx]} km/s")

Reentry Check

Filter out solutions that would reenter Earth's atmosphere:

import numpy as np
import lambert_rs

# Positions that might result in a reentry trajectory
r1 = np.array([7000., 0., 0.])  # km
r2 = np.array([0., -8000., 0.])  # km
tof = 1800.0  # Short time of flight might cause reentry

# Solve without reentry check
sol_no_check = lambert_rs.lambert_izzo_single(
    r1, r2, tof,
    max_rev=0,
    retrograde=False,
    mu=398600.4418,
    reentry_check=False  # Default: don't check for reentry
)

# Solve with reentry check enabled
sol_with_check = lambert_rs.lambert_izzo_single(
    r1, r2, tof,
    max_rev=0,
    retrograde=False,
    mu=398600.4418,
    reentry_check=True  # Filter out solutions that would reenter
)

# Count valid solutions (non-NaN velocities)
valid_no_check = ~np.isnan(sol_no_check.v1[:, 0])
valid_with_check = ~np.isnan(sol_with_check.v1[:, 0])
print(f"Solutions without check: {np.sum(valid_no_check)} valid")
print(f"Solutions with check: {np.sum(valid_with_check)} valid")

# The reentry check marks solutions as invalid if the perigee radius
# would be below Earth's surface (6378.137 km)

Find Best Delta-V Solution

Automatically find the best solution (prograde or retrograde):

import numpy as np
import lambert_rs

# States: shape (..., 6) for batch processing
s1 = np.array([7000., 0., 0., 0., 7.5, 0.])
s2 = np.array([0., -8000., 0., 8.0, 0., 0.]) 

# Time of flight (can be scalar or array)
tof = np.array([3600., 7200., 10800.])

# Full rendezvous: minimize total delta-v (dv1 + dv2)
# Use this when you need to match both position AND velocity at the target
dv1_rend, dv2_rend, dv_total_rend, valid_rend, dv1_vec_rend, dv2_vec_rend = lambert_rs.lambert_izzo_best_dv(
    s1, s2, tof,
    max_rev=0,
    mu=398600.4418,
    tol=1e-8,
    maxiter=50,
    rendezvous=True  # Default: minimize dv1 + dv2
)

# Transfer-only: minimize first impulse (dv1 only)
# Use this when you only need to reach the target position (not matching velocity)
dv1_trans, dv2_trans, dv_total_trans, valid_trans, dv1_vec_trans, dv2_vec_trans = lambert_rs.lambert_izzo_best_dv(
    s1, s2, tof,
    max_rev=0,
    mu=398600.4418,
    tol=1e-8,
    maxiter=50,
    rendezvous=False  # Minimize dv1 only
)

# Automatically selects best solution (prograde or retrograde)
print(f"Rendezvous total delta-v: {dv_total_rend} km/s")
print(f"Transfer-only total delta-v: {dv1_trans} km/s")

Two-Body Orbit Propagation

Propagate an orbit using universal variable formulation:

import numpy as np
import lambert_rs

# Initial state: [x, y, z, vx, vy, vz] (km, km/s)
initial_state = np.array([7000., 0., 0., 0., 7.5, 0.])

# Propagation time (seconds)
dt = 3600.0

# Propagate
r_final, v_final = lambert_rs.kepler_universal_2body_py(initial_state, dt)

print(f"Final position: {r_final} km")
print(f"Final velocity: {v_final} km/s")

Transfer Optimization (Grid Search)

Find optimal departure and arrival times using grid search:

import numpy as np
import lambert_rs

# Source and target initial states
source_state = np.array([7000., 0., 0., 0., 7.5, 0.])
target_state = np.array([0., -8000., 0., 8.0, 0., 0.]) 

# Time bounds: (min_departure, max_arrival) in seconds
time_bounds = (0., 86400.)  # 24 hours
dt_step = 3600.0  # 1 hour steps

# Find optimal transfer
sol = lambert_rs.optimize_lambert_transfer_py(
    source_state, target_state,
    time_bounds, dt_step,
    # mu=398600.4418,  # Optional: default is 398600.4418 (Earth)
    # retrograde=False  # Optional: default is False
)

# LambertBurnArray returns arrays - access first element for single solution
print(f"Optimal departure time: {sol.t0[0]:.1f} s")
print(f"Optimal arrival time: {sol.tf[0]:.1f} s")
print(f"Time of flight: {sol.tof[0]:.1f} s")
print(f"Minimum delta-v: {sol.cost[0]:.3f} km/s")
print(f"First burn: {sol.dv1[0]} km/s")
print(f"Second burn: {sol.dv2[0]} km/s")
# Keplerian elements of the transfer orbit
print(f"Semi-major axis: {sol.sma[0]:.2f} km")
print(f"Eccentricity: {sol.ecc[0]:.4f}")
print(f"Inclination: {sol.inc[0]:.2f} deg")
print(f"RAAN: {sol.raan[0]:.2f} deg")
print(f"Argument of Perigee: {sol.argp[0]:.2f} deg")
print(f"True Anomaly (departure): {sol.ta_depart[0]:.2f} deg")
print(f"True Anomaly (arrival): {sol.ta_arrive[0]:.2f} deg")
print(f"Radius of perigee: {sol.rp[0]:.2f} km")
print(f"Radius of apogee: {sol.ra[0]:.2f} km")

Transfer Grid (All Solutions)

Get all transfer solutions in a time window:

import numpy as np
import lambert_rs

source_state = np.array([7000., 0., 0., 0., 7.5, 0.])
target_state = np.array([0., -8000., 0., 8.0, 0., 0.]) 

time_bounds = (0., 86400.)
dt_step = 3600.0

# Get all solutions
sols = lambert_rs.optimize_lambert_transfer_grid_py(
    source_state, target_state,
    time_bounds, dt_step,
    # mu=398600.4418,  # Optional: default is 398600.4418 (Earth)
    # retrograde=False  # Optional: default is False
)

# Plot or analyze all solutions
import matplotlib.pyplot as plt
plt.scatter(sols.t0, sols.tf, c=sols.cost, cmap='viridis')
plt.colorbar(label='Delta-V (km/s)')
plt.xlabel('Departure Time (s)')
plt.ylabel('Arrival Time (s)')
plt.show()

Grid Search with Best Delta-V Selection

Use lambert_grid_search for a grid search that automatically selects the best solution (prograde or retrograde):

import numpy as np
import lambert_rs

source_state = np.array([7000., 0., 0., 0., 7.5, 0.])
target_state = np.array([0., -8000., 0., 8.0, 0., 0.]) 

time_bounds = (0., 86400.)

# Example 1: Full rendezvous (both burns, default)
# Optimizes for total delta-v (dv1 + dv2)
sols_rendezvous = lambert_rs.lambert_grid_search(
    source_state, target_state,
    time_bounds,
    dt_step=3600.0,  # 1 hour steps
    max_rev=0,        # Maximum revolutions (default: 0)
    rendezvous=True,  # Full rendezvous: optimize for dv1 + dv2 (default)
    # mu=398600.4418,  # Optional: default is 398600.4418 (Earth)
    # tol=1e-8,        # Optional: tolerance for iterative solver (default: 1e-8)
    # maxiter=50,      # Optional: maximum iterations (default: 50)
)

# Analyze all solutions
print(f"Found {len(sols_rendezvous.cost)} rendezvous solutions")
print(f"Best cost: {np.min(sols_rendezvous.cost):.3f} km/s")

# Get the best solution
best_idx = np.argmin(sols_rendezvous.cost)
print(f"\nBest rendezvous solution:")
print(f"  Times: t0={sols_rendezvous.t0[best_idx]:.1f} s, tf={sols_rendezvous.tf[best_idx]:.1f} s")
print(f"  Total cost (dv1 + dv2): {sols_rendezvous.cost[best_idx]:.3f} km/s")
print(f"  First burn: {sols_rendezvous.dv1[best_idx]} km/s")
print(f"  Second burn: {sols_rendezvous.dv2[best_idx]} km/s")

# Example 2: Transfer only (first burn)
# Optimizes for first burn delta-v only (dv1)
sols_transfer = lambert_rs.lambert_grid_search(
    source_state, target_state,
    time_bounds,
    dt_step=3600.0,  # 1 hour steps
    max_rev=0,
    rendezvous=False,  # Transfer only: optimize for dv1 only
)

# Analyze transfer solutions
print(f"\nFound {len(sols_transfer.cost)} transfer solutions")
print(f"Best cost: {np.min(sols_transfer.cost):.3f} km/s")

best_transfer_idx = np.argmin(sols_transfer.cost)
print(f"\nBest transfer solution:")
print(f"  Times: t0={sols_transfer.t0[best_transfer_idx]:.1f} s, tf={sols_transfer.tf[best_transfer_idx]:.1f} s")
print(f"  Cost (dv1 only): {sols_transfer.cost[best_transfer_idx]:.3f} km/s")
print(f"  First burn: {sols_transfer.dv1[best_transfer_idx]} km/s")
print(f"  Second burn: {sols_transfer.dv2[best_transfer_idx]} (NaN for transfer-only)")

# Compare to optimize_lambert_transfer_grid_py:
# - lambert_grid_search: automatically selects prograde/retrograde, can compute both impulses or transfer only
# - optimize_lambert_transfer_grid_py: requires retrograde flag, only computes transfer (dv1)

Nelder-Mead Optimization

Use gradient-free optimization to find optimal transfer:

import numpy as np
import lambert_rs

source_state = np.array([7000., 0., 0., 0., 7.5, 0.])
target_state = np.array([0., -8000., 0., 8.0, 0., 0.]) 

time_bounds = (0., 86400.)

# Optimize for transfer (single burn)
# Automatically selects best solution from both prograde and retrograde
sol = lambert_rs.optimize_lambert_nm_py(
    source_state, target_state,
    time_bounds,
    max_rev=0,
    optimize_rendezvous=False  # Transfer only (one burn)
    # mu=398600.4418,  # Optional: default is 398600.4418 (Earth)
    # max_iters=None,  # Optional: default is None (uses 1000 internally)
    # initial_guess=None  # Optional: (t0, tf) guess
)

# LambertBurnArray returns arrays - access first element for single solution
print(f"Optimal times: t0={sol.t0[0]:.1f} s, tf={sol.tf[0]:.1f} s")
print(f"Time of flight: {sol.tof[0]:.1f} s")
print(f"Cost (delta-v): {sol.cost[0]:.3f} km/s")
print(f"First burn: {sol.dv1[0]} km/s")

# Optimize for rendezvous (two burns)
# Automatically selects best solution from both prograde and retrograde
sol = lambert_rs.optimize_lambert_nm_py(
    source_state, target_state,
    time_bounds,
    max_rev=0,
    optimize_rendezvous=True,  # Full rendezvous (two burns)
    # mu=398600.4418,  # Optional: default is 398600.4418 (Earth)
    # max_iters=None,  # Optional: default is None (uses 1000 internally)
)

# LambertBurnArray returns arrays - access first element for single solution
print(f"Rendezvous cost: {sol.cost[0]:.3f} km/s")
print(f"Time of flight: {sol.tof[0]:.1f} s")
print(f"First burn: {sol.dv1[0]} km/s")
print(f"Second burn: {sol.dv2[0]} km/s")
# Access Keplerian elements
print(f"Transfer orbit SMA: {sol.sma[0]:.2f} km, Ecc: {sol.ecc[0]:.4f}")

Multi-Start Optimization

Find multiple local optima using parallel multi-start optimization:

import numpy as np
import lambert_rs

source_state = np.array([7000., 0., 0., 0., 7.5, 0.])
target_state = np.array([0., -8000., 0., 8.0, 0., 0.]) 

time_bounds = (0., 86400.)

# Run multiple optimizers in parallel
# Automatically selects best solution from both prograde and retrograde for each optimizer
sols = lambert_rs.optimize_lambert_nm_multi_py(
    source_state, target_state,
    time_bounds,
    max_rev=0,
    optimize_rendezvous=True,
    num_solvers=10,        # Number of parallel optimizers
    return_best_only=False, # Return all solutions
    remove_duplicates=True,  # Remove duplicate solutions
    # mu=398600.4418,  # Optional: default is 398600.4418 (Earth)
    # max_iters=None,  # Optional: default is None (uses 1000 internally)
)

# Analyze all solutions
print(f"Found {len(sols.cost)} solutions")
# Access Keplerian elements for all solutions
print(f"Transfer orbit SMAs: {sols.sma} km")
print(f"Eccentricities: {sols.ecc}")
print(f"Inclinations: {sols.inc} deg")
print(f"Best cost: {np.min(sols.cost):.3f} km/s")
print(f"Worst cost: {np.max(sols.cost):.3f} km/s")

# Get the best solution
best_idx = np.argmin(sols.cost)
print(f"\nBest solution:")
print(f"  Times: t0={sols.t0[best_idx]:.1f} s, tf={sols.tf[best_idx]:.1f} s")
print(f"  Time of flight: {sols.tof[best_idx]:.1f} s")
print(f"  Cost: {sols.cost[best_idx]:.3f} km/s")
print(f"  Transfer orbit SMA: {sols.sma[best_idx]:.2f} km, Ecc: {sols.ecc[best_idx]:.4f}")
print(f"  First burn: {sols.dv1[best_idx]} km/s")
print(f"  Second burn: {sols.dv2[best_idx]} km/s")

Array-based Multi-Start Optimization

For analyzing multiple source and target states simultaneously, use optimize_lambert_nm_multi_array:

import numpy as np
import lambert_rs

# Generate arrays of source and target states
# Example: 100 source states and 20 target states
num_sources = 100
num_targets = 20

# Create random orbital states (in practice, these would be your actual states)
source_states = np.random.randn(num_sources, 6) * 1000  # [N, 6] array
target_states = np.random.randn(num_targets, 6) * 1000  # [M, 6] array

# Time bounds for optimization
time_bounds = (0.0, 86400.0)  # 0 to 1 day in seconds

# Run optimization for all combinations
# This will compute 100 * 20 = 2000 combinations
result = lambert_rs.optimize_lambert_nm_multi_array(
    source_states, target_states,
    time_bounds,
    max_rev=0,
    optimize_rendezvous=True,
    num_solvers=10,        # Number of parallel optimizers per combination
    return_best_only=True, # Return only best solution per [i,j] combination
    remove_duplicates=True,
    # mu=398600.4418,  # Optional: default is 398600.4418 (Earth)
    # max_iters=None,  # Optional: default is None (uses 1000 internally)
)

# Result is a pandas DataFrame if pandas is available, otherwise a dictionary
# Columns/keys:
# - source_idx: index of source state (0 to num_sources-1)
# - target_idx: index of target state (0 to num_targets-1)
# - t0: departure time [s]
# - tf: arrival time [s]
# - cost: total delta-v [km/s]
# - dv1_x, dv1_y, dv1_z: first burn impulse components [km/s]
# - dv2_x, dv2_y, dv2_z: second burn impulse components [km/s]

# If pandas is available, result is a DataFrame
try:
    import pandas as pd
    if isinstance(result, pd.DataFrame):
        print(f"Total solutions: {len(result)}")
        print(f"\nBest solution:")
        best = result.sort_values("cost").head(1)
        print(best)
        
        # Find best solution for each source state
        best_per_source = result.sort_values("cost").groupby("source_idx").first()
        print(f"\nBest solution per source state:")
        print(best_per_source)
        
        # Filter solutions with cost < 5 km/s
        low_cost = result[result["cost"] < 5.0]
        print(f"\nSolutions with cost < 5 km/s: {len(low_cost)}")
except ImportError:
    pass

# If pandas is not available, result is a dictionary
if isinstance(result, dict):
    import numpy as np
    print(f"Total solutions: {len(result['cost'])}")
    
    # Find best solution
    best_idx = np.argmin(result['cost'])
    print(f"\nBest solution:")
    print(f"  source_idx: {result['source_idx'][best_idx]}")
    print(f"  target_idx: {result['target_idx'][best_idx]}")
    print(f"  cost: {result['cost'][best_idx]:.3f} km/s")
    
    # Filter solutions with cost < 5 km/s
    cost_array = np.array(result['cost'])
    low_cost_mask = cost_array < 5.0
    print(f"\nSolutions with cost < 5 km/s: {np.sum(low_cost_mask)}")

# When return_best_only=False, you can get multiple solutions per [i,j] combination
result_all = lambert_rs.optimize_lambert_nm_multi_array(
    source_states, target_states,
    time_bounds,
    max_rev=0,
    optimize_rendezvous=True,
    num_solvers=10,
    return_best_only=False, # Return all solutions per [i,j] combination
    remove_duplicates=True,
    # mu=398600.4418,  # Optional: default is 398600.4418 (Earth)
    # max_iters=None,  # Optional: default is None (uses 1000 internally)
)

# Analyze all solutions for a specific source-target pair
if isinstance(result_all, dict):
    import numpy as np
    source_idx = 0
    target_idx = 5
    mask = (np.array(result_all['source_idx']) == source_idx) & \
           (np.array(result_all['target_idx']) == target_idx)
    indices = np.where(mask)[0]
    costs = np.array(result_all['cost'])[indices]
    sorted_indices = indices[np.argsort(costs)]
    print(f"\nAll solutions for source[{source_idx}] -> target[{target_idx}]:")
    for idx in sorted_indices:
        print(f"  cost: {result_all['cost'][idx]:.3f} km/s, "
              f"t0: {result_all['t0'][idx]:.1f} s, tf: {result_all['tf'][idx]:.1f} s")
else:
    # pandas DataFrame
    source_0_target_5 = result_all[
        (result_all["source_idx"] == 0) & (result_all["target_idx"] == 5)
    ]
    print(f"\nAll solutions for source[0] -> target[5]:")
    print(source_0_target_5.sort_values("cost"))

Coordinate Frame Conversions

Convert between Earth-Centered Inertial (ECI) and Radial-In-Track-Cross (RIC) frames:

ECI to RIC State Conversion

Convert deputy spacecraft state from ECI to RIC relative state:

import numpy as np
import lambert_rs

# Deputy spacecraft state in ECI: [x, y, z, vx, vy, vz] (km, km/s)
dep_state_eci = np.array([7000., 0., 0., 0., 7.5, 0.])

# Reference spacecraft state in ECI: [x, y, z, vx, vy, vz] (km, km/s)
ref_state_eci = np.array([6800., 100., 0., 0., 7.6, 0.])

# Convert to RIC relative state
rel_state_ric = lambert_rs.eci2ric_state(dep_state_eci, ref_state_eci)

print(f"Relative position in RIC: {rel_state_ric[:3]} km")
print(f"Relative velocity in RIC: {rel_state_ric[3:]} km/s")

RIC to ECI State Conversion

Convert deputy spacecraft state from RIC back to ECI:

import numpy as np
import lambert_rs

# Relative state in RIC: [dr_r, dr_i, dr_c, dv_r, dv_i, dv_c] (km, km/s)
rel_state_ric = np.array([200., 0., 0., 0.1, 0., 0.])

# Reference spacecraft state in ECI
ref_state_eci = np.array([6800., 100., 0., 0., 7.6, 0.])

# Convert back to ECI
dep_state_eci = lambert_rs.ric2eci_state(rel_state_ric, ref_state_eci)

print(f"Deputy position in ECI: {dep_state_eci[:3]} km")
print(f"Deputy velocity in ECI: {dep_state_eci[3:]} km/s")

Batch Processing

Process multiple states at once:

import numpy as np
import lambert_rs

# Multiple deputy states: shape (N, 6)
dep_states = np.array([
    [7000., 0., 0., 0., 7.5, 0.],
    [7100., 0., 0., 0., 7.4, 0.],
    [7200., 0., 0., 0., 7.3, 0.],
])

# Single reference state (broadcasted to all deputies)
ref_state = np.array([6800., 100., 0., 0., 7.6, 0.])

# Convert all to RIC
rel_states_ric = lambert_rs.eci2ric_state(dep_states, ref_state)

print(f"Output shape: {rel_states_ric.shape}")  # (N, 6)

Direction Cosine Matrices

Get the rotation matrices for coordinate transformations:

import numpy as np
import lambert_rs

# Reference state in ECI
ref_state_eci = np.array([7000., 0., 0., 0., 7.5, 0.])

# Get DCM from ECI to RIC
dcm_eci2ric = lambert_rs.dcm_eci2ric(ref_state_eci)
print(f"DCM shape: {dcm_eci2ric.shape}")  # (3, 3)

# Get DCM from RIC to ECI (transpose of above)
dcm_ric2eci = lambert_rs.dcm_ric2eci(ref_state_eci)

# Verify they are transposes
assert np.allclose(dcm_ric2eci, dcm_eci2ric.T)

# Rotate a vector from ECI to RIC
vec_eci = np.array([1., 0., 0.])
vec_ric = dcm_eci2ric @ vec_eci
print(f"Vector in RIC: {vec_ric}")

Vector Rotation (ECI ↔ RIC)

Rotate position vectors between ECI and RIC frames:

import numpy as np
import lambert_rs

# Reference state defining the RIC frame
ref_state_eci = np.array([7000., 0., 0., 0., 7.5, 0.])

# ECI position vector to rotate
pos_eci = np.array([1000., 500., 200.])

# Rotate to RIC
pos_ric = lambert_rs.eci2ric_vec(pos_eci, ref_state_eci)
print(f"Position in RIC: {pos_ric}")

# Rotate back to ECI
pos_eci_recovered = lambert_rs.ric2eci_vec(pos_ric, ref_state_eci)
assert np.allclose(pos_eci, pos_eci_recovered, rtol=1e-10)

# Batch processing: multiple positions with single reference
positions_eci = np.array([
    [1000., 500., 200.],
    [2000., 1000., 400.],
    [3000., 1500., 600.],
])
positions_ric = lambert_rs.eci2ric_vec(positions_eci, ref_state_eci)
print(f"Positions in RIC shape: {positions_ric.shape}")  # (3, 3)

Note: These functions rotate raw position vectors. For converting relative states (which account for angular velocity effects), use eci2ric_state and ric2eci_state instead.

Angular Velocity Computation

Compute angular velocity vector from orbital state:

import numpy as np
import lambert_rs

# State vector: [x, y, z, vx, vy, vz] (km, km/s)
state = np.array([7000., 0., 0., 0., 7.5, 0.])

# Compute angular velocity: omega = cross(r, v) / ||r||^2
omega = lambert_rs.compute_angular_vel_rad_p_sec(state)

print(f"Angular velocity: {omega} rad/s")
print(f"Magnitude: {np.linalg.norm(omega):.6e} rad/s")

# For multiple states
states = np.array([
    [7000., 0., 0., 0., 7.5, 0.],
    [8000., 0., 0., 0., 7.0, 0.],
])
omegas = lambert_rs.compute_angular_vel_rad_p_sec(states)
print(f"Angular velocities shape: {omegas.shape}")  # (N, 3)

Round-Trip Conversion Test

Verify conversions are reversible:

import numpy as np
import lambert_rs

# Original deputy state
dep_state_original = np.array([7000., 0., 0., 0., 7.5, 0.])
ref_state = np.array([6800., 100., 0., 0., 7.6, 0.])

# Convert to RIC and back
rel_state_ric = lambert_rs.eci2ric_state(dep_state_original, ref_state)
dep_state_recovered = lambert_rs.ric2eci_state(rel_state_ric, ref_state)

# Should match original (within numerical precision)
assert np.allclose(dep_state_original, dep_state_recovered, rtol=1e-10)
print("Round-trip conversion successful!")

Sun Position Computation

Compute the Sun's ECI position for given Julian date(s) or POSIX time(s):

Using Julian Dates

import numpy as np
import lambert_rs

# Single Julian date (J2000.0 epoch)
jd = 2451545.0
sun_state = lambert_rs.compute_sun_state_jd(jd)

# Result is (1, 6) array: [X, Y, Z, NaN, NaN, NaN]
# Position in km, velocity is NaN (position-only approximation)
print(f"Sun position: {sun_state[0, :3]} km")
print(f"Distance: {np.linalg.norm(sun_state[0, :3]):.2f} km")

# Multiple Julian dates
jd_array = np.array([2451545.0, 2451546.0, 2451547.0])
sun_states = lambert_rs.compute_sun_state_jd(jd_array)

# Result is (3, 6) array
print(f"Sun positions shape: {sun_states.shape}")
print(f"Sun positions:\n{sun_states[:, :3]}")

Using POSIX Time

import numpy as np
import lambert_rs

# Single POSIX time (seconds since 1970-01-01 00:00:00 UTC)
# Example: 2000-01-01 12:00:00 UTC
posix = 946728000.0
sun_state = lambert_rs.compute_sun_state_posix(posix)

# Result is (1, 6) array: [X, Y, Z, NaN, NaN, NaN]
print(f"Sun position: {sun_state[0, :3]} km")

# Multiple POSIX times
posix_array = np.array([946728000.0, 946814400.0, 946900800.0])
sun_states = lambert_rs.compute_sun_state_posix(posix_array)

# Result is (3, 6) array
print(f"Sun positions shape: {sun_states.shape}")

Note: This is an approximation based on simplified orbital mechanics. The velocity components are set to NaN as this function only computes position. Both functions produce identical results when given equivalent times.

RIC Waypoint Offsets

Apply RIC waypoint offsets to target states for Lambert targeting. This allows you to target positions relative to a reference state in the RIC frame, which is useful for formation flying, rendezvous planning, and relative navigation scenarios.

Applying RIC Waypoint Offsets

Apply a waypoint offset to a target state:

import numpy as np
import lambert_rs

# Target state in ECI
target_state = np.array([7000., 0., 0., 0., 7.5, 0.])

# RIC waypoint offset: [dr_r, dr_i, dr_c] (position only, velocity offset is zero)
ric_waypoint = np.array([100.0, 50.0, 25.0])  # 100 km radial, 50 km in-track, 25 km cross-track

# Apply waypoint offset
modified_target = lambert_rs.apply_ric_waypoint_offset(target_state, ric_waypoint)

# Full RIC waypoint offset: [dr_r, dr_i, dr_c, dv_r, dv_i, dv_c] (position + velocity)
ric_waypoint_full = np.array([100.0, 50.0, 25.0, 0.01, 0.02, 0.005])
modified_target_full = lambert_rs.apply_ric_waypoint_offset(target_state, ric_waypoint_full)

Using RIC Waypoints with Optimization

Target a waypoint relative to the target state using optimization:

import numpy as np
import lambert_rs

source_state = np.array([7000., 0., 0., 0., 7.5, 0.])
target_state = np.array([0., -8000., 0., 8.0, 0., 0.])
time_bounds = (0., 86400.)

# RIC waypoint offset: 100 km ahead in-track, 50 km radial
ric_waypoint = np.array([50.0, 100.0, 0.0])

# Optimize with waypoint offset
sol = lambert_rs.optimize_lambert_nm_py(
    source_state,
    target_state,
    time_bounds,
    max_rev=0,
    optimize_rendezvous=True,
    mu=398600.4418,
    ric_waypoint=ric_waypoint,  # Target waypoint relative to target_state
)

# LambertBurnArray returns arrays - access first element for single solution
print(f"Optimal cost: {sol.cost[0]:.3f} km/s")
print(f"Departure time: {sol.t0[0]:.1f} s")
print(f"Arrival time: {sol.tf[0]:.1f} s")
print(f"Time of flight: {sol.tof[0]:.1f} s")

Batch Processing with Waypoints

Apply single waypoint to multiple target states:

import numpy as np
import lambert_rs

# Multiple target states
target_states = np.array([
    [7000., 0., 0., 0., 7.5, 0.],
    [8000., 0., 0., 0., 7.0, 0.],
    [9000., 0., 0., 0., 6.5, 0.],
])

# Single waypoint (applied to all targets)
ric_waypoint = np.array([100.0, 50.0, 25.0])

modified_targets = lambert_rs.apply_ric_waypoint_offset(target_states, ric_waypoint)
print(f"Modified targets shape: {modified_targets.shape}")  # (3, 6)

Apply array of waypoints (one per target state):

import numpy as np
import lambert_rs

# Multiple target states
target_states = np.array([
    [7000., 0., 0., 0., 7.5, 0.],
    [8000., 0., 0., 0., 7.0, 0.],
])

# Array of waypoints (one per target)
ric_waypoints = np.array([
    [100.0, 50.0, 25.0],  # Waypoint for first target
    [200.0, 100.0, 50.0], # Waypoint for second target
])

modified_targets = lambert_rs.apply_ric_waypoint_offset(target_states, ric_waypoints)

Array Optimization with Waypoints

Use waypoints with array-based optimization:

import numpy as np
import lambert_rs

# Multiple source and target states
source_states = np.array([
    [7000., 0., 0., 0., 7.5, 0.],
    [7100., 0., 0., 0., 7.4, 0.],
])
target_states = np.array([
    [0., -8000., 0., 8.0, 0., 0.],
    [0., -9000., 0., 7.5, 0., 0.],
])

time_bounds = (0., 86400.)

# Single waypoint (applied to all target states)
ric_waypoint = np.array([100.0, 50.0, 25.0])

result = lambert_rs.optimize_lambert_nm_multi_array(
    source_states,
    target_states,
    time_bounds,
    max_rev=0,
    optimize_rendezvous=True,
    num_solvers=10,
    return_best_only=True,
    mu=398600.4418,
    ric_waypoint=ric_waypoint,  # Single waypoint for all targets
)

# Or use array of waypoints (one per target)
ric_waypoints = np.array([
    [100.0, 50.0, 25.0],
    [200.0, 100.0, 50.0],
])

result = lambert_rs.optimize_lambert_nm_multi_array(
    source_states,
    target_states,
    time_bounds,
    max_rev=0,
    optimize_rendezvous=True,
    num_solvers=10,
    return_best_only=True,
    mu=398600.4418,
    ric_waypoint=ric_waypoints,  # Array of waypoints
)

Note: The waypoint offset is applied relative to the target state at the arrival time. The RIC frame is defined by the target state itself, so the waypoint represents an offset from the target's position and velocity in its own RIC frame.

Classical Orbital Elements

Convert between ECI state vectors (position and velocity) and classical orbital elements (COEs).

State Vector to Classical Orbital Elements

Convert ECI state vector(s) to classical orbital elements:

import numpy as np
import lambert_rs

# Single state vector: [x, y, z, vx, vy, vz] (km, km/s)
state = np.array([7000., 0., 0., 0., 7.5, 0.])

# Convert to COEs (angles in degrees by default)
coes = lambert_rs.rv2coe(state, deg=True, mu=398600.4418, tol=1e-8, special=False)

# Output: [sma, ecc, inc, raan, aop, ta]
# sma: semi-major axis (km)
# ecc: eccentricity
# inc: inclination (degrees)
# raan: right ascension of ascending node (degrees)
# aop: argument of perigee (degrees)
# ta: true anomaly (degrees)
print(f"Semi-major axis: {coes[0]:.2f} km")
print(f"Eccentricity: {coes[1]:.6f}")
print(f"Inclination: {coes[2]:.2f} deg")
print(f"RAAN: {coes[3]:.2f} deg")
print(f"AOP: {coes[4]:.2f} deg")
print(f"True Anomaly: {coes[5]:.2f} deg")

# Convert with angles in radians
coes_rad = lambert_rs.rv2coe(state, deg=False)
print(f"Inclination: {coes_rad[2]:.6f} rad")

# Get special elements (argument of latitude, true longitude, longitude of periapsis)
coes_special = lambert_rs.rv2coe(state, special=True)
# Output: [sma, ecc, inc, raan, aop, ta, arglat, truelon, lonper]
print(f"Special elements shape: {coes_special.shape}")  # (9,)

Batch Processing

Process multiple states at once:

import numpy as np
import lambert_rs

# Multiple state vectors: shape (N, 6)
states = np.array([
    [7000., 0., 0., 0., 7.5, 0.],
    [8000., 0., 0., 0., 7.0, 0.],
    [9000., 0., 0., 0., 6.5, 0.],
])

# Convert all to COEs
coes = lambert_rs.rv2coe(states, deg=True)
print(f"COEs shape: {coes.shape}")  # (N, 6)

# Access individual elements
print(f"All semi-major axes: {coes[:, 0]} km")
print(f"All eccentricities: {coes[:, 1]}")

Classical Orbital Elements to State Vector

Convert classical orbital elements to ECI state vector(s):

import numpy as np
import lambert_rs

# Classical orbital elements: [sma, ecc, inc, raan, aop, ta]
# Angles in degrees by default
coes = np.array([7000.0, 0.0, 45.0, 30.0, 60.0, 90.0])

# Convert to ECI state vector
state = lambert_rs.coe2rv(coes, deg=True, mu=398600.4418)

# Output: [x, y, z, vx, vy, vz] (km, km/s)
print(f"Position: {state[:3]} km")
print(f"Velocity: {state[3:]} km/s")

# Convert with angles in radians
coes_rad = np.array([7000.0, 0.0, np.deg2rad(45.0), np.deg2rad(30.0), 
                     np.deg2rad(60.0), np.deg2rad(90.0)])
state_rad = lambert_rs.coe2rv(coes_rad, deg=False)

Batch Processing

Process multiple COE sets at once:

import numpy as np
import lambert_rs

# Multiple COE sets: shape (N, 6)
coes_batch = np.array([
    [7000.0, 0.0, 45.0, 30.0, 60.0, 90.0],
    [8000.0, 0.1, 50.0, 35.0, 65.0, 100.0],
    [9000.0, 0.2, 55.0, 40.0, 70.0, 110.0],
])

# Convert all to state vectors
states = lambert_rs.coe2rv(coes_batch, deg=True)
print(f"States shape: {states.shape}")  # (N, 6)

Round-Trip Conversion Test

Verify conversions are reversible:

import numpy as np
import lambert_rs

# Original state vector
state_original = np.array([7000., 0., 0., 0., 7.5, 0.])

# Convert to COEs and back
coes = lambert_rs.rv2coe(state_original, deg=True)
state_recovered = lambert_rs.coe2rv(coes, deg=True)

# Should match original (within numerical precision)
# Note: For circular/equatorial orbits, some angles may be set to zero
# by convention, so exact match may not be possible
print(f"Original: {state_original}")
print(f"Recovered: {state_recovered}")
print(f"Difference: {np.abs(state_original - state_recovered)}")

Note:

• For circular orbits (eccentricity ≈ 0), argument of perigee and true anomaly are set to zero by convention
• For equatorial orbits (inclination ≈ 0° or 180°), right ascension of ascending node is set to zero by convention
• The special parameter in rv2coe returns additional elements (argument of latitude, true longitude, longitude of periapsis) that are used for special cases
• All angle conversions handle quadrant checks automatically using arctangent functions for numerical stability

Sunline Targeting

The optimize_lambert_nm_multi_sunline_py function targets a position offset from the target state along the target-to-sun vector. This is useful for missions that need to approach a target while maintaining a specific offset along the sunline direction.

import numpy as np
import lambert_rs
from datetime import datetime

# Define source and target states
source_state = np.array([6778.137, 0.0, 0.0, 0.0, 7.66855818, 0.0])  # LEO
target_state = np.array([42164.0, 0.0, 0.0, 0.0, 3.07466628, 0.0])  # GEO

# Time bounds (relative to t0_ref)
time_bounds = (0.0, 86400.0)  # 0 to 1 day in seconds

# Reference epoch: 2024-01-01 00:00:00 UTC
# Convert to POSIX timestamp
t0_ref = datetime(2024, 1, 1).timestamp()

# Sunline offset: 200 km along target-to-sun vector
sun_vector_offset_km = 200.0

# Optimize
sols = lambert_rs.optimize_lambert_nm_multi_sunline_py(
    source_state,
    target_state,
    time_bounds,
    max_rev=0,
    optimize_rendezvous=True,
    num_solvers=20,
    return_best_only=True,
    sun_vector_offset_km=sun_vector_offset_km,
    t0_ref=t0_ref,
    mu=398600.4418,
    remove_duplicates=True,
    max_iters=100,
)

# Get best solution
best_idx = np.argmin(sols.cost)
print(f"Best solution:")
print(f"  Departure time: {sols.t0[best_idx]:.1f} s (relative to t0_ref)")
print(f"  Arrival time: {sols.tf[best_idx]:.1f} s (relative to t0_ref)")
print(f"  Total delta-v: {sols.cost[best_idx]:.6f} km/s")
print(f"  First burn: {np.linalg.norm(sols.dv1[best_idx]):.6f} km/s")
print(f"  Second burn: {np.linalg.norm(sols.dv2[best_idx]):.6f} km/s")

Note: The time_bounds are relative to t0_ref. For example, if t0_ref is the POSIX timestamp for 2024-01-01 00:00:00 UTC and time_bounds is (0, 86400), then the optimization window is from 2024-01-01 00:00:00 UTC to 2024-01-02 00:00:00 UTC. The function computes the sun's position at the absolute arrival time (t0_ref + tf) to determine the target-to-sun vector.

Solution Types

The library uses two main solution types:

LambertSolutionArray

• Used for raw Lambert problem solutions (rarely used directly)
• Contains velocities (v1, v2), time of flight (tof), and orbital parameters
• Includes Keplerian elements of the transfer orbit:
  • sma, ecc, inc, raan, argp (semi-major axis, eccentricity, inclination, RAAN, argument of perigee)
  • ta_depart, ta_arrive (true anomaly at departure and arrival)
  • rp, ra (radius of perigee and apogee)

LambertBurnArray

• Used for optimized transfer solutions with computed impulses
• Replaces the previous OptimizedLambertSolutionPy and OptimizedLambertSolutionArray types
• Contains all fields from LambertSolutionArray plus:
  • t0, tf, tof (departure time, arrival time, time of flight)
  • cost (total delta-v)
  • dv1, dv2 (first and second burn impulses)
  • All Keplerian elements of the transfer orbit
• Always returns arrays, even for single solutions (access with [0] for single results)

API Reference

Core Functions

• lambert_izzo_single(r1, r2, tof, max_rev, retrograde, mu, tol, maxiter)

  • Solve a single Lambert problem
  • Returns: LambertSolutionArray object with fields:
    • v1: Velocity at r1 [km/s] - shape (n, 3)
    • v2: Velocity at r2 [km/s] - shape (n, 3)
    • tof: Time of flight [s] - shape (n,)
    • num_revs: Number of revolutions - shape (n,)
    • path_flag: Path flag - shape (n,)
    • retrograde: Whether solution uses retrograde motion - shape (n,)
    • sma: Semi-major axis [km] - shape (n,)
    • ecc: Eccentricity - shape (n,)
    • inc: Inclination [deg] - shape (n,)
    • raan: Right Ascension of Ascending Node [deg] - shape (n,)
    • argp: Argument of Perigee [deg] - shape (n,)
    • ta_depart: True anomaly at departure [deg] - shape (n,)
    • ta_arrive: True anomaly at arrival [deg] - shape (n,)
    • rp: Radius of perigee [km] - shape (n,)
    • ra: Radius of apogee [km] - shape (n,)

• lambert_izzo_vec(r1, r2, tof, max_rev, retrograde, mu, tol, maxiter)

  • Vectorized batch processing
  • Returns: LambertSolutionArray object with same fields as above

• lambert_izzo_vec_dv(s1, s2, tof, max_rev, retrograde, mu, tol, maxiter)

  • Calculate delta-v for transfers
  • Returns: (dv1_mag, dv2_mag, dv_total, valid, dv1_vec, dv2_vec)

• lambert_izzo_best_dv(s1, s2, tof, max_rev, mu, tol, maxiter, rendezvous)

  • Automatically find best solution (prograde or retrograde)
  • rendezvous=True (default): Minimize total delta-v (dv1 + dv2) for full rendezvous
  • rendezvous=False: Minimize first impulse (dv1) for transfer-only missions
  • Returns: (dv1_mag, dv2_mag, dv_total, valid, dv1_vec, dv2_vec)

Propagation

• kepler_universal_2body_py(init_state_vec, dt_sec)
  • Two-body orbit propagation using universal variables
  • Returns: (r_final, v_final)

Celestial Mechanics

• compute_sun_state_jd(jd)

  • Compute the Sun's ECI XYZ position vector for given Julian date(s)
  • Parameters:
    • jd: Julian date(s) as scalar f64 or array of f64 values
  • Returns: Array with shape (N, 6) where N is the number of input times
    • Columns are X, Y, Z positions (km) and velocity (NaN, NaN, NaN)
    • Velocity is set to NaN as this is a position-only approximation
  • Based on: https://www.mathworks.com/matlabcentral/fileexchange/78766-approxecisunposition

• compute_sun_state_posix(posix)

  • Compute the Sun's ECI XYZ position vector for given POSIX time(s)
  • Parameters:
    • posix: POSIX time(s) in seconds since 1970-01-01 00:00:00 UTC as scalar f64 or array of f64 values
  • Returns: Array with shape (N, 6) where N is the number of input times
    • Columns are X, Y, Z positions (km) and velocity (NaN, NaN, NaN)
    • Velocity is set to NaN as this is a position-only approximation
  • Internally converts POSIX time to Julian date and calls compute_sun_state_jd
  • Based on: https://www.mathworks.com/matlabcentral/fileexchange/78766-approxecisunposition

Coordinate Frame Conversions

• eci2ric_state(dep_state_eci, ref_state_eci)

  • Convert deputy ECI state to RIC relative state
  • Parameters:
    • dep_state_eci: Deputy state in ECI [..., 6] (position + velocity)
    • ref_state_eci: Reference state in ECI [..., 6] (position + velocity)
  • Returns: Relative state in RIC [..., 6] (relative position + relative velocity)
  • Supports broadcasting: if one input is scalar (shape [6]) and the other is array (shape [N, 6]), the scalar is broadcasted

• ric2eci_state(rel_state_ric, ref_state_eci)

  • Convert deputy RIC relative state to ECI state
  • Parameters:
    • rel_state_ric: Relative state in RIC [..., 6] (relative position + relative velocity)
    • ref_state_eci: Reference state in ECI [..., 6] (position + velocity)
  • Returns: Deputy state in ECI [..., 6] (position + velocity)
  • Supports broadcasting: if one input is scalar (shape [6]) and the other is array (shape [N, 6]), the scalar is broadcasted

• dcm_eci2ric(ref_state_eci)

  • Build direction cosine matrix (DCM) that rotates ECI vectors into RIC
  • Parameters:
    • ref_state_eci: Reference state in ECI [..., 6] (position + velocity)
  • Returns: Rotation matrix [..., 3, 3] such that v_ric = dcm @ v_eci

• dcm_ric2eci(ref_state_eci)

  • Build direction cosine matrix (DCM) that rotates RIC vectors into ECI
  • Parameters:
    • ref_state_eci: Reference state in ECI [..., 6] (position + velocity)
  • Returns: Rotation matrix [..., 3, 3] such that v_eci = dcm @ v_ric
  • Note: This is the transpose of dcm_eci2ric

• compute_angular_vel_rad_p_sec(input_state)

  • Compute angular velocity vector: omega = cross(r, v) / ||r||^2
  • Parameters:
    • input_state: State vector [..., 6] (position + velocity)
  • Returns: Angular velocity vector [..., 3] in radians per second

• eci2ric_vec(pos_eci, ref_state_eci)

  • Rotate ECI position vector(s) into RIC frame
  • Note: These are raw ECI points being rotated, not points relative to the reference state. For relative state conversion, use eci2ric_state.
  • Parameters:
    • pos_eci: Position vector(s) in ECI [..., 3]
    • ref_state_eci: Reference state in ECI [..., 6] (position + velocity) defining the RIC frame
  • Returns: Position vector(s) in RIC [..., 3]
  • Supports broadcasting: if one input is scalar and the other is array, the scalar is broadcasted

• ric2eci_vec(pos_ric, ref_state_eci)

  • Rotate RIC position vector(s) into ECI frame
  • Note: These are raw RIC points being rotated, not points relative to the reference state. For relative state conversion, use ric2eci_state.
  • Parameters:
    • pos_ric: Position vector(s) in RIC [..., 3]
    • ref_state_eci: Reference state in ECI [..., 6] (position + velocity) defining the RIC frame
  • Returns: Position vector(s) in ECI [..., 3]
  • Supports broadcasting: if one input is scalar and the other is array, the scalar is broadcasted

• apply_ric_waypoint_offset(target_state_eci, ric_waypoint)

  • Apply RIC waypoint offset to target state(s)
  • Parameters:
    • target_state_eci: Target state(s) in ECI [..., 6] (position + velocity)
    • ric_waypoint: RIC waypoint offset [..., 3] or [..., 6] - [dr_r, dr_i, dr_c] or [dr_r, dr_i, dr_c, dv_r, dv_i, dv_c]
      • If 3 elements: only position offset (velocity offset is zero)
      • If 6 elements: full state offset (position + velocity)
  • Returns: Modified target state(s) in ECI [..., 6]
  • Supports broadcasting: if one input is scalar and the other is array, the scalar is broadcasted
  • If both are arrays, they must have matching batch sizes (unless one is scalar)

Classical Orbital Elements

• rv2coe(eci_state_rv, deg=True, mu=398600.4418, tol=1e-8, special=False)

  • Convert ECI state vector(s) to classical orbital elements
  • Parameters:
    • eci_state_rv: State vector(s) in ECI [..., 6] (position + velocity)
    • deg: If True, return angles in degrees; if False, return in radians (default: True)
    • mu: Gravitational parameter [km³/s²] (default: 398600.4418 for Earth)
    • tol: Tolerance for numerical stability (default: 1e-8)
    • special: If True, return special elements in addition to the six classical (default: False)
  • Returns: Array with shape [..., 6] or [..., 9] if special=True
    • Columns: [sma, ecc, inc, raan, aop, ta] (or with [arglat, truelon, lonper] if special=True)
    • sma: Semi-major axis [km]
    • ecc: Eccentricity
    • inc: Inclination [deg or rad]
    • raan: Right Ascension of the Ascending Node [deg or rad]
    • aop: Argument of Perigee [deg or rad]
    • ta: True Anomaly [deg or rad]
    • Special elements (if special=True):
      • arglat: Argument of Latitude [deg or rad]
      • truelon: True Longitude [deg or rad]
      • lonper: Longitude of Periapsis [deg or rad]
  • Supports vectorized operations: input can be single state [6] or batch [N, 6]
  • Handles special cases:
    • Circular orbits: AOP and TA set to zero by convention
    • Equatorial orbits: RAAN set to zero by convention
    • Circular and equatorial: uses true longitude instead of TA
    • Circular and non-equatorial: uses argument of latitude instead of TA
    • Elliptical and equatorial: uses longitude of perigee instead of AOP

• coe2rv(coe, deg=True, mu=398600.4418)

  • Convert classical orbital elements to ECI state vector(s)
  • Parameters:
    • coe: Classical orbital elements [..., 6] - [sma, ecc, inc, raan, aop, ta]
      • sma: Semi-major axis [km]
      • ecc: Eccentricity
      • inc: Inclination [deg or rad]
      • raan: Right Ascension of the Ascending Node [deg or rad]
      • aop: Argument of Perigee [deg or rad]
      • ta: True Anomaly [deg or rad]
    • deg: If True, input angles are in degrees; if False, input angles are in radians (default: True)
    • mu: Gravitational parameter [km³/s²] (default: 398600.4418 for Earth)
  • Returns: State vector(s) in ECI [..., 6] - [x, y, z, vx, vy, vz] (km, km/s)
  • Supports vectorized operations: input can be single COE set [6] or batch [N, 6]
  • Ref: Vallado, Fundamentals of Astrodynamics and Applications, 4th ed., Algorithm 10, p. 118

Optimization

• optimize_lambert_transfer_py(source_state, target_state, time_bounds, dt_step, mu=398600.4418, retrograde=False)

  • Grid search for optimal transfer
  • Parameters:
    • source_state: Initial state vector [6] (position + velocity)
    • target_state: Target state vector [6] (position + velocity)
    • time_bounds: Tuple (min_departure_time, max_arrival_time) [s]
    • dt_step: Time step for grid search [s]
    • mu: Gravitational parameter [km³/s²] (default: 398600.4418 for Earth)
    • retrograde: Whether to allow retrograde transfers (default: False)
  • Returns: LambertBurnArray object with fields (arrays, shape (1,) for single solution):
    • t0: Departure time [s] - shape (n,)
    • tf: Arrival time [s] - shape (n,)
    • tof: Time of flight [s] - shape (n,)
    • cost: Total delta-v [km/s] - shape (n,)
    • dv1: First burn impulse [km/s] - shape (n, 3)
    • dv2: Second burn impulse [km/s] - shape (n, 3)
    • sma: Semi-major axis [km] - shape (n,)
    • ecc: Eccentricity - shape (n,)
    • inc: Inclination [deg] - shape (n,)
    • raan: Right Ascension of Ascending Node [deg] - shape (n,)
    • argp: Argument of Perigee [deg] - shape (n,)
    • ta_depart: True anomaly at departure [deg] - shape (n,)
    • ta_arrive: True anomaly at arrival [deg] - shape (n,)
    • rp: Radius of perigee [km] - shape (n,)
    • ra: Radius of apogee [km] - shape (n,)

• optimize_lambert_transfer_grid_py(source_state, target_state, time_bounds, dt_step, mu=398600.4418, retrograde=False)

  • Get all grid search solutions
  • Parameters: Same as optimize_lambert_transfer_py
  • Returns: LambertBurnArray object with fields:
    • t0: Departure times [s] - shape (n,)
    • tf: Arrival times [s] - shape (n,)
    • tof: Time of flight [s] - shape (n,)
    • cost: Total delta-v [km/s] - shape (n,)
    • dv1: First burn impulses [km/s] - shape (n, 3)
    • dv2: Second burn impulses [km/s] - shape (n, 3)
    • sma: Semi-major axis [km] - shape (n,)
    • ecc: Eccentricity - shape (n,)
    • inc: Inclination [deg] - shape (n,)
    • raan: Right Ascension of Ascending Node [deg] - shape (n,)
    • argp: Argument of Perigee [deg] - shape (n,)
    • ta_depart: True anomaly at departure [deg] - shape (n,)
    • ta_arrive: True anomaly at arrival [deg] - shape (n,)
    • rp: Radius of perigee [km] - shape (n,)
    • ra: Radius of apogee [km] - shape (n,)

• lambert_grid_search_py(source_state, target_state, time_bounds, dt_step, max_rev=0, mu=398600.4418, tol=1e-8, maxiter=50, rendezvous=true)

  • Grid search using best_dv logic with automatic prograde/retrograde selection
  • Parameters:
    • source_state: Initial state vector [6] (position + velocity)
    • target_state: Target state vector [6] (position + velocity)
    • time_bounds: Tuple (min_departure_time, max_arrival_time) [s]
    • dt_step: Time step for grid search [s]
    • max_rev: Maximum number of revolutions to consider (default: 0)
    • mu: Gravitational parameter [km³/s²] (default: 398600.4418 for Earth)
    • tol: Tolerance for iterative solver (default: 1e-8)
    • maxiter: Maximum iterations for iterative solver (default: 50)
    • rendezvous: If True, optimize for full rendezvous (both burns, cost = dv1 + dv2); if False, optimize for transfer only (first burn, cost = dv1, dv2 = NaN) (default: True)
  • Returns: LambertBurnArray object with fields:
    • t0: Departure times [s] - shape (n,)
    • tf: Arrival times [s] - shape (n,)
    • tof: Time of flight [s] - shape (n,)
    • cost: Total delta-v [km/s] - shape (n,)
      • If rendezvous=True: cost = dv1 + dv2
      • If rendezvous=False: cost = dv1 only
    • dv1: First burn impulses [km/s] - shape (n, 3)
    • dv2: Second burn impulses [km/s] - shape (n, 3)
      • If rendezvous=True: computed second burn
      • If rendezvous=False: set to NaN
    • sma, ecc, inc, raan, argp: Keplerian elements [km, deg] - shape (n,)
    • ta_depart, ta_arrive: True anomaly at departure/arrival [deg] - shape (n,)
    • rp, ra: Radius of perigee/apogee [km] - shape (n,)
  • Note: Automatically selects the best solution (prograde or retrograde) for each grid point, similar to lambert_izzo_best_dv

• optimize_lambert_nm_py(source_state, target_state, time_bounds, max_rev, optimize_rendezvous, mu=398600.4418, max_iters=None, initial_guess=None)

  • Nelder-Mead optimization
  • Automatically selects best solution from both prograde and retrograde transfers
  • Parameters:
    • source_state: Initial state vector [6] (position + velocity)
    • target_state: Target state vector [6] (position + velocity)
    • time_bounds: Tuple (min_departure_time, max_arrival_time) [s]
    • max_rev: Maximum number of revolutions to consider
    • optimize_rendezvous: If True, optimize for full rendezvous (both burns); if False, optimize for transfer only (first burn)
    • mu: Gravitational parameter [km³/s²] (default: 398600.4418 for Earth)
    • max_iters: Maximum iterations (default: None, uses 1000 internally)
    • initial_guess: Optional tuple (t0, tf) for initial guess [s]
  • Returns: LambertBurnArray object (same fields as optimize_lambert_transfer_py)

• optimize_lambert_nm_multi_py(source_state, target_state, time_bounds, max_rev, optimize_rendezvous, num_solvers, return_best_only, remove_duplicates, mu=398600.4418, max_iters=None, ric_waypoint=None)

  • Multi-start parallel optimization
  • Automatically selects best solution from both prograde and retrograde transfers for each optimizer
  • Parameters:
    • source_state: Initial state vector [6] (position + velocity)
    • target_state: Target state vector [6] (position + velocity)
    • time_bounds: Tuple (min_departure_time, max_arrival_time) [s]
    • max_rev: Maximum number of revolutions to consider
    • optimize_rendezvous: If True, optimize for full rendezvous (both burns); if False, optimize for transfer only (first burn)
    • num_solvers: Number of parallel optimizers to spawn
    • return_best_only: If True, return only the best solution; if False, return all solutions
    • remove_duplicates: If True, remove duplicate solutions that converged to the same point
    • mu: Gravitational parameter [km³/s²] (default: 398600.4418 for Earth)
    • max_iters: Maximum iterations (default: None, uses 1000 internally)
    • ric_waypoint: Optional RIC waypoint offset [3] or [6] - [dr_r, dr_i, dr_c] or [dr_r, dr_i, dr_c, dv_r, dv_i, dv_c] (default: None)
  • Returns: LambertBurnArray object (same fields as optimize_lambert_transfer_grid_py)

• optimize_lambert_nm_multi_sunline_py(source_state, target_state, time_bounds, max_rev, optimize_rendezvous, num_solvers, return_best_only, sun_vector_offset_km, t0_ref, mu=398600.4418, remove_duplicates=true, max_iters=None)

  • Multi-start parallel optimization targeting a position offset from target_state along the target-to-sun vector
  • Automatically selects best solution from both prograde and retrograde transfers for each optimizer
  • Parameters:
    • source_state: Initial state vector [6] (position + velocity)
    • target_state: Target state vector [6] (position + velocity)
    • time_bounds: Tuple (min_departure_time, max_arrival_time) [s] (relative to t0_ref)
    • max_rev: Maximum number of revolutions to consider
    • optimize_rendezvous: If True, optimize for full rendezvous (both burns); if False, optimize for transfer only (first burn)
    • num_solvers: Number of parallel optimizers to spawn
    • return_best_only: If True, return only the best solution; if False, return all solutions
    • sun_vector_offset_km: Offset distance along target-to-sun vector [km]
    • t0_ref: Reference epoch in UTC POSIX seconds (absolute start time of time_bounds[0])
    • mu: Gravitational parameter [km³/s²] (default: 398600.4418 for Earth)
    • remove_duplicates: If True, remove duplicate solutions that converged to the same point (default: True)
    • max_iters: Maximum iterations (default: None, uses 1000 internally)
  • Returns: LambertBurnArray object (same fields as optimize_lambert_transfer_grid_py)
  • Note: The function computes the sun's position at the arrival time using the compute_sun_state function, then offsets the target position along the target-to-sun vector by the specified distance.

• optimize_lambert_nm_multi_array(source_state, target_state, time_bounds, max_rev, optimize_rendezvous, num_solvers, return_best_only, remove_duplicates, mu=398600.4418, max_iters=None, ric_waypoint=None)

  • Array-based multi-start optimization for multiple source and target states
  • Parameters:
    • source_state: Array of shape [N, 6] where N is any number of source states
    • target_state: Array of shape [M, 6] where M is any number of target states
    • time_bounds: Tuple (min_departure_time, max_arrival_time) [s]
    • max_rev: Maximum number of revolutions to consider
    • optimize_rendezvous: If True, optimize for full rendezvous (both burns); if False, optimize for transfer only (first burn)
    • num_solvers: Number of parallel optimizers to spawn per combination
    • return_best_only: If True, return only the best solution per [i,j] combination; if False, return all solutions
    • remove_duplicates: If True, remove duplicate solutions that converged to the same point
    • mu: Gravitational parameter [km³/s²] (default: 398600.4418 for Earth)
    • max_iters: Maximum iterations (default: None, uses 1000 internally)
    • ric_waypoint: Optional RIC waypoint offset(s) - single [3] or [6] (applied to all targets), or array [M, 3] or [M, 6] (one per target state) (default: None)
  • Computes all N×M combinations in parallel
  • Automatically selects best solution from both prograde and retrograde transfers for each optimizer
  • Returns:
    • pandas DataFrame (if pandas is available) with columns:
      • source_idx: Index of source state (0 to N-1)
      • target_idx: Index of target state (0 to M-1)
      • t0: Departure time [s]
      • tf: Arrival time [s]
      • cost: Total delta-v [km/s]
      • dv1_x, dv1_y, dv1_z: First burn impulse components [km/s]
      • dv2_x, dv2_y, dv2_z: Second burn impulse components [km/s]
    • Dictionary (if pandas is not available) with the same keys as DataFrame columns, each containing a list of values
  • When return_best_only=True: One row/entry per [i,j] source-target combination
  • When return_best_only=False: Multiple rows/entries per [i,j] combination (all solutions)

Performance

• Parallel Processing: Automatically uses N-1 CPU cores for batch operations
• Vectorized Operations: Efficient NumPy array operations
• Rust Performance: Core algorithms implemented in Rust for maximum speed

Units

• Positions: kilometers (km)
• Velocities: kilometers per second (km/s)
• Time: seconds (s)
• Gravitational Parameter (mu): km³/s²
  • Earth: 398600.4418 km³/s²