Three astronauts; Goofus, Gallant and Zippy the Pinhead play a game. All three take their cold gas thruster-powered jet packs and quickly move 100 meters away from the ISS. For some reason after accelerating and decelerating to a stop, all three ran out of propellant.

Goofus moved 100 meters out in front of the ISS while Gallant moved to the side; and Zippy the Pinhead moved 100 meters farther from Earth so that he could be seen from ISS' purported zenith-facing window (1, 2, 3).

The four position vectors at the beginning for ISS, Goofy, Gallant and Zippy are:

``````x0 = [a, 0, 0.] + [a, 100, 0.] + [a, 0, 100.] + [a+100, 0, 0.]
``````

and all four starting velocity vectors are `[0, speed, 0.]`

After one orbit around the Earth had elapsed, how far did each of our three heroes end up from the ISS?

This is what I got, but I'd prefer to accept an answer based more on math and principles than brute-force numerical integration.

J2 is important but for the purposes of other answers it can be ignored. ``````import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

def deriv(X, t):
x, v = X.reshape(2, -1)
accs = []
for xx in x.reshape(-1, 3):
acc = -xx * GMe * ((xx**2).sum())**-1.5
if use_J2:
x, y, z = xx
xsq, ysq, zsq = xx**2
rm7 = (xsq + ysq + zsq)**-3.5
accJ2x = x * rm7 * (6*zsq - 1.5*(xsq + ysq))
accJ2y = y * rm7 * (6*zsq - 1.5*(xsq + ysq))
accJ2z = z * rm7 * (3*zsq - 4.5*(xsq + ysq))
acc -= J2 * np.hstack((accJ2x, accJ2y, accJ2z))
accs.append(acc)
return np.hstack((v, np.hstack((accs))))

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
GMe = 3.98600435436E+14 # earth (from DE430)
Re = 6378136.3 # meters
J2_earth = -1.08262545E-03  # unitless
J2 = J2_earth * (GMe * Re**2)

inc = 51*pi/180.

a = Re + 400*1000.
speed = np.sqrt(GMe/a) # close enough for government work
x0 = [a, 0, 0.] + [a, 100, 0.] + [a, 0, 100.] + [a+100, 0, 0.]
v0 = 4*[0, speed, 0.]
X0 = np.array(x0 + v0, dtype=float).reshape(-1, 3)
s, c = [f(inc) for f in (np.sin, np.cos)]
R = np.array([[1, 0, 0], [0, c, -s], [0, s, c]])
X0 = (X0[:, None, :] * R).sum(axis=2).flatten()

times = np.arange(0, 92*60, 10)

use_J2 = True
answer, info = ODEint(deriv, X0, times, full_output=True)
use_J2 = False
answer_noJ2, info = ODEint(deriv, X0, times, full_output=True)