A field guide
Orbital mechanics.
Throw a ball. It falls and hits the ground. Throw it harder, and it travels further before it falls. Now imagine throwing it hard enough that, as it falls toward the ground, the ground curves away underneath it just as fast. The ball never lands. It just keeps falling, forever, in a circle.
That is an orbit. Once you have that picture in your head, almost everything that follows is bookkeeping.
A satellite going around in a circle. The arrow toward the planet is gravity. The arrow tangent to the path is its velocity. They never quite line up, and that is the whole reason it stays up there.
II · What an orbit really is
Newton imagined a cannon on a very tall mountain. Fire a ball horizontally with a small charge and it arcs out a few kilometres before it lands. Pack in more gunpowder and the arc gets longer. Keep packing. At some point the ball's curving path matches the curve of the planet itself. The ball is now in orbit. It is still falling. It is just falling and missing.
Now I want you to notice something. The shape of the path depends on exactly one number: how fast you launched the thing. We measure that speed in units of the local circular speed $v_c = \sqrt{GM/r}$, which is a fancy way of saying "the speed at which falling and missing balance out." Below $v_c$, the ball arcs down and hits the surface. At exactly $v_c$, you get a perfect circle. Speed it up a bit and the orbit stretches into an ellipse. At $\sqrt{2}\,v_c$ something special happens: the ball has enough energy to climb away and never come back. We call this the escape speed. Above it, the trajectory opens out into a hyperbola.
Circle, ellipse, parabola, hyperbola. These are the same four curves the Greeks got by slicing a cone. Under one inverse-square pull, those are all the trajectories there are. Drag the speed slider and watch the orbit hop between them. The simulator is honestly integrating Newton's equations every frame, so what you see is the math.
III · Three rules from a stack of notebooks
Between 1609 and 1619, Johannes Kepler took twenty years of eyestrain (most of it Tycho Brahe's) and squeezed it into three sentences. He had no calculus. He had no theory of gravity. He had columns of numbers and stubbornness.
First, planets do not move in circles. They move on ellipses, with the Sun parked at one of the two foci. Second, draw a line from the Sun to the planet and let it sweep along as the planet moves. The area that line covers per second is the same whether the planet is near the Sun or far away. Third, square the orbital period, divide by the cube of the semi-major axis, and you get the same number for every planet. In symbols: $T^2 \propto a^3$.
The second law sounds magical until you realise it is just angular momentum hiding in geometric clothes. The little triangle the line sweeps in a moment of time has area $\tfrac12 |\mathbf r \times \mathbf v|\,dt$. For any radial force, $\mathbf r \times \mathbf v$ is constant. So the sweep rate is constant. The planet hustles through a wide angle when it is close to the Sun and crawls through a narrow one when it is far. The area comes out the same.
Try it on the left. Drag the planet around. The little box reports the sweep rate, and it does not flinch. On the right, every dot is a real planet on a log-log plot of $T$ against $a$. They all sit on one straight line of slope $3/2$. Slide your imaginary planet out to forty AU and the line tells you exactly how long its year would be.
IV · One equation that pays for itself
Energy is conserved. That is the whole derivation. Write down the kinetic energy plus the gravitational potential, set it equal to its value at any other point of the orbit, and rearrange. You get
$$v^2 = \mu\left(\frac{2}{r} - \frac{1}{a}\right).$$
Leibniz called kinetic energy vis viva, the "living force," and the name has stuck to this little equation ever since. It tells you how fast a body is moving at radius $r$ on an orbit of semi-major axis $a$. That is it. No time. No angle. Two numbers in, one speed out.
Pin $r = a$ and you recover the circular speed, $\sqrt{\mu/a}$. At the closest point, $r = a(1-e)$, the body is moving fastest; at the farthest point, $r = a(1+e)$, it is at its slowest. And here is the practical thing. Every interplanetary mission has a fuel budget, and that fuel budget is just a sum of differences of vis-viva speeds at the various $(r, a)$ pairs you pass through. That is what mission designers are doing all day.
V · The naive integrator, and why it lies
Newton's equations almost never have a tidy formula for the answer. The two-body problem does. Almost nothing else does. So if you want to know where the planet is at time $t$, you have to march forward in little steps and add things up. This is called numerical integration, and the simplest version of it is called forward Euler.
The recipe is two lines. Update the position using the current velocity. Update the velocity using the current acceleration. That is the whole algorithm.
$$\mathbf r_{n+1} = \mathbf r_n + \Delta t\, \mathbf v_n, \qquad \mathbf v_{n+1} = \mathbf v_n + \Delta t\, \mathbf a(\mathbf r_n).$$
Run it on a perfect circle of radius one with a time step of $0.01$, and you get a slow disaster. After five trips around, the radius has grown by 45%. The energy has grown too. Halve the time step and the drift halves, but it never reaches zero. This is the strange part. The problem is not that the steps are too big. The problem is the algorithm. Forward Euler reads the force at the old position and applies it to make the new velocity. The two are slightly out of phase, and the result is that every step pumps a tiny bit of extra energy into the orbit. Backward Euler, the obvious cousin, makes the opposite mistake and spirals the planet inward instead. There is no time step at which Euler simply works.
VI · The same orbit, done right
Suppose you do not push the velocity forward in one big step. Suppose, instead, you push it forward by half a step, use that half-step velocity to move the position the full way, and then push the velocity the other half-step using the force at the new place. Three lines instead of two. Same number of multiplies. Same number of additions.
$$\mathbf v_{n+\tfrac12} = \mathbf v_n + \tfrac{\Delta t}{2}\,\mathbf a_n, \quad \mathbf r_{n+1} = \mathbf r_n + \Delta t\,\mathbf v_{n+\tfrac12}, \quad \mathbf v_{n+1} = \mathbf v_{n+\tfrac12} + \tfrac{\Delta t}{2}\,\mathbf a_{n+1}.$$
This is called the velocity Verlet method, or leapfrog. It looks like a small rearrangement, and it is. But something quietly important happens. Verlet preserves a structure of physics called a symplectic form. I do not have a perfect intuition for what that means, but the consequence is sharp. The energy does not drift. It wobbles a tiny amount, back and forth, around the true value, and it stays in that wobble for as long as you want to run the simulation.
Below, the same orbit as the last section. Same time step. Same starting point. With Verlet, the radius stays at one. The energy chart hums in a band of width about $10^{-9}$ centred on $-\tfrac12$. Toggle to Euler and watch the same time step turn into a slow shipwreck.
VII · A solar system in fifty lines
Once you have the Verlet step, going from two bodies to many is almost free. Every body pulls on every other body with an inverse-square force, you sum those pulls to get each body's acceleration, and then you take the same little three-line step. Use astronomical units, years, and solar masses, and the gravitational constant becomes $G = 4\pi^2$, which is a pleasing little gift of unit conversion.
Below is the Sun plus the inner five planets, running live. Mercury laps Earth about four times each year. Jupiter takes twelve Earth years to drag itself once around the loop. The energy of the whole arrangement does not stay exactly constant. It wobbles, the way Verlet's energy always does. But it does not drift. You could leave this thing running overnight and the wobble would be the same width in the morning.
VIII · Three bodies, and the place predictions die
Two bodies under gravity have a closed-form solution. Newton found it. Add a third body and the whole thing falls apart. There is no formula. There is no series. There is only the integrator, and the integrator gives different answers if you nudge the initial conditions by a fraction of a percent. Henri Poincaré proved this around 1887 and accidentally invented chaos theory while he was trying to win a prize from the King of Sweden.
There is one strange exception. Three equal masses can chase each other along a single figure-8 curve, like cars on a track. Cris Moore stumbled into it numerically in 1993. Alain Chenciner and Richard Montgomery proved in 2000 that it really exists. With exactly the right starting conditions, this dance repeats forever.
Now I want to show you what chaos looks like. Take those exact starting conditions, copy them, and nudge the first body's $x$-coordinate by one part in a million. Run both versions side-by-side. Watch what happens to the distance between the two "first bodies."
A stable system would let that distance grow steadily, predictably. Chaos does something different. The separation jumps, plateaus, shrinks back partway, jumps again. The number you should look at on the chart is the count of steps where the two trajectories briefly move toward each other instead of apart. In a chaotic system that number is huge. In a non-chaotic one it is zero. That is the fingerprint.
IX · Five quiet places
Stand in a frame that rotates along with the Earth as it goes around the Sun. In this frame the Earth holds still. Now ask: where could a spacecraft sit so that gravity from the Sun, gravity from the Earth, and the centrifugal pull of the rotating frame all cancel each other out? Solve $\nabla U = 0$ for the effective potential
$$U(x,y) = \tfrac{1}{2}(x^2 + y^2) + \frac{1-\mu}{r_1} + \frac{\mu}{r_2}$$
and you get five answers. Three of them, L1, L2, and L3, sit on the line through the Sun and the Earth. Two of them, L4 and L5, sit at the corners of equilateral triangles drawn with the Sun and the Earth as the other two vertices. Five quiet places where the forces balance.
Now you might think that "quiet" means "stable." It does not. L1, L2, and L3 are saddle points. A spacecraft placed there is balanced the way a pencil is balanced on its tip. The James Webb Space Telescope sits at L2 and has to fire its thrusters every few weeks to stay there. L4 and L5 are even stranger. They are local maxima of the potential, which sounds disastrous, but the Coriolis force from the rotating frame steers things into orbits around them as long as the mass ratio of the two big bodies is below a critical value of about $0.0385$. That is why the Sun-Jupiter L4 and L5 hold the Trojan asteroids in clouds that have been stable for billions of years.
Click an L-point on the map to see what spacecraft live there. Click anywhere else to drop a test particle and watch the rotating-frame physics carry it around.
A small honesty note. This widget uses $\mu \approx 0.012$, which is more like the Earth-Moon ratio. The real Sun-Earth $\mu$ is about three parts in a million, so L1 and L2 would sit right on top of the Earth and you would not be able to see anything.
X · How you actually get to Mars
With vis-viva from chapter IV and the integrators from chapters V and VI, you have enough to design a real mission. Two tricks do most of the work.
The first trick is the Hohmann transfer, written down by Walter Hohmann in 1925. Suppose you are in a circular orbit at one radius and want to be in a circular orbit at a bigger radius. Fire your engine forward, briefly. This puts you on an ellipse whose nearest point is your old orbit and whose farthest point touches the new one. Coast for half a period. When you arrive at the far end, fire forward again. Now you are circular at the target radius. Two burns. The total fuel cost is just the vis-viva difference at the two ends. For Earth to Mars, that adds up to about $5.6$ km/s, and the coast takes about eight and a half months.
The second trick is the gravity assist. This one is harder to believe at first. A spacecraft flies past a planet on a hyperbola. In the planet's own frame, it comes in fast, swings around, and leaves at the same speed in a new direction. No energy gained, no energy lost. But the planet itself is moving through the solar system. In the Sun's frame, you have to add the planet's velocity to the spacecraft's. And because the new direction is different, that addition lands you somewhere different in the Sun's frame too. You can gain up to twice $v_\infty$ of heliocentric speed per encounter, depending on which side of the planet you fly. Voyager 2 used Jupiter, then Saturn, then Uranus this way, and reached Neptune on a fuel budget that should have stranded it just past the asteroid belt.
XI · A workshop Four experiments. Each one settles a question the prose above can only describe.
Up to here we have been watching. Now we will build the thing that does the watching. Not because anyone needs to learn Python, but because four facts about orbits really only land when you watch the arithmetic produce them for yourself.
First: how badly does the simplest possible orbit-integrator go wrong, and can you fix it by making the time step smaller? Second: what is the tiniest possible change that fixes it, and why does that change fix it so completely? Third: if you hand a working integrator the real masses and distances of the solar system, do Kepler's laws come back out, without ever telling the program about them? Fourth: take that same well-behaved integrator and feed it a three-body orbit that is famously stable. Nudge one body by a millionth of an inch. Watch what happens to that nudge.
You can read about these four facts. It is much harder to actually believe them without watching a circle become a spiral on your own screen, and then refusing to.
The naive integrator, in two lines
The simplest thing that could possibly work: update the position from the velocity, then update the velocity from the acceleration. Two lines of arithmetic per step. The mathematical content is what you already know — Newton's $\ddot{\vec r} = -GM\,\vec r/r^3$, evaluated and rolled forward in tiny ticks of time.
Run it on a perfect circle. The circle should stay a circle. The question is whether it does. The answer, when you watch the integrator actually run, is the point of this whole workshop — if the simplest reasonable thing is wrong, the rest of orbital mechanics has a real problem.
fragmentimport numpy as np
GM = 1.0
def acceleration(r):
return -GM * r / np.linalg.norm(r)**3
def step_euler(r, v, dt):
a = acceleration(r)
r_new = r + dt * v # position uses OLD velocity
v_new = v + dt * a # velocity uses OLD acceleration
return r_new, v_new
# Circular orbit: r=(1,0), v=(0,1), dt=0.01, 5 periods.
r = np.array([1.0, 0.0]); v = np.array([0.0, 1.0])
for period in range(5 * int(2*np.pi / 0.01)):
r, v = step_euler(r, v, 0.01)
# Expect |r|=1.0 forever. Watch what Euler actually does.
Five trips, and the radius has gained almost half its size. The integrator is broken. Halving the time step halves the drift, so it never goes away.
The same orbit, fixed
Same orbit. Same time step. Same arithmetic precision. The only change is the order of operations: instead of one velocity push per step, take a half-push, then the position update, then another half-push using the acceleration evaluated at the new position. One extra line of code.
The claim, which is hard to fully accept the first time you hear it, is that this small reordering keeps the orbit honest forever. Not just for a few more periods than Fragment I lasted. Forever. The energy number on the screen should refuse to drift, no matter how long you let the loop run.
fragmentdef step_verlet(r, v, dt, a):
v_half = v + 0.5 * dt * a
r_new = r + dt * v_half
a_new = acceleration(r_new)
v_new = v_half + 0.5 * dt * a_new
return r_new, v_new, a_new
r = np.array([1.0, 0.0]); v = np.array([0.0, 1.0]); a = acceleration(r)
for step in range(100 * int(2*np.pi / 0.01)):
r, v, a = step_verlet(r, v, 0.01, a)
# Print E every 10 periods. The number does not move.
One hundred orbits in, the energy has barely moved. The wobble is down at $10^{-9}$. Two extra lines of arithmetic bought you that.
Six bodies, real masses
Now apply the same Verlet step to six bodies. The only new piece of code is a small inner loop that sums up every body's pull on every other body. Twelve more lines of NumPy.
The integrator is not told about Kepler. It does not know that $T^2$ should equal $a^3$. It just computes accelerations from Newton's law and steps forward in time. And yet, when you watch Mercury go around forty-nine times while Jupiter completes one circuit, you are seeing Kepler's third law happen on the screen, derived in real time from the inverse-square. The law is a consequence, not a postulate.
fragmentG = 4 * np.pi**2 # AU^3 / (M_sun * year^2)
def accelerations(rs, ms):
accs = [np.zeros(2) for _ in ms]
for i in range(len(ms)):
for j in range(len(ms)):
if i == j: continue
d = rs[j] - rs[i]
accs[i] += G * ms[j] * d / np.linalg.norm(d)**3
return accs
bodies = [
("Sun", 1.0, [0.0, 0.0], [0.0, 0.0]),
("Mercury", 1.66e-7, [0.39, 0.0], [0.0, 10.07]),
("Venus", 2.45e-6, [0.72, 0.0], [0.0, 7.39]),
("Earth", 3.00e-6, [1.00, 0.0], [0.0, 6.28]),
("Mars", 3.21e-7, [1.52, 0.0], [0.0, 5.08]),
("Jupiter", 9.55e-4, [5.20, 0.0], [0.0, 2.76]),
]
# 12 years, dt = 0.001 yr. Count each body's loops.
Mercury makes about 49 trips, Earth makes 12, Jupiter barely finishes one. The 3:2 slope of Kepler's third law is right there in the numbers, with no theorem in sight.
The fingerprint of chaos
The figure-eight is a real, beautiful, exactly periodic three-body orbit, discovered numerically in 1993 and proved to exist rigorously in 2000. Three equal masses chase each other along one curve, forever. Energy is conserved. Angular momentum is conserved. Nothing is approximate.
Now run it twice in parallel. In one copy, nudge one body's starting position by $10^{-6}$. The integrator below is the same Verlet that held the two-body orbit steady to a part in a billion, so any difference between the two runs is not a numerical artifact — it is what the physics itself does to that nudge.
Count the timesteps where the separation between the two runs is shrinking rather than growing. In a stable system that count is zero; the nudge just spreads, like a yardstick stretching. In a chaotic one the count is huge, because the separation bounces back and forth as it grows. The percentage you'll see is the actual signature of chaos, more honest than any single number for the size of the divergence.
fragment# Chenciner-Montgomery figure-eight (2000):
rs0 = [
np.array([ 0.97000436, -0.24308753]),
np.array([-0.97000436, 0.24308753]),
np.array([ 0.0, 0.0]),
]
v3 = np.array([-0.93240737, -0.86473146])
vs0 = [-0.5*v3, -0.5*v3, v3]
# Perturb body 0 by +1e-6 in x, run both, compare.
rs0_pert = [r.copy() for r in rs0]
rs0_pert[0] += np.array([1e-6, 0.0])
traj_a = run(rs0, vs0, ms=[1,1,1], dt=0.001, n_steps=60000)
traj_b = run(rs0_pert, vs0, ms=[1,1,1], dt=0.001, n_steps=60000)
sep = np.linalg.norm(traj_a - traj_b, axis=1)
print(f"3-body decreased in {(np.diff(sep) < 0).sum()} of {len(sep)-1} steps")
# Compare to a 2-body control: predicted to decrease in 0 steps.
A chaotic system does not just diverge faster. It diverges in a different shape. Watch the percentage and feel the difference.
What the four experiments add up to
You ran the four experiments. Here is what each one actually showed, and where its conclusion stops being safe.
- A wrong choice in two lines of arithmetic can ruin a perfect circle in five revolutions. Smaller time steps slow the disease without curing it.
- A two-line reordering of the same arithmetic holds the orbit steady to a part in a billion for as long as you let the loop run. The structure of the algorithm matters more than the size of $\Delta t$.
- Feed the fixed integrator the actual masses and distances of the planets, and Kepler's third law shows up in the orbit counts on its own. Newton's gravity contains Kepler; it does not need to be told.
- The same integrator on a three-body orbit reveals chaos as a different shape of divergence, not just a bigger one. Two-body separations grow monotonically; three-body separations oscillate.
- None of this handles a close approach. The first tight flyby would eat the time step. Production codes use adaptive methods like IAS15 or Bulirsch-Stoer for exactly this reason.
- There are no general-relativistic corrections here. Mercury's famous 43 arcsec-per-century perihelion advance is not in this arithmetic.
- No atmospheric drag, no radiation pressure, no $J_2$ oblateness. A real low-Earth-orbit satellite would need all three of those before the model means anything.
If you want to keep poking: halve $\Delta t$ in the first experiment and watch the drift halve right back. Run the second one for ten thousand periods. Replace the figure-eight in the fourth experiment with the Pythagorean problem (three masses at the corners of a 3-4-5 triangle, all at rest) and watch the same Verlet fail at the first close encounter. That failure is the reason production codes look the way they do.
The integrator is the foundation. Everything else — ephemerides, mission design, long-term stability arguments — rides on top of getting this one little step right.
XII · Anatomy — how the parts fit
Three plates of the same machinery, taken apart on the table.