Lambert Targeting¶

This tutorial demonstrates how to use satkit's Lambert solver to design orbit transfers. We'll compute the delta-v required for a satellite to move between orbits and visualize the transfer geometry.

In [1]:

import matplotlib.pyplot as plt
import scienceplots  # noqa: F401
plt.style.use(["science", "no-latex", "../satkit.mplstyle"])
%config InlineBackend.figure_formats = ['svg']

import satkit as sk
import numpy as np

import matplotlib.pyplot as plt import scienceplots # noqa: F401 plt.style.use(["science", "no-latex", "../satkit.mplstyle"]) %config InlineBackend.figure_formats = ['svg'] import satkit as sk import numpy as np

Scenario: LEO to MEO Transfer¶

A satellite is in a 400 km circular orbit (LEO) and needs to transfer to a 2000 km circular orbit (MEO). We'll compute the optimal transfer for a range of transfer angles and find the one that minimizes total delta-v.

In [2]:

# Orbit parameters
alt_1 = 400e3    # departure altitude (m)
alt_2 = 2000e3   # arrival altitude (m)
r_earth = sk.consts.earth_radius
mu = sk.consts.mu_earth

r1_mag = r_earth + alt_1
r2_mag = r_earth + alt_2

# Circular velocities at each orbit
v_circ_1 = np.sqrt(mu / r1_mag)
v_circ_2 = np.sqrt(mu / r2_mag)

print(f"Departure orbit: {alt_1/1e3:.0f} km altitude, v = {v_circ_1:.1f} m/s")
print(f"Arrival orbit:   {alt_2/1e3:.0f} km altitude, v = {v_circ_2:.1f} m/s")

# Orbit parameters alt_1 = 400e3 # departure altitude (m) alt_2 = 2000e3 # arrival altitude (m) r_earth = sk.consts.earth_radius mu = sk.consts.mu_earth r1_mag = r_earth + alt_1 r2_mag = r_earth + alt_2 # Circular velocities at each orbit v_circ_1 = np.sqrt(mu / r1_mag) v_circ_2 = np.sqrt(mu / r2_mag) print(f"Departure orbit: {alt_1/1e3:.0f} km altitude, v = {v_circ_1:.1f} m/s") print(f"Arrival orbit: {alt_2/1e3:.0f} km altitude, v = {v_circ_2:.1f} m/s")

Departure orbit: 400 km altitude, v = 7668.6 m/s
Arrival orbit:   2000 km altitude, v = 6897.6 m/s

Delta-v vs Transfer Angle¶

We place the satellite at a fixed departure point and sweep the arrival point around the higher orbit, computing the Lambert solution and delta-v for each transfer angle.

In [3]:

# Departure position (on x-axis)
r1 = np.array([r1_mag, 0.0, 0.0])

# Sweep transfer angles from 30 to 170 degrees
angles_deg = np.linspace(30, 180, 200)
angles_rad = np.radians(angles_deg)

# Hohmann transfer time (reference)
a_transfer = (r1_mag + r2_mag) / 2
t_hohmann = np.pi * np.sqrt(a_transfer**3 / mu)

dv1_arr = []
dv2_arr = []
dv_total_arr = []

for theta in angles_rad:
    r2 = np.array([r2_mag * np.cos(theta), r2_mag * np.sin(theta), 0.0])

    # Scale TOF with transfer angle: fraction of Hohmann time
    tof = t_hohmann * (theta / np.pi)

    solutions = sk.lambert(r1, r2, tof)
    v1_transfer, v2_transfer = solutions[0]

    # Departure: satellite in circular orbit, moving in +y direction
    v1_circular = np.array([0.0, v_circ_1, 0.0])

    # Arrival: circular velocity tangent to orbit at arrival point
    v2_circular = np.array([
        -v_circ_2 * np.sin(theta),
        v_circ_2 * np.cos(theta),
        0.0,
    ])

    dv1 = np.linalg.norm(v1_transfer - v1_circular)
    dv2 = np.linalg.norm(v2_transfer - v2_circular)

    dv1_arr.append(dv1)
    dv2_arr.append(dv2)
    dv_total_arr.append(dv1 + dv2)

