Keplerian Orbits¶
Two-body orbital mechanics using np_cvode (SUNDIALS) and np_odeint (ODEPACK). We simulate orbits under inverse-square gravity, verify conservation laws, and visualise orbital trajectories.
Physics¶
A body orbiting a central mass $M$ obeys Newton's law of gravitation:
$$\ddot{\vec{r}} = -\frac{\mu}{|\vec{r}|^3}\vec{r}$$
where $\mu = GM$ is the gravitational parameter. In 2D Cartesian coordinates $(x, y)$ with velocity $(v_x, v_y)$:
$$\dot{x} = v_x, \quad \dot{y} = v_y, \quad \dot{v_x} = -\frac{\mu x}{r^3}, \quad \dot{v_y} = -\frac{\mu y}{r^3}$$
where $r = \sqrt{x^2 + y^2}$.
In [15]:
load("numerics")$
load("numerics-sundials")$
load("ax-plots")$
1. Circular Orbit¶
For a circular orbit of radius $r_0$, the required velocity is $v_c = \sqrt{\mu/r_0}$. The period is $T = 2\pi\sqrt{r_0^3/\mu}$.
We use $\mu = 1$ (normalised units) and $r_0 = 1$, giving $v_c = 1$ and $T = 2\pi$.
In [16]:
/* Gravitational parameter (normalised) */
mu : 1.0$
/* 2D two-body ODE: state = [x, y, vx, vy] */
/* r^3 = (x^2 + y^2)^(3/2) */
r3 : (x^2 + y^2)^(3/2)$
rhs : [vx, vy, -mu*x/r3, -mu*y/r3]$
vars : [t, x, y, vx, vy]$
/* Circular orbit: r0=1, v_circular = sqrt(mu/r0) = 1 */
y0_circ : [1.0, 0.0, 0.0, 1.0]$
T_period : 2 * float(%pi)$
tspan : np_linspace(0.0, T_period, 500)$
orbit : np_cvode(rhs, vars, y0_circ, tspan, 1e-12, 1e-12)$
print("Orbit shape:", np_shape(orbit))$
Orbit shape: [500,5]
In [17]:
/* Plot the circular orbit — should close perfectly */
xc : np_slice(orbit, all, 1)$
yc : np_slice(orbit, all, 2)$
ax_draw2d(
aspect_ratio = true,
line_width = 2,
color = royalblue,
lines(xc, yc),
color = gold,
marker_size = 8,
points([0], [0]), /* Central body */
title = "Circular Orbit (r = 1, T = 2pi)"
)$
In [18]:
/* Verify: position returns to start after one period */
xf : np_ref(orbit, 499, 1)$
yf : np_ref(orbit, 499, 2)$
print("Final position:", xf, yf)$
print("Closure error:", sqrt((xf - 1.0)^2 + yf^2))$
Final position: 1.0000000000630314 -2.559723405459735e-10 Closure error: 2.6361864730073745e-10
2. Elliptical Orbit¶
An elliptical orbit occurs when the velocity at periapsis exceeds circular velocity but is below escape velocity. For a body at $r_0 = 1$ with $v_0 = 0.8$, the orbit is elliptical.
The specific orbital energy $\varepsilon = v^2/2 - \mu/r$ is conserved. The semi-major axis is $a = -\mu/(2\varepsilon)$.
In [19]:
/* Elliptical orbit: less than circular velocity */
v0_ellip : 0.8$
y0_ellip : [1.0, 0.0, 0.0, v0_ellip]$
/* Orbital energy and semi-major axis */
epsilon : v0_ellip^2/2 - mu/1.0$
a_sma : -mu/(2*epsilon)$
T_ellip : 2 * float(%pi) * sqrt(a_sma^3/mu)$
print("Semi-major axis:", a_sma)$
print("Period:", T_ellip)$
tspan_e : np_linspace(0.0, T_ellip, 500)$
orbit_e : np_cvode(rhs, vars, y0_ellip, tspan_e, 1e-12, 1e-12, bdf)$
xe : np_slice(orbit_e, all, 1)$
ye : np_slice(orbit_e, all, 2)$
ax_draw2d(
aspect_ratio = true,
line_width = 2,
color = tomato,
lines(xe, ye),
color = gold,
marker_size = 8,
points([0], [0]),
title = sconcat("Elliptical Orbit (a = ", a_sma, ")")
)$
Semi-major axis: 0.7352941176470589 Period: 3.9616080528290403
3. Conservation Laws¶
In a two-body problem, three quantities are conserved:
- Specific energy: $\varepsilon = \frac{1}{2}|\vec{v}|^2 - \frac{\mu}{|\vec{r}|}$
- Specific angular momentum: $h = \vec{r} \times \vec{v}$ (scalar in 2D: $h = x v_y - y v_x$)
- Laplace-Runge-Lenz vector (direction of periapsis)
Let's verify energy and angular momentum conservation for the elliptical orbit.
In [20]:
/* Extract state columns */
t_col : np_slice(orbit_e, all, 0)$
xe_col : np_slice(orbit_e, all, 1)$
ye_col : np_slice(orbit_e, all, 2)$
vxe_col : np_slice(orbit_e, all, 3)$
vye_col : np_slice(orbit_e, all, 4)$
/* Compute r = sqrt(x^2 + y^2) */
r_col : np_sqrt(np_add(np_mul(xe_col, xe_col), np_mul(ye_col, ye_col)))$
/* Specific kinetic energy: v^2/2 */
KE : np_scale(0.5, np_add(np_mul(vxe_col, vxe_col), np_mul(vye_col, vye_col)))$
/* Specific potential energy: -mu/r */
PE : np_scale(-mu, np_div(np_ones(np_shape(r_col)), r_col))$
/* Total energy */
energy : np_add(KE, PE)$
/* Angular momentum (2D): h = x*vy - y*vx */
h_col : np_sub(np_mul(xe_col, vye_col), np_mul(ye_col, vxe_col))$
print("Energy range:", np_min(energy), "to", np_max(energy))$
print("Energy drift:", np_max(energy) - np_min(energy))$
print("Angular momentum range:", np_min(h_col), "to", np_max(h_col))$
print("h drift:", np_max(h_col) - np_min(h_col))$
Energy range: -0.6800000002161121 to -0.6799999994748769 Energy drift: 7.4123518434277e-10 Angular momentum range: 0.7999999998419117 to 0.8000000001747377 h drift: 3.328259889912033e-10
In [21]:
/* Plot energy and angular momentum over time */
ax_draw2d(
line_width = 2,
color = royalblue,
name = "Energy",
lines(t_col, energy),
color = tomato,
name = "Ang. momentum",
lines(t_col, h_col),
title = "Conservation Laws (Elliptical Orbit)",
xlabel = "t",
ylabel = "Value"
)$
4. Kepler's Equation¶
Kepler's equation relates the mean anomaly $M$ (which grows linearly with time) to the eccentric anomaly $E$:
$$M = E - e \sin E$$
where $e$ is the eccentricity. This transcendental equation must be solved numerically for $E$ given $M$. We use np_fsolve to solve it.
In [22]:
load("numerics-optimize")$
/* Eccentricity of our elliptical orbit */
/* e = 1 - r_periapsis/a = 1 - r0/a (launched from periapsis) */
e_ecc : 1.0 - 1.0/a_sma$
print("Eccentricity:", e_ecc)$
/* Solve Kepler's equation M = E - e*sin(E) for a range of M values */
M_vals : np_linspace(0.0, 2*float(%pi), 100)$
E_vals : np_zeros(100)$
for i : 0 thru 99 do (
M_i : np_ref(M_vals, i),
/* np_fsolve: solve E - e*sin(E) - M = 0 */
sol : first(np_fsolve([E_var - e_ecc*sin(E_var) - M_i], [E_var], [M_i])),
np_set(E_vals, i, np_ref(sol, 0))
)$
/* Convert E to true anomaly: tan(theta/2) = sqrt((1+e)/(1-e)) * tan(E/2) */
half_theta : np_atan2(
np_scale(sqrt((1+e_ecc)/(1-e_ecc)), np_sin(np_scale(0.5, E_vals))),
np_cos(np_scale(0.5, E_vals))
)$
theta_vals : np_scale(2.0, half_theta)$
/* Radius from true anomaly: r = a(1-e^2)/(1 + e*cos(theta)) */
p_param : a_sma * (1 - e_ecc^2)$
r_kepler : np_div(
np_full(100, p_param),
np_add(np_ones(100), np_scale(e_ecc, np_cos(theta_vals)))
)$
/* Convert polar to Cartesian */
xk : np_mul(r_kepler, np_cos(theta_vals))$
yk : np_mul(r_kepler, np_sin(theta_vals))$
Eccentricity: -0.3599999999999999
In [23]:
/* Overlay Kepler solution (analytic) on numerical integration */
ax_draw2d(
aspect_ratio = true,
line_width = 2,
color = tomato,
name = "Numerical (CVODE)",
lines(xe, ye),
line_width = 1,
dash = "dash",
color = black,
name = "Kepler equation",
lines(xk, yk),
color = gold,
marker_size = 8,
name = "",
points([0], [0]),
title = "Numerical vs Kepler Equation"
)$
5. Multiple Orbits — Energy Comparison¶
Varying the initial velocity produces circular, elliptical, parabolic, and hyperbolic trajectories. The escape velocity is $v_{esc} = \sqrt{2\mu/r_0}$.
In [24]:
v_esc : sqrt(2*mu/1.0)$
print("Escape velocity:", v_esc)$
/* Sub-circular (elliptical, low energy) */
v_list : [0.6, 0.8, 1.0, 1.2, v_esc + 0.01]$
colors : [steelblue, tomato, seagreen, darkorange, purple]$
orbit_names : ["v=0.6 (ellip)", "v=0.8 (ellip)", "v=1.0 (circ)", "v=1.2 (ellip)", "v>v_esc (hyper)"]$
objects : []$
for k : 1 thru 5 do (
v0k : v_list[k],
y0k : [1.0, 0.0, 0.0, v0k],
/* Use shorter integration for hyperbolic */
tf : if k = 5 then 6.0 else 4*float(%pi),
tsp : np_linspace(0.0, tf, 800),
orbk : np_cvode(rhs, vars, y0k, tsp, 1e-12, 1e-12, bdf),
xk : np_slice(orbk, all, 1),
yk : np_slice(orbk, all, 2),
objects : append(objects, [
color = colors[k],
name = orbit_names[k],
line_width = 2,
lines(xk, yk)
])
)$
apply(ax_draw2d, append([
aspect_ratio = true,
color = gold, marker_size = 8, name = "",
points([0], [0]),
title = "Orbits at Different Energies"
], objects))$
Escape velocity: 1.4142135623730951
6. 3D Orbit — Inclined Plane¶
Extend to 3D with an orbit inclined at 30 degrees to the reference plane. The equations are the same but in $(x, y, z)$ with initial velocity rotated out of the xy-plane.
In [25]:
/* 3D two-body ODE: state = [x, y, z, vx, vy, vz] */
r3_3d : (x3^2 + y3^2 + z3^2)^(3/2)$
rhs3d : [vx3, vy3, vz3, -mu*x3/r3_3d, -mu*y3/r3_3d, -mu*z3/r3_3d]$
vars3d : [t, x3, y3, z3, vx3, vy3, vz3]$
/* Inclined circular orbit: rotate velocity vector by 30 deg about x-axis */
inc : float(%pi/6)$ /* 30 degrees */
vy0 : cos(inc)$
vz0 : sin(inc)$
y0_3d : [1.0, 0.0, 0.0, 0.0, vy0, vz0]$
tspan3d : np_linspace(0.0, 2*float(%pi), 500)$
orbit3d : np_cvode(rhs3d, vars3d, y0_3d, tspan3d, 1e-12, 1e-12)$
x3c : np_slice(orbit3d, all, 1)$
y3c : np_slice(orbit3d, all, 2)$
z3c : np_slice(orbit3d, all, 3)$
/* Angular momentum vector L = r x v (should be constant) */
r_vec0 : ndarray([np_ref(orbit3d, 0, 1), np_ref(orbit3d, 0, 2), np_ref(orbit3d, 0, 3)], [3])$
v_vec0 : ndarray([np_ref(orbit3d, 0, 4), np_ref(orbit3d, 0, 5), np_ref(orbit3d, 0, 6)], [3])$
L_vec : np_cross(r_vec0, v_vec0)$
print("Angular momentum vector:", np_to_list(L_vec))$
print("|L| =", np_norm(L_vec))$
Angular momentum vector: [0.0,-0.49999999999999994,0.8660254037844387] |L| = 1.0
In [26]:
/* Project 3D orbit onto xy and xz planes */
ax_draw2d(
aspect_ratio = true,
line_width = 2,
color = royalblue,
name = "xy projection",
lines(x3c, y3c),
color = tomato,
name = "xz projection",
lines(x3c, z3c),
color = gold,
marker_size = 8,
name = "",
points([0], [0]),
title = "3D Orbit Projections (i = 30 deg)"
)$
7. Hohmann Transfer¶
The Hohmann transfer is the most fuel-efficient two-impulse manoeuvre between coplanar circular orbits. Two burns:
- At radius $r_1$: increase velocity to enter an elliptical transfer orbit
- At radius $r_2$: increase velocity again to circularise
Transfer orbit semi-major axis: $a_t = (r_1 + r_2)/2$
In [27]:
/* Transfer from r1=1 to r2=3 */
r1 : 1.0$ r2 : 3.0$
/* Circular velocities */
vc1 : sqrt(mu/r1)$
vc2 : sqrt(mu/r2)$
/* Transfer orbit */
a_t : (r1 + r2) / 2$
v_depart : sqrt(mu * (2/r1 - 1/a_t))$ /* velocity at r1 on transfer ellipse */
v_arrive : sqrt(mu * (2/r2 - 1/a_t))$ /* velocity at r2 on transfer ellipse */
dv1 : v_depart - vc1$
dv2 : vc2 - v_arrive$
print("Delta-v 1 (departure):", dv1)$
print("Delta-v 2 (arrival):", dv2)$
print("Total delta-v:", dv1 + dv2)$
/* Simulate: initial orbit, transfer, final orbit */
T_transfer : float(%pi) * sqrt(a_t^3/mu)$
/* Initial circular orbit (partial) */
tspan1 : np_linspace(0.0, float(%pi), 200)$
orb1 : np_cvode(rhs, vars, [r1, 0.0, 0.0, vc1], tspan1, 1e-12, 1e-12)$
/* Transfer ellipse */
tspan_t : np_linspace(0.0, T_transfer, 200)$
orb_t : np_cvode(rhs, vars, [r1, 0.0, 0.0, v_depart], tspan_t, 1e-12, 1e-12)$
/* Final circular orbit (partial) */
tspan2 : np_linspace(0.0, float(%pi), 200)$
orb2 : np_cvode(rhs, vars, [-r2, 0.0, 0.0, -vc2], tspan2, 1e-12, 1e-12)$
Delta-v 1 (departure): 0.22474487139158894 Delta-v 2 (arrival): 0.16910197872576277 Total delta-v: 0.3938468501173517
In [28]:
ax_draw2d(
aspect_ratio = true,
line_width = 2,
color = steelblue,
name = "Initial orbit (r=1)",
lines(np_slice(orb1, all, 1), np_slice(orb1, all, 2)),
color = tomato,
dash = "dash",
name = "Hohmann transfer",
lines(np_slice(orb_t, all, 1), np_slice(orb_t, all, 2)),
dash = "solid",
color = seagreen,
name = "Final orbit (r=3)",
lines(np_slice(orb2, all, 1), np_slice(orb2, all, 2)),
color = gold,
marker_size = 8,
name = "",
points([0], [0]),
color = red,
marker_size = 5,
name = "Burns",
points([r1, -r2], [0, 0]),
title = "Hohmann Transfer Orbit (r=1 to r=3)"
)$
Summary¶
| Concept | Technique |
|---|---|
| Two-body gravity | np_cvode with inverse-square force law |
| Conservation verification | np_dot, np_cross, energy/momentum computation |
| Kepler's equation | np_fsolve for transcendental root finding |
| Orbital classification | Varying initial velocity: elliptic, circular, hyperbolic |
| 3D orbits | Extend state vector, use np_cross for angular momentum |
| Hohmann transfer | Two-impulse delta-v calculation + trajectory simulation |