Невозможно определить ошибку в вычислении вектора эксцентриситета орбиты; величина равна единице вместо нуля (с кодом Python)

У меня есть гравитационное моделирование тела, для которого я хотел бы определить различные параметры орбиты. Для каждого тела у меня есть трехмерные векторы (x,y,z-пространство) для положения, скорости и ускорения. Я пытаюсь выполнить шаги, изложенные в этом посте, чтобы получить эксцентриситет каждой орбиты. Прежде чем добавить в симуляцию n тел, я проверяю алгоритм на более простых системах, таких как система из двух тел, в которой орбитальный путь Земли вокруг Солнца представляет собой почти идеальный круг. Поскольку орбита круговая, я ожидаю, что эксцентриситет будет равен нулю; это не то, что я получаю, поэтому я надеюсь, что кто-нибудь поможет мне определить мои ошибки (либо в понимании, либо в коде). В частности, я хотел бы знать, что я делаю неправильно, пытаясь вычислить эксцентриситет.

Заранее извините за длину этого поста; большая часть кода ниже предназначена для демонстрации того, что методология работает для получения векторов положения и скорости; последняя часть кода (пропустить ПРОБЛЕМУ ) должна «показать мою работу» по использованию этих параметров для вычисления векторов эксцентриситета. Помимо визуального осмотра, методы с этого поста использовались для обеспечения того, чтобы орбита была круговой.

Создать круговую орбиту через систему Солнце-Земля

Во-первых, мы инициализируем начальные условия наших связанных ОДУ и соответствующие параметры моделирования.

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)

Затем мы решаем связанные ОДУ, используя простой метод Эйлера.

## 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, :, :]

Чтобы убедиться, что симуляция прошла так, как ожидалось, мы строим график.

## 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)

ПРОБЛЕМА

Я больше знаком с угловым моментом, выраженным как$L = \vec{r} x \vec{p}$, где$\vec{p} = m\vec{v}$, хотя я полагаю, что можно интерпретировать приведенный ниже угловой момент, выраженный в единицах углового момента на единицу массы. В декартовых координатах,$\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)

Насколько я понимаю, эксцентриситет меняется в зависимости от$0 ≤ e < 1$для эллиптических орбит (круговые орбиты$e=0$),$e=1$для параболических орбит и$e>1$для гиперболических орбит. Значит что-то должно быть не так. Нужно ли считать координаты из конкретной системы отсчета? Или, может быть, я пропустил предположение для уравнений, которые раньше выполнялись? Может ли кто-нибудь указать на причину этой ошибки? Что менее важно, можно ли обобщить уравнение, используемое для расчета эксцентриситета, на все орбиты или только на эллиптические?