Himmelsmechanik: Bahnelemente und Ephemeriden

11. Mai 2015   Allgemein, Astronomie, Himmelsmechanik, Mathematik, Technik   No comments

Himmelsmechanik_4

In den vergangenen drei Teilen dieser Serie haben wir uns mit den Koordinatensystemen und mit der Herleitung der Keplergesetze beschäftigt. Dies ist nun das Rüstzeug um die Bahnelemente und die Ephemeriden von Planeten oder auch von anderen Objekten zu bestimmen. Dazu schauen wir uns mal konkret die Stellung der Venus am 11.06.2016 um 13:30 Uhr an. Dazu wurde eigens ein Matlab-Skript zum besseren Nachvollziehen der Ergebnisse erstellt, welches heruntergeladen werden kann.

Die keplerschen Bahnelemente

Definition aller Bahnelemente

Definition aller Bahnelemente

In der Abbildung sind folgende Bahnelemente definiert:

A Aphel (sonnenfernste Punkt)
P Perihel (sonnennächster Punkt)
M Mittelpunkt der Ellipse
a große Halbachse
e numerische Exzentrizität
i Inklination der Bahnebene zur Referenzebene
r Radiusvektor
\varphi  Winkel zwischen Perihel und Radiusvektor (wahre Anomalie)
☊ aufsteigender Knotenpunkt (Schnittpunkt zwischen Bahn- und Referenzebene)
☋ absteigender Knotenpunkt (Schnittpunkt zwischen Bahn- und Referenzebene)
\Omega  Winkel zwischen Frühlingspunkt und Schnittfläche der Referenzebene
\varphi Perihellänge
\vernal  Frühlingspunkt
\omega Argument der Periapsis

Die Referenzebene kann jede Ebene sein, wir gehen in den folgenden Überlegungen von der Ekliptikebene aus.

Bestimmung des Ortes

Um nun die Position eines Planeten oder eines Objektes auf der elliptischen Bahn zu bestimmen, sind einige geometrische Überlegungen erforderlich.

Definition der mittleren und wahren Anomalie

Definition der mittleren und wahren Anomalie

Zuerst zerlegen wir nun den Radiusvekor \vec{r} in seine x- und y-Komponenten. Es ergibt sich:

(1)   \begin{align*} \vec{r}_{x}&=r\cos \varphi\\ \vec{r}_{y}&=r\sin \varphi\\ \end{align*}

Wenden wir nun einen Trick an. Um die Ellipsenbahn denken wir uns einen Hilfskreis mit dem Radius a. So können wir alle Berechnungen als Kreisbewegungen annehmen. Das fiktive Objekt auf dieser Kreisbahn wird auch mittleres Objekt O' genannt.  Die exzentrische Anomalie ist der Winkel vom Mittelpunkt der Ellipse bzw. des Kreises zu O'. Betrachtet man die Ellipsengleichung als Funktion der exzentrischen Anomalie, so ergibt sich folgender Zusammenhang:

(2)   \begin{equation*} \vec{r}=(a\cos E-e)\cdot \vec{I}_{x}+\dfrac{b}{a} \cdot a\sin E \cdot \vec{I}_{y} \end{equation*}

Den Stauchungsfaktor \frac{b}{a} vereinfachen wir:

(3)   \begin{equation*} b=a \cdot \sqrt{1-e^{2}}\\ \end{equation*}

Es lässt sich nun der Radiusvektor bestimmen:

(4)   \begin{equation*} |\vec{r}|=r=a (1-e \cos E) \end{equation*}

Einfacher ist es, wenn man den Radiusvektor in seine Komponente zerlegt. Dann resultiert durch den Satz des Pythagoras der Betrag des Radiusvektors:

(5)   \begin{align*} \vec {r}_{x}&=(a\cdot \cos E-e)\cdot \vec{I_{x}}\\ \vec {r}_{y}&=a \cdot \sqrt{1-e^{2}} \cdot \sin E \cdot \vec{I_{y}}\\ r&=\sqrt{{r}_{x}^{2}+{r}_{y}^{2}}\\ \end{align*}

