Orbits and Kepler's third law

Script by Ordinate 2D ✓ · Gravitation and orbits

Gravitation between bodies (s.G, body.gravitates): circular planets give T²/r³ constant, a slow comet gives an ellipse and conserves r·v.

Use in Ordinate 2DOpens in the Script tab: it runs when you press Run. Free, no account needed.

"""Orbits: gravitation between bodies, and Kepler's third law.

Three planets circle a fixed sun under F = G·M·m/r²: the sun and the planets *gravitate* (body.gravitates, the constant s.G), and the
engine works the pull out in every substep, so the orbits close (the scene's own gravity is off).
Each planet is started on a circular orbit, v = √(GM/r), and its period is timed from the angle it sweeps. T²/r³ should be the
same for all of them, 4π²/GM. A fourth planet starts too slowly: its orbit is an ellipse, and r·v is the same at its closest
and its furthest point (angular momentum is conserved).

Runs in the studio (Script tab → Run) and locally:  python3 orbit_kepler.py
"""
import math
import studio as s

SUN = 100.0                                 # kg
s.clear()
s.gravity = (0, 0)
s.G = 0.4                                   # m³/(kg s²): far stronger than the real 6.674e-11, for orbits of metres and seconds
GM = s.G * SUN                              # 40 m³/s², the sun's strength
sun = s.circle(0, 0, 0.25, mass=SUN, static=True, name="Sun", color="#f5a524", collide=False)
sun.gravitates = True
planets = []                                # [body, radius, angle swept, last angle, period]
for name, r, hex_ in [("Mercury", 1.0, "#8b939b"), ("Earth", 1.8, "#4f8ff7"), ("Mars", 2.8, "#e5484d")]:
    # light planets: they pull each other too, as every gravitating body does; at 0.01 kg that is 10⁻⁴ of the sun's pull
    p = s.circle(r, 0, 0.1, mass=0.01, name=name, color=hex_, collide=False, velocity=(0, math.sqrt(GM / r)))
    p.gravitates = True
    planets.append([p, r, 0.0, 0.0, None])
comet = s.circle(3.4, 0, 0.07, mass=0.01, name="Comet", color="#30a46c", collide=False, velocity=(0, 0.55 * math.sqrt(GM / 3.4)))
comet.gravitates = True
r_min = r_max = 3.4
v_at_min = v_at_max = comet.speed


@s.on_step
def watch(t, dt):
    global r_min, r_max, v_at_min, v_at_max
    for entry in planets:                                    # time one revolution from the angle swept
        p, r, swept, last, period = entry
        a = math.atan2(p.y, p.x)
        entry[2] += (a - last + math.pi) % (2 * math.pi) - math.pi
        entry[3] = a
        if period is None and entry[2] >= 2 * math.pi:
            entry[4] = t
            s.record(planet=r, radius=r, period=t, T2_over_r3=t ** 2 / r ** 3, theory=4 * math.pi ** 2 / GM)
    rc = math.hypot(comet.x, comet.y)
    if rc < r_min:
        r_min, v_at_min = rc, comet.speed
    if rc > r_max:
        r_max, v_at_max = rc, comet.speed


@s.on_end
def report():
    print(f"Kepler III: 4π²/GM = {4 * math.pi ** 2 / GM:.4f} s²/m³")
    for p, r, _, _, period in planets:
        if period:
            print(f"{p.name:8s} r = {r:.1f} m  T = {period:.3f} s  T²/r³ = {period ** 2 / r ** 3:.4f}")
    print(f"Comet: closest {r_min:.2f} m at {v_at_min:.2f} m/s, furthest {r_max:.2f} m at {v_at_max:.2f} m/s → r·v = {r_min * v_at_min:.2f} and {r_max * v_at_max:.2f}")


s.run(14)

table gravitation