Frage:
Wie kann man Orbitalelemente mithilfe von Positions- / Geschwindigkeitsvektoren programmgesteuert berechnen?
Stu
2013-09-11 16:49:49 UTC
view on stackexchange narkive permalink

Ich möchte eine mechanische Orbital-Software von Grund auf neu erstellen. Ich bin der Meinung, dass dies eine großartige Möglichkeit wäre, die Schritte zu lernen, die erforderlich sind, um verschiedene Kepler-Orbitalelemente eines Objekts zu berechnen, Umlaufbahnen zu zeichnen und vorherzusagen, wo sich das Objekt zu einem späteren Zeitpunkt befinden wird.

Insbesondere möchte ich mit der Berechnung der Keplarian-Elemente beginnen. Die Eingaben, die ich dem Programm geben würde, wären die Positions- und Geschwindigkeitsvektoren zusammen mit einer Zeit. Diese Eingabevektoren sind relativ zum Erdmittelpunkt, daher muss ich möglicherweise auch eine Koordinatenübertragung durchführen, wenn ich einen bestimmten Ort auf der Oberfläche als Referenzpunkt verwenden möchte.

Ich habe den gesehen Mathematik zur Berechnung von Kepler-Orbitalelementen aus diesem Buch, und ich weiß, dass im Laufe der Jahre eine Menge Software entwickelt wurde, um sie zu berechnen, aber es fällt mir schwer, die beiden zu verbinden. Die Mathematik in dem Buch ist etwas verwirrend, und ich denke, es wäre für mich einfacher zu verstehen, wenn ich die Schritte in einer Programmiersprache "ausgeschrieben" sehen würde.

Ich kann R oder Python verwenden. Vektorisierter Code aus den meisten Sprachen Ich sollte in eine dieser beiden Sprachen übersetzen können. Vielen Dank! @Chris Ich werde die Frage noch einmal mit ein wenig mehr Informationen über meine Probleme bei der Übersetzung der Mathematik in Code aktualisieren.
Überprüfen Sie http://orsa.sourceforge.net/ für ihre Lösungen / Methoden. Lange Diskussion hier http://www.physicsforums.com/showthread.php?t=232778. Und für die grundlegende Python-F = ma-Simulation: https://fiftyexamples.readthedocs.org/en/latest/gravity.html
FWIW, wenn Ihr Ziel die Berechnung der zukünftigen Position ist, gibt es normalerweise wenig Grund, von Position und Geschwindigkeitsvektor in Kepler-Elemente umzuwandeln. Berechnen und wenden Sie einfach den Luftwiderstand und die Schwerkraft an Ihrem aktuellen Standort und Ihrer Geschwindigkeit an und integrieren Sie sie in kleinen Schritten nach vorne.
Drei antworten:
#1
+40
user29
2013-09-12 19:02:44 UTC
view on stackexchange narkive permalink

Angesichts der erdzentrierten Positions- und Geschwindigkeitsvektoren $ \ vec {r} $ und $ \ vec {v} $ für die Trägheit (ECI) können Sie direkt nach den klassischen Orbitalelementen $ (a, e, i, \) suchen Omega, \ omega, \ nu) $ wie folgt (Algorithmen zuerst, gefolgt von Pseudocode unten):

Lösen Sie zuerst den Drehimpuls $$ \ vec {h} = \ vec {r} \ times \ vec {v} $$

dann der Knotenvektor $$ \ hat {n} = \ hat {K} \ times \ vec {h} $$, der später verwendet wird.

Der Exzentrizitätsvektor ist dann $$ \ vec {e} = \ frac {(v ^ 2- \ mu / r) \ vec {r} - (\ vec {r} \ cdot \ vec {v}) ) \ vec {v}} {\ mu} $$

und $ e = | \ vec {e} | $.

Die spezifische mechanische Energie ist $$ E = \ frac {v ^ 2} {2} - \ frac {\ mu} {r} $$

Wenn $ e \ neq 1 $, dann $$ a = - \ frac {\ mu} {2E} $$$$ p = a (1-e ^ 2) $$ Andernfalls ist $$ p = \ frac {h ^ 2} {\ mu} $$ $$ a = \ infty $$

