Hufeisenbahnen und Integration in C

Ich werde weder Ihre Gleichungen noch Ihren Code sorgfältig durchgehen, aber da Sie die genauen Startpositionen, die Sie verwendet haben, nicht erwähnt haben, vermute ich, dass Sie nicht zuerst nach einem stabilen Hufeisenbahn-Zustandsvektor gelöst haben .

Die Positionen von Sonne und Erde sind im synodischen (rotierenden) Rahmen fixiert. Aber wo hast du deinen Astroiden angefangen? Für eine bestimmte Startposition müssen Sie dem Asteroiden auch die richtige Startgeschwindigkeit (Geschwindigkeit und Richtung) geben, sonst ist er möglicherweise nicht stabil und wandert einfach ab oder dreht sich um, wie Sie erwähnt haben.

Schauen Sie sich diese Antwort an , in der ich die richtigen Gleichungen zeige und dann nach stabilen Hufeisenbahnen löse, bevor ich sie zeichne.

Ich beginne mit einer Position auf der x-Achse, gegenüber der Erde und etwas weiter oder näher von der Sonne entfernt als die Erde. Dann starte ich es mit geringer Geschwindigkeit (im rotierenden Rahmen) genau in Richtung der y-Achse. Ich beobachte, wie es zur Erde driftet und dann als Bumerang zurück in die Region des Startgebiets. Wenn es die x-Achse wieder kreuzt, überprüfe ich, um zu sehen, wie nahe die Geschwindigkeit an y ist. Ich verwende einen Optimierungs-Nulllöser, brentqum die Anfangsgeschwindigkeit anzupassen, bis die Rückkehrgeschwindigkeit genau in die entgegengesetzte Richtung zeigt.

Wenn das stimmt, dann wird die Umlaufbahn periodisch und sich wiederholend sein, und wir können sie unter den Einschränkungen des CR3BP als abgestandene Umlaufbahn bezeichnen.

Aus dieser Antwort (dort gibt es noch viel mehr zu lesen, einschließlich mehrerer Referenzen!):

$$\ddot{x} = x + 2\dot{y} - \frac{(1-\mu)(x+\mu)}{r_1^3} - \frac{\mu(x-1+\mu)}{r_2^3}$$ $$\ddot{y} = y - 2\dot{x} - \frac{(1-\mu)y}{r_1^3} - \frac{\mu y}{r_2^3}$$ $$\ddot{z} = -\frac{(1-\mu) z}{r_1^3} - \frac{\mu z}{r_2^3}$$

oben: Halbzyklen einiger wackeliger Hufeisenbahnen

oben: Zeiten bis zum ersten Kreuzen der x-Achse der gleichen wackeligen Hufeisenbahnen, die zur Berechnung der Halbzykluszeiten verwendet werden.

oben: Zykluszeiten aus dieser Berechnung (schwarze Punkte) versus aus der synodischen Periodenschätzungsmethode (rote Punkte). Gute qualitative Übereinstimmung. Auch die Startgeschwindigkeiten y an jedem Startpunkt in x.

unten: Python-Skript für diese Plots.

def x_acc(x, ydot):
    r1    = np.abs(x-x1)
    r2    = np.abs(x-x2)
    xddot = x + 2*ydot  -  ((1-mu)/r1**3)*(x+mu) - (mu/r2**3)*(x-(1-mu))
    return xddot

def C_calc(x, y, z, xdot, ydot, zdot):
    r1 = np.sqrt((x-x1)**2 + y**2 + z**2)
    r2 = np.sqrt((x-x2)**2 + y**2 + z**2)
    C = (x**2 + y**2 + 2.*(1-mu)/r1 + 2.*mu/r2 - (xdot**2 + ydot**2 + zdot**2))
    return C

