N-Body Orbits, Hohmann Transfers, and Galactic Rotation Curves#
This tutorial covers three scales of gravitational dynamics with
physicskit.astro.nbody, physicskit.astro.orbital_mechanics,
and physicskit.astro.galactic_dynamics: a symplectically-integrated
two-body orbit, an interplanetary transfer maneuver, and the flat
rotation curve that first pointed astronomers toward dark matter.
A symplectically-integrated circular orbit#
NBodySystem advances any number of
gravitating bodies with a leapfrog (symplectic) integrator – the same
family of integrator prized elsewhere in this library for long-term
orbital stability. Two equal masses on a circular orbit around their
common center of mass is the simplest nontrivial test:
import numpy as np
from physicskit.astro.nbody import NBodySystem
G, m, d = 1.0, 1.0, 2.0
v = np.sqrt(G * m / (2 * d)) # circular-orbit speed for this separation
positions = np.array([[d / 2, 0, 0], [-d / 2, 0, 0]])
velocities = np.array([[0, v, 0], [0, -v, 0]])
system = NBodySystem(positions, velocities, masses=np.array([m, m]), G=G)
E0 = system.total_energy()
period = 2 * np.pi * (d / 2) / v
system.simulate(dt=period / 200, n_steps=2000) # 10 full orbits
E1 = system.total_energy()
print((E1 - E0) / E0) # ~1.6e-12 -- energy conserved to near machine precision
Ten complete orbits leave the total energy essentially unchanged – exactly the long-term stability a symplectic integrator is chosen for, where a naive (non-symplectic) integrator would accumulate a steady energy drift over the same run.
A Hohmann transfer between two orbits#
hohmann_transfer() computes the
two burns of the most fuel-efficient two-impulse transfer between
circular orbits – the maneuver actually used to send spacecraft between
planets. In units where 1 AU is the length unit, 1 year is the time unit,
and the Sun’s gravitational parameter is \(\mu_\odot=4\pi^2\), an
Earth-to-Mars transfer (\(r_1=1\,{\rm AU}\),
\(r_2=1.524\,{\rm AU}\)) looks like:
from physicskit.astro.orbital_mechanics import hohmann_transfer
mu_sun = 4 * np.pi ** 2
dv1, dv2, transfer_time = hohmann_transfer(r1=1.0, r2=1.524, mu=mu_sun)
print(dv1, dv2, transfer_time)
# 0.6215 0.5590 0.7089 (speed changes in AU/yr, time in years)
A transfer time of roughly 0.71 years – about 8.5 months – matches the real-world figure quoted for Hohmann-transfer trajectories to Mars.
Galactic rotation curves and dark matter#
circular_velocity(), given any
enclosed-mass profile, gives the orbital speed a test particle needs to
stay on a circular orbit at each radius. For an ordinary (luminous, point
or centrally-concentrated) mass distribution this speed falls off as
\(1/\sqrt r\) at large radius – but for a Navarro-Frenk-White dark
matter halo, nfw_enclosed_mass()
keeps growing (logarithmically) with radius, and the two effects very
nearly cancel:
from physicskit.astro.galactic_dynamics import nfw_enclosed_mass, circular_velocity
rho_s, r_s = 1.0, 1.0
r = np.array([2, 5, 10, 20, 30]) * r_s
v = circular_velocity(r, lambda rr: nfw_enclosed_mass(rr, rho_s, r_s))
print(v)
# [1.6474 1.552 1.3678 1.1465 1.0164]
Beyond a few scale radii, the rotation curve is nearly flat (the speed at 30 \(r_s\) is still within about 16% of its value at 5 \(r_s\), rather than falling off by a large factor as Keplerian rotation around a fixed point mass would) – exactly the puzzle Rubin and Ford’s measurements presented in the 1970s, resolved by exactly this kind of extended, mostly-invisible mass distribution wrapped around every galaxy.