Die Vektoren \vec{I}_{x} und \vec{I}_{y} sind jeweils die Einheitsvektoren.

Die hervorgehobene Fläche ist gleich dem Stauchungsfaktor der Fläche SO'P

Die hervorgehobene Fläche ist gleich dem Stauchungsfaktor der Fläche SO’P

Aus der  Gleichung (2) entsteht ein neuer Winkel E. Dieser beschreibt den Winkel von der x-Achse zu dem fiktiven Objekt O', welche man auch die exzentrische Anomalie nennt. Da wir dieses als Kreis annehmen, vereinfachen sich die Berechnungsgänge wesentlich. Betrachten wir uns nochmal den Betrag des Richtungsvektores |\vec{r}|:

(6)   \begin{align*} |\vec{r}|=a (1-e \cos E) \end{align*}

Um nun für E zu einem bestimmten Zeitpunkt zu bestimmen, wenden wir das zweite Keplersche Gesetz an, welches besagt, das zu gleichen Zeiten gleiche Flächen überschritten werden. In der obigen Abbildung ist die hervorgehobene Fläche:

(7)   \begin{align*} A= \dfrac{\pi ab (t-\tau)}{t_{\mathrm{U}}} \end{align*}

wobei t-\tau die verstrichene Zeit seit dem Periheldurchgang  und t_{\mathrm{U}} die Umlaufzeit darstellt. Die Fläche SO'P ist um den Faktor \frac{b}{a} gestaucht. Subtrahieren mit der Dreiecksfläche SO'P bringt:

(8)   \begin{align*} A&=\dfrac{b}{a} \times A_{SO'P} - A_{\bigtriangleup, SO'P }\\ A&=\dfrac{1}{2} ab (E - e \sin E) \end{align*}

Nach Gleichsetzten von (7) und (8)  entsteht die Keplergleichung:

(9)   \begin{align*} E-e\sin E = M(t) \end{align*}

 wobei M(t)=\frac{2\pi (t-\tau) }{t_\mathrm{{U}}} die mittlere Anomalie ist.

Die Lösung der Keplergleichung

Die Gleichung (9) ist nicht analytisch lösbar, da diese Gleichung transzendent ist. Doch wie löst man nun effektiv diese Gleichung? Hier soll uns die numerische Mathematik einen Lösungsansatz liefern. Im wesentlichen sind hier das Iterationsverfahren und die Reihenentwicklung zu nennen. Das Iterationsverfahren eignet sich für kleine numerischen Exzentrizitäten bis zu einem Wert 0,5 sehr gut. Ist die numerische Exzentrizität größer als 0,5, ist das Konvergenzverhalten sehr langsam. Im folgenden wollen wir nun das Iterationsverfahren mal genauer ansehen:

(10)   \begin{align*} E_{0}&=M(t)\\ E_{1}&=M(t)+ e\cdot \sin E_{0};\, |E_{1}-E_{0}|=\Delta_{1}\\ E_{2}&=M(t)+ e\cdot \sin E_{1};\, |E_{2}-E_{1}|=\Delta_{2}\\ ...\\ E_{n+1}&=M(t)+ e\cdot \sin E_{n};\, |E_{n+1}-E_{n}|=\Delta_{n+1}\\ \end{align*}

Die Iteration wird abgebrochen, wenn \Delta_{n+1} kleiner als der vorgebende Wert ist.

Für e>0,5 wählt man einen ähnlichen Ansatz:

(11)   \begin{align*} E_{0}&=M(t)\\ E_{1}&=E_{0}+ \dfrac{M(t)+e\cdot \sin E_{0} -E_{0}}{1-e\cdot \cos E_{0}};\, |E_{1}-E_{0}|=\Delta_{1}\\ E_{2}&=E_{1}+\dfrac{M(t)+e\cdot \sin E_{1} -E_{1}}{1-e\cdot \cos E_{1}};\, |E_{1}-E_{0}|=\Delta_{2}\\ ...\\ E_{n+1}&=E_{n}+ \dfrac{M(t)+e\cdot \sin E_{n} -E_{n}}{1-e\cdot \cos E_{n}};\, |E_{n+1}-E_{n}|=\Delta_{n+1}\\ \end{align*}

