# Solution to two-body problem in orbital mechanics for r(t) and theta(t), rather than r(theta)?

I have written a simple numerical integration code to calculate the orbits of two planetary bodies orbiting a star, in order to calculate the transit-timing variation for one body due to the gravitational perturbations of the other. I am stuck on two issues:

(1) I want to be able to check my numerical code for the one-planet case by deriving an analytical expression for x(t) and y(t), or equivalently theta(t) and r(t), as an initial value problem using initial positions and velocities. I can find plenty of derivations online for r(theta), but none that explicitly solve for position as a function of time. My understanding is that this should be possible, but I can't find any references that work it out. I am ok assuming the unperturbed orbit is elliptical and stable, but hopefully not necessarily circular. How would I go about doing this? Note: I have seen the answers to this question, which state that there are analytic solutions but they are elliptic functions. However, none of the answers show how to derive this, or whether there might be nicer solutions under certain simplifications/assumptions.

(2) I would also like to be able to check that the values I calculate for the transit-timing variations are reasonable. It seems like this is a very involved calculation (ie. as discussed in this paper). Is there any simplified calculation that could be made (perhaps making some restrictive assumptions on the orbital configuration) that would allow me to check my numerical results analytically?

Thank you in advance for any pointers!

• Great question for the astronomy stack exchange site. Not so good for this one. Nov 18, 2022 at 18:24
• Analytical solutions would be of the form $t(r)$ rather than $r(t)$, see answers to What is the analytical closed-form solution of the two-body problem to verify its numerical integration results? Also, answers to [Can the radial oscillations of an elliptical orbit be solved using a fictitious centrifugal potential?]() in Physics SE may shed some light.
– uhoh
Nov 18, 2022 at 18:34
• What uhoh said. To check your integration analytically, it's much easier to calculate mean anomaly (time) as a function of true anomaly (angle), rather than the other way around, which involves solving Kepler's equation. Nov 18, 2022 at 18:37
• @PM2Ring Asking is Physics SE was wrong. I will answer here where techniques and challenges related to accurate numerical simulation of orbits is discussed, questioned and answered all the time.
– uhoh
Nov 18, 2022 at 22:14
• @PM2Ring cf. space.stackexchange.com/a/37335/12102
– uhoh
Nov 18, 2022 at 22:20

The OP has

written a simple numerical integration code to calculate the orbits of two planetary bodies orbiting a star, in order to calculate the transit-timing variation for one body due to the gravitational perturbations of the other.

That means they need to carefully check the accuracy of the integration. They want

to be able to check my numerical code for the one-planet case by deriving an analytical expression for x(t) and y(t), or equivalently theta(t) and r(t), as an initial value problem using initial positions and velocities...

I will write an answer based on my best understanding of the question and the OP's comment:

I am looking for something that can be solved as an initial value problem (ie. I specify the initial location and velocity vectors of the planet, write down the differential equations for acceleration, and go from there), rather than something that depends on knowing (and having well-defined) orbital parameters like eccentricity from the start.

which will benefit from @PM2Ring's comment:

...And to get the eccentricity & the orientation of the semi-major axis, calculate the eccentricity vector.

We can search this site for eccentricity vector and also review all the great answers to Converting Orbital Elements to Cartesian State Vectors and have fun inverting them in 3D. See also answers to:

But I'll answer in 2D assuming the two inter-perturbing planets will be coplanar.

Let's start with this answer to Why does the eccentricity vector equation always equal -1?. If you have a position vector $$\mathbf{x}$$ and velocity vector $$v$$ of an object in orbit around a central mass $$\mu$$, then

$$\mathbf{e} = {v^2 \mathbf{r} \over {\mu}} - {(\mathbf{r} \cdot \mathbf{v} ) \mathbf{v} \over{\mu}} - {\mathbf{r}\over{r}}$$

where $$\mathbf{e}$$ is called the eccentricity vector. Bold like $$\mathbf{e, r, v}$$ indicates vector quantities and italics like $$e, r, v$$ are their scalar magnitude.

With $$\mu=1$$ for now, let's look at $$\mathbf{r} = [0.42, 1.42]$$ and $$\mathbf{v} = [0.86, 0.45]$$ that gives $$e=$$0.75069393.

We have the scalar vis-viva equation $$v^2 = \mu(2/r - 1/a)$$ and we know $$r_{peri}, r_{apo} = a(1-e), a(1+e)$$ and we know that a vector in the direction of $$\mathbf{e}$$ with length $$r_{peri}$$ will point to the location of the periapsis.

First let's check by numerical integration that our derived orbital elements are correct. Big dots indicate periapsis, starting point and apoapsis, integration is for one period, and yes it returns to the same spot to within about 1 part per million. You can search for different/"better" integrators for better accuracy. If we switch to polar coordinates for the analytical models, the periapsis is at $$\theta_0 = \text{arctan2}(\mathbf{e}_y, \mathbf{e}_x)$$ or -174.23915 degrees.

