Kepler’s second law

Why does a planet move fastest when it is closest to the sun? Kepler’s second law says that the line from the sun to the planet sweeps out equal areas in equal times. Near the sun the line is short, so in the same time the planet has to cover a longer arc.

Where the planet is at a given time follows from Kepler’s equation,

\[M = E - e \sin E.\]

The mean anomaly \(M = 2\pi t / T\) grows evenly over the period \(T\). The eccentric anomaly \(E\) gives the planet’s point \((a \cos E, \; b \sin E)\) on the ellipse with semi-axes \(a\) and \(b = a \sqrt{1 - e^2}\), measured from its center, and the sun sits at the focus \((a e, 0)\).

The example is an orbit with \(a = 10\) and eccentricity \(e = 0.6\), more eccentric than any planet’s. It solves Kepler’s equation for twelve equal steps of time, cuts the sectors swept in them from planex2d’s exact ellipse and compares each area with \(\pi a b / 12\).

This page is part of a collection of geometry and design examples using planex2d.

Open this example in the playground, the interactive code editor in the browser, to change the code and run it.

import math

import numpy as np
import planex as px

The orbit

px.Ellipse draws the ellipse as one exact edge, so the orbit’s area is exactly \(\pi a b\):

a, e = 10, 0.6
b = a * math.sqrt(1 - e**2)
sun = (a * e, 0)
orbit = px.Face([px.Ellipse((0, 0), a, b)])
print("area", orbit.mass, "| π a b", math.pi * a * b)

area 251.32741228718317 | π a b 251.32741228718345

Positions at equal times

Kepler’s equation has no closed-form solution for \(E\), but Newton’s method finds it in a few steps from \(E = M\). The function returns the planet’s positions at \(N + 1\) equal steps of time. The first and the last are both at the point closest to the sun:

def positions(N):
    """the planet at N + 1 equal steps of time over one orbit, from Kepler's equation"""
    M = 2 * math.pi * np.arange(N + 1) / N
    E = M.copy()
    for _ in range(30):
        E -= (E - e * np.sin(E) - M) / (1 - e * np.cos(E))
    return np.column_stack([a * np.cos(E), b * np.sin(E)])

The swept sectors

The sector swept between two positions is the part of the orbit inside the wedge from the sun through them. The wedge is a polygon from the sun out along five rays to far beyond the orbit, so that a wedge wider than a right angle near the sun keeps its shape. The Boolean intersection * cuts the sector out of the exact ellipse (Figure 1):

def sector(p, q):
    """the part of the orbit swept from p to q, seen from the sun"""
    t0 = math.atan2(p[1], p[0] - sun[0])
    t1 = math.atan2(q[1], q[0] - sun[0])
    if t1 <= t0:
        t1 += 2 * math.pi
    rays = [(sun[0] + 4 * a * math.cos(t), 4 * a * math.sin(t)) for t in np.linspace(t0, t1, 5)]
    return orbit * px.FromPolyline([sun] + rays, closed=True).Face()


points = positions(12)
sectors = [sector(points[k], points[k + 1]) for k in range(12)]
areas = [s.mass for s in sectors]
print("areas from", min(areas), "to", max(areas), "| π a b / 12", math.pi * a * b / 12)
# px.Plot(px.Compound(sectors))

areas from 20.94395102391355 to 20.943951023948113 | π a b / 12 20.943951023931955

Figure 1: The orbit, the sun (cross) and the planet at twelve equal steps of time

Faster near the sun

The planet’s path in each step is the sector’s curved edge, whose length is the sector’s perimeter less its two straight edges from the sun. Over equal times, the longest path lies next to the sun and the shortest opposite it. As the steps get shorter, the ratio of the two tends to the ratio of the speeds at the closest and the farthest point, \((1 + e)/(1 - e) = 4\), which follows from the planet’s angular momentum \(r v\) being the same at both:

for N in (12, 60, 240):
    points = positions(N)
    paths = [sector(p, q).wires[0].mass - math.dist(sun, p) - math.dist(sun, q) for p, q in zip(points, points[1:])]
    print("N", N, "| longest path / shortest", max(paths) / min(paths), "| (1 + e) / (1 - e)", (1 + e) / (1 - e))

N 12 | longest path / shortest 3.371095640112987 | (1 + e) / (1 - e) 4.0
N 60 | longest path / shortest 3.9563377402297473 | (1 + e) / (1 - e) 4.0
N 240 | longest path / shortest 3.9971375230643336 | (1 + e) / (1 - e) 4.0

See Two faces for * and Splines for px.Ellipse.