Beispielprogramm mit Matlab

Das folgende Listing wurde mit Matlab geschrieben und als Funktion kepler deklariert. Die Eingangsvariablen ist die mittlere Anomalie M(t)=M und die numerische Exzentrizität e.  Ist \Delta_{n+1}\le 10^{-20}  oder 150 Iterationsschritte erledigt, so wird die Schleife beendet. Die Funktion gibt als Ausgangsvariable den Winkel E aus. Am Ende dieses Beitrages können alle Matlab-Files heruntergeladen werden.

%% Die Funktion ‚kepler‘ berechnet die
%% Keplergleichung E=M+sin(E0) mit Hilfe
%% des Newton-Näherungsverfahrens.
%% skywatch-blog.de

function E=kepler(M,e)
E0=M;
    for i=1:150;
        E=M+e.*sin(E0);
            if abs(E-E0)<=10^-20
                break    
            end
        E0=E;
    end
return

Berechnung der wahren Anomalie

Um von den Winkel \varphi, also den Winkel zwischen Perihel und den Aufenthaltsort auf der Umlaufbahn berechnen zu können, müssen wir vom Hilfskreis wieder in die die Ellipse transformieren.

(12)   \begin{align*} r\cos\varphi&=a\cos E - ae\\ \cos\varphi&=\dfrac{a \cos E-e}{a (1-e \cos E)} \end{align*}

Eine alternative Formel  welche wir im folgenden Beispiel verwenden werden, lautet:

(13)   \begin{align*} \tan \dfrac{\varphi}{2}=\sqrt{\dfrac{1+e}{1-e}}\tan \dfrac {E}{2} \end{align*}

Nun haben wir das Rüstzeug zusammen um z.B. die Ephemeriden eines Planeten zu berechnen. Das folgende Verfahren kann auf beliebige Himmelsobjekte angewendet werden.

Ephemeridenrechnung: Beispiel am Planeten Venus

Welche Position und welchen Abstand wird die Venus am 11.06.2016 um 13:30 Uhr haben? Bild: NASA

Welche Position und welchen Abstand wird die Venus am 11.06.2016 um 13:30 Uhr haben? Bild: NASA

Im folgenden wollen wir die Stellung der Venus am 11.06.2016 um 13:30 Uhr sowie den Abstand von der Erde berechnen. Als Referenz dienen die Ergebnisse  aus Stellarium. Die Bahnelemente können aus astronomischen Jahrbüchern oder alternativ aus der unten verlinkte PDF-Datei entnommen werden. Die folgenden Elemente stammen aus dem Buch von Oliver Montenbruck: Grundlagen der Ephemeridenrechnung. Es gilt die Epoche J2000.0. Die Bahnelemente der Venus lauten:

(14)   \begin{align*} a_{\Venus}&=0,723332\,\mathrm{AE}\\ i_{\Venus}&=3,3946^\circ-0.000048^\circ\cdot T\\ M_{\Venus}&=50,4071^\circ+58517,8039^\circ\cdot T\\ \varpi_{\Venus}&=131,5718^\circ+0,0110^\circ \cdot T\\ e_{\Venus}&=0,006773-0,000048 \cdot T\\ \varOmega_{\Venus}&=76,680^\circ-0,278^\circ\cdot T \end{align*}

Die Variable \varpi ist die Länge des Perihels ausgehend vom Frühlingspunkt bis zum Perihel. Es gilt daher:

(15)   \begin{align*} \varpi=\omega + \Omega \end{align*}

Die Variable T ist das julianische Jahrhundert. Dieses berechnet sich aus:

(16)   \begin{align*} T=\dfrac{t}{36525} \end{align*}