dv1_arr = np.array(dv1_arr)
dv2_arr = np.array(dv2_arr)
dv_total_arr = np.array(dv_total_arr)

# Find minimum total delta-v
idx_min = np.argmin(dv_total_arr)
print(f"Minimum total delta-v: {dv_total_arr[idx_min]:.1f} m/s at {angles_deg[idx_min]:.1f} degrees")
print(f"  Departure burn: {dv1_arr[idx_min]:.1f} m/s")
print(f"  Arrival burn:   {dv2_arr[idx_min]:.1f} m/s")

# Departure position (on x-axis) r1 = np.array([r1_mag, 0.0, 0.0]) # Sweep transfer angles from 30 to 170 degrees angles_deg = np.linspace(30, 180, 200) angles_rad = np.radians(angles_deg) # Hohmann transfer time (reference) a_transfer = (r1_mag + r2_mag) / 2 t_hohmann = np.pi * np.sqrt(a_transfer**3 / mu) dv1_arr = [] dv2_arr = [] dv_total_arr = [] for theta in angles_rad: r2 = np.array([r2_mag * np.cos(theta), r2_mag * np.sin(theta), 0.0]) # Scale TOF with transfer angle: fraction of Hohmann time tof = t_hohmann * (theta / np.pi) solutions = sk.lambert(r1, r2, tof) v1_transfer, v2_transfer = solutions[0] # Departure: satellite in circular orbit, moving in +y direction v1_circular = np.array([0.0, v_circ_1, 0.0]) # Arrival: circular velocity tangent to orbit at arrival point v2_circular = np.array([ -v_circ_2 * np.sin(theta), v_circ_2 * np.cos(theta), 0.0, ]) dv1 = np.linalg.norm(v1_transfer - v1_circular) dv2 = np.linalg.norm(v2_transfer - v2_circular) dv1_arr.append(dv1) dv2_arr.append(dv2) dv_total_arr.append(dv1 + dv2) dv1_arr = np.array(dv1_arr) dv2_arr = np.array(dv2_arr) dv_total_arr = np.array(dv_total_arr) # Find minimum total delta-v idx_min = np.argmin(dv_total_arr) print(f"Minimum total delta-v: {dv_total_arr[idx_min]:.1f} m/s at {angles_deg[idx_min]:.1f} degrees") print(f" Departure burn: {dv1_arr[idx_min]:.1f} m/s") print(f" Arrival burn: {dv2_arr[idx_min]:.1f} m/s")

Minimum total delta-v: 768.8 m/s at 180.0 degrees
  Departure burn: 394.6 m/s
  Arrival burn:   374.2 m/s

In [4]:

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(angles_deg, dv1_arr, label=r"$\Delta v_1$ (departure)")
ax.plot(angles_deg, dv2_arr, label=r"$\Delta v_2$ (arrival)")
ax.plot(angles_deg, dv_total_arr, label=r"$\Delta v_{total}$", linewidth=2)
ax.axvline(angles_deg[idx_min], color="#BBBBBB", linestyle=":", linewidth=1)
ax.set_xlabel("Transfer Angle (degrees)")
ax.set_ylabel(r"$\Delta v$ (m/s)")
ax.set_title("LEO to MEO Transfer")
ax.legend()
plt.tight_layout()
plt.show()

fig, ax = plt.subplots(figsize=(8, 5)) ax.plot(angles_deg, dv1_arr, label=r"$\Delta v_1$ (departure)") ax.plot(angles_deg, dv2_arr, label=r"$\Delta v_2$ (arrival)") ax.plot(angles_deg, dv_total_arr, label=r"$\Delta v_{total}$", linewidth=2) ax.axvline(angles_deg[idx_min], color="#BBBBBB", linestyle=":", linewidth=1) ax.set_xlabel("Transfer Angle (degrees)") ax.set_ylabel(r"$\Delta v$ (m/s)") ax.set_title("LEO to MEO Transfer") ax.legend() plt.tight_layout() plt.show()

Pork-Chop Plot: Delta-v vs Transfer Angle and Time¶

A pork-chop plot shows total delta-v as a function of both transfer angle and time of flight, revealing the trade-off between fuel cost and transfer duration.

In [5]:

