Keplerian Orbital Elements¶

This tutorial demonstrates how to work with Keplerian orbital elements in satkit. Keplerian elements are a set of six parameters that uniquely describe the shape, orientation, and current position of a satellite orbit.

The Six Keplerian Elements¶

Element	Symbol	Description
Semi-major axis	$a$	Half the longest diameter of the orbital ellipse; sets the orbit size
Eccentricity	$e$	Shape of the ellipse: 0 = circular, 0 < $e$ < 1 = elliptical
Inclination	$i$	Tilt of the orbit plane relative to the equatorial plane
Right Ascension of Ascending Node	$\Omega$ (RAAN)	Angle from the vernal equinox to where the orbit crosses the equator going north
Argument of Perigee	$\omega$	Angle from the ascending node to the closest point of the orbit (perigee)
True Anomaly	$\nu$	Current angular position of the satellite along the orbit, measured from perigee

The first two elements ($a$, $e$) define the orbit shape, the next two ($i$, $\Omega$) define the orbit plane orientation, $\omega$ orients the ellipse within the plane, and $\nu$ locates the satellite on the orbit.

In [1]:

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

Creating Keplerian Elements¶

The sk.kepler constructor takes six parameters: semi-major axis (meters), eccentricity, inclination (radians), RAAN (radians), argument of perigee (radians), and anomaly (radians).

By default, the sixth argument is the true anomaly. You can alternatively specify the anomaly using keyword arguments mean_anomaly or eccentric_anomaly.

In [2]:

# Define an ISS-like orbit: ~400 km altitude, 51.6 deg inclination, nearly circular
altitude = 400e3  # meters
a = sk.consts.earth_radius + altitude
e = 0.001
i = math.radians(51.6)
raan = math.radians(45.0)
w = math.radians(30.0)
nu = math.radians(0.0)

# Create with true anomaly (default)
kep = sk.kepler(a, e, i, raan, w, nu)
print("Created with true anomaly:")
print(kep)
print(f"  Semi-major axis: {kep.a / 1e3:.1f} km")
print(f"  Eccentricity:    {kep.eccen:.4f}")
print(f"  Inclination:     {math.degrees(kep.inclination):.1f} deg")
print(f"  RAAN:            {math.degrees(kep.raan):.1f} deg")
print(f"  Arg of Perigee:  {math.degrees(kep.w):.1f} deg")
print(f"  True Anomaly:    {math.degrees(kep.nu):.1f} deg")

# Define an ISS-like orbit: ~400 km altitude, 51.6 deg inclination, nearly circular altitude = 400e3 # meters a = sk.consts.earth_radius + altitude e = 0.001 i = math.radians(51.6) raan = math.radians(45.0) w = math.radians(30.0) nu = math.radians(0.0) # Create with true anomaly (default) kep = sk.kepler(a, e, i, raan, w, nu) print("Created with true anomaly:") print(kep) print(f" Semi-major axis: {kep.a / 1e3:.1f} km") print(f" Eccentricity: {kep.eccen:.4f}") print(f" Inclination: {math.degrees(kep.inclination):.1f} deg") print(f" RAAN: {math.degrees(kep.raan):.1f} deg") print(f" Arg of Perigee: {math.degrees(kep.w):.1f} deg") print(f" True Anomaly: {math.degrees(kep.nu):.1f} deg")

Created with true anomaly:
Keplerian Elements:
  a = 6778137 m
  e = 0.001
  i = 0.901 rad
  Ω = 0.785 rad
  ω = 0.524 rad
  ν = 0.000 rad

  Semi-major axis: 6778.1 km
  Eccentricity:    0.0010
  Inclination:     51.6 deg
  RAAN:            45.0 deg
  Arg of Perigee:  30.0 deg
  True Anomaly:    0.0 deg

In [3]:

# Create with mean anomaly (keyword argument)
kep_mean = sk.kepler(a, e, i, raan, w, mean_anomaly=math.radians(10.0))
print("Created with mean anomaly = 10 deg:")
print(f"  True anomaly:      {math.degrees(kep_mean.true_anomaly):.4f} deg")
print(f"  Mean anomaly:      {math.degrees(kep_mean.mean_anomaly):.4f} deg")
print(f"  Eccentric anomaly: {math.degrees(kep_mean.eccentric_anomaly):.4f} deg")

