Orbit Maneuvers¶

This tutorial demonstrates how to model orbit maneuvers in satkit, covering both impulsive (chemical) and continuous (electric) thrust. We'll work through three progressively complex examples:

1. Hohmann Transfer -- the classic two-burn orbit raise using impulsive maneuvers
2. Lambert-Targeted Transfer -- using the Lambert solver to compute transfer burns
3. Low-Thrust Orbit Raising -- continuous electric propulsion spiral-out

Maneuver Types in SatKit¶

Type	Class	Where	Use Case
Impulsive	satstate.add_maneuver()	Instantaneous dv at a specific time	Chemical propulsion, station-keeping
Continuous	satproperties with thrust arcs	Sustained acceleration over a time window	Electric propulsion, drag makeup

The RTN Coordinate Frame¶

Both maneuver types support the RTN (Radial / Tangential / Normal) coordinate frame, which is the CCSDS OEM/OMM standard for orbit data messages and the canonical name in satkit. The same frame is also known as RSW (Vallado's textbook) and RIC (older NASA / Clohessy-Wiltshire literature); sk.frame.RSW and sk.frame.RIC are Python-level aliases that resolve to sk.frame.RTN.

The RTN axes are defined relative to the satellite's current orbital state:

• R (Radial): Points away from the Earth center along the position vector
• T (Tangential / In-track): Perpendicular to the radial direction in the orbital plane, roughly aligned with the velocity vector. For circular orbits this is exactly the velocity direction; for eccentric orbits it differs from velocity by the flight-path angle (use sk.frame.NTW if you need "strictly along velocity")
• N (Normal / Cross-track): Completes the right-handed triad; aligned with the orbit angular momentum vector (orbit normal)

Thrust or delta-v vectors specified as [R, T, N] in sk.frame.RTN are automatically rotated to GCRF by the propagator using the satellite's instantaneous position and velocity. See the Theory: Maneuver Coordinate Frames guide for a side-by-side comparison of RTN, NTW, LVLH, and GCRF.

In [1]:

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']

mu = sk.consts.mu_earth
Re = sk.consts.earth_radius

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'] mu = sk.consts.mu_earth Re = sk.consts.earth_radius

1. Hohmann Transfer (Impulsive)¶

A Hohmann transfer is the most fuel-efficient two-impulse transfer between coplanar circular orbits. It consists of:

1. Burn 1 at perigee: accelerate in-track to enter an elliptical transfer orbit
2. Coast along the transfer ellipse for half an orbit
3. Burn 2 at apogee: accelerate in-track to circularize at the target altitude

The Δv for each burn is computed from the vis-viva equation:

$$v = \sqrt{\mu \left(\frac{2}{r} - \frac{1}{a}\right)}$$

In [2]:

# Orbit altitudes
alt_initial = 400e3   # 400 km LEO
alt_final   = 800e3   # 800 km

r1 = Re + alt_initial
r2 = Re + alt_final
a_transfer = (r1 + r2) / 2  # semi-major axis of transfer ellipse

# Circular velocities
v_circ_1 = np.sqrt(mu / r1)
v_circ_2 = np.sqrt(mu / r2)

# Transfer orbit velocities (vis-viva)
v_transfer_perigee = np.sqrt(mu * (2/r1 - 1/a_transfer))
v_transfer_apogee  = np.sqrt(mu * (2/r2 - 1/a_transfer))

# Delta-v for each burn (both in-track / prograde)
dv1 = v_transfer_perigee - v_circ_1
dv2 = v_circ_2 - v_transfer_apogee
transfer_time = np.pi * np.sqrt(a_transfer**3 / mu)  # half orbit

print(f"Initial orbit: {alt_initial/1e3:.0f} km, v = {v_circ_1:.1f} m/s")
print(f"Final orbit:   {alt_final/1e3:.0f} km, v = {v_circ_2:.1f} m/s")
print("")
print(f"Burn 1 (perigee):  Δv = {dv1:.2f} m/s in-track")
print(f"Burn 2 (apogee):   Δv = {dv2:.2f} m/s in-track")
print(f"Total Δv:          {dv1 + dv2:.2f} m/s")
print(f"Transfer time:     {transfer_time/60:.1f} min")

# Orbit altitudes alt_initial = 400e3 # 400 km LEO alt_final = 800e3 # 800 km r1 = Re + alt_initial r2 = Re + alt_final a_transfer = (r1 + r2) / 2 # semi-major axis of transfer ellipse # Circular velocities v_circ_1 = np.sqrt(mu / r1) v_circ_2 = np.sqrt(mu / r2) # Transfer orbit velocities (vis-viva) v_transfer_perigee = np.sqrt(mu * (2/r1 - 1/a_transfer)) v_transfer_apogee = np.sqrt(mu * (2/r2 - 1/a_transfer)) # Delta-v for each burn (both in-track / prograde) dv1 = v_transfer_perigee - v_circ_1 dv2 = v_circ_2 - v_transfer_apogee transfer_time = np.pi * np.sqrt(a_transfer**3 / mu) # half orbit print(f"Initial orbit: {alt_initial/1e3:.0f} km, v = {v_circ_1:.1f} m/s") print(f"Final orbit: {alt_final/1e3:.0f} km, v = {v_circ_2:.1f} m/s") print("") print(f"Burn 1 (perigee): Δv = {dv1:.2f} m/s in-track") print(f"Burn 2 (apogee): Δv = {dv2:.2f} m/s in-track") print(f"Total Δv: {dv1 + dv2:.2f} m/s") print(f"Transfer time: {transfer_time/60:.1f} min")

Initial orbit: 400 km, v = 7668.6 m/s
Final orbit:   800 km, v = 7451.8 m/s

Burn 1 (perigee):  Δv = 109.12 m/s in-track
Burn 2 (apogee):   Δv = 107.56 m/s in-track
Total Δv:          216.68 m/s
Transfer time:     48.3 min

Setting Up the Maneuvers¶

We use satstate.add_maneuver() to schedule the two burns. The Δv is specified in the RTN frame where the components are [radial, tangential, normal]. Both burns are purely tangential (prograde for a circular orbit).

In [3]:

# Set up initial state: circular orbit at 400 km, equatorial
t0 = sk.time(2025, 6, 15, 0, 0, 0)
pos0 = np.array([r1, 0, 0])
vel0 = np.array([0, v_circ_1, 0])

sat = sk.satstate(time=t0, pos=pos0, vel=vel0)

# Schedule the two Hohmann burns
t_burn1 = t0  # burn immediately
t_burn2 = t0 + sk.duration(seconds=transfer_time)  # at apogee

# RTN frame: [radial, tangential, normal]
sat.add_maneuver(t_burn1, [0, dv1, 0], frame=sk.frame.RTN)
sat.add_maneuver(t_burn2, [0, dv2, 0], frame=sk.frame.RTN)

print(f"Maneuvers scheduled: {sat.num_maneuvers}")
print(sat)

# Set up initial state: circular orbit at 400 km, equatorial t0 = sk.time(2025, 6, 15, 0, 0, 0) pos0 = np.array([r1, 0, 0]) vel0 = np.array([0, v_circ_1, 0]) sat = sk.satstate(time=t0, pos=pos0, vel=vel0) # Schedule the two Hohmann burns t_burn1 = t0 # burn immediately t_burn2 = t0 + sk.duration(seconds=transfer_time) # at apogee # RTN frame: [radial, tangential, normal] sat.add_maneuver(t_burn1, [0, dv1, 0], frame=sk.frame.RTN) sat.add_maneuver(t_burn2, [0, dv2, 0], frame=sk.frame.RTN) print(f"Maneuvers scheduled: {sat.num_maneuvers}") print(sat)

Maneuvers scheduled: 2
Satellite State
                       Time: 2025-06-15T00:00:00.000000Z
              GCRF Position: [+6778137,       +0,       +0] m,
              GCRF Velocity: [  +0.000, +7668.558,   +0.000] m/s
       Maneuvers: 2

Propagating Through the Transfer¶

The propagator automatically segments at each maneuver time — no manual bookkeeping required. We propagate past both burns and sample the trajectory.

In [4]:

# Propagate well past the second burn to show the circularized orbit
t_end = t_burn2 + sk.duration(minutes=120)

# Use simple gravity for clean Hohmann geometry
settings = sk.propsettings(gravity_degree=1)

# Sample trajectory
npts = 500
times = [t0 + sk.duration(seconds=float(s)) for s in np.linspace(0, (t_end - t0).seconds, npts)]
positions = np.zeros((npts, 3))
altitudes = np.zeros(npts)

for i, t in enumerate(times):
    state = sat.propagate(t, propsettings=settings)
    positions[i] = state.pos / 1e3  # km
    altitudes[i] = np.linalg.norm(state.pos) - Re

# Also propagate without maneuvers for reference
sat_ref = sk.satstate(time=t0, pos=pos0, vel=vel0)
pos_ref = np.zeros((npts, 3))
for i, t in enumerate(times):
    state = sat_ref.propagate(t, propsettings=settings)
    pos_ref[i] = state.pos / 1e3

# Final state check
final = sat.propagate(t_end, propsettings=settings)
r_final = np.linalg.norm(final.pos)
v_final = np.linalg.norm(final.vel)
print(f"Final altitude: {(r_final - Re)/1e3:.1f} km (target: {alt_final/1e3:.0f} km)")
print(f"Final velocity: {v_final:.1f} m/s (circular: {v_circ_2:.1f} m/s)")

# Propagate well past the second burn to show the circularized orbit t_end = t_burn2 + sk.duration(minutes=120) # Use simple gravity for clean Hohmann geometry settings = sk.propsettings(gravity_degree=1) # Sample trajectory npts = 500 times = [t0 + sk.duration(seconds=float(s)) for s in np.linspace(0, (t_end - t0).seconds, npts)] positions = np.zeros((npts, 3)) altitudes = np.zeros(npts) for i, t in enumerate(times): state = sat.propagate(t, propsettings=settings) positions[i] = state.pos / 1e3 # km altitudes[i] = np.linalg.norm(state.pos) - Re # Also propagate without maneuvers for reference sat_ref = sk.satstate(time=t0, pos=pos0, vel=vel0) pos_ref = np.zeros((npts, 3)) for i, t in enumerate(times): state = sat_ref.propagate(t, propsettings=settings) pos_ref[i] = state.pos / 1e3 # Final state check final = sat.propagate(t_end, propsettings=settings) r_final = np.linalg.norm(final.pos) v_final = np.linalg.norm(final.vel) print(f"Final altitude: {(r_final - Re)/1e3:.1f} km (target: {alt_final/1e3:.0f} km)") print(f"Final velocity: {v_final:.1f} m/s (circular: {v_circ_2:.1f} m/s)")

Final altitude: 800.0 km (target: 800 km)
Final velocity: 7451.8 m/s (circular: 7451.8 m/s)

In [5]:

# Plot the transfer
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Left: 2D orbit plot
ax = axes[0]
earth = plt.Circle((0, 0), Re/1e3, color='skyblue', alpha=0.3)
ax.add_patch(earth)
ax.plot(pos_ref[:, 0], pos_ref[:, 1], '--', alpha=0.4, label=f'Initial orbit ({alt_initial/1e3:.0f} km)')
ax.plot(positions[:, 0], positions[:, 1], '-', linewidth=1.5, label='Hohmann transfer')

# Mark burn locations with arrows showing dv direction
burn1_state = sat.propagate(t_burn1 + sk.duration(seconds=1), propsettings=settings)
burn2_state = sat.propagate(t_burn2 + sk.duration(seconds=1), propsettings=settings)

# Burn 1: prograde at perigee
b1p = burn1_state.pos / 1e3
b1v = burn1_state.vel / np.linalg.norm(burn1_state.vel)  # unit velocity direction
ax.plot(b1p[0], b1p[1], 'o', color='green', markersize=10, zorder=5)
ax.annotate(f'Burn 1: \u0394v={dv1:.0f} m/s', xy=(b1p[0], b1p[1]),
            xytext=(b1p[0]+500, b1p[1]+800), fontsize=10,
            arrowprops=dict(arrowstyle='->', color='green', lw=1.5))

# Burn 2: prograde at apogee
b2p = burn2_state.pos / 1e3
ax.plot(b2p[0], b2p[1], 'o', color='red', markersize=10, zorder=5)
ax.annotate(f'Burn 2: \u0394v={dv2:.0f} m/s', xy=(b2p[0], b2p[1]),
            xytext=(b2p[0]-2500, b2p[1]+800), fontsize=10,
            arrowprops=dict(arrowstyle='->', color='red', lw=1.5))

# Draw target orbit circle
theta = np.linspace(0, 2*np.pi, 200)
ax.plot(r2/1e3*np.cos(theta), r2/1e3*np.sin(theta), ':', alpha=0.4, color='gray', label=f'Target orbit ({alt_final/1e3:.0f} km)')

ax.set_xlabel('X (km)')
ax.set_ylabel('Y (km)')
ax.set_title('Hohmann Transfer')
ax.set_aspect('equal')
ax.legend(loc='lower left', fontsize=10)

# Right: altitude vs time
ax = axes[1]
dt_min = np.array([(t - t0).seconds / 60 for t in times])
ax.plot(dt_min, altitudes / 1e3, linewidth=1.5)
ax.axhline(alt_initial/1e3, ls='--', alpha=0.4, label=f'Initial: {alt_initial/1e3:.0f} km')
ax.axhline(alt_final/1e3, ls='--', alpha=0.4, label=f'Target: {alt_final/1e3:.0f} km')
ax.axvline(transfer_time/60, ls=':', alpha=0.3, color='red', label='Burn 2')
ax.set_xlabel('Time (minutes)')
ax.set_ylabel('Altitude (km)')
ax.set_title('Altitude During Transfer')
ax.legend()

plt.tight_layout()
plt.show()

# Plot the transfer fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # Left: 2D orbit plot ax = axes[0] earth = plt.Circle((0, 0), Re/1e3, color='skyblue', alpha=0.3) ax.add_patch(earth) ax.plot(pos_ref[:, 0], pos_ref[:, 1], '--', alpha=0.4, label=f'Initial orbit ({alt_initial/1e3:.0f} km)') ax.plot(positions[:, 0], positions[:, 1], '-', linewidth=1.5, label='Hohmann transfer') # Mark burn locations with arrows showing dv direction burn1_state = sat.propagate(t_burn1 + sk.duration(seconds=1), propsettings=settings) burn2_state = sat.propagate(t_burn2 + sk.duration(seconds=1), propsettings=settings) # Burn 1: prograde at perigee b1p = burn1_state.pos / 1e3 b1v = burn1_state.vel / np.linalg.norm(burn1_state.vel) # unit velocity direction ax.plot(b1p[0], b1p[1], 'o', color='green', markersize=10, zorder=5) ax.annotate(f'Burn 1: \u0394v={dv1:.0f} m/s', xy=(b1p[0], b1p[1]), xytext=(b1p[0]+500, b1p[1]+800), fontsize=10, arrowprops=dict(arrowstyle='->', color='green', lw=1.5)) # Burn 2: prograde at apogee b2p = burn2_state.pos / 1e3 ax.plot(b2p[0], b2p[1], 'o', color='red', markersize=10, zorder=5) ax.annotate(f'Burn 2: \u0394v={dv2:.0f} m/s', xy=(b2p[0], b2p[1]), xytext=(b2p[0]-2500, b2p[1]+800), fontsize=10, arrowprops=dict(arrowstyle='->', color='red', lw=1.5)) # Draw target orbit circle theta = np.linspace(0, 2*np.pi, 200) ax.plot(r2/1e3*np.cos(theta), r2/1e3*np.sin(theta), ':', alpha=0.4, color='gray', label=f'Target orbit ({alt_final/1e3:.0f} km)') ax.set_xlabel('X (km)') ax.set_ylabel('Y (km)') ax.set_title('Hohmann Transfer') ax.set_aspect('equal') ax.legend(loc='lower left', fontsize=10) # Right: altitude vs time ax = axes[1] dt_min = np.array([(t - t0).seconds / 60 for t in times]) ax.plot(dt_min, altitudes / 1e3, linewidth=1.5) ax.axhline(alt_initial/1e3, ls='--', alpha=0.4, label=f'Initial: {alt_initial/1e3:.0f} km') ax.axhline(alt_final/1e3, ls='--', alpha=0.4, label=f'Target: {alt_final/1e3:.0f} km') ax.axvline(transfer_time/60, ls=':', alpha=0.3, color='red', label='Burn 2') ax.set_xlabel('Time (minutes)') ax.set_ylabel('Altitude (km)') ax.set_title('Altitude During Transfer') ax.legend() plt.tight_layout() plt.show()

2. Lambert-Targeted Transfer (Impulsive)¶

The Hohmann transfer works for coplanar circular orbits, but real mission design often requires transferring between arbitrary positions. Lambert's problem solves this: given two position vectors and a time of flight (TOF), find the connecting orbit.

Here we use Lambert to design a LEO-to-GEO transfer — the classic geostationary transfer orbit. We depart from a 400 km parking orbit and arrive at geostationary altitude (35,786 km), then compare how the total Δv varies with time of flight. The Hohmann transfer emerges as the minimum-Δv solution.

In [6]:

# LEO to GEO transfer via Lambert
r_leo = Re + 400e3
r_geo = 42164e3  # GEO radius
v_leo = np.sqrt(mu / r_leo)
v_geo = np.sqrt(mu / r_geo)

# Departure: 400 km equatorial
r_depart = np.array([r_leo, 0, 0])
v_depart = np.array([0, v_leo, 0])

# Arrival: GEO, 180 deg ahead (opposite side of Earth)
r_arrive = np.array([-r_geo, 0, 0])
v_arrive_circ = np.array([0, -v_geo, 0])

# Hohmann reference for comparison
a_hohmann = (r_leo + r_geo) / 2
hohmann_tof = np.pi * np.sqrt(a_hohmann**3 / mu)
dv1_hohmann = np.sqrt(mu * (2/r_leo - 1/a_hohmann)) - v_leo
dv2_hohmann = v_geo - np.sqrt(mu * (2/r_geo - 1/a_hohmann))

# Solve Lambert for a range of TOFs
tof_hours = np.linspace(3.5, 10, 50)
dv_total = np.zeros_like(tof_hours)

for i, hours in enumerate(tof_hours):
    tof = hours * 3600
    solutions = sk.lambert(r_depart, r_arrive, tof)
    v1, v2 = solutions[0]
    dv1 = np.linalg.norm(v1 - v_depart)
    dv2 = np.linalg.norm(v_arrive_circ - v2)
    dv_total[i] = dv1 + dv2

# Use optimal TOF (near Hohmann)
optimal_idx = np.argmin(dv_total)
optimal_tof = tof_hours[optimal_idx] * 3600
solutions = sk.lambert(r_depart, r_arrive, optimal_tof)
v1_transfer, v2_transfer = solutions[0]
dv_depart = v1_transfer - v_depart
dv_arrive = v_arrive_circ - v2_transfer

print(f"Hohmann reference: \u0394v = {dv1_hohmann + dv2_hohmann:.0f} m/s, TOF = {hohmann_tof/3600:.2f} hours")
print(f"Lambert optimal:   \u0394v = {dv_total[optimal_idx]:.0f} m/s, TOF = {tof_hours[optimal_idx]:.2f} hours")
print(f"  Departure \u0394v: {np.linalg.norm(dv_depart):.0f} m/s")
print(f"  Arrival \u0394v:   {np.linalg.norm(dv_arrive):.0f} m/s")

# Plot dv vs TOF
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(tof_hours, dv_total / 1e3, linewidth=1.5)
ax.axvline(hohmann_tof/3600, ls='--', color='red', alpha=0.5, label=f'Hohmann TOF ({hohmann_tof/3600:.1f} h)')
ax.axhline((dv1_hohmann + dv2_hohmann)/1e3, ls=':', color='gray', alpha=0.4, label=f'Hohmann \u0394v ({(dv1_hohmann+dv2_hohmann)/1e3:.2f} km/s)')
ax.set_xlabel('Time of Flight (hours)')
ax.set_ylabel('Total \u0394v (km/s)')
ax.set_title('LEO -> GEO Transfer: \u0394v vs. Time of Flight')
ax.legend()
plt.tight_layout()
plt.show()

# LEO to GEO transfer via Lambert r_leo = Re + 400e3 r_geo = 42164e3 # GEO radius v_leo = np.sqrt(mu / r_leo) v_geo = np.sqrt(mu / r_geo) # Departure: 400 km equatorial r_depart = np.array([r_leo, 0, 0]) v_depart = np.array([0, v_leo, 0]) # Arrival: GEO, 180 deg ahead (opposite side of Earth) r_arrive = np.array([-r_geo, 0, 0]) v_arrive_circ = np.array([0, -v_geo, 0]) # Hohmann reference for comparison a_hohmann = (r_leo + r_geo) / 2 hohmann_tof = np.pi * np.sqrt(a_hohmann**3 / mu) dv1_hohmann = np.sqrt(mu * (2/r_leo - 1/a_hohmann)) - v_leo dv2_hohmann = v_geo - np.sqrt(mu * (2/r_geo - 1/a_hohmann)) # Solve Lambert for a range of TOFs tof_hours = np.linspace(3.5, 10, 50) dv_total = np.zeros_like(tof_hours) for i, hours in enumerate(tof_hours): tof = hours * 3600 solutions = sk.lambert(r_depart, r_arrive, tof) v1, v2 = solutions[0] dv1 = np.linalg.norm(v1 - v_depart) dv2 = np.linalg.norm(v_arrive_circ - v2) dv_total[i] = dv1 + dv2 # Use optimal TOF (near Hohmann) optimal_idx = np.argmin(dv_total) optimal_tof = tof_hours[optimal_idx] * 3600 solutions = sk.lambert(r_depart, r_arrive, optimal_tof) v1_transfer, v2_transfer = solutions[0] dv_depart = v1_transfer - v_depart dv_arrive = v_arrive_circ - v2_transfer print(f"Hohmann reference: \u0394v = {dv1_hohmann + dv2_hohmann:.0f} m/s, TOF = {hohmann_tof/3600:.2f} hours") print(f"Lambert optimal: \u0394v = {dv_total[optimal_idx]:.0f} m/s, TOF = {tof_hours[optimal_idx]:.2f} hours") print(f" Departure \u0394v: {np.linalg.norm(dv_depart):.0f} m/s") print(f" Arrival \u0394v: {np.linalg.norm(dv_arrive):.0f} m/s") # Plot dv vs TOF fig, ax = plt.subplots(figsize=(8, 5)) ax.plot(tof_hours, dv_total / 1e3, linewidth=1.5) ax.axvline(hohmann_tof/3600, ls='--', color='red', alpha=0.5, label=f'Hohmann TOF ({hohmann_tof/3600:.1f} h)') ax.axhline((dv1_hohmann + dv2_hohmann)/1e3, ls=':', color='gray', alpha=0.4, label=f'Hohmann \u0394v ({(dv1_hohmann+dv2_hohmann)/1e3:.2f} km/s)') ax.set_xlabel('Time of Flight (hours)') ax.set_ylabel('Total \u0394v (km/s)') ax.set_title('LEO -> GEO Transfer: \u0394v vs. Time of Flight') ax.legend() plt.tight_layout() plt.show()

Hohmann reference: Δv = 3854 m/s, TOF = 5.29 hours
Lambert optimal:   Δv = 3855 m/s, TOF = 5.36 hours
  Departure Δv: 2398 m/s
  Arrival Δv:   1457 m/s

In [7]:

# Apply the optimal Lambert transfer as impulsive maneuvers
t0 = sk.time(2025, 6, 15, 0, 0, 0)
t_depart = t0
t_arrive = t0 + sk.duration(seconds=optimal_tof)

sat_lambert = sk.satstate(time=t0, pos=r_depart, vel=v_depart)
sat_lambert.add_maneuver(t_depart, dv_depart, frame=sk.frame.GCRF)
sat_lambert.add_maneuver(t_arrive, dv_arrive, frame=sk.frame.GCRF)

# Propagate and sample trajectory
settings = sk.propsettings(gravity_degree=6)
t_end = t_arrive + sk.duration(hours=6)
npts = 400
dt_arr = np.linspace(0, (t_end - t0).seconds, npts)
pos_lambert = np.zeros((npts, 3))

for i, dt in enumerate(dt_arr):
    t = t0 + sk.duration(seconds=float(dt))
    state = sat_lambert.propagate(t, propsettings=settings)
    pos_lambert[i] = state.pos / 1e3

# Verify arrival accuracy
state_transfer = np.array([*r_depart, *v1_transfer])
result_check = sk.propagate(state_transfer, t0, end=t_arrive, propsettings=settings)
pos_err = np.linalg.norm(result_check.pos - r_arrive)
print(f"Lambert arrival position error: {pos_err:.0f} m")

final = sat_lambert.propagate(t_end, propsettings=settings)
r_final = np.linalg.norm(final.pos)
print(f"Final orbit radius: {r_final/1e3:.0f} km (GEO = {r_geo/1e3:.0f} km)")

# Apply the optimal Lambert transfer as impulsive maneuvers t0 = sk.time(2025, 6, 15, 0, 0, 0) t_depart = t0 t_arrive = t0 + sk.duration(seconds=optimal_tof) sat_lambert = sk.satstate(time=t0, pos=r_depart, vel=v_depart) sat_lambert.add_maneuver(t_depart, dv_depart, frame=sk.frame.GCRF) sat_lambert.add_maneuver(t_arrive, dv_arrive, frame=sk.frame.GCRF) # Propagate and sample trajectory settings = sk.propsettings(gravity_degree=6) t_end = t_arrive + sk.duration(hours=6) npts = 400 dt_arr = np.linspace(0, (t_end - t0).seconds, npts) pos_lambert = np.zeros((npts, 3)) for i, dt in enumerate(dt_arr): t = t0 + sk.duration(seconds=float(dt)) state = sat_lambert.propagate(t, propsettings=settings) pos_lambert[i] = state.pos / 1e3 # Verify arrival accuracy state_transfer = np.array([*r_depart, *v1_transfer]) result_check = sk.propagate(state_transfer, t0, end=t_arrive, propsettings=settings) pos_err = np.linalg.norm(result_check.pos - r_arrive) print(f"Lambert arrival position error: {pos_err:.0f} m") final = sat_lambert.propagate(t_end, propsettings=settings) r_final = np.linalg.norm(final.pos) print(f"Final orbit radius: {r_final/1e3:.0f} km (GEO = {r_geo/1e3:.0f} km)")

Lambert arrival position error: 285593 m
Final orbit radius: 41876 km (GEO = 42164 km)

In [8]:

# Plot the GTO transfer
fig, ax = plt.subplots(figsize=(9, 9))

# Draw Earth
earth = plt.Circle((0, 0), Re/1e3, color='skyblue', alpha=0.3)
ax.add_patch(earth)

# Draw LEO and GEO circles
theta = np.linspace(0, 2*np.pi, 200)
ax.plot(r_leo/1e3*np.cos(theta), r_leo/1e3*np.sin(theta), '--', alpha=0.3, color='gray', label=f'LEO ({400} km)')
ax.plot(r_geo/1e3*np.cos(theta), r_geo/1e3*np.sin(theta), '--', alpha=0.3, color='gray', label=f'GEO ({35786} km)')

# Transfer trajectory
ax.plot(pos_lambert[:, 0], pos_lambert[:, 1], '-', linewidth=1.5, label='GTO transfer')

# Mark burns with annotations
ax.plot(r_depart[0]/1e3, r_depart[1]/1e3, 'o', color='green', markersize=12, zorder=5)
ax.annotate(f'TLI: \u0394v={np.linalg.norm(dv_depart)/1e3:.2f} km/s',
            xy=(r_depart[0]/1e3, r_depart[1]/1e3),
            xytext=(r_depart[0]/1e3 + 5000, 8000), fontsize=11,
            arrowprops=dict(arrowstyle='->', color='green', lw=1.5))

ax.plot(r_arrive[0]/1e3, r_arrive[1]/1e3, 'o', color='red', markersize=12, zorder=5)
ax.annotate(f'Circularization: \u0394v={np.linalg.norm(dv_arrive)/1e3:.2f} km/s',
            xy=(r_arrive[0]/1e3, r_arrive[1]/1e3),
            xytext=(r_arrive[0]/1e3 + 5000, -8000), fontsize=11,
            arrowprops=dict(arrowstyle='->', color='red', lw=1.5))

ax.set_xlabel('X (km)')
ax.set_ylabel('Y (km)')
ax.set_title('LEO -> GEO Transfer (Lambert Solution)')
ax.set_aspect('equal')
ax.legend(loc='lower right', fontsize=10)
plt.tight_layout()
plt.show()

# Plot the GTO transfer fig, ax = plt.subplots(figsize=(9, 9)) # Draw Earth earth = plt.Circle((0, 0), Re/1e3, color='skyblue', alpha=0.3) ax.add_patch(earth) # Draw LEO and GEO circles theta = np.linspace(0, 2*np.pi, 200) ax.plot(r_leo/1e3*np.cos(theta), r_leo/1e3*np.sin(theta), '--', alpha=0.3, color='gray', label=f'LEO ({400} km)') ax.plot(r_geo/1e3*np.cos(theta), r_geo/1e3*np.sin(theta), '--', alpha=0.3, color='gray', label=f'GEO ({35786} km)') # Transfer trajectory ax.plot(pos_lambert[:, 0], pos_lambert[:, 1], '-', linewidth=1.5, label='GTO transfer') # Mark burns with annotations ax.plot(r_depart[0]/1e3, r_depart[1]/1e3, 'o', color='green', markersize=12, zorder=5) ax.annotate(f'TLI: \u0394v={np.linalg.norm(dv_depart)/1e3:.2f} km/s', xy=(r_depart[0]/1e3, r_depart[1]/1e3), xytext=(r_depart[0]/1e3 + 5000, 8000), fontsize=11, arrowprops=dict(arrowstyle='->', color='green', lw=1.5)) ax.plot(r_arrive[0]/1e3, r_arrive[1]/1e3, 'o', color='red', markersize=12, zorder=5) ax.annotate(f'Circularization: \u0394v={np.linalg.norm(dv_arrive)/1e3:.2f} km/s', xy=(r_arrive[0]/1e3, r_arrive[1]/1e3), xytext=(r_arrive[0]/1e3 + 5000, -8000), fontsize=11, arrowprops=dict(arrowstyle='->', color='red', lw=1.5)) ax.set_xlabel('X (km)') ax.set_ylabel('Y (km)') ax.set_title('LEO -> GEO Transfer (Lambert Solution)') ax.set_aspect('equal') ax.legend(loc='lower right', fontsize=10) plt.tight_layout() plt.show()

Multi-Revolution Solutions¶

For longer times of flight, Lambert's problem admits multi-revolution solutions where the satellite completes one or more full orbits before arriving. sk.lambert() returns all physically possible solutions automatically.

At long TOFs, the highest-revolution solution can be significantly cheaper than the direct (zero-revolution) transfer — it follows a tighter ellipse closer to the original orbit. The minimum-Δv multi-rev solution approaches the Hohmann limit as the number of revolutions increases.

In [9]:

# Multi-revolution Lambert solutions at various TOFs
print(f"Hohmann reference: {dv1_hohmann + dv2_hohmann:.0f} m/s\n")

for hours in [20, 30, 40, 50]:
    tof = hours * 3600
    solutions = sk.lambert(r_depart, r_arrive, tof)
    print(f"TOF = {hours}h: {len(solutions)} solution(s)")
    best_dv = np.inf
    for i, (v1, v2) in enumerate(solutions):
        dv = np.linalg.norm(v1 - v_depart) + np.linalg.norm(v_arrive_circ - v2)
        if i == 0:
            label = "0-rev"
        else:
            label = f"{(i+1)//2}-rev {'short' if i%2==1 else 'long'}"
        marker = " <-- best" if dv < best_dv else ""
        if dv < best_dv:
            best_dv = dv
        print(f"  {label:12s}: {dv:7.0f} m/s{marker}")
    print()

# Multi-revolution Lambert solutions at various TOFs print(f"Hohmann reference: {dv1_hohmann + dv2_hohmann:.0f} m/s\n") for hours in [20, 30, 40, 50]: tof = hours * 3600 solutions = sk.lambert(r_depart, r_arrive, tof) print(f"TOF = {hours}h: {len(solutions)} solution(s)") best_dv = np.inf for i, (v1, v2) in enumerate(solutions): dv = np.linalg.norm(v1 - v_depart) + np.linalg.norm(v_arrive_circ - v2) if i == 0: label = "0-rev" else: label = f"{(i+1)//2}-rev {'short' if i%2==1 else 'long'}" marker = " <-- best" if dv < best_dv else "" if dv < best_dv: best_dv = dv print(f" {label:12s}: {dv:7.0f} m/s{marker}") print()

Hohmann reference: 3854 m/s

TOF = 20h: 3 solution(s)
  0-rev       :    6450 m/s <-- best
  1-rev short :    4486 m/s <-- best
  1-rev long  :    5720 m/s

TOF = 30h: 5 solution(s)
  0-rev       :    7064 m/s <-- best
  1-rev short :    5673 m/s <-- best
  1-rev long  :    6776 m/s
  2-rev short :    4252 m/s <-- best
  2-rev long  :    4917 m/s

TOF = 40h: 7 solution(s)
  0-rev       :    7410 m/s <-- best
  1-rev short :    6316 m/s <-- best
  1-rev long  :    7249 m/s
  2-rev short :    5243 m/s <-- best
  2-rev long  :    5966 m/s
  3-rev short :    4114 m/s <-- best
  3-rev long  :    4528 m/s

TOF = 50h: 9 solution(s)
  0-rev       :    7634 m/s <-- best
  1-rev short :    6721 m/s <-- best
  1-rev long  :    7530 m/s
  2-rev short :    5853 m/s <-- best
  2-rev long  :    6510 m/s
  3-rev short :    4964 m/s <-- best
  3-rev long  :    5488 m/s
  4-rev short :    4023 m/s <-- best
  4-rev long  :    4298 m/s

3. Low-Thrust Orbit Raising (Continuous)¶

Electric propulsion produces very low thrust over long durations, spiraling outward gradually rather than making discrete burns. This is modeled using sk.thrust.constant() with the satproperties object.

We'll raise the same 400 km → 800 km orbit as the Hohmann example, but using continuous in-track thrust. The trajectory will be a slow spiral rather than an elliptical transfer, and requires significantly more time.

In [10]:

# Low-thrust parameters
# For a 100 kg satellite with a 50 mN thruster:
thrust_accel = 5e-4  # m/s^2
thrust_duration_hours = 120  # ~5 days of continuous thrust

dv_applied = thrust_accel * thrust_duration_hours * 3600

print(f"Thrust acceleration: {thrust_accel*1e3:.1f} mm/s\u00b2")
print(f"Thrust duration: {thrust_duration_hours} hours ({thrust_duration_hours/24:.1f} days)")
print(f"Total \u0394v applied: {dv_applied:.0f} m/s")
print(f"(Hohmann \u0394v for same altitude change: {dv1 + dv2:.0f} m/s)")
print("Low-thrust requires more \u0394v due to gravity losses")

# Low-thrust parameters # For a 100 kg satellite with a 50 mN thruster: thrust_accel = 5e-4 # m/s^2 thrust_duration_hours = 120 # ~5 days of continuous thrust dv_applied = thrust_accel * thrust_duration_hours * 3600 print(f"Thrust acceleration: {thrust_accel*1e3:.1f} mm/s\u00b2") print(f"Thrust duration: {thrust_duration_hours} hours ({thrust_duration_hours/24:.1f} days)") print(f"Total \u0394v applied: {dv_applied:.0f} m/s") print(f"(Hohmann \u0394v for same altitude change: {dv1 + dv2:.0f} m/s)") print("Low-thrust requires more \u0394v due to gravity losses")

Thrust acceleration: 0.5 mm/s²
Thrust duration: 120 hours (5.0 days)
Total Δv applied: 216 m/s
(Hohmann Δv for same altitude change: 4999 m/s)
Low-thrust requires more Δv due to gravity losses

In [11]:

# Set up continuous tangential thrust in RTN frame
t0 = sk.time(2025, 6, 15, 0, 0, 0)
t_thrust_end = t0 + sk.duration(hours=thrust_duration_hours)
t_coast_end = t_thrust_end + sk.duration(hours=3)

# Constant tangential acceleration in RTN: [radial, tangential, normal]
thrust_arc = sk.thrust.constant(
    [0, thrust_accel, 0],
    t0,
    t_thrust_end,
    frame=sk.frame.RTN,
)

props = sk.satproperties(thrusts=[thrust_arc])

# Initial state: 400 km circular (same as Hohmann example)
state0 = np.array([r1, 0, 0, 0, v_circ_1, 0])

# Use lower-order integrator for robustness with long-duration thrust
settings = sk.propsettings(
    gravity_degree=1,
    integrator=sk.integrator.rkts54,
)

result = sk.propagate(
    state0, t0,
    end=t_coast_end,
    propsettings=settings,
    satproperties=props,
)

# Sample the trajectory
npts = 1000
dt_arr = np.linspace(0, (t_coast_end - t0).seconds, npts)
time_arr = [t0 + sk.duration(seconds=float(dt)) for dt in dt_arr]

pos_lt = np.zeros((npts, 3))
alt_lt = np.zeros(npts)

for i, t in enumerate(time_arr):
    s = result.interp(t)
    pos_lt[i] = s[0:3] / 1e3
    alt_lt[i] = np.linalg.norm(s[0:3]) - Re

final_alt = alt_lt[-1]
print(f"Final altitude: {final_alt/1e3:.1f} km (target: {alt_final/1e3:.0f} km)")
print(f"Thrust duration: {thrust_duration_hours/24:.1f} days")
print(f"Hohmann transfer time: {transfer_time/60:.0f} minutes")

# Set up continuous tangential thrust in RTN frame t0 = sk.time(2025, 6, 15, 0, 0, 0) t_thrust_end = t0 + sk.duration(hours=thrust_duration_hours) t_coast_end = t_thrust_end + sk.duration(hours=3) # Constant tangential acceleration in RTN: [radial, tangential, normal] thrust_arc = sk.thrust.constant( [0, thrust_accel, 0], t0, t_thrust_end, frame=sk.frame.RTN, ) props = sk.satproperties(thrusts=[thrust_arc]) # Initial state: 400 km circular (same as Hohmann example) state0 = np.array([r1, 0, 0, 0, v_circ_1, 0]) # Use lower-order integrator for robustness with long-duration thrust settings = sk.propsettings( gravity_degree=1, integrator=sk.integrator.rkts54, ) result = sk.propagate( state0, t0, end=t_coast_end, propsettings=settings, satproperties=props, ) # Sample the trajectory npts = 1000 dt_arr = np.linspace(0, (t_coast_end - t0).seconds, npts) time_arr = [t0 + sk.duration(seconds=float(dt)) for dt in dt_arr] pos_lt = np.zeros((npts, 3)) alt_lt = np.zeros(npts) for i, t in enumerate(time_arr): s = result.interp(t) pos_lt[i] = s[0:3] / 1e3 alt_lt[i] = np.linalg.norm(s[0:3]) - Re final_alt = alt_lt[-1] print(f"Final altitude: {final_alt/1e3:.1f} km (target: {alt_final/1e3:.0f} km)") print(f"Thrust duration: {thrust_duration_hours/24:.1f} days") print(f"Hohmann transfer time: {transfer_time/60:.0f} minutes")

Final altitude: 797.0 km (target: 800 km)
Thrust duration: 5.0 days
Hohmann transfer time: 48 minutes

In [12]:

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Left: spiral trajectory
ax = axes[0]
earth = plt.Circle((0, 0), Re/1e3, color='skyblue', alpha=0.3)
ax.add_patch(earth)
ax.plot(pos_lt[:, 0], pos_lt[:, 1], '-', linewidth=0.5, alpha=0.7, label='Low-thrust spiral')

# Mark initial and target orbits
theta = np.linspace(0, 2*np.pi, 200)
ax.plot(r1/1e3*np.cos(theta), r1/1e3*np.sin(theta), '--', alpha=0.3, color='gray', label=f'{alt_initial/1e3:.0f} km')
ax.plot(r2/1e3*np.cos(theta), r2/1e3*np.sin(theta), '--', alpha=0.3, color='gray', label=f'{alt_final/1e3:.0f} km')

ax.set_xlabel('X (km)')
ax.set_ylabel('Y (km)')
ax.set_title('Low-Thrust Orbit Raising')
ax.set_aspect('equal')
ax.legend(fontsize=10)

# Right: altitude comparison
ax = axes[1]
hours = dt_arr / 3600
ax.plot(hours, alt_lt / 1e3, linewidth=1.5, label='Low-thrust')
ax.axhline(alt_initial/1e3, ls='--', alpha=0.4, color='gray')
ax.axhline(alt_final/1e3, ls='--', alpha=0.4, color='gray')
ax.axvline(thrust_duration_hours, ls=':', alpha=0.3, color='red', label='Thrust cutoff')
ax.set_xlabel('Time (hours)')
ax.set_ylabel('Altitude (km)')
ax.set_title('Altitude vs Time')
ax.legend()

plt.tight_layout()
plt.show()

fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # Left: spiral trajectory ax = axes[0] earth = plt.Circle((0, 0), Re/1e3, color='skyblue', alpha=0.3) ax.add_patch(earth) ax.plot(pos_lt[:, 0], pos_lt[:, 1], '-', linewidth=0.5, alpha=0.7, label='Low-thrust spiral') # Mark initial and target orbits theta = np.linspace(0, 2*np.pi, 200) ax.plot(r1/1e3*np.cos(theta), r1/1e3*np.sin(theta), '--', alpha=0.3, color='gray', label=f'{alt_initial/1e3:.0f} km') ax.plot(r2/1e3*np.cos(theta), r2/1e3*np.sin(theta), '--', alpha=0.3, color='gray', label=f'{alt_final/1e3:.0f} km') ax.set_xlabel('X (km)') ax.set_ylabel('Y (km)') ax.set_title('Low-Thrust Orbit Raising') ax.set_aspect('equal') ax.legend(fontsize=10) # Right: altitude comparison ax = axes[1] hours = dt_arr / 3600 ax.plot(hours, alt_lt / 1e3, linewidth=1.5, label='Low-thrust') ax.axhline(alt_initial/1e3, ls='--', alpha=0.4, color='gray') ax.axhline(alt_final/1e3, ls='--', alpha=0.4, color='gray') ax.axvline(thrust_duration_hours, ls=':', alpha=0.3, color='red', label='Thrust cutoff') ax.set_xlabel('Time (hours)') ax.set_ylabel('Altitude (km)') ax.set_title('Altitude vs Time') ax.legend() plt.tight_layout() plt.show()

Comparison: Impulsive vs. Low-Thrust¶

	Hohmann (Impulsive)	Low-Thrust (Continuous)
Mechanism	Two discrete burns	Sustained tangential acceleration
Δv	Minimum for 2-impulse	Higher due to gravity losses
Transfer time	~Half orbit period	Hours to days
Trajectory	Elliptical coast arc	Spiral
API	satstate.add_maneuver()	satproperties(thrusts=[...])
Frame	sk.frame.RTN or sk.frame.GCRF	sk.frame.RTN or sk.frame.GCRF

The choice depends on the propulsion system: chemical rockets deliver high thrust for short durations (impulsive), while electric propulsion provides low thrust over long arcs (continuous).