No se puede identificar el error al calcular el vector de excentricidad orbital; magnitud es igual a uno en lugar de cero (con código python)

Tengo una simulación de un cuerpo gravitacional, para la cual me gustaría determinar varios parámetros orbitales. Para cada cuerpo, tengo vectores tridimensionales (espacio x, y, z) para la posición, la velocidad y la aceleración. Estoy tratando de seguir los pasos establecidos en esta publicación para obtener la excentricidad de cada órbita. Antes de incluir n cuerpos en la simulación, estoy probando el algoritmo en sistemas más simples, como un sistema de 2 cuerpos en el que la trayectoria orbital de la Tierra alrededor del Sol es casi un círculo perfecto. Como la órbita es circular, espero que la excentricidad sea cero; este no es el resultado que obtengo, por lo que espero que alguien pueda ayudarme a identificar mis errores (ya sea en la comprensión o en el código). Específicamente, me gustaría saber qué estoy haciendo mal al tratar de calcular la excentricidad.

Lo siento de antemano por la longitud de esta publicación; la mayor parte del código a continuación es para mostrar que la metodología funciona para obtener vectores de posición y velocidad; la última parte del código (saltar al PROBLEMA ) es "mostrar mi trabajo" al usar estos parámetros para calcular los vectores de excentricidad. Además de la inspección visual, se utilizaron métodos de esta publicación para garantizar que la órbita sea circular.

Crea una órbita circular a través del sistema Sol-Tierra

Primero, inicializaremos las condiciones iniciales de nuestras ODE acopladas y los parámetros de simulación relevantes.

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)

Luego, resolvemos las EDO acopladas utilizando un método de Euler simple.

## 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.append(dacc)
    ak = np.array(ak)
    vk = vi + ak * dt
    xk = xi + vk * dt
    X[:, :, ti] = xk
    V[:, :, ti] = vk
    xi, vi = xk, vk

## GET POSITION VECTORS PER BODY
Xs = X[0, :, :]
Xe = X[1, :, :]

## GET VELOCITY VECTORS PER BODY
Vs = V[0, :, :]
Ve = V[1, :, :]

Para verificar que la simulación se ejecutó como se esperaba, graficamos.

## VERIFY -- SHOW POSITION VECTORS
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_aspect('equal')
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)
plt.show()
plt.close(fig)

PROBLEMA

Estoy más familiarizado con ver el momento angular expresado como$L = \vec{r} x \vec{p}$, dónde$\vec{p} = m\vec{v}$, aunque supongo que uno puede interpretar el momento angular a continuación expresado en unidades de momento angular por unidad de masa. En coordenadas cartesianas,$\vec{r} = \vec{x} + \vec{y} + \vec{z} = x\hat{x} + y\hat{y} + z\hat{z}$.

## GET ANGULAR MOMENTUM VECTORS PER BODY
Le = np.cross(Xe, Ve, axis=0)
Ls = np.cross(Xs, Vs, axis=0)

## GET ORBITAL ECCENTRICITY PER BODY
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))

## VERIFY -- SHOW ORBITAL ECCENTRICITY VECTORS PER BODY
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)
plt.show()
plt.close(fig)

## VERIFY -- SHOW ORBITAL ECCENTRICITY MAGNITUDES PER BODY
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)
plt.show()
plt.close(fig)

Tengo entendido que la excentricidad varía según$0 ≤ e < 1$para órbitas elípticas (siendo las órbitas circulares$e=0$),$e=1$para órbitas parabólicas, y$e>1$para órbitas hiperbólicas. Así que algo debe estar mal. ¿Necesito considerar las coordenadas de un marco de referencia específico? ¿O tal vez me perdí una suposición para las ecuaciones utilizadas? ¿Alguien puede señalar la causa de este error? Menos importante, ¿la ecuación utilizada para calcular la excentricidad es generalizable a todas las órbitas o solo a las elípticas?