# Create with eccentric anomaly
kep_ecc = sk.kepler(a, e, i, raan, w, eccentric_anomaly=math.radians(10.0))
print("\nCreated with eccentric anomaly = 10 deg:")
print(f"  True anomaly:      {math.degrees(kep_ecc.true_anomaly):.4f} deg")
print(f"  Mean anomaly:      {math.degrees(kep_ecc.mean_anomaly):.4f} deg")
print(f"  Eccentric anomaly: {math.degrees(kep_ecc.eccentric_anomaly):.4f} deg")

# Create with mean anomaly (keyword argument) kep_mean = sk.kepler(a, e, i, raan, w, mean_anomaly=math.radians(10.0)) print("Created with mean anomaly = 10 deg:") print(f" True anomaly: {math.degrees(kep_mean.true_anomaly):.4f} deg") print(f" Mean anomaly: {math.degrees(kep_mean.mean_anomaly):.4f} deg") print(f" Eccentric anomaly: {math.degrees(kep_mean.eccentric_anomaly):.4f} deg") # Create with eccentric anomaly kep_ecc = sk.kepler(a, e, i, raan, w, eccentric_anomaly=math.radians(10.0)) print("\nCreated with eccentric anomaly = 10 deg:") print(f" True anomaly: {math.degrees(kep_ecc.true_anomaly):.4f} deg") print(f" Mean anomaly: {math.degrees(kep_ecc.mean_anomaly):.4f} deg") print(f" Eccentric anomaly: {math.degrees(kep_ecc.eccentric_anomaly):.4f} deg")

Created with mean anomaly = 10 deg:
  True anomaly:      10.0199 deg
  Mean anomaly:      10.0000 deg
  Eccentric anomaly: 10.0100 deg

Created with eccentric anomaly = 10 deg:
  True anomaly:      10.0100 deg
  Mean anomaly:      9.9901 deg
  Eccentric anomaly: 10.0000 deg

Converting to Cartesian Coordinates¶

Use to_pv() to convert Keplerian elements to position and velocity vectors in the perifocal/inertial frame. The position is in meters and velocity in meters per second.

In [4]:

pos, vel = kep.to_pv()

print(f"Position (km): [{pos[0]/1e3:.3f}, {pos[1]/1e3:.3f}, {pos[2]/1e3:.3f}]")
print(f"Velocity (km/s): [{vel[0]/1e3:.4f}, {vel[1]/1e3:.4f}, {vel[2]/1e3:.4f}]")
print(f"\nPosition magnitude: {np.linalg.norm(pos)/1e3:.1f} km")
print(f"Velocity magnitude: {np.linalg.norm(vel)/1e3:.4f} km/s")

pos, vel = kep.to_pv() print(f"Position (km): [{pos[0]/1e3:.3f}, {pos[1]/1e3:.3f}, {pos[2]/1e3:.3f}]") print(f"Velocity (km/s): [{vel[0]/1e3:.4f}, {vel[1]/1e3:.4f}, {vel[2]/1e3:.4f}]") print(f"\nPosition magnitude: {np.linalg.norm(pos)/1e3:.1f} km") print(f"Velocity magnitude: {np.linalg.norm(vel)/1e3:.4f} km/s")

Position (km): [2659.543, 5633.644, 2653.335]
Velocity (km/s): [-5.6338, 0.2059, 5.2098]

Position magnitude: 6771.4 km
Velocity magnitude: 7.6762 km/s

Converting from Cartesian Coordinates¶

Use sk.kepler.from_pv() to recover Keplerian elements from position and velocity vectors. This is the inverse of to_pv().

The example below uses a test case from Vallado's Fundamentals of Astrodynamics and Applications, Example 2-6.

In [5]:

# Vallado Example 2-6
r = np.array([6524.834, 6862.875, 6448.296]) * 1e3  # meters
v = np.array([4.901327, 5.533756, -1.976341]) * 1e3  # m/s

kep_from_pv = sk.kepler.from_pv(r, v)

print(f"Semi-major axis: {kep_from_pv.a / 1e3:.1f} km")
print(f"Eccentricity:    {kep_from_pv.eccen:.5f}")
print(f"Inclination:     {math.degrees(kep_from_pv.inclination):.2f} deg")
print(f"RAAN:            {math.degrees(kep_from_pv.raan):.2f} deg")
print(f"Arg of Perigee:  {math.degrees(kep_from_pv.w):.2f} deg")
print(f"True Anomaly:    {math.degrees(kep_from_pv.nu):.3f} deg")

