I have a gravitational nbody simulation, for which I would like to determine various orbital parameters. For each body, I have 3-D vectors (x,y,z -space) for position, velocity, and acceleration. I am trying to follow the steps laid out in this post to obtain the eccentricity of each orbit. Before throwing n bodies into the simulation, I am testing the algorithm on simpler systems, such as a 2-body system in which the orbital path of the Earth around the Sun is nearly a perfect circle. Because the orbit is circular, I expect the eccentricity to be zero; this is not the output I get, so I am hoping someone can help me identify my errors (either in understanding or in code). Specifically, I would like to know what I am doing wrong in trying to calculate the eccentricity.

Sorry ahead of time for the length of this post; most of the code below is to show that the methodology works to obtain vectors of position and velocity; the last part of code (skip down to PROBLEM) is to "show my work" in using these parameters to calculate the eccentricity vectors. Aside from visual inspection, methods from this post were used to ensure that the orbit is circular.

Create circular orbit via Sun-Earth system

First, we'll initialize the initial conditions of our coupled ODEs and relevant simulation parameters.

import numpy as np
import matplotlib.pyplot as plt

## simulation parameters
ndim = 3 ## x,y,z
gravitational_constant = 6.67e-11 ## SI units
nbodies = 2 ## sun, earth
duration = 365*24*60*60 ## duration; 1 years --> seconds; day/yr * hr/day * min/hr * sec/min
dt = 2 * 24 * 60 * 60 ## time-step; 2 days --> seconds
t = np.arange(duration/dt)

meters_to_au = 1.496e11 ## 1.496e11 meters = 1 AU

## BODY 1 (sun)
m_sun = 1.989e30 ## kilograms
x_sun = np.zeros(ndim) ## position (x,y,z); meters
v_sun = np.zeros(ndim) ## velocity (x,y,z); m/s

## BODY 2 (earth)
m_earth = 5.972e24 ## kilograms
x_earth = np.array([meters_to_au, 0, 0]) ##
_v = np.sqrt(gravitational_constant * m_sun / meters_to_au)
v_earth = np.array([0, _v, 0])

## standard gravitational parameters and reduced mass
mu = np.array([m_sun, m_earth]) * gravitational_constant
mred = (m_sun * m_earth) / (m_sun + m_earth)

Then, we solve the coupled ODEs using a simple Euler method.

## initialize SOLUTION SPACE
X = np.zeros((nbodies, ndim, t.size))
V = np.zeros((nbodies, ndim, t.size))
xi = np.array([x_sun, x_earth])
X[:, :, 0] = xi ## position of bodies at time t=0
vi = np.array([v_sun, v_earth])
V[:, :, 0] = vi ## velocity of bodies at time t=0

## ITERATE (i --> k=i+1)
for ti in range(1, t.size): ## t=1, ..., t=end
    ak = []
    for j in range(nbodies):
        dacc = 0
        for k in range(nbodies):
            if j != k:
                dpos = xi[j, :] - xi[k, :]
                r = np.sum(np.square(dpos))
                dacc -= mu[k] * dpos / np.sqrt(r**3)
    ak = np.array(ak)
    vk = vi + ak * dt
    xk = xi + vk * dt
    X[:, :, ti] = xk
    V[:, :, ti] = vk
    xi, vi = xk, vk

Xs = X[0, :, :]
Xe = X[1, :, :]

Vs = V[0, :, :]
Ve = V[1, :, :]

To verify that the simulation ran as expected, we plot.

fig, ax = plt.subplots(figsize=(7,7))
ax.scatter(Xe[0, :] / meters_to_au, Xe[1, :] / meters_to_au, marker='.', color='steelblue', s=2, label='Earth')
ax.scatter(Xs[0, :] / meters_to_au, Xs[1, :] / meters_to_au, marker='*', color='darkorange', s=5, label='Sun')
ax.set_xlabel('X (AU)', fontsize=8)
ax.set_ylabel('Y (AU)', fontsize=9)
fig.legend(mode='expand', loc='lower center', ncol=2, fontsize=8)

Position vectors


I am more familiar with seeing angular momentum expressed as $L = \vec{r} x \vec{p}$, where $\vec{p} = m\vec{v}$, though I suppose one can interpret the angular momentum below expressed in units of angular momentum per unit mass. In Cartesian coordinates, $\vec{r} = \vec{x} + \vec{y} + \vec{z} = x\hat{x} + y\hat{y} + z\hat{z}$.

Le = np.cross(Xe, Ve, axis=0)
Ls = np.cross(Xs, Vs, axis=0)

Ee = np.cross(Ve, Le, axis=0) / mred - Xe / np.sqrt(np.sum(np.square(Xe), axis=0))
Es = np.cross(Vs, Ls, axis=0) / mred - Xs / np.sqrt(np.sum(np.square(Xs), axis=0))
mag_Ee = np.sqrt(np.sum(np.square(Ee), axis=0))
mag_Es = np.sqrt(np.sum(np.square(Es), axis=0))

fig, ax = plt.subplots(figsize=(7,7))
ax.scatter(Ee[0, :], Ee[1, :], marker='.', color='steelblue', s=2, label='Earth')
ax.scatter(Es[0, :], Es[1, :], marker='*', color='darkorange', s=5, label='Sun')
ax.set_aspect('equal') ## x- and y- scales are equal; nearly perfect circle
ax.set_xlabel(r'eccentricity $\hat{x}$', fontsize=8)
ax.set_ylabel(r'eccentricity $\hat{y}$', fontsize=8)
fig.legend(mode='expand', loc='lower center', ncol=2, fontsize=8)

