Sun Almanac selbst berechnet

Wo ist die Sonne?

Ohne ein Nautisches Jahrbuch bzw. einen Sun Almanac ist eine Navigation mit der Sonne unmöglich. Man muss schließlich wissen auf welcher Position über der Erdoberfläche die Sonne in einer vorgegebenen Sekunde eines Jahres gerade im Zenit steht. Wird die Entfernung zu diesem Punkt in der gegebenen Sekunde mit einem Sextanten gemessen, dann kann man sich auf der Erdkugel eine Kreislinie mit gleicher Entfernung zum Bildpunkt vorstellen, auf der man sich gerade befinden könnte. Wenn die Sonne dann ein paar Stunden später einen anderen Bildpunkt hat, kann auch der Radius dieses zweiten Kreises gemessen werden. Beide Kreise überlappen sich und kreuzen sich dadurch an zwei Stellen und einer dieser Kreuzungspunkte ist dann der eigene Standort.

Nach diesem Prinzip arbeiten fast alle Methoden zur Standortbestimmung auf hoher See. Weil sich aber die jeweiligen Positionen der Bildpunkte der Sonne an kein festes Datum und schon gar nicht an eine feste Zeit halten, wird für jedes Jahr ein neu aufgelegtes nautisches Jahrbuch benötigt. Das muss dann auch immer wieder neu beschafft werden, sonst kann man eine astronomische Standortbestimmung vergessen. Besser wäre deshalb eine Möglichkeit, die Positionen des Sonnenbildpunktes mit ausreichender Genauigkeit selbst berechnen zu können. Die heutigen Mobiltelefone oder auch ein Bordcomputer sind dazu durchaus in der Lage und dann wäre man unabhängig. Man sollte vielleicht meinen, dass diese Berechnungen äußerst kompliziert sind. Das ist aber nicht so, wie die es nachfolgenden Ausführungen zeigen werden.

Bahnparameter

Zur Berechnung der Position der Sonne, wo sie also in einer vorgegebenen Sekunde an einem vorgegebenen Datum im Zenit steht, werden zunächst die Bahnparameter der Erde benötigt. Diese sind:

  1.  L = mittlere Länge in ab Perihel in Grad (_A1, _B1)
  2.  𝜛 = Länge des Perihels in Grad (_A2, _B2)
  3.   e = numerische Exzentrizität (_A3, _B3)
  4.   𝜖 = Schiefe der Ekliptik in Grad (_A4, _B4)

Die folgende Tabelle zeigt mit den Daten _A1 bis _A4 den Stand am letzten Referenzdatum ref = 1.1.2000 um 12:00:00 UT1. Die Daten _B1 bis _B4 zeigen die Veränderungen, wie sie sich innerhalb eines Tages ergeben.

Nr. Name Parameter
1 _A1 = 280,4656°
2 _A2 = 282,94°
3 _A3 = 0,016709
4 _A4 = 23,4391667°
5 _B1 = 0,98564733744011°
6 _B2 = 0,0000470691307323751°
7 _B3 = 0,00000000114989733059548
8 _B4 = 0,000000356072705148681°
9 ref = 1.1.2000

Die Berechnung der konkreten Werte für L, 𝜛, e und 𝜖 für einen Zeitpunkt nach dem Referenzdatum 1.1.2000 erfolgt mit Hilfe einer linearen Gleichung in der Form

    \[\text{Wert}=A_x +T\cdot B_x\]

Hierin sind x die Indizes 1 bis 4 aus der vorstehenden Tabelle und „Wert“ entspricht den Parametern in der obigen Aufzählung von 1. bis 4.

Um die Werte ausrechnen zu können, muss die Zeit T bestimmt werden. Diese Zeitangabe erfolgt in Tagen mit einer genügenden Anzahl von Nachkommastellen. Die Zeit T wird nach folgender Formel berechnet:

    \[T=\mathbb{Z}\Bigr(\big(JJJJ-2000\big)\cdot 365,25\Bigr)+N-0,5+\frac{UT1}{24}-SJ\cdot\]

Mit ℤ ist der ganzzahlige Teil des Wertes gemeint, der in der nachfolgenden Klammer steht. Die darin zu subtrahierende Jahreszahl 2000 ist das Jahr des Referenzdatums. N bezeichnet die im betrachteten Jahr vergangenen Tage einschließlich des aktuellen Tages. SJ ist die für das 21. Jahrhundert geltende Schaltjahresregel:

    \[SJ=\text{wenn}\;4|(JJJJ -2000)\;\text{dann}\;1\;\text{sonst}\;0\cdot\]