# Vallado Example 2-6 r = np.array([6524.834, 6862.875, 6448.296]) * 1e3 # meters v = np.array([4.901327, 5.533756, -1.976341]) * 1e3 # m/s kep_from_pv = sk.kepler.from_pv(r, v) print(f"Semi-major axis: {kep_from_pv.a / 1e3:.1f} km") print(f"Eccentricity: {kep_from_pv.eccen:.5f}") print(f"Inclination: {math.degrees(kep_from_pv.inclination):.2f} deg") print(f"RAAN: {math.degrees(kep_from_pv.raan):.2f} deg") print(f"Arg of Perigee: {math.degrees(kep_from_pv.w):.2f} deg") print(f"True Anomaly: {math.degrees(kep_from_pv.nu):.3f} deg")

Semi-major axis: 36127.3 km
Eccentricity:    0.83285
Inclination:     87.87 deg
RAAN:            227.90 deg
Arg of Perigee:  53.38 deg
True Anomaly:    92.335 deg

In [6]:

# Round-trip test: convert to PV and back
pos2, vel2 = kep_from_pv.to_pv()
kep_roundtrip = sk.kepler.from_pv(pos2, vel2)

print("Round-trip consistency check:")
print(f"  a difference:     {abs(kep_from_pv.a - kep_roundtrip.a):.6e} m")
print(f"  e difference:     {abs(kep_from_pv.eccen - kep_roundtrip.eccen):.6e}")
print(f"  i difference:     {abs(kep_from_pv.inclination - kep_roundtrip.inclination):.6e} rad")
print(f"  RAAN difference:  {abs(kep_from_pv.raan - kep_roundtrip.raan):.6e} rad")
print(f"  w difference:     {abs(kep_from_pv.w - kep_roundtrip.w):.6e} rad")
print(f"  nu difference:    {abs(kep_from_pv.nu - kep_roundtrip.nu):.6e} rad")

# Round-trip test: convert to PV and back pos2, vel2 = kep_from_pv.to_pv() kep_roundtrip = sk.kepler.from_pv(pos2, vel2) print("Round-trip consistency check:") print(f" a difference: {abs(kep_from_pv.a - kep_roundtrip.a):.6e} m") print(f" e difference: {abs(kep_from_pv.eccen - kep_roundtrip.eccen):.6e}") print(f" i difference: {abs(kep_from_pv.inclination - kep_roundtrip.inclination):.6e} rad") print(f" RAAN difference: {abs(kep_from_pv.raan - kep_roundtrip.raan):.6e} rad") print(f" w difference: {abs(kep_from_pv.w - kep_roundtrip.w):.6e} rad") print(f" nu difference: {abs(kep_from_pv.nu - kep_roundtrip.nu):.6e} rad")

Round-trip consistency check:
  a difference:     5.215406e-08 m
  e difference:     1.110223e-16
  i difference:     2.220446e-16 rad
  RAAN difference:  0.000000e+00 rad
  w difference:     2.220446e-16 rad
  nu difference:    2.220446e-16 rad

Orbital Period and Mean Motion¶

The period property returns the orbital period in seconds, and mean_motion returns the mean angular rate in radians per second.

In [7]:

print(f"Orbital period:  {kep.period:.1f} seconds = {kep.period / 60:.2f} minutes")
print(f"Mean motion:     {kep.mean_motion:.6e} rad/s")
print(f"Mean motion:     {kep.mean_motion * 86400 / (2 * math.pi):.4f} revs/day")

# Verify the relationship: period = 2*pi / mean_motion
print(f"\n2*pi / mean_motion = {2 * math.pi / kep.mean_motion:.1f} seconds (should match period)")

print(f"Orbital period: {kep.period:.1f} seconds = {kep.period / 60:.2f} minutes") print(f"Mean motion: {kep.mean_motion:.6e} rad/s") print(f"Mean motion: {kep.mean_motion * 86400 / (2 * math.pi):.4f} revs/day") # Verify the relationship: period = 2*pi / mean_motion print(f"\n2*pi / mean_motion = {2 * math.pi / kep.mean_motion:.1f} seconds (should match period)")

Orbital period:  5553.6 seconds = 92.56 minutes
Mean motion:     1.131367e-03 rad/s
Mean motion:     15.5574 revs/day