mit t  die Anzahl der Tage, die seit dem 01.01.2000 vergangen sind.

Schritt 1: Berechnung des julianischen Jahrhundert

Nach (16) berechnet sich für den 11.06.2016 13:30 Uhr  ein  Wert von T=0,1644.

Schritt 2: Berechnung des Winkel \omega

Der Winkel \omega beschreibt denjenigen Winkel zwischen dem aufsteigenden Knotens und Perihels. Stellt man Gleichung (15) nach \omega um, so ergibt sich:

(17)   \begin{align*} \omega_{\Venus}&=\varpi_{\Venus}-\varOmega_{\Venus}\\ \omega_{\Venus}&=54,9393^\circ \end{align*}

Schritt 3: Lösen der Keplergleichung

Wir lösen nun die Keplergleichung nach dem Verfahren wie in (10) beschrieben. Wichtig ist es hierbei, dass alle Werte ins Bogenmaß transformiert werden müssen. Es ergibt sich hier einen Wert für E= 9672,90^\circ. Der hohe Winkel resultiert aus den Bahnelementen. Möchte man den wirklichen Winkel berechnen, so kann dieser mit modulo 360 ermittelt werden.

Schritt 4:  Nun können wir den Radiusvektor nach Gleichung (5) und die wahre Anomalie nach Gleichung (13) berechnen. Es ergeben sich folgende Werte:

(18)   \begin{align*} r_{\Venus, x}&=0,4875\,\mathrm{AE}\\ r_{\Venus, y}&=-0,5299\,\mathrm{AE}\\ r_{\Venus}&=0,7200\,\mathrm{AE}\\ \varphi_{\Venus}&= - 47,3835^\circ\\ \end{align*}

Schritt 5: Nun haben wir die Koordinaten der Venus. An den Vorzeichen erkennt man, im welchen Quadranten wir uns befinden. Als nächstes müssen wir die heliozentrischen ekliptikalen Koordinaten berechnen. Dabei soll uns der Ansatz von Standish (Link, siehe unten) weiterhelfen:

(19)   \begin{align*} \begin{pmatrix}x_{p}\\y_{p}\\z_{p}\end{pmatrix}&=\begin{pmatrix}\cos \omega \cos \varOmega - \sin \omega \sin \varOmega \cos i\\ \cos \omega \sin \varOmega +\sin \omega \cos \varOmega \cos i\\ \sin \omega \sin i \end{pmatrix}r_{x}\\&+\begin{pmatrix}-\sin \omega \cos \varOmega - \cos \omega \sin \varOmega \cos i\\ - \sin \omega \sin \varOmega + \cos \omega \cos \varOmega \cos i\\ \cos \omega \sin i \end{pmatrix}r_{y} \end{align*}

Es errechnen sich für die Venus folgende Koordinaten:

(20)   \begin{align*} \begin{pmatrix}x_{\Venus}\\ y_{\Venus}\\ z_{\Venus} \end {pmatrix}=\begin{pmatrix}0,0730\\ 0,7163\\ 0,0056 \end{pmatrix}\,\mathrm{AE} \end{align*}

Nun haben wir den ersten Teil der Berechnung erledigt.

Schritt 6: Nun muss nach dem gleichen Verfahren, die Erdkoordinaten berechnet werden. Es ergeben sich mit den Bahnelementen der Erde:

(21)   \begin{align*} a_{\Earth}&=1,000000031\,\mathrm{AE}\\ i_{\Earth}&=0^\circ+0,0131^\circ\cdot T\\ M_{\Earth}&=357,5256^\circ+35999,0498^\circ \cdot T\\ \varpi_{\Earth}&=102,9400^\circ+0,3222^\circ\cdot T\\ e_{\Earth}&=0,016709^\circ-0,000042^\circ\cdot T\\ \varOmega_{\Earth}&=174,876^\circ-0,242^\circ\cdot T \end{align*}

folgende Ergebnisse:

(22)   \begin{align*} \omega_{\mathrm{\Earth}}&=-71,8432^\circ\\ r_{\Earth, x}&= -0,9379\,\mathrm{AE}\\ r_{\Earth, y}&= 0,3890\,\mathrm{AE}\\ r_{\Earth}&=1,0154\,\mathrm{AE}\\ \varphi_{\Earth}&=157,4722^\circ \end{align*}

und folgende Koordinaten:

(23)   \begin{align*} \begin{pmatrix}x_{\Earth}\\y_{\Earth}\\z_{\Earth} \end{pmatrix}=\begin{pmatrix}-0,1682\\-1,0014\\3,8064\cdot10^{-5} \end{pmatrix}\,\mathrm{AE} \end{align*}

Schritt 7: Nun müssen die geozentrisch orthogonale Koordinaten in geozentrisch orthogonale sphärische Ekliptikalkoordinaten transformieren:

(24)   \begin{align*} \begin{pmatrix}X\\Y\\Z\end{pmatrix}=\begin{pmatrix}x_{\Venus}-x_{\Earth}\\y_{\Venus}-y_{\Earth}\\z_{\Venus}-z_{\Earth} \end{pmatrix} \end{align*}

Der Abstand l zwischen der Erde und des Planeten ergibt sich wieder aus dem Satz des Pythagoras:

(25)   \begin{align*} l=\sqrt{X^2+Y^2+Z^2} \end{align*}

Für unser Beispiel errechnen sich folgende Werte:

(26)   \begin{align*} \begin{pmatrix}X\\ Y\\ Z \end{pmatrix}&=\begin{pmatrix}0,2412\\ 1,7176\\ 0,0056 \end{pmatrix}\,\mathrm{AE}\\ l_{\Venus,\Earth}&=1,7345\,\mathrm{AE} \end{align*}

→ Die Venus ist also am 11.06.2016 um 13:30 Uhr genau 1,7345 AE von der Erde entfernt.

Schritt 8: Nun wollen wir diese Fragestellung beantworten:

Welche Rektaszension \alpha und Deklination \delta besitzt die Venus am 11.06.2016 13:30 Uhr?

Die geozentrische orthogonale Ekliptikalkoordinaten werden nun in geozentrische spährische Ekliptikalkoordinaten transformiert:

(27)   \begin{align*} \beta&=\mathrm{atan2} (Z, \sqrt{X^2+Y^2})\\ \lambda&=\mathrm{atan2}(Y, X) \end{align*}

Es wird mit Absicht auf die Funktion \mathrm{atan2}(y, x) (Arkustangens mit zwei Argumenten) zurückgegriffen. Somit braucht sich nicht die Frage zu stellen, in welchen Quadranten man sich aktuell befindet.

Achtung! Die Winkel sind im Bogenmaß. Für unser Beispiel errechnen sich folgende Zahlenwerte in Grad:

(28)   \begin{align*} \beta&=0,1839^\circ\\ \lambda&= 82,0050^\circ\\ \end{align*}

Schritt 9: Das Ziel ist erreicht: Die Rektaszension und die Deklination können wir nun berechnen:

Mit dem Erdneigungswinkel \epsilon=23,439291-0,013004\cdot T ergeben sich die äquatoriale Koordinaten:

(29)   \begin{align*} \delta&=\mathrm{asin}(\sin\beta\cdot\cos\epsilon+\cos \beta\cdot \sin\epsilon\cdot \sin \lambda)\\ \alpha&=\mathrm{atan2}(\sin \lambda\cdot \cos \epsilon-\tan\beta\cdot \sin\epsilon, \cos \lambda)\\ \end{align*}

Auch hier ist die Rektaszension im Bogenmaß angegben. Für unser Beispiel errechnet sich nun folgende Winkel in Grad:

(30)   \begin{align*} \delta&=23,3795^\circ\\ \alpha&=81,2846^\circ\\ \end{align*}

Ergebnisbewertung