Für diese aufwändige Formel gibt es in den Rechnerprogrammen oft einfachere Ausdrücke. So lässt sich die Zeit T in Excel einfach mit dem Ausdruck „=TAGE(dat;ref) +ot-0,5“ berechnen, wobei ot = observation time bedeutet und der Abzug von 0,5 berücksichtigt, dass sich das Referenzdatum auf 12:00:00 UT1 und nicht auf 00:00:00 UT1 bezieht. Nehmen wir beispielsweise den Zeitpunkt

  • Datum            30.05.2023            Name dat,
  • Zeit                 07:36:07 UT1        Name ot,

dann berechnet sich dafür die Zeit T mit genau T = 8549,81674768518 Tagen nach dem Referenzdatum 1.1.2000 | 12:00:00 UT1.

Rechentabelle

Damit können nun L, 𝜛, e und 𝜖 ausgerechnet werden. Rechenprogramme arbeiten mit dem Bogenmaß. Dieses ist das Produkt von Winkelmaß in Grad und dem Faktor 𝜋/180. Die Ergebnisse stehen in der nachstehenden Tabelle in den Zeilen 2, 3, 5 und 8 in der Spalte „Observation“. Die großen Zahlen in den Zeilen 2, 4 und später auch 7 mögen zunächst irritieren. Man könnte hier eine ganzzahlige Menge von 2𝜋 (360°) bis auf einen verbliebenen Rest, den sogenannten Hauptwinkel, entfernen. Das muss man aber nicht, denn die Winkelfunktionen haben nach einem um jeweils 2𝜋 vergrößerten Winkel immer wieder den gleichen Wert.

In der folgenden Tabelle ist die gesamte Berechnung für den Greenwichwinkel Grt und die Deklination 𝛿 für den genannten Beobachtungszeitpunkt 30. Mai 23, 7:36:07 UT1  angegeben.

Nr. Name Symbol Observation
1 T T = 8549,817
2 L Lrad = 151,9758
3 Pi 𝜛rad = 4,9453
4 M M = 147,0305
5 e e = 0,0167
6 _C C = 0,0192
7 Lam Λ = 151,9950
8 eps εrad = 0,4090
9 g grad = 0,0109
10 Grt0 Grt0 = 5,1427
11 Grt Grt = 5,1427
12 d δ = 0,3795
13 Grt° Grt° = 294,65°
14 δ° = 21,75°

In beiden Tabellen enthält die Spalte „Name“ einen Vorschlag für die in den Rechenprogrammen, z. B. Excel, zu verwendeten Namen, weil das Formelsymbol selbst dafür meist nicht verwendet werden kann. So beziehen sich alle Variablen immer auf die Namen und nicht mehr auf die Felder. Wer also die Berechnung in Excel auf seinem eigenen Rechner nachvollziehen möchte, der muss den Feldern in den Spalten „Observation“ bzw. „Parameter“ in beiden Tabellen die in der Spalte „Name“ angegebenen Namen zuweisen. Auch dem Beobachtungsdatum (dat) und der Beobachtungszeit (ot) müssen ihre Namen zugewiesen werden. Es sind somit insgesamt 25 Namen anzulegen.

Die Zeitgleichung

Mit der Differenz von L und 𝜛, der mittleren Anomalie M, kann jetzt die Mittelpunktsgleichung C berechnet werden. Damit wird die Elliptizität der Umlaufbahn berücksichtigt. Die Mittelpunksgleichung lautet:

    \[C=(2e-\frac{e^3}{4})\cdot\sin M+\frac{5}{4}e^2\cdot\sin (2M)+\frac{13}{12}e^3\cdot\sin (3M)+ ...\cdot\]

Es ist eine Reihenentwicklung, die nach dem dritten Glied abgebrochen werden kann, weil weitere Glieder keine bessere Genauigkeit liefern würden. Die Summe von C und dem Hauptwinkel L = 𝛼m liefert nun den Weg

    \[\Lambda=\alpha_m+C\]

der wahren Sonne S auf der Ekliptik, während 𝛼m der zeitgleich zurückgelegte Weg einer Vergleichssonne Z ist, die mit derselben Geschwindigkeit wie die der wahren Sonne auf der Himmelsäquator-Ebene ihre Bahn zieht. Bild 1 zeigt das in einem grafischen Modell. Das Bild entspricht einem geozentrischen Verständnis, wie es der Navigator an Bord eines Schiffes auch wahrnimmt.

