Kann Fehler bei der Berechnung des orbitalen Exzentrizitätsvektors nicht identifizieren; Magnitude gleich eins statt null (mit Python-Code)

Ich habe eine Gravitations-Körper-Simulation, für die ich verschiedene Bahnparameter bestimmen möchte. Für jeden Körper habe ich 3-D-Vektoren (x,y,z-Raum) für Position, Geschwindigkeit und Beschleunigung. Ich versuche, die in diesem Beitrag beschriebenen Schritte zu befolgen, um die Exzentrizität jeder Umlaufbahn zu ermitteln. Bevor ich n Körper in die Simulation werfe, teste ich den Algorithmus an einfacheren Systemen, wie z. B. einem 2-Körper-System, in dem die Umlaufbahn der Erde um die Sonne nahezu ein perfekter Kreis ist. Da die Umlaufbahn kreisförmig ist, erwarte ich, dass die Exzentrizität Null ist; Dies ist nicht die Ausgabe, die ich bekomme, also hoffe ich, dass mir jemand helfen kann, meine Fehler zu identifizieren (entweder beim Verständnis oder im Code). Insbesondere würde ich gerne wissen, was ich falsch mache, wenn ich versuche, die Exzentrizität zu berechnen.

Entschuldigung im Voraus für die Länge dieses Beitrags; Der größte Teil des folgenden Codes soll zeigen, dass die Methodik funktioniert, um Positions- und Geschwindigkeitsvektoren zu erhalten. Der letzte Teil des Codes (springen Sie nach unten zu PROBLEM ) besteht darin, "meine Arbeit zu zeigen", indem Sie diese Parameter verwenden, um die Exzentrizitätsvektoren zu berechnen. Abgesehen von der visuellen Inspektion wurden Methoden aus diesem Beitrag verwendet, um sicherzustellen, dass die Umlaufbahn kreisförmig ist.

Erstellen Sie eine kreisförmige Umlaufbahn über das Sonne-Erde-System

Zunächst initialisieren wir die Anfangsbedingungen unserer gekoppelten ODEs und relevanten Simulationsparameter.

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)

Dann lösen wir die gekoppelten ODEs unter Verwendung eines einfachen Euler-Verfahrens.

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

Um zu überprüfen, ob die Simulation wie erwartet gelaufen ist, zeichnen wir.

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

PROBLEM

Ich bin eher damit vertraut, den Drehimpuls ausgedrückt zu sehen als$L = \vec{r} x \vec{p}$, Wo$\vec{p} = m\vec{v}$, obwohl ich annehme, dass man den Drehimpuls unten interpretieren kann, ausgedrückt in Einheiten des Drehimpulses pro Masseneinheit. In kartesischen Koordinaten,$\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)

Nach meinem Verständnis variiert die Exzentrizität$0 ≤ e < 1$für elliptische Bahnen (kreisförmige Bahnen sind$e=0$),$e=1$für Parabelbahnen und$e>1$für hyperbolische Bahnen. Irgendwas muss also stimmen. Muss ich die Koordinaten aus einem bestimmten Bezugssystem berücksichtigen? Oder habe ich vielleicht eine Annahme für die verwendeten Gleichungen übersehen? Kann jemand auf die Ursache dieses Fehlers hinweisen? Weniger wichtig, ist die Gleichung, die zur Berechnung der Exzentrizität verwendet wird, auf alle Umlaufbahnen oder nur auf elliptische Umlaufbahnen verallgemeinerbar?