# Sweep both transfer angle and time of flight
angles = np.linspace(60, 300, 120)
tofs = np.linspace(0.4 * t_hohmann, 3.0 * t_hohmann, 100)

DV = np.full((len(tofs), len(angles)), np.nan)

for i, t in enumerate(tofs):
    for j, theta_deg in enumerate(angles):
        theta = np.radians(theta_deg)
        r2 = np.array([r2_mag * np.cos(theta), r2_mag * np.sin(theta), 0.0])

        try:
            solutions = sk.lambert(r1, r2, t)
            v1_t, v2_t = solutions[0]

            v1_c = np.array([0.0, v_circ_1, 0.0])
            v2_c = np.array([
                -v_circ_2 * np.sin(theta),
                v_circ_2 * np.cos(theta),
                0.0,
            ])

            dv = np.linalg.norm(v1_t - v1_c) + np.linalg.norm(v2_t - v2_c)
            if dv < 15000:
                DV[i, j] = dv
        except Exception:
            pass

fig, ax = plt.subplots(figsize=(9, 6))
levels = np.linspace(500, 8000, 20)
cs = ax.contourf(angles, tofs / 60, DV, levels=levels, cmap="viridis", extend="max")
ax.contour(angles, tofs / 60, DV, levels=levels, colors="white", linewidths=0.3)
cbar = fig.colorbar(cs, ax=ax, label=r"Total $\Delta v$ (m/s)")
ax.set_xlabel("Transfer Angle (degrees)")
ax.set_ylabel("Time of Flight (minutes)")
ax.set_title("LEO to MEO Transfer: Pork-Chop Plot")
plt.tight_layout()
plt.show()

# Sweep both transfer angle and time of flight angles = np.linspace(60, 300, 120) tofs = np.linspace(0.4 * t_hohmann, 3.0 * t_hohmann, 100) DV = np.full((len(tofs), len(angles)), np.nan) for i, t in enumerate(tofs): for j, theta_deg in enumerate(angles): theta = np.radians(theta_deg) r2 = np.array([r2_mag * np.cos(theta), r2_mag * np.sin(theta), 0.0]) try: solutions = sk.lambert(r1, r2, t) v1_t, v2_t = solutions[0] v1_c = np.array([0.0, v_circ_1, 0.0]) v2_c = np.array([ -v_circ_2 * np.sin(theta), v_circ_2 * np.cos(theta), 0.0, ]) dv = np.linalg.norm(v1_t - v1_c) + np.linalg.norm(v2_t - v2_c) if dv < 15000: DV[i, j] = dv except Exception: pass fig, ax = plt.subplots(figsize=(9, 6)) levels = np.linspace(500, 8000, 20) cs = ax.contourf(angles, tofs / 60, DV, levels=levels, cmap="viridis", extend="max") ax.contour(angles, tofs / 60, DV, levels=levels, colors="white", linewidths=0.3) cbar = fig.colorbar(cs, ax=ax, label=r"Total $\Delta v$ (m/s)") ax.set_xlabel("Transfer Angle (degrees)") ax.set_ylabel("Time of Flight (minutes)") ax.set_title("LEO to MEO Transfer: Pork-Chop Plot") plt.tight_layout() plt.show()

Transfer Orbit Visualization¶

Let's visualize the optimal transfer orbit alongside the departure and arrival orbits.

In [6]:

# Optimal transfer
theta_opt = np.radians(angles_deg[idx_min])
r2_opt = np.array([r2_mag * np.cos(theta_opt), r2_mag * np.sin(theta_opt), 0.0])

tof_opt = t_hohmann * (theta_opt / np.pi)
solutions = sk.lambert(r1, r2_opt, tof_opt)
v1_transfer, v2_transfer = solutions[0]

# Propagate transfer orbit via two-body integration
from scipy.integrate import solve_ivp

def twobody(t, state):
    r = state[:3]
    v = state[3:]
    r_norm = np.linalg.norm(r)
    a = -mu / r_norm**3 * r
    return np.concatenate([v, a])

state0 = np.concatenate([r1, v1_transfer])
sol = solve_ivp(twobody, [0, tof_opt], state0, rtol=1e-10, atol=1e-12,
                dense_output=True)
t_plot = np.linspace(0, tof_opt, 500)
positions = sol.sol(t_plot)[:3].T