Eccentricity Vectors

rescaled_t = t * dt
fig, ax = plt.subplots(figsize=(7,7))
ax.scatter(rescaled_t, mag_Ee, marker='.', color='steelblue', s=2, label='Earth', alpha=0.5)
ax.scatter(rescaled_t, mag_Es, marker='*', color='darkorange', s=5, label='Sun', alpha=0.5)
ax.set_xlabel('Time', fontsize=8)
ax.set_ylabel('Eccentricity', fontsize=8)
ax.set_ylim(bottom=-0.1, top=1.2)
fig.legend(mode='expand', loc='lower center', ncol=2, fontsize=8)

Magnitude of Eccentricity over Time

It is my understanding that eccentricity varies as $0 ≤ e < 1$ for elliptical orbits (circular orbits being $e=0$), $e=1$ for parabolic orbits, and $e>1$ for hyperbolic orbits. So something must be off. Do I need to consider the coordinates from a specific reference frame? Or maybe I missed an assumption for the equations used to hold? Can someone point to the cause of this error? Less importantly, is the equation used to calculate eccentricity generalizable to all orbits or just elliptical ones?


1 Answer 1


You are doing many things wrong.

  1. You are computing the eccentricity of one body with respect to the center of mass. You need to compute the eccentricity of one body with respect to the other.

  2. You are using reduced mass in np.cross(Ve, Le, axis=0) / mred - Xe / np.sqrt(np.sum(np.square(Xe), axis=0)) This is wrong for multiple reasons. First off, look at the units! The first term, np.cross(Ve, Le, axis=0) / mred, has units of length^3/time^2/mass. The second term, np.sqrt(np.sum(np.square(Xe), axis=0)), is unitless. And you should not be using reduced mass at all. You should be using the combined gravitational parameter (not the reduced gravitational parameter). A gravitational parameter has units of length^3/time^2.

  3. To calculate the eccentricity correct, compute the position of the Earth with respect to the Sun (Xrel = Xe - Xs and the velocity of the Earth with respect to the Sun (Vrel = Ve - Vs). Then compute the cross product of these two (Lrel = np.cross(Xrel, Vrel) to yield the specific angular momentum of the Sun-Earth system. Finally, compute the eccentricity vector via np.cross(Vrel, Lrel) / mu_combined - Xrel / np.sqrt(np.sum(np.square(XRel))), where mu_combined is the sum of the gravitational parameters of the Sun and the Earth.

Finally, as a comment rather than a critique, it is best not use mass and the universal gravitational constant. It is much better to use gravitational parameters. You can find a fairly accurate list of solar system gravitational parameters in the wikipedia standard gravitational parameter article. Conceptually, a body's gravitational parameter is equal to the product of its mass and the gravitational constant. Another way of looking at it is that a body's mass is the body's gravitational parameter divided by the gravitational constant. The problem is that the gravitational constant is only known to four or five decimal places, while a body's gravitational parameter is observable and is known to six or more decimal places.

  • $\begingroup$ Thanks for the detailed response; very helpful! I do have one related follow-up question, regarding how to treat (1) and (3) in a system with more than two bodies. If I consider a system with a few planets (say, 3) in circular orbits about a star, then is it best to calculate the eccentricity of each planetary orbit wrt only dominant central mass (ie, not account for perturbations from other bodies)? Also, if we consider an exotic orbit (perhaps a stable horseshoe or an unstable figure-eight), then does one consider orbits piece-wise (wrt to dominant central mass) or in some other way? $\endgroup$
    – user33354
    Commented Sep 30, 2020 at 7:32
  • $\begingroup$ @allthemikeysaretaken - You essentially are asking about osculating orbital elements, or perhaps some time-averaged version of them. When one does this, one finds that the five Keplerian elements that are supposed to be constant are not constant. Mercury's orbit precesses by about 575 arc seconds per Julian century. Most of this precession, about 532 arc seconds per century, is attributable to Newtonian influences by other planets. The discrepancy, 43 arc seconds per century, is due to general relativity. $\endgroup$ Commented Oct 1, 2020 at 10:53
  • $\begingroup$ If I understood the linked post correctly, one chooses $t$ as the time elapsed since periapsis (for a 2-body system); one calculates the 3 anomalies for any given time $t$. If so, then the time-dependence of these parameters makes sense (ie, perturbations). I realize this simulation won't be accurate to relativistic order. (I read somewhere on stackexchange about a relativistic correction to the gravitational potential that I wanted to play with; will search for link tomorrow). That said, I'd like to understand this better and every example I've seen is a classic ~circular orbit. $\endgroup$
    – user33354
    Commented Oct 1, 2020 at 11:34
  • $\begingroup$ @allthemikeysaretaken - Osculating elements are relatively easy to calculate. (1) Compute the specific angular momentum vector, $\vec h = \vec r \times \vec v$. The $z$ component yields the inclination while the $x$ and $y$ components yield the right ascension of ascending node. (2) Compute the eccentricity vector, $\vec v \times \vec h / \mu - \hat r$. The magnitude yields the eccentricity and the direction points from the central mass to the periapsis point, hence yielding the argument of periapsis. (continued) $\endgroup$ Commented Oct 2, 2020 at 5:59
  • $\begingroup$ (3) Compute the semi-major axis length via the vis viva equation, $v^2 = \mu\left(\frac2r-\frac1a\right)$. (4) The true anomaly is the angle between the eccentricity vector and current location. $\endgroup$ Commented Oct 2, 2020 at 5:59

You must log in to answer this question.