Die Beantwortung der Fragestellung lässt sich folgendermaßen beantworten:

  • Die Venus wird am 11.06.2016 um 13:30 Uhr einen Abstand von l=1,7345\,\mathrm{AE} besitzen. Damit liegt sie in der Umlaufbahn zwischen Sonne und Erde.
  • Die Deklination beträgt \delta=23,3795^\circ und die Rektaszension beträgt \alpha=81,2846^\circ.
Planetenstellung am 11.06.2016 13:30, Blick oberhalb von der z-Achse

Planetenstellung am 11.06.2016 13:30, Blick oberhalb von der z-Achse

Vergleicht man die Werte mit Stellarium, so ergibt sich eine sehr gute Genauigkeit. Für die Deklination zeigt Stellarium \delta=23,0394^\circ  und für die Rektaszension \alpha=81,27083^\circ an. Für den Abstand zeigt Stellarium l=1,735\,\mathrm{AE} Die Abweichungen resultieren deshalb, weil wir die Aufgabenstellung nach wie vor als Zweikörperproblem (ZKP) betrachtet haben. Stellarium wählt hier einen anderen Weg. Das Programm wendet hier die VSOP87-Planetentheorie an. In dieser Theorie sind die Bahnelemente zeitlich veränderlich und beschreiben die Keplerbahn, die mit der tatsächlichen Bahn oskuliert und werden über Potenzreihen berechnet.  Als Argument dient hier die Zeit und wird bis zur 5. Potenz und jeweiligen vorkommende Faktoren werden mittels einer Fourierranalyse ermittelt. Somit kann eine höhere Genauigkeit im Ergebnis erzielt werden.

Vergleichswerte in Stellarium

Vergleichswerte in Stellarium

Fazit und Danksagung

Dieser Beitrag war der umfangreichste Teil dieser Serie, denn die Ephemeridenrechnung ist nicht einfach.  Das Durchrechnen und die mathematische Modellierung in Matlab hat sehr viel Zeit und Nerven gekostet. Natürlich berechnet heute jedes Astronomieprogramm die Planetenstellung. Dennoch habt Ihr jetzt einen Einblick hinter den Kulissen und wisst nun, was für einen Aufwand hinter diesen Berechnungen steckt.

An dieser Stelle möchte ich auch zwei Danksagungen loswerden, denn mit Ihrer Hilfe habe ich die mathematische Modellierung und die Implementierung in Matlab erst geschafft. Zum einem sind es die Nutzer des Astronomieforums astrotreff.de und zum anderen, Herr Prof. i. R. Dr.-Ing. Helmut Haase, Professor an der Leibnitz-Universtät Hannover. Vielen Dank für die tolle Unterstützung 🙂

Zum Nachvollziehen der Ergebnisse kann man die Matlab-Files herunterladen und selber ausprobieren. Im letzten Teil dieser Serie möchte ich mit den Bahnelementen, die Zeitgleichung besprechen. Dies wird auch der letzte Beitrag zu dem Thema Himmelsmechanik sein. Falls etwas unklar ist oder Ihr Fehler entdeckt habt, lasst es mich wissen.

Literatur

Jean Meeus: Astronomical Algorithmus, Second Edition, Atlantic Books 1998
Karttunen, Kröger u.a. : Astronomie. Eine Einführung, Springer Verlag 1990
Oliver Montenbruck: Einführung in die Ephemeridenrechnung, Springer Verlag 2010

Links

Himmelsmechanik Teil 1: Koordinatensysteme
Himmelsmechanik Teil 2: Störung von Koordinatensystemen
Himmelsmechanik Teil 3: Die Keplerschen Gesetze
Helmut Haase Scilab Lernseite
Thread im Astrotreff.de-Forum zur Fragestellung

Downloads

Matlab-Dateien zum Nachvollziehen der Ergebnisse
E.M. Standish: Keplerian Elements for Approximate Positions of the Major Planets, (pdf-Datei)

, , , , ,

Schreibe einen Kommentar

Deine E-Mail-Adresse wird nicht veröffentlicht. Erforderliche Felder sind mit * markiert.