# Draw orbits
fig, ax = plt.subplots(figsize=(8, 8))
theta_circ = np.linspace(0, 2 * np.pi, 500)

# Earth
earth = plt.Circle((0, 0), r_earth / 1e3, color="lightblue", ec="steelblue", linewidth=1)
ax.add_patch(earth)

# Departure orbit
ax.plot(r1_mag * np.cos(theta_circ) / 1e3, r1_mag * np.sin(theta_circ) / 1e3,
        "--", color="#BBBBBB", linewidth=1, label=f"LEO ({alt_1/1e3:.0f} km)")

# Arrival orbit
ax.plot(r2_mag * np.cos(theta_circ) / 1e3, r2_mag * np.sin(theta_circ) / 1e3,
        "--", color="#BBBBBB", linewidth=1, label=f"MEO ({alt_2/1e3:.0f} km)")

# Transfer arc
ax.plot(positions[:, 0] / 1e3, positions[:, 1] / 1e3,
        linewidth=2, color="#CC3311", label="Transfer orbit")

# Endpoints
ax.plot(r1[0] / 1e3, r1[1] / 1e3, "o", color="#0077BB", markersize=8, zorder=5)
ax.annotate("Departure", xy=(r1[0] / 1e3, r1[1] / 1e3),
            xytext=(15, 15), textcoords="offset points", fontsize=11)
ax.plot(r2_opt[0] / 1e3, r2_opt[1] / 1e3, "s", color="#009988", markersize=8, zorder=5)
ax.annotate("Arrival", xy=(r2_opt[0] / 1e3, r2_opt[1] / 1e3),
            xytext=(15, -20), textcoords="offset points", fontsize=11)

ax.set_xlabel("x (km)")
ax.set_ylabel("y (km)")
ax.set_title(f"Transfer Orbit  ({angles_deg[idx_min]:.0f}°, " + r"$\Delta v$ = " + f"{dv_total_arr[idx_min]:.0f} m/s)")
ax.set_aspect("equal")
ax.legend(loc="lower left")
plt.tight_layout()
plt.show()

# Optimal transfer theta_opt = np.radians(angles_deg[idx_min]) r2_opt = np.array([r2_mag * np.cos(theta_opt), r2_mag * np.sin(theta_opt), 0.0]) tof_opt = t_hohmann * (theta_opt / np.pi) solutions = sk.lambert(r1, r2_opt, tof_opt) v1_transfer, v2_transfer = solutions[0] # Propagate transfer orbit via two-body integration from scipy.integrate import solve_ivp def twobody(t, state): r = state[:3] v = state[3:] r_norm = np.linalg.norm(r) a = -mu / r_norm**3 * r return np.concatenate([v, a]) state0 = np.concatenate([r1, v1_transfer]) sol = solve_ivp(twobody, [0, tof_opt], state0, rtol=1e-10, atol=1e-12, dense_output=True) t_plot = np.linspace(0, tof_opt, 500) positions = sol.sol(t_plot)[:3].T # Draw orbits fig, ax = plt.subplots(figsize=(8, 8)) theta_circ = np.linspace(0, 2 * np.pi, 500) # Earth earth = plt.Circle((0, 0), r_earth / 1e3, color="lightblue", ec="steelblue", linewidth=1) ax.add_patch(earth) # Departure orbit ax.plot(r1_mag * np.cos(theta_circ) / 1e3, r1_mag * np.sin(theta_circ) / 1e3, "--", color="#BBBBBB", linewidth=1, label=f"LEO ({alt_1/1e3:.0f} km)") # Arrival orbit ax.plot(r2_mag * np.cos(theta_circ) / 1e3, r2_mag * np.sin(theta_circ) / 1e3, "--", color="#BBBBBB", linewidth=1, label=f"MEO ({alt_2/1e3:.0f} km)") # Transfer arc ax.plot(positions[:, 0] / 1e3, positions[:, 1] / 1e3, linewidth=2, color="#CC3311", label="Transfer orbit") # Endpoints ax.plot(r1[0] / 1e3, r1[1] / 1e3, "o", color="#0077BB", markersize=8, zorder=5) ax.annotate("Departure", xy=(r1[0] / 1e3, r1[1] / 1e3), xytext=(15, 15), textcoords="offset points", fontsize=11) ax.plot(r2_opt[0] / 1e3, r2_opt[1] / 1e3, "s", color="#009988", markersize=8, zorder=5) ax.annotate("Arrival", xy=(r2_opt[0] / 1e3, r2_opt[1] / 1e3), xytext=(15, -20), textcoords="offset points", fontsize=11) ax.set_xlabel("x (km)") ax.set_ylabel("y (km)") ax.set_title(f"Transfer Orbit ({angles_deg[idx_min]:.0f}°, " + r"$\Delta v$ = " + f"{dv_total_arr[idx_min]:.0f} m/s)") ax.set_aspect("equal") ax.legend(loc="lower left") plt.tight_layout() plt.show()

