A very good place to start analyzing any motion is choosing a frame of reference and setting up coordinates. Orbital mechanics systems are often considered with respect to the fixed stars.

We usually try to find a reference frame that is not accelerating. This is hard! Everything turns out to be wagging a little, even what we can’t see or feel. But supposing we do find an inertial frame, we would see that Newton’s laws of motion apply:

  1. A body remains at rest, or in motion at a constant speed in a straight line, unless acted upon by a force.
  2. When a body is acted upon by a force, the time rate of change of its momentum equals the force.
  3. If two bodies exert forces on each other, these forces have the same magnitude but opposite directions.

In this way we avoid having to think through fictitious forces such as the centrifugal or Coriolis.

A system of two bodies

We first look at the mutual gravitation of two masses, \(m_0\) and \(m_1\) with positions \({\bf r}_0\) and \({\bf r}_1\) described with their Cartesian coordinates:

\[ {\bf r}_0 = x_0 {\bf \hat{x}} + y_0 {\bf \hat{y}} + z_0 {\bf \hat{z}} \\ {\bf r}_1 = x_1 {\bf \hat{x}} + y_1 {\bf \hat{y}} + z_1 {\bf \hat{z}} \]

The relative position, \({\bf r}\) and its unit vector, \(\hat{\bf u}_\mathrm{r}\) can be calculated.

\[ {\bf r} = {\bf r}_1 - {\bf r}_0 = (x_1 - x_0) {\bf \hat{i}} + (y_1 - y_0) {\bf \hat{j}} + (z_1 - z_0) {\bf \hat{k}} \\ \hat{\bf u}_\mathrm{r} = \frac{\bf r}{r} \]

Example

df      <- data.frame(m = 1.0e26, x = 100, y = 120)
df[2, ] <- data.frame(m = 5.0e25, x = 150, y = 50)

Here we are imagining the only two bodies in the universe! (Or that everyone else is really, really, far away.) So there are no other forces acting and the force on \(m_0\) is opposite that on \(m_1\) by Newtons’ third law, i.e. \({\bf F}_0 = -{\bf F}_1\). Using the law of gravitation (\(F = G \frac{m_0 m_1}{r^2}\)) we can find the force acting on the line between the masses.

\[ \begin{align} {\bf F}_0 &= \frac{G m_0 m_1}{r^2} \hat{\bf u}_\mathrm{r} \\ {\bf F}_1 &= - \frac{G m_0 m_1}{r^2} \hat{\bf u}_\mathrm{r} \end{align} \]

Now applying Newton’s second law (\({\bf F}_i = m {\bf \ddot{r}}_i\)) and the definition of \(\hat{\bf u}_\mathrm{r}\) above we get the accelerations of the two masses, which depend only on the relative position.

\[ \begin{align} \ddot{\bf r}_0 &= G m_1 \frac{\bf r}{r^3} \\ \ddot{\bf r}_1 &= - G m_0 \frac{\bf r}{r^3} \end{align} \]

We can integrate the accelerations to get equations for the velocities and positions of each mass \(m_\mathrm{i}\). Keeping in mind that in general the positions, \({\bf r}_\mathrm{i}\), velocities, \(\dot{\bf r}_\mathrm{i}\), and accelerations, \(\ddot{\bf r}_\mathrm{i}\), depend on time.

\[ \dot{\bf r}(t)_\mathrm{i} = \dot{\bf r}(0)_\mathrm{i} + \int_{0}^{t} \ddot{\bf r}(t)_\mathrm{i} \mathrm{d}t \\ {\bf r}(t)_\mathrm{i} = {\bf r}(0)_\mathrm{i} + \int_{0}^{t} \dot{\bf r}(t)_\mathrm{i} \mathrm{d}t \]

We see that we will need the initial positions, \({\bf r}(0)_\mathrm{i}\), and velocities, \(\dot{\bf r}(0)_\mathrm{i}\), to solve the motion. Together we call these 12 quantities the state vector of the system: \(S = [{\bf r}_0, {\bf r}_1, {\bf \dot{r}}_0, {\bf \dot{r}}_1]\).

Numerically in R

Initial conditions

We input the gravitational constant, \(G = 6.674 \times 10^{-2} \text{km}^3 \text{kg}^{-1} \text{s}^{-2}\), and an equal mass for the particles, \(m_0 = m_1 = 10^{26} \text{kg}\). For reference Earth’s mass is about \(6 \times 10^{24} \text{kg}\), so these are somewhat larger.

G <- 6.674e-2
m_0 <- 1.0e26
m_1 <- m_0

We represent the state vector as a data.frame where each record contains the state vector of a single mass and a timestamp.

s_vector      <- data.frame(n_mass=0, ts=0, x=0, y=0, z=0, dot.x=0, dot.y=0, dot.z=0)
s_vector[2, ] <- data.frame(n_mass=1, ts=0, x=2000, y=0, z=0, dot.x=-2000, dot.y=500, dot.z=0)

We have \(m_0\) starting out stationary at \((0, 0, 0)\) and \(m_1\) starting out at \((2000, 0, 0)\) with some initially counterclockwise velocity.

Stepwise solutions

