{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Keplerian Orbital Elements\n",
"\n",
"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.\n",
"\n",
"## The Six Keplerian Elements\n",
"\n",
"| Element | Symbol | Description |\n",
"|---|---|---|\n",
"| Semi-major axis | $a$ | Half the longest diameter of the orbital ellipse; sets the orbit size |\n",
"| Eccentricity | $e$ | Shape of the ellipse: 0 = circular, 0 < $e$ < 1 = elliptical |\n",
"| Inclination | $i$ | Tilt of the orbit plane relative to the equatorial plane |\n",
"| Right Ascension of Ascending Node | $\\Omega$ (RAAN) | Angle from the vernal equinox to where the orbit crosses the equator going north |\n",
"| Argument of Perigee | $\\omega$ | Angle from the ascending node to the closest point of the orbit (perigee) |\n",
"| True Anomaly | $\\nu$ | Current angular position of the satellite along the orbit, measured from perigee |\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import satkit as sk\n",
"import numpy as np\n",
"import math\n",
"import matplotlib.pyplot as plt\n",
"import scienceplots # noqa: F401\n",
"plt.style.use([\"science\", \"no-latex\", \"../satkit.mplstyle\"])\n",
"%config InlineBackend.figure_formats = ['svg']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Creating Keplerian Elements\n",
"\n",
"The `sk.kepler` constructor takes six parameters: semi-major axis (meters), eccentricity, inclination (radians), RAAN (radians), argument of perigee (radians), and anomaly (radians).\n",
"\n",
"By default, the sixth argument is the true anomaly. You can alternatively specify the anomaly using keyword arguments `mean_anomaly` or `eccentric_anomaly`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Define an ISS-like orbit: ~400 km altitude, 51.6 deg inclination, nearly circular\n",
"altitude = 400e3 # meters\n",
"a = sk.consts.earth_radius + altitude\n",
"e = 0.001\n",
"i = math.radians(51.6)\n",
"raan = math.radians(45.0)\n",
"w = math.radians(30.0)\n",
"nu = math.radians(0.0)\n",
"\n",
"# Create with true anomaly (default)\n",
"kep = sk.kepler(a, e, i, raan, w, nu)\n",
"print(\"Created with true anomaly:\")\n",
"print(kep)\n",
"print(f\" Semi-major axis: {kep.a / 1e3:.1f} km\")\n",
"print(f\" Eccentricity: {kep.eccen:.4f}\")\n",
"print(f\" Inclination: {math.degrees(kep.inclination):.1f} deg\")\n",
"print(f\" RAAN: {math.degrees(kep.raan):.1f} deg\")\n",
"print(f\" Arg of Perigee: {math.degrees(kep.w):.1f} deg\")\n",
"print(f\" True Anomaly: {math.degrees(kep.nu):.1f} deg\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Create with mean anomaly (keyword argument)\n",
"kep_mean = sk.kepler(a, e, i, raan, w, mean_anomaly=math.radians(10.0))\n",
"print(\"Created with mean anomaly = 10 deg:\")\n",
"print(f\" True anomaly: {math.degrees(kep_mean.true_anomaly):.4f} deg\")\n",
"print(f\" Mean anomaly: {math.degrees(kep_mean.mean_anomaly):.4f} deg\")\n",
"print(f\" Eccentric anomaly: {math.degrees(kep_mean.eccentric_anomaly):.4f} deg\")\n",
"\n",
"# Create with eccentric anomaly\n",
"kep_ecc = sk.kepler(a, e, i, raan, w, eccentric_anomaly=math.radians(10.0))\n",
"print(\"\\nCreated with eccentric anomaly = 10 deg:\")\n",
"print(f\" True anomaly: {math.degrees(kep_ecc.true_anomaly):.4f} deg\")\n",
"print(f\" Mean anomaly: {math.degrees(kep_ecc.mean_anomaly):.4f} deg\")\n",
"print(f\" Eccentric anomaly: {math.degrees(kep_ecc.eccentric_anomaly):.4f} deg\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Converting to Cartesian Coordinates\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"pos, vel = kep.to_pv()\n",
"\n",
"print(f\"Position (km): [{pos[0]/1e3:.3f}, {pos[1]/1e3:.3f}, {pos[2]/1e3:.3f}]\")\n",
"print(f\"Velocity (km/s): [{vel[0]/1e3:.4f}, {vel[1]/1e3:.4f}, {vel[2]/1e3:.4f}]\")\n",
"print(f\"\\nPosition magnitude: {np.linalg.norm(pos)/1e3:.1f} km\")\n",
"print(f\"Velocity magnitude: {np.linalg.norm(vel)/1e3:.4f} km/s\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Converting from Cartesian Coordinates\n",
"\n",
"Use `sk.kepler.from_pv()` to recover Keplerian elements from position and velocity vectors. This is the inverse of `to_pv()`.\n",
"\n",
"The example below uses a test case from Vallado's *Fundamentals of Astrodynamics and Applications*, Example 2-6."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Vallado Example 2-6\n",
"r = np.array([6524.834, 6862.875, 6448.296]) * 1e3 # meters\n",
"v = np.array([4.901327, 5.533756, -1.976341]) * 1e3 # m/s\n",
"\n",
"kep_from_pv = sk.kepler.from_pv(r, v)\n",
"\n",
"print(f\"Semi-major axis: {kep_from_pv.a / 1e3:.1f} km\")\n",
"print(f\"Eccentricity: {kep_from_pv.eccen:.5f}\")\n",
"print(f\"Inclination: {math.degrees(kep_from_pv.inclination):.2f} deg\")\n",
"print(f\"RAAN: {math.degrees(kep_from_pv.raan):.2f} deg\")\n",
"print(f\"Arg of Perigee: {math.degrees(kep_from_pv.w):.2f} deg\")\n",
"print(f\"True Anomaly: {math.degrees(kep_from_pv.nu):.3f} deg\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Round-trip test: convert to PV and back\n",
"pos2, vel2 = kep_from_pv.to_pv()\n",
"kep_roundtrip = sk.kepler.from_pv(pos2, vel2)\n",
"\n",
"print(\"Round-trip consistency check:\")\n",
"print(f\" a difference: {abs(kep_from_pv.a - kep_roundtrip.a):.6e} m\")\n",
"print(f\" e difference: {abs(kep_from_pv.eccen - kep_roundtrip.eccen):.6e}\")\n",
"print(f\" i difference: {abs(kep_from_pv.inclination - kep_roundtrip.inclination):.6e} rad\")\n",
"print(f\" RAAN difference: {abs(kep_from_pv.raan - kep_roundtrip.raan):.6e} rad\")\n",
"print(f\" w difference: {abs(kep_from_pv.w - kep_roundtrip.w):.6e} rad\")\n",
"print(f\" nu difference: {abs(kep_from_pv.nu - kep_roundtrip.nu):.6e} rad\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Orbital Period and Mean Motion\n",
"\n",
"The `period` property returns the orbital period in seconds, and `mean_motion` returns the mean angular rate in radians per second."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(f\"Orbital period: {kep.period:.1f} seconds = {kep.period / 60:.2f} minutes\")\n",
"print(f\"Mean motion: {kep.mean_motion:.6e} rad/s\")\n",
"print(f\"Mean motion: {kep.mean_motion * 86400 / (2 * math.pi):.4f} revs/day\")\n",
"\n",
"# Verify the relationship: period = 2*pi / mean_motion\n",
"print(f\"\\n2*pi / mean_motion = {2 * math.pi / kep.mean_motion:.1f} seconds (should match period)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Keplerian Propagation\n",
"\n",
"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.\n",
"\n",
"Below we propagate an orbit over one full period, sampling positions to visualize the orbit."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Propagate over one orbit, sampling 200 points\n",
"npts = 200\n",
"T = kep.period\n",
"times = np.linspace(0, T, npts)\n",
"\n",
"positions = np.zeros((npts, 3))\n",
"for idx, dt in enumerate(times):\n",
" kep_t = kep.propagate(dt)\n",
" pos_t, _ = kep_t.to_pv()\n",
" positions[idx, :] = pos_t / 1e3 # convert to km\n",
"\n",
"# 3D orbit plot\n",
"fig = plt.figure(figsize=(9, 9))\n",
"ax = fig.add_subplot(111, projection='3d')\n",
"ax.plot(positions[:, 0], positions[:, 1], positions[:, 2], 'b-', linewidth=1.5)\n",
"ax.scatter(*positions[0, :], color='green', s=80, label='Start (perigee)', zorder=5)\n",
"\n",
"# Draw Earth sphere\n",
"u = np.linspace(0, 2 * np.pi, 40)\n",
"v = np.linspace(0, np.pi, 20)\n",
"Re = sk.consts.earth_radius / 1e3\n",
"x_e = Re * np.outer(np.cos(u), np.sin(v))\n",
"y_e = Re * np.outer(np.sin(u), np.sin(v))\n",
"z_e = Re * np.outer(np.ones(np.size(u)), np.cos(v))\n",
"ax.plot_surface(x_e, y_e, z_e, alpha=0.2, color='skyblue')\n",
"\n",
"ax.set_xlabel('X (km)')\n",
"ax.set_ylabel('Y (km)')\n",
"ax.set_zlabel('Z (km)')\n",
"ax.set_title('ISS-like Orbit (One Period)')\n",
"ax.legend()\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Comparing Kepler vs. Numerical Propagation\n",
"\n",
"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.\n",
"\n",
"Here we propagate the same initial orbit using both:\n",
"- **Kepler** (two-body analytical)\n",
"- **High-precision numerical** propagation with full force models\n",
"\n",
"and plot the position difference over time to visualize the perturbation effects."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Define orbit: 550 km circular, 45 deg inclination\n",
"alt = 550e3\n",
"a_circ = sk.consts.earth_radius + alt\n",
"kep0 = sk.kepler(a_circ, 0.001, math.radians(45.0), 0.0, 0.0, 0.0)\n",
"\n",
"# Get initial state vector for numerical propagation\n",
"state0 = np.concatenate(kep0.to_pv())\n",
"\n",
"# Propagation epoch and duration (6 hours)\n",
"epoch = sk.time(2025, 6, 15, 0, 0, 0)\n",
"prop_duration = sk.duration(hours=6)\n",
"\n",
"# Numerical propagation with high-precision force model\n",
"settings = sk.propsettings(\n",
" abs_error=1e-12,\n",
" rel_error=1e-12,\n",
" gravity_degree=10,\n",
")\n",
"settings.precompute_terms(epoch, epoch + prop_duration)\n",
"\n",
"result = sk.propagate(\n",
" state0,\n",
" epoch,\n",
" epoch + prop_duration,\n",
" propsettings=settings,\n",
")\n",
"\n",
"# Sample at 1-minute intervals\n",
"npts = 361\n",
"dt_seconds = np.linspace(0, prop_duration.seconds, npts)\n",
"time_array = [epoch + sk.duration(seconds=float(dt)) for dt in dt_seconds]\n",
"\n",
"# Get positions from both propagators\n",
"pos_kepler = np.zeros((npts, 3))\n",
"pos_numerical = np.zeros((npts, 3))\n",
"\n",
"for idx, (dt, t) in enumerate(zip(dt_seconds, time_array)):\n",
" # Kepler propagation\n",
" kep_t = kep0.propagate(dt)\n",
" pos_k, _ = kep_t.to_pv()\n",
" pos_kepler[idx, :] = pos_k\n",
"\n",
" # Numerical propagation (interpolated)\n",
" state_num = result.interp(t)\n",
" pos_numerical[idx, :] = state_num[0:3]\n",
"\n",
"# Compute position difference\n",
"pos_diff = pos_numerical - pos_kepler\n",
"pos_diff_mag = np.linalg.norm(pos_diff, axis=1)\n",
"\n",
"# Plot\n",
"hours = dt_seconds / 3600\n",
"\n",
"fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8), sharex=True)\n",
"\n",
"ax1.plot(hours, pos_diff[:, 0] / 1e3, label='X')\n",
"ax1.plot(hours, pos_diff[:, 1] / 1e3, label='Y')\n",
"ax1.plot(hours, pos_diff[:, 2] / 1e3, label='Z')\n",
"ax1.set_ylabel('Position Difference (km)')\n",
"ax1.set_title('Numerical - Keplerian Position Difference (per component)')\n",
"ax1.legend()\n",
"ax1.grid(True, alpha=0.3)\n",
"\n",
"ax2.plot(hours, pos_diff_mag / 1e3, 'k-', linewidth=1.5)\n",
"ax2.set_xlabel('Time (hours)')\n",
"ax2.set_ylabel('Position Difference (km)')\n",
"ax2.set_title('Total Position Difference Magnitude')\n",
"ax2.grid(True, alpha=0.3)\n",
"\n",
"plt.tight_layout()\n",
"plt.show()\n",
"\n",
"print(f\"Maximum position difference over 6 hours: {pos_diff_mag[-1]/1e3:.2f} km\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"name": "python",
"version": "3.12.0"
}
},
"nbformat": 4,
"nbformat_minor": 4
}