Die ausgezeichnete Antwort von @Julio beschreibt einen Flugbahnwinkel und erklärt, dass es sich um den Winkel zwischen der tangentialen Richtung (senkrecht zum radialen Vektor zum Zentralkörper) und dem aktuellen Geschwindigkeitsvektor handelt.
Ich habe zuerst versucht, den Winkel aus diesem Ausdruck herauszubekommen, aber seitdem ist er offensichtlich falsch ist eine gerade Funktion und der Winkel kann aus gehen zu :
Ich habe Orbits für GM integriert ( ) und SMA ( ) der Einheit und Startabstände von 0,2 bis 1,8. Das macht die Periode immer . Wenn ich das Ergebnis meiner Funktion zeichne, erhalte ich zu viele Wackelbewegungen.
Welchen Ausdruck kann ich verwenden, um ausgehend von Zustandsvektoren den korrekten Flugbahnwinkel Gamma zu erhalten?
Überarbeitetes Python für den fehlerhaften Teil wäre wünschenswert, aber für eine Antwort sicherlich nicht erforderlich.
def deriv(X, t):
x, v = X.reshape(2, -1)
acc = -x * ((x**2).sum())**-1.5
return np.hstack((v, acc))
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint
halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
T = twopi
time = np.linspace(0, twopi, 201)
a = 1.0
rstarts = 0.2 * np.arange(1, 10)
vstarts = np.sqrt(2./rstarts - 1./a) # from vis-viva equation
answers = []
for r, v in zip(rstarts, vstarts):
X0 = np.array([r, 0, 0, v])
answer, info = ODEint(deriv, X0, time, full_output= True)
answers.append(answer.T)
gammas = []
for a in answers:
xx, vv = a.reshape(2, 2, -1)
dotted = ((xx*vv)**2).sum(axis=0)
rabs, vabs = [np.sqrt((thing**2).sum(axis=0)) for thing in (xx, vv)]
gamma = np.arccos(dotted/(rabs*vabs)) - halfpi
gammas.append(gamma)
if True:
plt.figure()
plt.subplot(4, 1, 1)
for x, y, vx, vy in answers:
plt.plot(x, y)
plt.plot(x[:1], y[:1], '.k')
plt.plot([0], [0], 'ok')
plt.title('y vs x')
plt.subplot(4, 1, 2)
for x, y, vx, vy in answers:
plt.plot(time, x, '-b')
plt.plot(time, y, '--r')
plt.title('x (blue) y (red, dashed)')
plt.xlim(0, twopi)
plt.subplot(4, 1, 3)
for x, y, vx, vy in answers:
plt.plot(time, vx, '-b')
plt.plot(time, vy, '--r')
plt.title('vx (blue) vy (red), dashed')
plt.xlim(0, twopi)
plt.subplot(4, 1, 4)
for gamma in gammas:
plt.plot(time, gamma)
plt.title('gamma?')
plt.xlim(0, twopi)
plt.show()
Dies ist ein Problem, das Gruppen von Menschen geplagt hat, die sich mit Orbitaldynamik sehr gut auskennen, aber aus verschiedenen Lehrbüchern gelernt haben: Es gibt zwei verschiedene Definitionen von "Flugbahnwinkel"!!
Zusätzlich zu , dem Winkel zwischen der Tangentialrichtung und dem Geschwindigkeitsvektor, besteht , dem Winkel zwischen der radialen Richtung und dem Geschwindigkeitsvektor. Die Leute sagen oft "Flugbahnwinkel", ohne zu sagen, welche Definition sie verwenden . Verwirrend! (Mir ist gerade aufgefallen, dass das Diagramm in Julios Antwort auch zeigt )
Wenn Sie mit arbeiten anstatt , wird von gegeben
die von 0 ("gerade nach oben") bis geht ("gerade nach unten"). Verwenden , "gerade nach oben" ist und "gerade nach unten" ist , also konvertieren zu du subtrahierst einfach von :
Dies entspricht
Ich bin mit der Sprache, die Sie für Ihre Berechnungen und Diagramme verwendet haben, nicht vertraut, daher habe ich mir Ihren Algorithmus nicht angesehen, um zu sehen, warum es "zu viele Wackeln" gibt.
Ich habe den Fehler im Skript gefunden, er lag an meinem "Homebrew" -Punktprodukt. Ich hatte eine zusätzliche Quadrierung:
dotted = ((xx*vv)**2).sum(axis=0) # WRONG
dotted = (xx*vv).sum(axis=0) # Correct
Unter Verwendung dieser und der hervorragenden Erläuterungen von @TomSpilker habe ich die folgenden zwei Methoden zur Berechnung von Gamma verwendet:
Methode 1:
Methode 2:
Eine Brute-Force-Alternativmethode zur doppelten Überprüfung:
Die Modulo-Operation wird nur im Computerprogramm wirklich benötigt, da jedes Theta aus einer separaten arctan2-Operation stammt:
gammas_1, gammas_2 = [], []
for a in answers:
xx, vv = a.reshape(2, 2, -1)
dotted = (xx*vv).sum(axis=0)
rabs, vabs = [np.sqrt((thing**2).sum(axis=0)) for thing in (xx, vv)]
gamma_1 = np.arcsin(dotted/(rabs*vabs)) # Per Tom Spilker's answer Eq. 3
theta_r = np.arctan2(xx[1], xx[0])
theta_v = np.arctan2(vv[1], vv[0])
theta_tanj = theta_r + halfpi
gamma_2 = theta_tanj - theta_v
gamma_2 = np.mod(gamma_2 + pi, twopi) - pi
gammas_1.append(gamma_1)
gammas_2.append(gamma_2)
plt.figure()
plt.subplot(2, 1, 1)
for gamma_1 in gammas_1:
plt.plot(time, gamma_1)
plt.title('gammas_1', fontsize=16)
plt.subplot(2, 1, 2)
for gamma_2 in gammas_2:
plt.plot(time, gamma_2)
plt.title('gammas_2', fontsize=16)
Benutzer20636