Bestimmung der Orbitalposition zu einem zukünftigen Zeitpunkt

Problemstellung

Das Problem, das Sie lösen wollen, heißt Kepler-Problem . In Ihrer Formulierung des Problems gehen Sie von den kartesischen Bahnzustandsvektoren (auch kartesische Elemente genannt ) aus, also der Anfangsposition und -geschwindigkeit.

Wie Sie festgestellt haben, ist die einzige Möglichkeit, die kartesischen Elemente zeitlich vorwärts zu propagieren, die numerische Integration. Dies funktioniert in Ordnung, kann jedoch langsam sein, wenn Sie eine hohe Genauigkeit wünschen, und es gibt einige numerische Probleme (durch Rundung verursachte Fehler [sammeln sich langsam und viele Integratoren verursachen eine Energiedrift ). Sie können einige dieser Probleme umgehen, indem Sie einen Integrator höherer Ordnung verwenden ( Runge-Kutta ist ein beliebter), der es Ihnen ermöglicht, größere Schritte für die gleiche Genauigkeit zu machen oder eine bessere Genauigkeit für die gleiche Schrittgröße zu erhalten. Dies ist jedoch für eine einfache Simulation etwas übertrieben.

Wenn Ihre Simulation als Zwei-Körper-Problem behandelt werden kann , vereinfachen sich die Dinge dramatisch. Das Zwei-Körper-Problem ist eine gute Vereinfachung, wenn die Simulationsobjekte hauptsächlich von einem einzigen, großen Objekt beeinflusst werden. Zum Beispiel werden die Erde, die um die Sonne wandert, oder ein Raumschiff, das sich in einer niedrigen Erdumlaufbahn bewegt, gut als Zwei-Körper-Problem modelliert; Ein Raumschiff, das von der Erde zum Mond reist, ist es jedoch nicht (dazu später mehr).

Da Sie versuchen, die Positionen der Planeten mit mittlerer Genauigkeit zu modellieren, sollte die Reduktion auf das Zwei-Körper-Problem für Sie funktionieren.

Begriffsdefinitionen

Die traditionelle Lösung des Zwei-Körper-Problems beinhaltet eine andere Art, die Position des umlaufenden Körpers darzustellen, die sogenannten Kepler-Orbitalelemente (auch nur Orbitalelemente genannt ). Anstatt Position und Geschwindigkeit anzugeben, geben sie sechs verschiedene Parameter der Umlaufbahn an (wenn Sie nur zum Code gelangen möchten, können Sie diesen Teil überspringen):

1. Große Halbachse,$a$: Halber maximaler Durchmesser der Ellipsenbahn, ( = Kreisradius bei kreisförmiger Bahn). Die Energie und Periode der Umlaufbahn hängen nur davon ab$a$. Das Semi-Latus-Rektum $\ell$, die "Breite" der Umlaufbahn, kann eine bessere Wahl für Umlaufbahnen sein, die nahezu parabolisch sind (wie bei Asteroiden) oder die von elliptisch zu hyperbolisch wechseln (wie bei interplanetaren Raumfahrzeugen). Die beiden sind verwandt durch$\ell=a(1-e^2)$.

2. Exzentrizität,$e$: Die "Spitzigkeit* der Umlaufbahn. Reicht von$e=0$für eine perfekt kreisförmige Umlaufbahn, zu$e=1$für eine parabolische Umlaufbahn, zu$e>1$für hyperbolische Bahnen. Merkur ist der exzentrischste Planet mit$e \approx 0.2$. Erdumkreisende Raumfahrzeuge haben normalerweise$e<0.01$.

Abgesehen von$e$und$a$wir können den am weitesten entfernten und den nächsten Punkt in der Umlaufbahn bestimmen, die Apoapsis und Periapsis (zusammen Apsiden ):

Die beiden Parameter$a$und$e$reichen aus, um die Form der Umlaufbahn zu bestimmen. Die nächsten drei Parameter definieren die Orientierung der Umlaufbahn relativ zu einem Koordinatensystem, das aus einer Referenzebene und einer Referenzrichtung (parallel zur Ebene) besteht.

Für fast alle Umlaufbahnen im Sonnensystem ist das verwendete Koordinatensystem das Ekliptik-Koordinatensystem . Die Bezugsebene ist die Ekliptikebene , die Ebene der Umlaufbahn der Erde um die Sonne. Die Bezugsrichtung ist der Frühlingsäquinoktiumspunkt , die Richtung von der Erde zur Sonne im Moment des Frühlingsäquinoktiums. Da diese beiden Referenzen im Laufe der Zeit langsam driften, müssen wir eine bestimmte Zeit angeben, zu der diese Referenzen definiert werden, die als Epoche bezeichnet wird . Die häufigste ist J2000 , Mittag am 1. Januar 2000 (UTC).

Erdzentrierte Umlaufbahnen verwenden oft das äquatoriale Koordinatensystem , dessen Bezugsebene der Äquator der Erde ist. Die Situation mit der Epoche ist etwas kompliziert, deshalb gehe ich hier nicht darauf ein.

Folgende Parameter lokalisieren die Umlaufbahn bezüglich der Erdumlaufbahn:

3. Neigung,$i$: der Winkel zwischen der Bahnebene und der Bezugsebene. Eine Neigung zwischen 90 und 180 Grad bezieht sich auf eine rückläufige Umlaufbahn, die von der üblichen Richtung "rückwärts" umkreist.

4. Längengrad des aufsteigenden Knotens,$\Omega$: Der aufsteigende Knoten ist der Ort, an dem die Umlaufbahn von unterhalb der Referenzebene nach oben kreuzt. (Es ist am Schnittpunkt zwischen der Orbitalebene und der Referenzebene)$\Omega$ist der Winkel zwischen diesem Punkt und der Bezugsrichtung, gemessen gegen den Uhrzeigersinn.

5. Argument der Periapsis,$\omega$: der Winkel zwischen dem aufsteigenden Knoten und der Periapsis (dem tiefsten Punkt der Umlaufbahn). Für Umlaufbahnen mit sehr geringer Neigung, bei denen die Position des aufsteigenden Knotens schwer zu bestimmen ist (da es sich um den Schnittpunkt zwischen zwei fast parallelen Ebenen handelt), verwenden wir stattdessen die Länge der Periapsis $\varpi = \Omega + \omega$.

Der sechste Parameter definiert die Position des Objekts in seiner Umlaufbahn. Es gibt ein paar verschiedene Möglichkeiten, aber die häufigste ist:

• Mittlere Anomalie,$M$: Ein "imaginärer" Winkel, der an der Periapsis Null ist und mit einer konstanten Rate von 360 Grad pro Umlaufbahn zunimmt.

Der Kurs, zu dem$M$Änderungen wird die mittlere Bewegung genannt ,$n$, gleicht$2\pi/T$. Normalerweise haben Sie eine Messung von$M$zu einer bestimmten Epoche$t_0$, genannt (nicht überraschend) die mittlere Anomalie in Epoche ,$M_0$.

Genau wie beim Argument der Periapsis verwenden wir für Umlaufbahnen mit geringer Neigung einen verwandten Wert, den mittleren Längengrad ,$L=\varpi + M$.

Der tatsächliche Winkel zwischen dem umlaufenden Körper und der Periapsis wird als wahre Anomalie bezeichnet.$\nu$. Das ist der Winkel, den wir brauchen, um die Position des Körpers zu berechnen. Eine direkte Berechnung ist leider nicht möglich$\nu$aus$M$. Stattdessen lösen wir zuerst nach der exzentrischen Anomalie $E$:

Dies wird als Kepler-Gleichung bezeichnet und kann nicht analytisch gelöst werden. Sobald wir haben$E$Es gibt jedoch einen relativ einfachen Ausdruck für$\nu$.

Berechnung der Position aus Orbitalelementen

Wir führen diese Berechnung in drei Schritten durch: Zuerst lösen wir die Kepler-Gleichung. Zweitens berechnen wir die 2. Position des Körpers in der Orbitalebene. Zuletzt drehen wir unsere 2D-Position in 3D-Koordinaten. Ich werde für die meisten dieser Aufgaben etwas "Pseudocode" in Javascript geben.

Ich gehe davon aus, dass Sie eine Reihe von Elementen wie diese von der JPL-Website verwenden . Diese verwenden$L$und$\varpi$Anstatt von$M$und$\omega$. Die Tabelle gibt zwei Werte für jedes der Elemente an; die zweite ist die zeitliche Ableitung. Wenn Sie die Werte in dieser Tabelle verwenden, sollten Sie auch die Derivate verwenden.

Berechne die Zeit$t$in Jahrhunderten ab J2000:

// month is zero-indexed, so 0 is January
var tMillisFromJ2000 = Date.now() - Date.UTC(2000, 0, 1, 12, 0, 0);
var tCenturiesFromJ2000 = tMillisFromJ2000 / (1000*60*60*24*365.25*100);

Jetzt berechnen wir die aktuellen Werte jedes Orbitalparameters. Zum Beispiel die große Halbachse der Erde mit den Werten aus Tabelle 1 (gültig von 1800–2500):

// a0 = 1.00000261; adot = 0.00000562
var a = a0 + adot * tCenturiesFromJ2000;

(Beachten Sie, dass die Werte tatsächlich für "EM Barycenter " angegeben sind, den Schwerpunkt des Erde-Mond-Systems. Die Erde befindet sich etwa 4600 Kilometer vom Barycenter entfernt in der entgegengesetzten Richtung zum Mond. Wenn Sie dies korrigieren möchten Ungenauigkeit müssen Sie auch die Bewegung des Mondes simulieren, aber das ist wahrscheinlich übertrieben.)

Tabelle 2a enthält Elemente, die von 3000 v. Chr. bis 3000 n. Chr. Genau sind; wenn Sie jedoch die Elemente aus Tabelle 2a verwenden, müssen Sie diese um Korrekturen ergänzen$L$aus Tabelle 2b! Hier wird zum Beispiel der Längengrad von Saturn berechnet:

// L0 = 34.33479152; Ldot = 3034.90371757
// b = -0.00012452          
// c =  0.06064060
// s = -0.35635438
// f = 38.35125000
var L = L0 + Ldot * tCenturiesFromJ2000
           + b * Math.pow(tCenturiesFromJ2000, 2)
           + c * Math.cos(f * tCenturiesFromJ2000)
           + s * Math.sin(f * tCenturiesFromJ2000);

Wir müssen die mittlere Bewegung nicht explizit berechnen und hinzufügen$L$, da es in beiden Tabellen enthalten ist$\dot L$.

Jetzt sind wir bereit zu berechnen$M$und$\omega$( w):

var M = L - p \\ p is the longitude of periapsis
var w = p - W \\ W is the longitude of the ascending node

Weiter zu Schritt 2: Wir müssen die Kepler-Gleichung lösen:

Wir können dies numerisch mit dem Newton-Verfahren lösen . Das Lösen der Kepler-Gleichung entspricht dem Finden der Wurzeln von$f(E) = E - e \sin E - M$. Gegeben$E_i$, eine Schätzung von$E$, können wir die Newton-Methode verwenden, um eine bessere Schätzung zu finden:

Da der nichtlineare Teil$e \sin E$sehr klein ist, können wir mit der Schätzung beginnen$E=M$. Unser Code sieht in etwa so aus:

E = M;
while(true) {
  var dE = (E - e * Math.sin(E) - M)/(1 - e * Math.cos(E));
  E -= dE;
  if( Math.abs(dE) < 1e-6 ) break;
}

Nun gibt es zwei Möglichkeiten, die Position aus der exzentrischen Anomalie zu berechnen. Wir können zuerst die wahre Anomalie und den Radius (die Position des Objekts in Polarkoordinaten) berechnen und dann in rechtwinklige Koordinaten umwandeln; Wenn wir jedoch ein wenig Geometrie anwenden, können wir stattdessen die Koordinaten direkt daraus berechnen$E$:

var P = a * (Math.cos(E) - e);
var Q = a * Math.sin(E) * Math.sqrt(1 - Math.pow(e, 2));

( Pund Qbilden ein 2D-Koordinatensystem in der Ebene der Umlaufbahn, +Pdas in Richtung Periapsis zeigt.)

Schließlich können wir diese Koordinaten in das vollständige 3D-Koordinatensystem drehen:

// rotate by argument of periapsis
var x = Math.cos(w) * P - Math.sin(w) * Q;
var y = Math.sin(w) * P + Math.cos(w) * Q;
// rotate by inclination
var z = Math.sin(i) * y;
    y = Math.cos(i) * y;
// rotate by longitude of ascending node
var xtemp = x;
x = Math.cos(W) * xtemp - Math.sin(W) * y;
y = Math.sin(W) * xtemp + Math.cos(W) * y;

( x, y, und zwerden in Einheiten von AE angegeben.)

Und du bist fertig!

Einige Tipps:

• Wenn Sie auch die Geschwindigkeit berechnen möchten, können Sie dies gleichzeitig mit der Berechnung tun$P$und$Q$, dann drehen Sie es auf die gleiche Weise.

  $$\dot M = n = \dot L \\ \dot M = \dot E - e (\cos E) \dot E \\ \dot E = \dot M / (1 - e \cos E) \\ \dot P = -a (\sin E) \dot E \qquad \dot Q = a (\cos E) \dot E \sqrt{1 - e^2}$$ Beachten Sie, dass ich keine der Derivate einschließe (außer$\dot L$) in dieser Berechnung, da sie das Ergebnis nicht stark beeinflussen. Sie könnten dies wie folgt codieren:

   var vP = - a * Math.sin(E) * Ldot / (1 - e * Math.cos(E));
   var vQ = a * Math.cos(E) * Math.sqrt(1 - e*e) * Ldot / (1 - e * Math.cos(E));
  

  Beachten Sie, dass die Geschwindigkeiten in AE pro Jahrhundert angegeben werden.

• Wenn Sie die Positionen sehr häufig aktualisieren, können Sie den vorherigen Wert von verwenden$E$um Newtons Methode zu impfen und eine feste Anzahl von Iterationen durchzuführen (wahrscheinlich würde nur eine ausreichen). Beachten Sie jedoch, dass Sie diesen Wert beibehalten müssen$E$lokal für jedes Objekt!

• Sie können auch einfach eine feste Anzahl von Iterationen für die anfängliche Lösung verwenden. Sogar für$e=0.2$, nach drei Iterationen der Fehler in$E$geht es nur um$10^{-13}$, und nach vier Iterationen ist der Fehler kleiner als der Rundungsfehler einer IEEE-Verdopplung$e=0.42$.

Wenn Sie weitere Informationen wünschen, können Sie online suchen, aber wenn Sie wirklich interessiert sind, sollten Sie einen Einführungstext zur Orbitalmechanik lesen. Ich persönlich empfehle Fundamentals of Astrodynamics von Bate, Mueller und White (pdf) . Mein Vater benutzte dieses Buch, als er auf dem College war, und ich fand es besser lesbar als mein College-Lehrbuch. Sie interessieren sich für Kapitel 4, Position und Geschwindigkeit als Funktion der Zeit.

Da es nur ein Spiel ist, wären Sie mit kreisförmigen Umlaufbahnen und den Umlaufbahnen der Planeten zufrieden, die nur vom Zentralkörper beeinflusst werden? In diesem Fall ist die Vermehrung recht einfach. In der Bahnebene mit dem Zentralkörper bei (0,0) ist die Position als Funktion der Zeit:

wo$a$ist die große Halbachse oder in diesem Fall wirklich nur der Umlaufbahnradius,$T$ist die Umlaufzeit, und$t_0$bestimmt die Phasenlage der Umlaufbahn, wo an$t=t_0$, der Planet ist auf der x-Achse auf der positiven Seite.

Um die Umlaufbahnen der verschiedenen Planeten miteinander in Einklang zu bringen, müssen Sie nur die definieren$GM$des zentralen Organs, das wir anrufen werden$\mu$. Dann für jeden Bahnradius$a$, auf die die Umlaufzeit bezogen ist$a$durch:

Während es bereits seit Jahren eine qualitativ hochwertige akzeptierte Antwort gibt, finden Sie hier einige zusätzliche Hintergrundinformationen, einige besonders hilfreiche Ressourcen und zusätzliche Tipps für die erstmalige Ausbreitung der Umlaufbahn.

Wenn Sie keine N-Körper-Physik betreiben, also die Planeten nicht interagieren, können Sie analytische Lösungen für das Kepler-Problem verwenden. Irgendwann werden Sie feststellen, dass Sie irgendwann auch hyperbolische Bahnen lösen müssen. Das führt Sie zu universellen Variablenformulierungen zur Lösung des Kepler-Problems.

Die besten Lösungen dafür werden wahrscheinlich die Methode von Goodyear sein:

W. Goodyear, „Vollständig allgemeine Lösung in geschlossener Form für Koordinaten und partielle Ableitungen des Zweikörperproblems“, The Astronomical Journal, Vol. 3, No. 70, Nr. 3, 1965, S. 189–192 (oder das NASA NTRS TD-Dokument zum gleichen Material )

Shepherds Methode:

Shepperd, SW Celestial Mechanics (1985) 35: 129. https://doi.org/10.1007/BF01227666

Oder Danby-Stumpff:

Danby, JMA Celestial Mechanics (1987) 40: 303. https://doi.org/10.1007/BF01235847

Es gibt hier einigen MATLAB-Code, der nützlich (und weitaus zugänglicher) sein könnte, obwohl zufällige Code-Schnipsel auf matlabcentral bei weitem nicht garantiert fehlerfrei sind und es so aussieht, als ob diesem Code eine nützliche Normalisierung seiner Eingaben fehlt (im Allgemeinen gehen Sie auf die Skala Ihres Problems normalisieren wollen, so dass Sie in Einheiten rechnen, in denen r0-bar = 1,0 und mu-bar = 1,0 und wobei v-bar = 1 die Geschwindigkeit in einer kreisförmigen Umlaufbahn bei r0 oder so ähnlich ist) .

Wenn Sie eine N-Körper-Integration der Planetenbewegung durchführen, müssen Sie meiner Meinung nach die numerische Integration verwenden. Runge-Kutta verstößt gegen die Energieerhaltung, sodass Sie wahrscheinlich Symplektische Integration verwenden möchten . Der symplektische Integrator 4. Ordnung in diesem Artikel ist nicht so schwer zu codieren - obwohl Sie dadurch die Schwierigkeit haben, den richtigen Zeitschritt zu erraten (wiederum hilft die Normalisierung, da eine kreisförmige Planetenbahn und ein kreisförmiger LEO dasselbe Problem sind, nur mit unterschiedlichen Entfernungsskalen ) und mit Interpolation der inneren Punkte (und Sie müssen auf das Runge-Phänomen achten , aber damit habe ich nicht gerungen, also weiß ich nicht, welchen Ansatz ich da nehmen soll).

Wenn Sie Runge-Kutta verwenden, dann ist Dormand-Prince mit dynamischer Schrittseite und seinem Interpolanten 3. Ordnung sehr praktisch und wird von Matlab in seinem ode45-Solver verwendet.

Ich würde wahrscheinlich raten, mit der einfachsten Runge-Kutta-Implementierung zu beginnen, basierend auf der einfachen Codierung, aber wenn Sie Runge-Kutta bei jedem Physik-Tick machen, um einen Schritt vorwärts zu kommen, dann ist das ziemlich brutal und die Fehler summieren sich schließlich. aber man könnte es so prototypisieren. Irgendwann möchten Sie zu einem System gehen, in dem Sie das Problem für viele Zeitschritte in die Zukunft lösen, und dann verwenden Sie eine Interpolationsfunktion, um die Lösung in dazwischenliegenden Zeitschritten abzugreifen (was der Punkt ist, an dem ich Dormand erwähne). Prince und seine interpolierende Funktion).