R_earth <- 6.371e6 # m
g <- 9.81 # m/s^2
T_moon <- 27.32 * seconds_per_day # sidereal month
g * (R_earth / distance_earth_moon)^2 # inverse-square prediction
#> [1] 0.002694744
4 * pi^2 * distance_earth_moon / T_moon^2 # centripetal, from the orbit
#> [1] 0.0027236683 Newton’s Law and the N-Body Problem
Everything orbitr does rests on one equation. This chapter writes it down carefully, in the vector form the code actually uses, and then asks the questions a careful reader should ask about it: why the mass of the body being pulled drops out, why a planet can be treated as a point, what the constant \(G\) is and how well it is known, and what the computer does with the equation once it has it. By the end of the chapter you will have written your own version of orbitr’s inner loop in plain R and checked it against the package.
3.1 From the apple to the Moon
Newton’s law of universal gravitation says that two point masses \(m_1\) and \(m_2\) separated by a distance \(r\) attract each other with a force of magnitude
\[ F = \frac{G\,m_1 m_2}{r^2}. \tag{3.1}\]
The force is along the line joining them, and it falls off as the inverse square of the distance.
3.2 The law in vector form
A simulation needs directions, not just magnitudes. Let \(\mathbf{r}_j\) be the position vector of body \(j\). The vector from body \(j\) to body \(k\) is \(\mathbf{r}_k - \mathbf{r}_j\), its length is \(r_{jk} = |\mathbf{r}_k - \mathbf{r}_j|\), and the unit vector pointing from \(j\) toward \(k\) is
\[ \hat{\mathbf{r}}_{jk} = \frac{\mathbf{r}_k - \mathbf{r}_j}{r_{jk}}. \]
The force on \(j\) due to \(k\) points toward \(k\) (gravity attracts), so
\[ \mathbf{F}_{jk} = \frac{G\,m_j m_k}{r_{jk}^2}\,\hat{\mathbf{r}}_{jk}. \tag{3.2}\]
Swapping \(j\) and \(k\) flips the unit vector and nothing else, so \(\mathbf{F}_{kj} = -\mathbf{F}_{jk}\): Newton’s third law is built in. This is the fact behind the drifting Earth in Chapter 2. Because every internal force comes in an equal and opposite pair, the total momentum of a system of gravitating bodies never changes, no matter how complicated the motion.
3.3 Superposition: the N-body acceleration
Newton’s second law gives the acceleration of body \(j\) from the net force on it, \(m_j \ddot{\mathbf{r}}_j = \sum_k \mathbf{F}_{jk}\). Gravitational forces add as vectors; the pull of the Sun on Earth does not care that the Moon is also pulling. Summing Equation 3.2 over every other body and dividing by \(m_j\),
\[ \ddot{\mathbf{r}}_j = \sum_{k \neq j} \frac{G\,m_k}{r_{jk}^2}\,\hat{\mathbf{r}}_{jk} = \sum_{k \neq j} G\,m_k\,\frac{\mathbf{r}_k - \mathbf{r}_j}{|\mathbf{r}_k - \mathbf{r}_j|^3}. \tag{3.3}\]
This is the equation. The first form is the one in the package documentation; the second is the one the code uses, because it avoids computing the unit vector separately. The cube in the denominator is one power for turning the difference vector into a unit vector and two for the inverse-square law.
Two things to notice.
The mass of body \(j\) has cancelled. The acceleration of a body in a gravitational field does not depend on its own mass. A feather and a hammer fall together; the International Space Station and an astronaut floating inside it follow the same orbit. This is only true because the \(m_j\) in Newton’s second law (inertial mass, the resistance to being accelerated) is the same quantity as the \(m_j\) in the law of gravitation (gravitational mass, the thing that gravity pulls on). Nothing in Newtonian mechanics requires those to be equal; they just are. The equality has been tested to about one part in \(10^{15}\) (Touboul et al. 2022), and Einstein took it as the starting point of general relativity. For the simulation it means something practical: a body’s mass matters for how it pulls on the others, but a test particle of negligible mass follows exactly the same path it would with any other negligible mass.
The sum is over every pair. With \(N\) bodies there are \(N(N-1)\) ordered pairs, and each needs a subtraction, a norm, and a division. Doubling the number of bodies quadruples the work. This \(O(N^2)\) scaling is the single fact that dominates the design of every N-body code ever written, and it is why orbitr’s inner loop is in C++ (Chapter 15). For the handful of bodies in this book it is not a problem.
3.4 Why point masses are enough
Earth is not a point. It is a ball 12,742 km across, and every piece of it pulls on the Moon from a slightly different direction and distance. Why is it legitimate to write Equation 3.3 with a single \(\mathbf{r}_{\mathrm{Earth}}\)?
Theorem 3.1 (Newton’s shell theorem) A spherically symmetric body attracts an external point mass exactly as if all of its mass were concentrated at its center. A spherically symmetric shell exerts no net gravitational force on a point mass anywhere inside it.
Newton proved this geometrically in the Principia (Book I, Propositions 70 and 71), and he is often said to have delayed publishing the theory of gravitation until he had the proof, since without it the Moon test above is not justified. The modern proof takes three lines if you know Gauss’s law for gravity, which says that the flux of the gravitational field \(\mathbf{g}\) through any closed surface is proportional to the mass enclosed:
\[ \oint \mathbf{g}\cdot d\mathbf{A} = -4\pi G\,M_{\mathrm{enc}}. \]
Take the surface to be a sphere of radius \(r\) centered on the body. By symmetry, \(\mathbf{g}\) points radially and has the same magnitude \(g(r)\) everywhere on the sphere, so the integral is \(-g(r)\cdot 4\pi r^2\) (negative because the field points inward). Hence \(g(r) = G M_{\mathrm{enc}}/r^2\), directed toward the center. Outside the body, \(M_{\mathrm{enc}}\) is the whole mass, and the field is that of a point mass. Inside a shell, \(M_{\mathrm{enc}} = 0\), and the field vanishes.
So for any body whose mass is arranged in spherical shells, which is a very good description of stars and planets, the point-mass equation is not an approximation at all. It is exact, as long as the bodies do not overlap.
Where it fails is instructive. Earth bulges at the equator by about 21 km because it spins; the extra mass there pulls satellites in low orbit slightly off the point-mass prediction, and the orbits of GPS satellites precess measurably because of it. Astronomers call the leading correction \(J_2\). Two stars in a close binary raise tides on each other and are not spherical either. orbitr ignores all of this, which is the right choice for orbital distances and the wrong one for anything skimming a surface. The package roadmap lists \(J_2\) as a possible addition; Chapter 16 discusses what it would take.
3.5 Units, \(G\), and how well we know it
Everything in orbitr is in SI: meters, kilograms, seconds. The gravitational constant is
gravitational_constant
#> [1] 6.6743e-11in \(\mathrm{m^3\,kg^{-1}\,s^{-2}}\). That is the CODATA 2018 recommended value, \(6.67430(15)\times10^{-11}\) (Tiesinga et al. 2021), and the digits in parentheses are the uncertainty in the last two places: about 2 parts in \(10^5\). By the standards of fundamental constants that is terrible. The speed of light is exact by definition; Planck’s constant is now exact by definition; the fine-structure constant is known to a few parts in \(10^{10}\). \(G\) is hard to measure because gravity is so weak that laboratory experiments have to detect the pull of a kilogram-scale mass on another, and because there is no way to shield the apparatus from everything else.
The constant that astronomy actually uses is not \(G\) but the product \(GM\) for each body, called the gravitational parameter and written \(\mu\). For the Sun, \(GM_\odot\) is known to about one part in \(10^{10}\), because it is what you get directly from fitting planetary orbits and spacecraft trajectories, no laboratory required. When the package reports the Sun’s mass as \(1.989\times10^{30}\) kg, that number is \(GM_\odot\) divided by a poorly known \(G\), and its uncertainty is essentially the uncertainty in \(G\). None of this affects a simulation, because the equations of motion only ever contain \(G\) and a mass together. It does mean that if you ever compare a simulated orbital period to a measured one at the fifth significant figure, the discrepancy may be in \(G\), not in the integrator.
3.6 The equations of motion as a system
Equation 3.3 is a second-order differential equation for each of \(N\) bodies. A simulation does not work with second-order equations; it works with a state and a rule for advancing it. The state of the system at time \(t\) is the list of all positions and velocities,
\[ \big(\mathbf{r}_1, \ldots, \mathbf{r}_N,\ \mathbf{v}_1, \ldots, \mathbf{v}_N\big), \]
\(6N\) numbers in all. The rule for how the state changes is a pair of first-order equations for each body,
\[ \dot{\mathbf{r}}_j = \mathbf{v}_j, \qquad \dot{\mathbf{v}}_j = \mathbf{a}_j(\mathbf{r}_1, \ldots, \mathbf{r}_N), \tag{3.4}\]
where \(\mathbf{a}_j\) is the right-hand side of Equation 3.3. The accelerations depend only on positions, not on velocities. That is a special feature of gravity (and of any force that comes from a potential), and Chapter 5 will show that the best integrators exploit it.
Here is the structure of every N-body code, including this one, in pseudo-code:
state <- initial positions and velocities
for each time step:
a <- accelerations(state) # the N-body sum, all pairs
state <- advance(state, a, dt) # the integrator
record stateThe accelerations() function is the physics. The advance() function is the numerics. orbitr keeps them separate, and so will this book: this chapter is about the first function, Chapter 5 is about the second.
3.7 Writing the acceleration yourself
The best way to be sure what Equation 3.3 says is to implement it. Here it is in plain R, with no attempt at speed:
accelerations <- function(bodies, G = gravitational_constant) {
n <- nrow(bodies)
pos <- as.matrix(bodies[, c("x", "y", "z")])
a <- matrix(0, nrow = n, ncol = 3)
for (j in seq_len(n)) {
for (k in seq_len(n)) {
if (j == k) next
d <- pos[k, ] - pos[j, ] # vector from j toward k
r <- sqrt(sum(d^2)) # distance
a[j, ] <- a[j, ] + G * bodies$mass[k] * d / r^3
}
}
tibble(id = bodies$id, ax = a[, 1], ay = a[, 2], az = a[, 3])
}The double loop runs over all ordered pairs \((j, k)\) with \(j \neq k\). For each pair it forms the difference vector, its length, and the term in the sum. bodies is any data frame with columns id, mass, x, y, z, which is exactly what get_bodies() returns:
em <- create_system() |>
add_body("Earth", mass = mass_earth) |>
add_body("Moon", mass = mass_moon, x = distance_earth_moon, vy = speed_moon)
accelerations(get_bodies(em))
#> # A tibble: 2 × 4
#> id ax ay az
#> <chr> <dbl> <dbl> <dbl>
#> 1 Earth 0.0000332 0 0
#> 2 Moon -0.00270 0 0The Moon’s acceleration is in the \(-x\) direction, toward Earth, with magnitude \(G M_{\mathrm{Earth}}/r^2\):
gravitational_constant * mass_earth / distance_earth_moon^2
#> [1] 0.002697483and Earth’s is in the \(+x\) direction, toward the Moon, smaller by the mass ratio \(m_{\mathrm{Moon}}/M_{\mathrm{Earth}} \approx 1/81\). Equal and opposite forces; unequal accelerations.
Now a check against the package. The package’s output does not include accelerations, but it does include velocities at every step, and over one hour the velocity change of the Moon should be close to \(\mathbf{a}\,\Delta t\):
sim <- simulate_system(em, time_step = seconds_per_hour,
duration = seconds_per_hour)
sim |>
filter(id == "Moon") |>
summarise(dvx_per_second = diff(vx) / seconds_per_hour)
#> # A tibble: 1 × 1
#> dvx_per_second
#> <dbl>
#> 1 -0.00270That is the \(x\) acceleration computed above, to about four significant figures. (It is not exact because the acceleration changes slightly during the hour as the Moon moves; the integrator uses the average of the accelerations at the start and end of the step, which is one of the things that makes it good.)
3.8 A three-body example
With three bodies the sums in Equation 3.3 have two terms each. Put the Sun at the origin, Earth at 1 AU, and the Moon 384,400 km beyond Earth, and ask how hard the Sun pulls on the Moon compared to how hard Earth does:
sem <- create_system() |>
add_sun() |>
add_body("Earth", mass = mass_earth, x = distance_earth_sun) |>
add_body("Moon", mass = mass_moon,
x = distance_earth_sun + distance_earth_moon)
a_sun_on_moon <- gravitational_constant * mass_sun /
(distance_earth_sun + distance_earth_moon)^2
a_earth_on_moon <- gravitational_constant * mass_earth /
distance_earth_moon^2
a_sun_on_moon / a_earth_on_moon
#> [1] 2.187709The Sun pulls on the Moon more than twice as hard as Earth does. This surprises almost everyone the first time. The Moon orbits Earth not because Earth’s pull is stronger but because Earth and the Moon are both falling toward the Sun together, at nearly the same rate; what keeps the Moon bound to Earth is the small difference between the Sun’s pull on the two of them. Seen from the Sun, the Moon’s path is a slightly wobbly ellipse that is everywhere concave toward the Sun. Chapter 10 draws it.
The full acceleration of the Moon is the vector sum of both terms, and our function gets it right with no changes:
accelerations(get_bodies(sem)) |>
filter(id == "Moon")
#> # A tibble: 1 × 4
#> id ax ay az
#> <chr> <dbl> <dbl> <dbl>
#> 1 Moon -0.00860 0 03.9 What the engine does differently
The function above is correct and slow. Each call does \(N(N-1)\) iterations of interpreted R, and a simulation calls it once or twice per step for hundreds of thousands of steps. The package’s pure-R fallback vectorizes the same computation with matrix outer products so that R does the loop in compiled code internally; its C++ engine writes the double loop directly in C++ and runs about as fast as the hardware allows. Both compute exactly Equation 3.3, with one small addition: an optional softening length \(\varepsilon\) that replaces \(r^3\) in the denominator with \((r^2 + \varepsilon^2)^{3/2}\). With \(\varepsilon = 0\), the default, there is no difference. Chapter 6 explains what softening is for.