Nun ist $$ i = \ cos ^ {- 1} {\ frac {h_K} {h}} $$$$ \ Omega = \ cos ^ {- 1} {\ frac {n_I} {n}} $$$ $ \ omega = \ cos ^ {- 1} {\ frac {\ vec {n} \ cdot \ vec {e}} {ne}} $$$$ \ nu = \ cos ^ {- 1} {\ frac { \ vec {e} \ cdot \ vec {r}} {er}} $$

Und Sie müssen die folgenden Überprüfungen durchführen: Wenn $ n_J<0 $, dann $ \ Omega = 360 ^ {\ circ} - \ Omega $,

Wenn $ e_K<0 $, dann $ \ omega = 360 ^ {\ circ} - \ omega $ und

Wenn $ \ vec {r} \ cdot \ vec {v} <0 $, dann $ \ nu = 360 ^ {\ circ} - \ nu $.

Beachten Sie, dass Sie mit Sicherheit auf Probleme (Singularitäten) stoßen werden Fälle: Kreisbahnen ($ e \ ca. 0 $) und äquatoriale Bahnen ($ i \ ca. 0 $), insbesondere. In diesen Fällen führen Sie normalerweise eine neue, weniger störende Variable ein, z. B. die mittlere Länge oder die wahre Länge des Perigäums.

  h = Kreuz (r, v) nhat = Kreuz ([0 0 1], h) evec = ((mag (v) ^ 2-mu / mag (r)) * r-Punkt (r, v) * v) / mue = mag (evec) Energie = mag (v) ^ 2/2- mu / mag (r) wenn abs (e-1.0) >eps a = -mu / (2 * Energie) p = a * (1-e ^ 2) sonst p = mag (h) ^ 2 / mu a = infi = acos (h (3) / mag (h)) Omega = acos (n (1) / mag (n)), wenn n (2) <0 Omega = 360-Omegaargp = acos (Punkt (n, evec) / (mag () n) * e)) wenn e (3) <0 argp = 360-argpnu = acos (Punkt (evec, r) / (e * mag (r)) wenn Punkt (r, v) <0 nu = 360 - nu  

Hinweis : Dies folgt aus der in Fundamentals of Astrodynamics and Applications von Vallado, 2007, beschriebenen Methode.

Endlich in Python neu erstellt. Danke für die Hilfe! Es war ziemlich einfach und jetzt verstehe ich die Mathematik viel besser.
Ich würde gerne den Beweis für diese Gleichungen sowie ihre Anwendung in einem echten Problem sehen, das eine Reihe von Orbitalelementen ableitet. Sind Orbitalelemente überhaupt notwendig, wenn Sie diese Positions- und Geschwindigkeitsvektoren haben? Es scheint, dass die Vektorrechnung damit einfach genug umgehen kann. Ich bin auch verwirrt von der Tatsache, dass du nicht wenig mu definiert hast. Auch der v-Vektor ist die erste Ableitung des r-Vektors. Der n-Hut ist ein Einheitsvektor, aber Sie haben die Mathematik für die Einheit nicht gezeigt. Was sind die tiefgestellten Werte n und h? Außerdem haben Sie nicht gezeigt, wie Sie eine mittlere Anomalie und eine mittlere Bewegung ableiten können
Was ist mit der Berechnung der mittleren Anomalie?
Gibt es einen Unterschied zwischen nhat und n? Ich bin nicht sicher, was n ist (vorausgesetzt, nhat ist der Knotenvektor), aber es wurde zur Berechnung des Arguments der Periapsis verwendet. Ich nehme an, n ist dasselbe wie nhat und n wurde versehentlich verwendet?
Für die Programmierung würde ich äquivalente Beziehungen mit der atan2-Funktion bevorzugen, da dies automatisch Quadrantenauflösungen ausführt und Sie davor schützt, den Bereich vor den meisten der oben genannten Anwendungen von acos überprüfen zu müssen. (Manchmal führt die numerische Genauigkeit dazu, dass das Argument für eine acos-Funktion geringfügig außerhalb des gültigen Bereichs von -1,00000 bis +1,00000 liegt. In diesem Fall stürzt das gezeigte Programm ab.) Es ist auch möglich, dass h Null ist. Erlauben einer Division durch Null Fehler.
Was für ein Khat in der zweiten Zeile? Ich sehe im Beispielcode, dass es [0 0 1] ist, dass der Normalenvektor der Erdäquatorebene ist?
Ein Problem besteht darin, dass ein Punkt r und v genommen und zu diesem Zeitpunkt in Elemente umgewandelt wird, was nicht die üblichen Mittelwerte in einem zweizeiligen Elementsatz sind. Um einen TLE aus Vektoren zu generieren, wird der Vektor mithilfe der numerischen Integration vorwärts propagiert und dann die TLE-Terme entsprechend angepasst. Wenn Sie einigermaßen hoch sind, können Sie den Drag-Begriff wahrscheinlich ignorieren. Wenn Sie keine besonders hohe Präzision benötigen, können Sie die Erde als Punktmasse modellieren.
Was ist ein? Ich schaue auf [dieses Referenzbild] (https://user-images.githubusercontent.com/1646875/47543388-f8cf1f00-d8af-11e8-8775-a97e44df246f.png).
Was ist der Knotenvektor `n ^`? Was ist K ^? Was ist "μ"? Was ist "r" und wie unterscheidet es sich von "r"? Was ist mit `v` und` v⃗`? Was ist "p"?
Was ist "h_K" und "n_I"?
@MattJessick Wie würden Sie dies mit Atan2 tun?
#2
+6
PearsonArtPhoto
2013-09-11 20:34:20 UTC
view on stackexchange narkive permalink

Der erste Schlüssel, um dies herauszufinden, ist die Richtigkeit Ihres Koordinatensystems. Es gibt zwei häufig verwendete Koordinatensysteme für solche Dinge. Dies sind die Rahmen Earth Earthed Earth Fixed (ECEF) und Earth Centered Interial (ECI). Um Mitternacht reihen sich diese beiden genau aneinander, aber sie unterscheiden sich zu anderen Zeiten, basierend auf der Rotation der Erde. ECEF funktioniert am besten für Dinge auf der Erde (Wenn Sie sich nicht bewegen, sollten Sie eine Geschwindigkeit von 0 haben. ECEF berücksichtigt dies, ECI-Geschwindigkeit lässt Sie sich mit der Erdrotation bewegen), ECI funktioniert am besten für Dinge im Orbit (Orbiting) Objekte interessieren sich nicht für die Rotation der Erde, zumindest ist es der Physik egal). Stellen Sie sicher, dass die Koordinatensysteme korrekt sind!

Okay, Sie haben also eine Position und Geschwindigkeit in ECI-Koordinaten. Was tun Sie? Es gibt ein ausgezeichnetes Papier, das den gesamten Prozess beschreibt. Ich werde die Endformeln hier kopieren. Es gibt auch einige gute Quellen hier, hier und hier. Ich empfehle dringend, sie sorgfältig zu lesen. Die Unsicherheit ist viel schwieriger. Nehmen wir also an, Sie haben ein perfektes Wissen über Geschwindigkeit und Position. Insbesondere sind die 6 klassischen keplarischen Elemente Exzentrizität (e), Neigung (i), rechter Aufstieg des aufsteigenden Knotens ($ \ Omega $), Argument des Perigäums ($ \ omega $), Halb-Hauptachse (a) und Zeit der Perigäumpassage ($ T_O $).

Ich sollte erwähnen, dass ich hauptsächlich der Laplace-Methode zur Orbitalbestimmung folge. Es gibt eine konkurrierende Methode, die als Gauß-Methode bekannt ist. Aber schließlich kam es darauf an, Matlab-Code zu entschlüsseln.

Semi-Major-Achse

$ W_s = \ frac {1} {2} * v ^ 2s - \ text {mus} ./ r; $

$ a = -mus / 2. / W_s $; % Semi-Major-Achse

Exzentrizität

  L = [rs (2, :). * vs (3, :) - rs ( 3, :). * Vs (2, :); ... rs (3, :). * Vs (1, :) - rs (1, :). * Vs (3, :); ... rs (1, :). * Vs (2, :) - rs (2, :). * Vs (1, :)]; % Drehimpuls  

$ p = \ sum {L ^ 2} ./ mus; $% semi-latus rectum

$ e = \ sqrt {1 - p / a}; % Exzentrizität $

Neigung

$ I = atan (\ frac {\ sqrt {L (1, :) ^ 2 + L (2,: ) ^ 2}} {L (3, :)}); $

Argumente des Perizentrums

$ \ omega = atan2 (\ frac {( vs (1, :). * L (2, :) - vs (2, :). * L (1, :)) ./ mus - rs (3, :) ./ r) ./ (e. * sin (I))} {((\ sqrt {L2s}. * vs (3, :)) ./ mus - (L (1, :). * rs (2, :) - L (2, :). * rs (1, :)) ./ (\ sqrt {L2s}. * r)) ./ (e. * sin (I)))} $