Bild 1: Die Erde im Zentrum der Himmelkugel.

Daraus kann jetzt mit 𝛼m die Zeitgleichung g = 𝛼m – 𝛼 berechnet werden:

    \[g=\alpha_m-\arctan_\Lambda(\tan\Lambda\cdot\cos\epsilon)\cdot\]

Nach Anwendung eines Tangens Additionstheorems folgt daraus eine zwar größere, aber besser zu händelnde Gleichung (tan 𝛼m = tan L):

    \[g=\arctan\frac{\tan L-\tan \Lambda\cdot\cos \epsilon}{1+\tan L\cdot\tan \Lambda\cdot\cos \epsilon}\cdot\]

Hat man die Zeitgleichung einer Sekunde des Tagesdatums, dann hat man auch den Greenwichwinkel Grt, denn es gilt:

    \[Grt=g+15\cdot\big(t+12^h\big);\;\;\;\;\;0\le\big(t+12^h\big)\le 24\cdot\]

Hierin ist t die Zeit UT1 (bzw. UTC) des Tages im 24-Stunden-Lauf, für die der Grt berechnet wird. Die Zeit 12 h wird von der Sonne gebraucht, um von 180° E kommend den Nullmeridian zu erreichen. Zu diesen 12 h muss die Zeit, für die der Grt bestimmt werden soll, addiert werden. Die Summe der Zeiten liefert dann nach Multiplikation mit der Winkelgeschwindigkeit der Sonne von 15°/h den mittleren Ortsstundenwinkel von Greenwich. Nach Addition der Zeitgleichung g erhält man daraus den aktuellen Ortsstundenwinkel Grt. Der Klammerinhalt darf nicht überlaufen, also größer als 24h werden. Für Nachmittags-Zeiten müssen deshalb immer 24 Stunden subtrahiert werden. Das Ergebnis kann auch negativ sein und dann müssen 2𝜋 addiert werden.

Die Deklination 𝛿 ist eine Kathete des rechtwinkligen sphärischen Dreiecks 𝛼 𝛿 𝛬 und somit einfach mit der Formel

    \[\delta=\arcsin(\sin\epsilon\cdot\sin\Lambda)\]

zu berechnen. In den Zeilen 13 und 14 sind die Ergebnisse dann in Gradangaben umgerechnet.

Berechnung in Excel

Es ist also kein Hexenwerk, die Ephemeriden der Sonne selbst zu berechnen. Das hier benutzte Beispiel wurde mit Hilfe von Excel ausgeführt, weil fast jeder diese Tabellenkalkulation auf seinem heimischen PC hat. Wer alles genau nachvollziehen will, der findet in der folgenden Aufführung die genauen Excel-Formeln. Sie können kopiert und in das entsprechende Feld in der Spalte Observation eingeschrieben werden.

1  =TAGE(dat;ref)+ot-0,5
2  =BOGENMASS(_A1+_B1*T)
3  =BOGENMASS(_A2+_B2*T)
4  =L-Pi
5  =_A3-_B3*T
6  =(2*e-e^3/4)*SIN(M)+5*e^2/4*SIN(2*M)+13*e^3/12*SIN(3*M)
7  =L+_C
8  =BOGENMASS(_A4-_B4*T)
9  =ARCTAN( (TAN(L)-TAN(Lam)*COS(eps)) /(1+TAN(L)*TAN(Lam)*COS(eps)))
10  =WENN((ot*24+12)>24;g+BOGENMASS(15)*(ot*24-12);g+BOGENMASS(15)*(ot*24+12))
11  =WENN(Grt0<0;Grt0+2*PI();Grt0)
12  =ARCSIN(SIN(Lam)*SIN(eps))
13  =GRAD(Grt)
14  =GRAD(d)

Eine fertige Excel Datei, die mit genau den vorstehenden Ausführungen arbeitet kann mit dem folgenden Link auch heruntergeladen werden. Ihre Gültigkeit ist bis 31.12.2099 gegeben.

Sun Almanac

Fehlerbetrachtung