2*pi / mean_motion = 5553.6 seconds (should match period)

Keplerian Propagation¶

The propagate() method advances the orbit forward in time using two-body (Keplerian) dynamics. Only the true anomaly changes; all other elements remain constant. This is an analytical solution, so it is exact (within the two-body assumption) and very fast.

Below we propagate an orbit over one full period, sampling positions to visualize the orbit.

In [8]:

# Propagate over one orbit, sampling 200 points
npts = 200
T = kep.period
times = np.linspace(0, T, npts)

positions = np.zeros((npts, 3))
for idx, dt in enumerate(times):
    kep_t = kep.propagate(dt)
    pos_t, _ = kep_t.to_pv()
    positions[idx, :] = pos_t / 1e3  # convert to km

# 3D orbit plot
fig = plt.figure(figsize=(9, 9))
ax = fig.add_subplot(111, projection='3d')
ax.plot(positions[:, 0], positions[:, 1], positions[:, 2], 'b-', linewidth=1.5)
ax.scatter(*positions[0, :], color='green', s=80, label='Start (perigee)', zorder=5)

# Draw Earth sphere
u = np.linspace(0, 2 * np.pi, 40)
v = np.linspace(0, np.pi, 20)
Re = sk.consts.earth_radius / 1e3
x_e = Re * np.outer(np.cos(u), np.sin(v))
y_e = Re * np.outer(np.sin(u), np.sin(v))
z_e = Re * np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(x_e, y_e, z_e, alpha=0.2, color='skyblue')

ax.set_xlabel('X (km)')
ax.set_ylabel('Y (km)')
ax.set_zlabel('Z (km)')
ax.set_title('ISS-like Orbit (One Period)')
ax.legend()
plt.tight_layout()
plt.show()

# Propagate over one orbit, sampling 200 points npts = 200 T = kep.period times = np.linspace(0, T, npts) positions = np.zeros((npts, 3)) for idx, dt in enumerate(times): kep_t = kep.propagate(dt) pos_t, _ = kep_t.to_pv() positions[idx, :] = pos_t / 1e3 # convert to km # 3D orbit plot fig = plt.figure(figsize=(9, 9)) ax = fig.add_subplot(111, projection='3d') ax.plot(positions[:, 0], positions[:, 1], positions[:, 2], 'b-', linewidth=1.5) ax.scatter(*positions[0, :], color='green', s=80, label='Start (perigee)', zorder=5) # Draw Earth sphere u = np.linspace(0, 2 * np.pi, 40) v = np.linspace(0, np.pi, 20) Re = sk.consts.earth_radius / 1e3 x_e = Re * np.outer(np.cos(u), np.sin(v)) y_e = Re * np.outer(np.sin(u), np.sin(v)) z_e = Re * np.outer(np.ones(np.size(u)), np.cos(v)) ax.plot_surface(x_e, y_e, z_e, alpha=0.2, color='skyblue') ax.set_xlabel('X (km)') ax.set_ylabel('Y (km)') ax.set_zlabel('Z (km)') ax.set_title('ISS-like Orbit (One Period)') ax.legend() plt.tight_layout() plt.show()

Comparing Kepler vs. Numerical Propagation¶

Keplerian propagation assumes a point-mass Earth with no perturbations. In reality, the Earth's oblateness ($J_2$ and higher-order gravity harmonics), atmospheric drag, third-body gravity from the Sun and Moon, and solar radiation pressure all perturb the orbit.

Here we propagate the same initial orbit using both:

• Kepler (two-body analytical)
• High-precision numerical propagation with full force models

and plot the position difference over time to visualize the perturbation effects.

In [9]:

# Define orbit: 550 km circular, 45 deg inclination
alt = 550e3
a_circ = sk.consts.earth_radius + alt
kep0 = sk.kepler(a_circ, 0.001, math.radians(45.0), 0.0, 0.0, 0.0)

# Get initial state vector for numerical propagation
state0 = np.concatenate(kep0.to_pv())

# Propagation epoch and duration (6 hours)
epoch = sk.time(2025, 6, 15, 0, 0, 0)
prop_duration = sk.duration(hours=6)

# Numerical propagation with high-precision force model
settings = sk.propsettings(
    abs_error=1e-12,
    rel_error=1e-12,
    gravity_degree=10,
)
settings.precompute_terms(epoch, epoch + prop_duration)

