lambert-rs 0.0.17
Rust + PyO3 implementation of the Izzo Lambert solver
License: MIT
Requires Python: >=3.8
Keywords: lambert, orbital-mechanics, astrodynamics, rust

# lambert_rs

High-performance Rust + PyO3 implementation of the Izzo Lambert solver with advanced optimization capabilities.

## Installation

### From PyPI

```bash
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:

```python
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:

```python
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:

```python
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:

```python
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):

```python
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:

```python
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:

```python
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:

```python
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):

```python
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:

```python
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:

```python
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`:

```python
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:

```python
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:

```python
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:

```python
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:

```python
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:

```python
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:

```python
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:

```python
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

```python
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

```python
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:

```python
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:

```python
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:

```python
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):

```python
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:

```python
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:

```python
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:

```python
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):

```python
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:

```python
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:

```python
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.

```python
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