Man könnte meinen, dass diese so simplen Berechnungen nicht präzise genug sind. Um einigermaßen Aufschluss darüber zu erhalten, wurden für einen im Jahr 2021 definierten Standort im Mittelmeergebiet die Abweichungen berechnet, die sich aus dem Unterschied zwischen den selbst berechneten Daten von Greenwichwinkel Grt und Deklination \delta und den Daten aus einem präzise berechneten Almanach ergeben. Konkret wurde ein vorgegebener Standort für jeden Tag des Jahres mit der Gauß Methode zweimal berechnet. Einmal unter Verwendung von selbst berechneten Grt und \delta und dann noch einmal mit Grt und \delta aus dem online verfügbaren Almanach https://www.nauticalalmanac.it/. Die sich daraus ergebenen Differenzen sind im Bild 2 dargestellt.

Bild2: Rechnerische Abweichungen eines Standortes mit selbst berechneten Grt und \delta gegenüber einer Berechnung, die den Einfluß von Mond und Planeten berücksichtigt.

Die Abweichungen kommen zustande, weil die Einflüsse des Mondes und der großen Planeten in den eigenen Berechnungen nicht erfasst sind, damit der Rechenumfang nicht ausufert. Sie sind jedoch niemals größer als eine halbe Meile. Was bedeutet das? Ein guter Navigator bringt es in der astronomischen Standortbestimmung auf einem nicht oder nur wenig schwankenden Deck auf eine Standortabweichung von 2 NM. Nach dem Fehlerfortpflanzungsgesetz addieren sich Einzelfehler jedoch geometrisch. Dafür ein Beispiel

  • F_M=\text{Summe der Messfehler}\quad\pm 1,94\,NM
  • \underline{F_C=\text{Rechenfehler}\quad\quad\quad\Quad\quad\,\;\;\;\;\;\pm0,50\,NM}
  • \text{maximal möglich}\quad\quad\quad\quad\quad\,\;\;\;\;\pm2,44\,NM

Man weiß aber nie, ob die Einzelfehler gerade positiv oder negativ ausfallen. Dass alle gerade positiv oder alle gerade negativ sind, ist unwahrscheinlich. In diesem Beispiel gilt deshalb eine statistisch zu erwartende Standortabweichung von

    \[F^2_{total}= F^2_M+F^2_C\,;\;\;F_{total}=\sqrt{ F^2_M+F^2_C}=2,00\, NM\cdot\]

Setzen wir den rechnerischen Fehler auf null, weil z. B. ein amtlicher Almanach benutzt wird, dann beträgt der Fehler

    \[F^2_{total}= F^2_M+0\,;\;\;F_{total}=\sqrt{ F^2_M+0}=1,94\,NM\cdot\]

Der Unterschied zwischen diesen Werten beträgt 0,063 NM und das sind 3,2%, die auf die im Bild 2 gezeigten Abweichungen zurückzuführen sind. Dieser geringe Anteil sinkt sogar drastisch, sobald andere Fehler zunehmen. Die im Bild dokumentierten Abweichungen infolge selbst berechneter Sonnenephemeriden haben somit keine praktische Bedeutung. Das Höhenverfahren von Saint Hilaire besitzt in der Regel größere Rechenfehler, nicht durch fehlerhafte Ephemeriden, sondern durch die Ungenauigkeit einer Standortschätzung und unpräzise ausgeführte grafische Arbeiten.

Interessant ist die Welligkeit der im Bild 2 gezeigten Kurven. Der höherfrequente Anteil der Welligkeit ist durch die Mondzyklen bedingt, denn es ist nicht der Erdmittelpunkt, der auf seiner elliptischen Bahn die Sonne umkreist. Vielmehr kreist der Schwerpunkt des Massensystems Erde-Mond um die Sonne, wobei dieser Schwerpunkt 4760 km vom Erdmittelpunkt entfernt ist. Dadurch „eiert“ der Erdmittelpunkt auf seiner Bahn um die Sonne hin und her, ist mal voran und dann wieder zurück. Eine Schwingungsperiode dauert dabei so lange wie ein Mondumlauf, also etwa einen Monat.
Der niederfrequente Anteil, zeigt die Annäherung der Erde auf ihrer jährlichen Umlaufbahn um die Sonne an die zwei massereichsten Planeten Jupiter und Saturn, deren Umlauf um die Sonne 12 bzw. 30 Jahre betragen kann. Diese großen Wellen, wie sie Bild 1 zeigt, werden sich in den nächsten Erdjahren also immer ein wenig verschieben.

Schreibe einen Kommentar

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