We can advance the state vector by an arbitrary time step \(\Delta t\) if we assume the velocities and accelerations are constant over that time slice. This can be a good approximation if the time slice is small compared to the characteristic time scale of the motion.

\[ \dot{\bf r}(t_0 + \Delta t)_\mathrm{i} = \dot{\bf r}(t_0)_\mathrm{i} + \ddot{\bf r}(t_0)_\mathrm{i} \Delta t \\ {\bf r}(t_0 + \Delta t)_\mathrm{i} = {\bf r}(t_0)_\mathrm{i} + \dot{\bf r}(t_0)_\mathrm{i} \Delta t. \]

We develop a function that takes this state vector and returns a new state vector including the current relative distance and accelerations as well as the next positions and velocities.

Show code 
#' Get relative distance and accelerations from initial positions
#' @return A data frame containing the relative distance and accelerations
get_derived_initial_values <- function(pos_0, pos_1) {
  # initial distance
  r0 <- sqrt(sum((pos_1 - pos_0)^2))
  # initial accelerations
  ddot_r0_0 <- m_1 * G * (pos_1 - pos_0) / (r0^3)
  ddot_r0_1 <- m_0 * G * (pos_0 - pos_1) / (r0^3)
  init_vals <- data.frame(r=r0, ddot_r0=ddot_r0_0, ddot_r1=ddot_r0_1)
  init_vals
}

#' Step forward linearly
step <- function(initial, dot_initial, delta) {
  final <- initial + dot_initial * delta
  final
}

#' Return new state vector after advancing by stepsize
#'
#' @return A data frame containing the original state vector with new values
#'         appended after stepsize as well as the derived initial quantities
get_next_s_vector <- function(s_vector, stepsize=1) {
  if (nrow(s_vector) < 2) {
    return(666)
  }
  # extract current positions (last rows)
  pos_0 <- s_vector[nrow(s_vector) - 1, 3:5]
  pos_1 <- s_vector[nrow(s_vector), 3:5]
  # extract current velocities (last rows)
  vel_0 <- s_vector[nrow(s_vector) - 1, 6:8]
  vel_1 <- s_vector[nrow(s_vector), 6:8]

  # calculate current distance and accelerations
  init_vals = get_derived_initial_values(pos_0, pos_1)

  # annotate current state vector with derived initial values
  s_vector[nrow(s_vector) - 1, 9:12] <- data.frame(init_vals[1], init_vals[2:4])
  s_vector[nrow(s_vector), 9:12] <- data.frame(init_vals[1], init_vals[5:7])

  # step forward
  new_pos_0 <- step(pos_0, vel_0, stepsize)
  new_pos_1 <- step(pos_1, vel_1, stepsize)
  new_vel_0 <- step(vel_0, init_vals[2:4], stepsize)
  new_vel_1 <- step(vel_1, init_vals[5:7], stepsize)
  s_vector[nrow(s_vector) + 1, 1:8] <- data.frame(n_mass=0, ts=s_vector[nrow(s_vector) - 1, 2] + stepsize,
                                                  new_pos_0, new_vel_0)
  s_vector[nrow(s_vector) + 1, 1:8] <- data.frame(n_mass=1, ts=s_vector[nrow(s_vector) - 1, 2] + stepsize,
                                                  new_pos_1, new_vel_1)
  s_vector
}

And now we can see how the motion of the two masses evolves in time.

stepsize <- 5e-10
new_s_vector <- s_vector
for(timestep in 1:58) { new_s_vector <- get_next_s_vector(new_s_vector, stepsize) }

Here are the first few time slices for \(m_0\):

n_mass ts x y z dot.x dot.y dot.z r ddot_r0.x ddot_r0.y ddot_r0.z
0 0e+00 0.000000 0 0 0 0.0000000 0 2000.000 1.668500e+18 0 0
0 0e+00 0.000000 0 0 834250000 0.0000000 0 2000.000 1.668500e+18 208562500 0
0 1e-09 0.417125 0 0 1668500001 0.1042813 0 1999.166 1.669893e+18 417647417 0
0 2e-09 1.251375 0 0 2503446411 0.3131050 0 1997.497 1.672684e+18 627954979 0

and for \(m_1\):

n_mass ts x y z dot.x dot.y dot.z r ddot_r0.x ddot_r0.y ddot_r0.z
1 0e+00 2000.000 0.0e+00 0 -2000 500.0000 0 2000.000 -1.668500e+18 0 0
1 0e+00 2000.000 2.5e-07 0 -834252000 500.0000 0 2000.000 -1.668500e+18 -208562500 0
1 1e-09 1999.583 5.0e-07 0 -1668502001 499.8957 0 1999.166 -1.669893e+18 -417647417 0
1 2e-09 1998.749 7.5e-07 0 -2503448411 499.6869 0 1997.497 -1.672684e+18 -627954979 0

The derived quantities \(r\), and \({\bf \ddot{r}}\) are not populated for the last record as they are generated during the next time step.

All the motion has been set in the \(x\)-\(y\) plane through the initial conditions.

We can see the mutual attraction pull \(m_0\) off of the origin as the two bodies approach each other, eventually getting within \(80 \text{km}\) where the maximum gravitational acceleration slingshots them back away from each other at high speed.

Center of Gravity ->

Orbits! | code