Wie berechnet man den Flugbahnwinkel γ aus einem Zustandsvektor?

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 arccos ist eine gerade Funktion und der Winkel kann aus gehen π / 2 zu π / 2 :

arccos ( r v | r |   | v | ) π 2        (falsch!)

Ich habe Orbits für GM integriert ( μ ) und SMA ( a ) der Einheit und Startabstände von 0,2 bis 1,8. Das macht die Periode immer 2 π . 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.

Orbit-Plots

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()
sollte diese Frage TLDRed sein, um anzuzeigen, dass es sich um einen Codierungsfehler handelt, da es immer noch zu fragen scheint, was mit der Formel nicht stimmt

Antworten (2)

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

(1) arccos ( r v | r |   | v | )

die von 0 ("gerade nach oben") bis geht π ("gerade nach unten"). Verwenden γ , "gerade nach oben" ist π / 2 und "gerade nach unten" ist π / 2 , also konvertieren β zu γ du subtrahierst einfach β von π / 2 :

(2) γ = π / 2 arccos ( r v | r |   | v | )

Dies entspricht

(3) γ = arcsin ( r v | r |   | v | )

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.

Vielen Dank! Ich habe Gleichungen Tags (Zahlen) hinzugefügt. Würden Sie sagen, dass es zu viele Wackelbewegungen gibt, oder ist dieses Wackelverhalten tatsächlich vernünftig? Da Ihr β (Gl. 1) ist das gleiche wie mein Fehler γ Abgesehen von einem Versatz von einem halben Pi sollten die Wackelbewegungen in meinem Plot die gleichen sein wie in einem richtigen Plot von Ihnen β (Gl. 1).
Sieht für mich nach zu vielen Wackeln aus. Ich werde das später überprüfen.
@uhoh, eigentlich ist meine Gleichung 1 nur das Negative Ihrer Gleichung. Etwas anderes ist falsch. Natürlich wissen Sie, dass etwas nicht stimmt, weil alle geplottet haben γ s sind negativ oder null, was nicht anders sein kann als eine nach innen gerichtete Spirale. Für eine keplersche exzentrische Umlaufbahn γ sollte genau zweimal Null kreuzen, bei Periapsis und Apoapsis, und zwischen den Extrema monoton sein, sowohl für die kurzen (Extremum durch Periapsis zum anderen Extremum) als auch für die langen (Extremum durch Apoapsis zum anderen Extremum) Segmente. Ich werde mal sehen, ob ich ein Beispiel dafür zeichnen kann γ Kurve sollte aussehen.
@uhoh , Oh, ja: in der Tat zu viele Wackeln. Darauf bezieht sich das "monoton". Irgendwie die γ die rechnung ist durcheinander gekommen. Stellen Sie sicher, dass das Argument der arccos das ist, was Sie beabsichtigen.
Okay danke. Ich gehe davon aus, dass jemand in der Lage sein wird, herauszufinden, was tatsächlich falsch ist. Es könnte sich herausstellen, dass ich es bin, ich werde später darauf zurückkommen.
Hoppla, ich hätte oben sagen sollen: "Meine Gl. 2 ist nur das Negativ von Ihnen." Ich sollte mich abmelden und ins Bett gehen!
Danke für all eure Hilfe! Ich habe den Fehler gefunden und habe ein besseres Verständnis dafür β und γ .
Ich habe Sie hier zitiert , sehen Sie sich an, ob ich es richtig gemacht habe oder ob Sie noch etwas hinzuzufügen haben. Vielen Dank!
@uhoh "tangential" zu einer Kugel, die auf der Primärkugel zentriert ist, nicht auf der Umlaufbahn. Persönlich würde ich lieber "Quergeschwindigkeit" sagen, aber mein erster Professor für Orbitaldynamik in Stanford verwendete "tangential".
Das habe ich schon vermutet .

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:

(3) γ 1 = arcsin ( r v | r |   | v | )

Methode 2:

Eine Brute-Force-Alternativmethode zur doppelten Überprüfung:

θ r = arctan 2 ( j , x )

θ v = arctan 2 ( v j , x )

θ t a n j = θ r + π 2

γ 2 = θ t a n j θ v

γ 2 m Ö d = Mod ( γ 2 + π , 2 π ) π

Die Modulo-Operation wird nur im Computerprogramm wirklich benötigt, da jedes Theta aus einer separaten arctan2-Operation stammt:

Geben Sie hier die Bildbeschreibung ein

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)
Tatsächlich das Neue γ Handlung ist, was ich erwartet hatte. Hurra! Guter Spürsinn.