def deriv(X, t): 
    x, y, z, xdot, ydot, zdot = X
    r1 = np.sqrt((x-x1)**2 + y**2 + z**2)
    r2 = np.sqrt((x-x2)**2 + y**2 + z**2)
    xddot = x + 2*ydot  -  ((1-mu)/r1**3)*(x+mu) - (mu/r2**3)*(x-(1-mu))
    yddot = y - 2*xdot  -  ((1-mu)/r1**3)*y      - (mu/r2**3)*y
    zddot =             -  ((1-mu)/r1**3)*z      - (mu/r2**3)*z
    return np.hstack((xdot, ydot, zdot, xddot, yddot, zddot))

# http://cosweb1.fau.edu/~jmirelesjames/hw4Notes.pdf

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint
from scipy.optimize import brentq

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]

mu = 0.001

x1 = -mu
x2 = 1. - mu

x = np.linspace(-1.4, 1.4, 1201)
y = np.linspace(-1.4, 1.4, 1201)

Y, X = np.meshgrid(y, x, indexing='ij')
Z    = np.zeros_like(X)

xdot, ydot, zdot = [np.zeros_like(X) for i in range(3)]

C = C_calc(X, Y, Z, xdot, ydot, zdot)
C[C>8] = np.nan

if True:
    plt.figure()
    plt.imshow(C)
    plt.colorbar()
    levels = np.arange(2.9, 3.2, 0.04) 
    CS = plt.contour(C, levels,
                 origin='lower',
                 linewidths=2) 
    plt.show()

ydot0s   = np.linspace(-0.08, 0.08, 20)
x0ydot0s = []
for ydot0 in ydot0s:
    x0, infob =  brentq(x_acc, -1.5, -0.5, args=(ydot0), xtol=1E-11, rtol=1E-11,
                           maxiter=100, full_output=True, disp=True)
    x0ydot0s.append((x0, ydot0))

states = [np.array([x0, 0, 0, 0, ydot0, 0]) for (x0, ydot0) in x0ydot0s]

times  = np.arange(0, 150, 0.01)

results = []
for X0 in states:
    answer, info = ODEint(deriv, X0, times, atol = 1E-11, full_output=True)
    results.append(answer.T.copy())

resultz = []
for x0ydot0, thing in zip(x0ydot0s, results):
    y     = thing[1]
    check = y[2:]*y[1:-1] < 0
    zc    = np.argmax(y[2:]*y[1:-1] < 0) + 1
    if zc > 10:
        resultz.append((thing, zc, x0ydot0))

if True:
    plt.figure()
    hw = 1.6
    for j, (thing, zc, x0ydot0) in enumerate(resultz):
        x, y = thing[:2,:zc]
        plt.plot(x, y)
    plt.xlim(-hw, hw)
    plt.ylim(-hw, hw)
    plt.plot([x1], [0], 'ok')
    plt.plot([x2], [0], 'ok')
    plt.show()

if True:
    plt.figure()
    for j, (thing, zc, x0ydot0) in enumerate(resultz):
        x, y = thing[:2]
        plt.plot(times[:zc], y[:zc])
    plt.show()

if True:
    plt.figure()
    for j, (thing, zc, x0ydot0) in enumerate(resultz):
        x0, ydot0 = x0ydot0
        cycle_time = 2. * times[zc] / twopi
        ratio = abs(x0/x2)
        T_simple_model = twopi * abs(x0/x2)**1.5
        T_synodic_simple_model = 1. / (1. - twopi/T_simple_model) # https://astronomy.stackexchange.com/a/25002/7982
        plt.subplot(2, 1, 1)
        plt.plot(x0, cycle_time, 'ok')
        plt.plot(x0, abs(T_synodic_simple_model), 'or')
        plt.subplot(2, 1, 2)
        plt.plot(x0, ydot0, 'ok')
    plt.subplot(2, 1, 1)
    plt.xlabel('x0', fontsize=16)
    plt.ylabel('cycle times (periods)', fontsize=16)
    plt.subplot(2, 1, 2)
    plt.xlabel('x0', fontsize=16)
    plt.ylabel('ydot0', fontsize=16)
    plt.show()