result = sk.propagate(
    state0,
    epoch,
    epoch + prop_duration,
    propsettings=settings,
)

# Sample at 1-minute intervals
npts = 361
dt_seconds = np.linspace(0, prop_duration.seconds, npts)
time_array = [epoch + sk.duration(seconds=float(dt)) for dt in dt_seconds]

# Get positions from both propagators
pos_kepler = np.zeros((npts, 3))
pos_numerical = np.zeros((npts, 3))

for idx, (dt, t) in enumerate(zip(dt_seconds, time_array)):
    # Kepler propagation
    kep_t = kep0.propagate(dt)
    pos_k, _ = kep_t.to_pv()
    pos_kepler[idx, :] = pos_k

    # Numerical propagation (interpolated)
    state_num = result.interp(t)
    pos_numerical[idx, :] = state_num[0:3]

# Compute position difference
pos_diff = pos_numerical - pos_kepler
pos_diff_mag = np.linalg.norm(pos_diff, axis=1)

# Plot
hours = dt_seconds / 3600

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8), sharex=True)

ax1.plot(hours, pos_diff[:, 0] / 1e3, label='X')
ax1.plot(hours, pos_diff[:, 1] / 1e3, label='Y')
ax1.plot(hours, pos_diff[:, 2] / 1e3, label='Z')
ax1.set_ylabel('Position Difference (km)')
ax1.set_title('Numerical - Keplerian Position Difference (per component)')
ax1.legend()
ax1.grid(True, alpha=0.3)

ax2.plot(hours, pos_diff_mag / 1e3, 'k-', linewidth=1.5)
ax2.set_xlabel('Time (hours)')
ax2.set_ylabel('Position Difference (km)')
ax2.set_title('Total Position Difference Magnitude')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"Maximum position difference over 6 hours: {pos_diff_mag[-1]/1e3:.2f} km")

# Define orbit: 550 km circular, 45 deg inclination alt = 550e3 a_circ = sk.consts.earth_radius + alt kep0 = sk.kepler(a_circ, 0.001, math.radians(45.0), 0.0, 0.0, 0.0) # Get initial state vector for numerical propagation state0 = np.concatenate(kep0.to_pv()) # Propagation epoch and duration (6 hours) epoch = sk.time(2025, 6, 15, 0, 0, 0) prop_duration = sk.duration(hours=6) # Numerical propagation with high-precision force model settings = sk.propsettings( abs_error=1e-12, rel_error=1e-12, gravity_degree=10, ) settings.precompute_terms(epoch, epoch + prop_duration) result = sk.propagate( state0, epoch, epoch + prop_duration, propsettings=settings, ) # Sample at 1-minute intervals npts = 361 dt_seconds = np.linspace(0, prop_duration.seconds, npts) time_array = [epoch + sk.duration(seconds=float(dt)) for dt in dt_seconds] # Get positions from both propagators pos_kepler = np.zeros((npts, 3)) pos_numerical = np.zeros((npts, 3)) for idx, (dt, t) in enumerate(zip(dt_seconds, time_array)): # Kepler propagation kep_t = kep0.propagate(dt) pos_k, _ = kep_t.to_pv() pos_kepler[idx, :] = pos_k # Numerical propagation (interpolated) state_num = result.interp(t) pos_numerical[idx, :] = state_num[0:3] # Compute position difference pos_diff = pos_numerical - pos_kepler pos_diff_mag = np.linalg.norm(pos_diff, axis=1) # Plot hours = dt_seconds / 3600 fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8), sharex=True) ax1.plot(hours, pos_diff[:, 0] / 1e3, label='X') ax1.plot(hours, pos_diff[:, 1] / 1e3, label='Y') ax1.plot(hours, pos_diff[:, 2] / 1e3, label='Z') ax1.set_ylabel('Position Difference (km)') ax1.set_title('Numerical - Keplerian Position Difference (per component)') ax1.legend() ax1.grid(True, alpha=0.3) ax2.plot(hours, pos_diff_mag / 1e3, 'k-', linewidth=1.5) ax2.set_xlabel('Time (hours)') ax2.set_ylabel('Position Difference (km)') ax2.set_title('Total Position Difference Magnitude') ax2.grid(True, alpha=0.3) plt.tight_layout() plt.show() print(f"Maximum position difference over 6 hours: {pos_diff_mag[-1]/1e3:.2f} km")

Maximum position difference over 6 hours: 292.06 km