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[0], t_eval[-1]
args = (mu, )
initial_state_vector = np.hstack([rvec, vvec])
t_start = time.process_time()
answer = solve_ivp(deriv, t_span=t_span, y0=initial_state_vector,
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)
state_vectors = answer['y']
print('state_vectors.shape: ', state_vectors.shape)
if True:
fig, ax = plt.subplots(1, 1)
plt.plot([0], [0], '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[0], evec[1], 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[1], evec[0])
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([0], [0], '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[0], evec[1], 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()