We'll use that as an offset to convert the analytical solution in polar coordinates back to cartesian. From this answer to What is the analytical closed-form solution of the two-body problem to verify its numerical integration results?:

The relationship between $$E$$ and the true anomaly $$\theta = \arctan2(y, x)$$ is

$$\tan \frac{\theta}{2} = \sqrt{ \frac{1+e}{1-e} } \tan \frac{E}{2}$$

and solving for $$E$$:

$$E(\theta) = 2 \arctan \sqrt{ \frac{1-e}{1+e} } \tan \frac{\theta}{2}.$$

plugging back in to the first equation (but not writing it all out):

$$t(\theta) = a \sqrt{\frac{a}{\mu}}\left(E(\theta) - e \sin E(\theta) \right)$$

I've done that in the Python script below, and superimposed the analytical solution recast in cartesian coordinates onto the same plot. The individual dots don't line up of course since the numerical starting point was arbitrary, but you can use some scheme to find points in time where the two coincide. I'll leave that as an exercise for the reader. Python script for the plots and math

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import time

def get_evec_and_a(rvec, vvec, mu):
vsq  = (vvec**2).sum()
r = np.sqrt((rvec**2).sum())
evec = (vsq*rvec - (rvec*vvec).sum() * vvec) / mu - rvec/r
a = 1 / (2/r - vsq/mu) # vis-viva equation
return evec, a

mu = 1.0
rvec = np.array([0.42, 1.414])  # https://en.wikipedia.org/wiki/42_(number)
vvec = np.array([0.86, 0.45])  # https://english.stackexchange.com/q/399555/217285

evec, a = get_evec_and_a(rvec, vvec, mu)
e = np.sqrt((evec**2).sum())

r_peri, r_apo = a*(1-e), a*(1+e)

# https://space.stackexchange.com/a/30127/12102
evec_norm = evec / e # points towards periapsis

rvec_peri = r_peri * evec_norm
rvec_apo = -r_apo * evec_norm

print('eccentricity e: ', e)
print('periapsis direction (normalized): ', evec_norm)
print('peri, apo: ', r_peri, r_apo)

# integrate for one period T to check closure

def deriv(t, state_vector, mu):
rvec, vvec = state_vector.reshape(2, -1)
rsq = (rvec**2).sum()
acc = -mu * rvec * rsq**-1.5
return np.hstack([vvec, acc])

T = 2 * np.pi * np.sqrt(a**3 / mu)

t_eval = np.linspace(0, T, 101) # interpolated output values, not internal timesteps
t_span = t_eval, t_eval[-1]
args = (mu, )

initial_state_vector = np.hstack([rvec, vvec])

t_start = time.process_time()

args=args, t_eval=t_eval, rtol=1E-12, method='DOP853',
dense_output=False, events=None)

process_time = time.process_time() - t_start
print('process time: ', process_time)

print('state_vectors.shape: ', state_vectors.shape)

if True:
fig, ax = plt.subplots(1, 1)
plt.plot(, , 'oy', ms=20) # Star
for vec, color in zip((rvec_peri, rvec, rvec_apo), 'rgb'):
ax.plot(vec[:1], vec[1:], 'o', color=color, ms=10)
ax.arrow(0, 0, evec, evec, color='r', width=0.01, head_width=0.07)
ax.set_aspect('equal')
x, y, vx, vy = state_vectors
ax.plot(x, y, '-k', linewidth=0.5)
ax.plot(x, y, '.k')
plt.show()

print('initial state vector: ', state_vectors[:, 0])
print('final state vector: ', state_vectors[:, -1])
print('difference: ', state_vectors[:, -1] - state_vectors[:, 0])

theta_peri = np.arctan2(evec, evec)
print('theta_peri = ', np.degrees(theta_peri))

# analytical solution for Kelper orbit
# https://math.stackexchange.com/a/3463417/284619
# https://space.stackexchange.com/a/43237/12102
theta = np.linspace(0, 2 * np.pi, 201)
r = a * (1 - e**2) / (1 + e * np.cos(theta))
E = 2. * np.arctan(np.sqrt((1.-e)/(1.+e)) * np.tan(theta/2))
t = a * np.sqrt(a/mu) * (E - e * np.sin(E))
x_analytic = r * np.cos(theta + theta_peri)
y_analytic = r * np.sin(theta + theta_peri)

if True:
fig, ax = plt.subplots(1, 1)
plt.plot(, , 'oy', ms=20) # Star
for vec, color in zip((rvec_peri, rvec, rvec_apo), 'rgb'):
ax.plot(vec[:1], vec[1:], 'o', color=color, ms=10)
ax.arrow(0, 0, evec, evec, color='r', width=0.01, head_width=0.07)
ax.set_aspect('equal')
x, y, vx, vy = state_vectors
ax.plot(x, y, '-k', linewidth=0.5)
ax.plot(x, y, '.k')
ax.plot(x_analytic, y_analytic, '-b')
ax.plot(x_analytic, y_analytic, '.b')
plt.show()