Länge des aufsteigenden Knotens

$ \ Omega = atan2 (-L (2, :), L (1, :)); $

Zeit, in der das Perigäum vergeht: stark>

$ T_0 = - (E - e. * sin (E)) ./ \ sqrt {mus. * a. ^ - 3} $

Die Methode von Laplace besteht in der Bestimmung der anfänglichen Umlaufbahn aus Winkelmessungen und (glaube ich) weit über den Rahmen dessen hinaus, wonach OP sucht. Wenn Sie über eine ECI-Position und -Geschwindigkeit verfügen, ist das Abrufen von Kepler-Elementen nur eine einfache [Koordinatentransformation] (http://ccar.colorado.edu/ASEN5070/primers/cart2kep/cart2kep.htm).
@Chris: Ich wusste, dass es einen einfacheren Weg geben musste, um die Transformation durchzuführen ... Seufz.
Was ist mit dem Koordinatensystem in Bezug auf Händigkeit und Orientierung? Soweit ich weiß, verwenden die meisten Guides Z-is-up. Wie würden diese Formeln in einem Y-is-up-System für Linkshänder wie Unity oder einem Y-is-up-System für Rechtshänder wie Godot implementiert?
#3
+5
alexamici
2016-04-08 11:50:18 UTC
view on stackexchange narkive permalink

OrbitalPy verfügt über eine praktische Funktion elements_from_state_vector , die genau das tut:

https://github.com/RazerM/orbital/ blob / 0.7.0 / orbital / utilities.py # L252

Sie können überprüfen, ob die Mathematik mit der Antwort von user29 übereinstimmt.

Hey, das ist ziemlich cool! Ich freue mich darauf, es zu versuchen. Kann OrbitalPy im zweiten Beispiel in den Dokumenten mit dem Titel "* molniya orbit * erstellen" Präzession implementieren? Gibt es einen Platz zum Hinzufügen eines J2? (und ich denke, Molniya sollte groß geschrieben werden - ich denke, es ist ein Eigenname).
Ich bin mit OrbitalPy selbst nicht vertraut, aber aus der begrenzten Analyse des Quellcodes geht hervor, dass es sich um eine reine Zwei-Körper-Kepler-Ausbreitung ohne Störung handelt, also keine Ellipsoidkorrekturen.
OK, danke für die Info - ich werde es versuchen.
Ist dies mit Blick auf ein "Z-is-up" -Koordinatensystem geschrieben? Wenn ich ein "Y-is-up" -Koordinatensystem habe, sollte ich "h.z" mit "h.y" und "[0, 0, 1]" mit "[0, 1, 0]" austauschen?
"Das Argument der Periapsis ist der Winkel zwischen dem Exzentrizitätsvektor und seiner x-Komponente." Warum?


Diese Fragen und Antworten wurden automatisch aus der englischen Sprache übersetzt.Der ursprüngliche Inhalt ist auf stackexchange verfügbar. Wir danken ihm für die cc by-sa 3.0-Lizenz, unter der er vertrieben wird.
Loading...