Interplanetary Example: Earth to Mars¶

Lambert targeting also works for interplanetary transfers by using the Sun's gravitational parameter and heliocentric positions.

In [7]:

mu_sun = sk.consts.mu_sun
au = sk.consts.au

# Simplified circular coplanar orbits
r_earth_orbit = 1.0 * au
r_mars_orbit = 1.524 * au

# Earth at departure
r1_helio = np.array([r_earth_orbit, 0, 0])

# Mars at arrival (~135 degree transfer)
transfer_angle = np.radians(135)
r2_helio = np.array([
    r_mars_orbit * np.cos(transfer_angle),
    r_mars_orbit * np.sin(transfer_angle),
    0,
])

# Transfer time: ~8 months
tof_days = 245
tof_interp = tof_days * 86400

solutions = sk.lambert(r1_helio, r2_helio, tof_interp, mu=mu_sun)
v1_t, v2_t = solutions[0]

# Circular velocities
v_earth = np.sqrt(mu_sun / r_earth_orbit)
v_mars = np.sqrt(mu_sun / r_mars_orbit)

v1_earth = np.array([0, v_earth, 0])
v2_mars = np.array([
    -v_mars * np.sin(transfer_angle),
    v_mars * np.cos(transfer_angle),
    0,
])

dv_depart = np.linalg.norm(v1_t - v1_earth)
dv_arrive = np.linalg.norm(v2_t - v2_mars)

print(f"Earth-Mars transfer ({tof_days} days):")
print(f"  Departure delta-v:  {dv_depart/1e3:.2f} km/s")
print(f"  Arrival delta-v:    {dv_arrive/1e3:.2f} km/s")
print(f"  Total delta-v:      {(dv_depart + dv_arrive)/1e3:.2f} km/s")

mu_sun = sk.consts.mu_sun au = sk.consts.au # Simplified circular coplanar orbits r_earth_orbit = 1.0 * au r_mars_orbit = 1.524 * au # Earth at departure r1_helio = np.array([r_earth_orbit, 0, 0]) # Mars at arrival (~135 degree transfer) transfer_angle = np.radians(135) r2_helio = np.array([ r_mars_orbit * np.cos(transfer_angle), r_mars_orbit * np.sin(transfer_angle), 0, ]) # Transfer time: ~8 months tof_days = 245 tof_interp = tof_days * 86400 solutions = sk.lambert(r1_helio, r2_helio, tof_interp, mu=mu_sun) v1_t, v2_t = solutions[0] # Circular velocities v_earth = np.sqrt(mu_sun / r_earth_orbit) v_mars = np.sqrt(mu_sun / r_mars_orbit) v1_earth = np.array([0, v_earth, 0]) v2_mars = np.array([ -v_mars * np.sin(transfer_angle), v_mars * np.cos(transfer_angle), 0, ]) dv_depart = np.linalg.norm(v1_t - v1_earth) dv_arrive = np.linalg.norm(v2_t - v2_mars) print(f"Earth-Mars transfer ({tof_days} days):") print(f" Departure delta-v: {dv_depart/1e3:.2f} km/s") print(f" Arrival delta-v: {dv_arrive/1e3:.2f} km/s") print(f" Total delta-v: {(dv_depart + dv_arrive)/1e3:.2f} km/s")

Earth-Mars transfer (245 days):
  Departure delta-v:  7.89 km/s
  Arrival delta-v:    4.81 km/s
  Total delta-v:      12.70 km/s