Astrosimulationen

Im Jahr 1990 kaufte mir mein Vater fürs Studium einen damals brandaktuellen PC, konkret einen 386-er mit ganzen 33 MHz Taktrate und 42 MB Festplatte! Damit konnte ich in den darauffolgenden Jahren das eine oder andere Programm (zumeist in TurboPascal) aus dem Bereich Physik bzw. Astronomie schreiben.

Im Alter von 17 Jahren stieß ich zufällig in der amerikanischen Zeitschrift Astronomy (welche ich mir immer in einem Laden vornehmlich für Sexzeitschriften besorgte 😉 ) auf ein unscheinbares BASIC-Programm, welches trotz seiner wenigen Zeilen in der Lage zu sein schien, komplexe kosmische Bewegungen zu simulieren. Konkret ging es um die gravitative Wechselwirkung eines massereichen schwarzen Lochs mit dem Zentralkörper und den diesen umrundenden Sternen einer Galaxie.

Ich bin bei meiner Recherche nach dieser Zeitschrift fündig geworden. Es handelt sich bei der gesuchten Ausgabe um jene vom Dezember 1988. Auf den Seiten 90-96 wird unter dem Titel „Galactic collisions on your computer“ von Michael C. Schroeder und Neil F. Comins diese Simulation präsentiert.

Natürlich habe ich den Programmcode gleich abgeschrieben und ich war wirklich verblüfft zu sehen, wie nun das schwarze Loch völlig realistisch meine Galaxie am Computermonitor sukzessive verformte. Den mathematischen/physikalischen Hintergrund verstand ich damals noch nicht wirklich. Heute weiß ich, dass hinter der Simulation von damals das sog. Euler-Verfahren zur Lösung der Differentialgleichungen steckt. Das eulersche Polygonzugverfahren ist das einfachste Verfahren zur numerischen Lösung eines Anfangswertproblems. Es wurde von Leonhard Euler 1768 in seinem Buch Institutiones Calculi Integralis vorgestellt.

Was steckt dahinter? Nehmen wir unsere Astrosimulation. Die wirkenden Kräfte sind rein gravitativer Natur. Das Newtonsche Gravitationsgesetz beinhaltet zur Berechnung der Kraft zwischen 2 Himmelskörper neben der Gravitationskonstante G nur deren Massen M1 und M2 und deren Abstand d:

Kenne ich alle Positionen der Sterne, des Zentralkörpers und des schwarzen Lochs, so kann ich einmal die auf jedes dieser Objekte wirkende Gesamtkraft F berechnen. Was fange ich aber nun mit diesen Kraftvektoren an?

Das zweite Newtonsche Axiom beinhaltet die sog. Bewegungsgleichung

F = m · a.

Diese wohl wichtigste Formel der ganzen Mechanik verknüpft also die auf einen Körper wirkende Kraft F mit dessen Beschleunigung a. Kenne ich F, so kenne ich auch a. Wir kennen aber auch den Zusammenhang zwischen der Beschleunigung a und der Geschwindigkeitsänderung dv bzw. Δv innerhalb des Zeitintervalls dt bzw. Δt. Er lautet:

a = dv/dt (infinitesimal) bzw. bei größeren Zeitintervallen a = Δv/Δt.

Formt man diese Gleichung nach Δv = v(t + Δt) v(t) um erhält man:

v(t + Δt) – v(t) = a · Δt      bzw.       v(t + Δt) = v(t) + a · Δt

Damit kennen wir also bereits die Geschwindigkeiten v(t + Δt), also zum etwas späteren Zeitpunkt t + Δt! Nun möchten wir nur noch die neuen Orte zum Zeitpunkt t + Δt kennen. Hierfür bedienen wir uns des Zusammenhangs zwischen v und der Ortsänderung ds bzw. Δs innerhalb des Zeitintervalls dt bzw. Δt. Er lautet:

v = ds/dt (infinitesimal) bzw. bei größeren Zeitintervallen v = Δs/Δt.

Mit dieser Beziehung sind wir [wie vorhin bei der Berechnung der Geschwindigkeit v(t + Δt)] in der Lage den aktualisierten Ort s(t + Δt) zum Zeitpunkt t + Δt zu berechnen. Es gilt:

Δs = s(t + Δt) – s(t) = v · Δt      bzw.       s(t + Δt) = s(t) + v · Δt

Welches v setzen wir aber nun hier ein? Wir könnten die Geschwindigkeit v(t) zu Beginn des betrachteten Zeitintervalls [t, t+Δt] einsetzen, was von der Herleitung her naheliegend erscheinen würde. Bezüglich der Geschwindigkeiten kennen wir aber jetzt nicht nur die Geschwindigkeit v(t) zu Beginn des betrachteten Zeitintervalls [t, t+Δt], sondern auch jene an dessen Ende in Form von v(t + Δt). Diesen Vorteil wollen wir natürlich nützen und setzen daher das arithmetische von v(t) und v(t + Δt) ein und erhalten

s(t + Δt) = s(t) + ½ · [v(t) + v(t + Δt)] · Δt

Wir kennen damit also nicht nur die „neuen“ Geschwindigkeiten zum Zeitpunkt t + Δt, sondern auch die aktualisierten Orte s(t + Δt) zum Zeitpunkt t + Δt. Und genau mit diesen neuen Werten werden abermals in Form einer Berechnungsschleife die aktuellen Kräfte F, Beschleunigungen a, Geschwindigkeiten v und Orte s berechnet. Dieses Spiel wiederholt sich immer wieder, was natürlich ein gefundenes Fressen für einen Blechtrottel ist. Hier kann er seine ganze Stärke ausspielen und vollführt tausende Berechnungen innerhalb einer einzelnen Sekunde.

Dieses Euler-Verfahren ist wie oben schon erwähnt das einfachste Verfahren zur Lösung dieses Anfangwertproblems. Es lässt sich aber im Gegensatz zu komplexeren Methoden (wie etwa das Runge-Kutta-Verfahren) auch für Schüler einigermaßen leicht nachvollziehen. Die Genauigkeit des Euler-Verfahrens hängt natürlich entscheidend von der gewählten Schrittweite Δt ab. Je feiner man diese wählt, umso genauer wird die Simulation bzw. umso später erst weicht dieses simple Verfahren deutlich von den realen Verhältnissen ab. Zum Abschluss also nochmals die grafische Zusammenfassung des hier zum Einsatz kommenden Euler-Verfahrens:

Mit diesem Verfahren habe ich etliche meiner astrophysikalischen Simulationen vom Sonnenanalemma bis hin zur Galaxienkollision berechnet. Einige davon werden nachfolgend kurz vorgestellt.


Sonnenanalemma

Bildquelle: https://de.wikipedia.org/wiki/Analemma

Fotografiert man die Sonne im Laufe einen ganzen Jahres immer zur selben Uhrzeit (z.B. 12 Uhr mittags), so ändert sich nicht nur die Höhe über dem Horizont (im Winter tiefer, im Sommer höher), sondern auch der sog. Azimut, also der Positionswinkel. Denn die Sonne befindet sich keinesfalls so wie man annehmen könnte um 12 Uhr mittags immer genau im Süden! Dies ist nur an 4 Tagen im Jahr genau der Fall. An allen anderen Tagen befindet sich die Sonne etwas in Richtung Osten (geht quasi nach) oder Westen (geht quasi vor).

Basis der Berechnung des Analemmas war die Simulation der Erdbewegung um die Sonne im Laufe eines Jahres. Die Erde bewegt sich ja nach Johannes Kepler auf einer Ellipse um die Sonne, wobei sich diese in einem der beiden Brennpunkte befindet. In der unteren Abbildung ist die Exzentrizität der Erdbahn deutlich vergrößert dargestellt. In Wirklichkeit unterscheiden sich die Abstände Erde-Sonne in den Punkten 1 (Perihel) und 2 (Aphel) bedeutend weniger als skizziert.


Satellitensimulation

Die Simulation von Satellitenbahnen um die Erde ist deutlich einfacher als die Herleitung des Analemmas. Ich griff bei der Positionsberechnung des Satelliten abermals auf das Euler-Verfahren zurück. Mit dieser Simulation lässt sich etwa die erste kosmische Geschwindigkeit vI simulieren bzw. wie die Bahn aussieht, wenn der Satellit eine Geschwindigkeit größer bzw. kleiner als die für seinen Abstand notwendige Kreisbahngeschwindigkeit besitzt. Die Schüler erkennen dann sofort, dass sich dabei Ellipsenbahnen ergeben. Einen besonderen Spaß bereitet es natürlich den Schülern, den Satelliten abstürzen zu lassen…


Galaxienkollisionen

Diese Simulation wurde inspiriert durch den Eingangs bereits erwähnten Artikel in der Zeitschrift Astronomy Ende der 80-er Jahre. In der ersten Variante stört ein massereiches schwarzes Loch (dargestellt durch einen roten Kreis) eine Galaxie (deren massereicher Zentralkörper ist der grüne Kreis) und ihre Sterne (weiße Punkte).

In der zweiten Variante ersetze ich das schwarze Loch durch eine zweite vollständige Galaxie (roter Kreis und rote Sterne). Ziel war es, solche engen kosmischen Begegnungen wie sie etwa bei den beiden Mäusegalaxien NGC 4676 A und NGC 4676 B im Sternbild Haar der Berenike oder bei NGC4038 und NGC4039 im Sternbild Rabe stattgefunden haben zumindest ansatzweise zu simulieren.


Variante 1: Galaxie vs. Schwarzes Loch

Bei dieser Simulation tritt ein schwarzes Loch gegen eine Galaxie an, die über einen massereichen Zentralkörper verfügt. Dessen Masse und die Masse des schwarzen Lochs können mittels Potentiometer eingegeben werden. Zudem noch die Startposition (x/y/z) und Startgeschwindigkeit (vx/vy/vz) des schwarzen Lochs. Die Anzahl der Sterne in der Galaxie ist auch variabel, wobei aufgrund des limitierten Speichers des Arduino Due maximal 80 Sterne gewählt werden können.

Nach der Werteeingabe gelangt man zur Simulation. In der linken Bildschirmhälfte wird die x-y Ansicht, rechts die x-z Ansicht ausgegeben. Mittels Taster sind die Zeitschrittweite dt und der Zoomfaktor einstellbar.

Arduino-Code:

Bilddatei der Galaxie M51:  M51

 


Variante 2: Galaxie vs. Galaxie

Bei der zweiten Variante wird das schwarze Loch gegen eine weitere Galaxie ausgetauscht. Die Masse der beiden Zentralkörper ist ebenso einstellbar wie die Startposition (x/y) und Startgeschwindigkeit (vx/vy) der zweiten Galaxie. Die Anzahl der Sterne pro Galaxie ist hier aufgrund des Arbeitsspeichers auf 40 Sterne beschränkt.

Da sich diese Simulation nur in der x-y-Ebene abspielt, ist der Bildschirm nicht zweigeteilt. Zu sehen ist die Draufsicht der beiden sich nähernden Galaxien.

Während der Simulation sind wieder wie bei der Variante 1 die Zeitschritte dt und der Zoomfaktor veränderbar.

Arduino-Code:

Bilddatei der Galaxie M51:  M51

Ein abschließendes Wort noch zur doch recht bescheidenen Anzahl an Sternen. Leider lässt der Speicher des Arduino Due keine allzu große Anzahl an Sternen zu. Im Fall Galaxy vs. black hole sind es max. 80 Sterne, im Fall Galaxy vs. Galaxy max. je 48 Sterne. Der Grund für diese Limitierung besteht darin, dass es bei n gravitativ wechselwirkenden Objekten insgesamt n·(n–1) Wechselwirkungen gibt. Man muss aber zum Glück nur die Hälfte dieser Wechselwirkungen ausrechnen, da erstens die Kraft auf sich selbst ja 0 ist für alle Objekte und zudem die Gravitationskraft ja symmetrisch ist. Die Kraft des i-ten Objekts auf das j-te Objekt ist die negative Kraft des j-ten Objekts auf das i-te. Durch diese Symmetrie verbleiben bei n Objekten also „nur“ noch ½·n·(n–1) zu berechnende Wechselwirkungen. Grob gesprochen wächst aber der notwendige Speicher und auch die Rechenzeit mit n² stark an.

Hier noch das Youtube-Video:


Venustransit und Merkurtransit

Bei dieser Simulation geht es um die Stellung der inneren Planeten Merkur und Venus bei einer sog. unteren Konjunktion. Dabei befindet sich der jeweilige Planet zwischen Sonne und Erde. Jetzt könnte man meinen, dass Merkur bzw. Venus bei einer jeden unteren Konjunktion von der Erde aus betrachtet vor der Sonnenscheibe vorbeiwandert. Dies ist aber nicht der Fall, da die Umlaufsebenen der Planeten „leicht“ unterschiedich sind. Daher ist es viel wahrscheinlicher, dass Merkur bzw. Venus bei einer unteren Konjunktion oberhalb oder unterhalb der Sonne vorbeizieht. Die Umlaufsebene des Planeten Merkur ist zum Beispiel um 7° gegenüber der Umlaufsebene der Erde (= Ekliptik) geneigt!

Ich habe weiters untersucht, ob die unteren Konjunktionen bei allen Positionswinkeln der Erde in etwa gleich wahrscheinlich sind. Bei 15000 unteren Konjunktionen erhielt ich folgende (Gleich)Verteilung:

Mit diesem Ergebnis lässt sich dann die durchschnittliche Häufigkeit eines Merkur- oder Venustransits berechnen. Die Sonne besitzt von der Erde aus betrachtet einen Winkeldurchmesser von ca. 0.5° = 30 Bogenminuten. Merkur wandert bei einer unteren Konjunktion maximal um ca. 4.4° oberhalb bzw. unterhalb der Sonne vorüber (siehe untere Abbildung). Wenn nun alle Positionswinkel bei der eine untere Konjunktion stattfindet gleich wahrscheinlich sind, so gilt für die Wahrscheinlichkeit eines Merkurtransits:

P_Merkurtransit = günstige Positionswinkel/ mögliche Positionswinkel = 4 · 3.257° / 360° = 0.0362 = 3.62%.

Eine untere Konjunktion des Planeten Merkur findet ca. alle 116 Tage statt. Wenn nur 3.59% aller unteren Konjunktionen zu einem Merkurtransit führen, muss man im Schnitt 116/0.0362 = 3204 Tage auf einen Merkurtransit warten. Dies sind 8.78 Jahre. Laut Internet kommt es insgesamt ca. 13-mal pro Jahrhundert zu einem Merkurdurchgang, was also einem durchschnittlichen Abstand zwischen 2 Merkurtransiten von 100/13 = 7.69 Jahren entspricht. Beide Werte liegen nicht allzuweit voneinander entfernt…


Venushelligkeit

Die Konjunktion der Planeten Venus und Mars (ca. in der Bildmitte) im Jahr 2015

Die scheinbare Helligkeit der Venus variiert recht stark und besitzt maximal einen Wert von -4.4 mag (lat. magnitudo für Größenklasse). Wovon hängt deren scheinbare Helligkeit, also von der Erde aus betrachtet, ab? Zum einen einmal sicher vom Abstand Erde-Venus. Die Lichtintensität einer punktförmigen Lichtquelle nimmt ja nach dem Abstandsgesetz 1/r² ab. Dies kann man auf die Venus übertragen. Demnach müsste deren Helligkeit in der unteren Konjunktion (Abstand Erde-Venus ist minimal) am größten sein. Dies ist aber nicht der Fall. Grund sind die unterschiedlichen Phasen der Venus. In der unteren Konjunktion beträgt der von der Erde aus sichtbare beleuchtete Teil der Venusoberfläche nahezu 0%. Daher wird/ist die Venus in der unteren Konjunktion zumeist unsichtbar. Die größtmögliche Phase von 100% besitzt die Venus in der oberen Konjunktion, also von der Erde aus betrachtet hinter der Sonne. In dieser Position besitzt aber die Venus zugleich auch den größten Abstand zur Erde. Daher wird das Helligkeitsmaximum irgendwo zwischen oberer und unterer Konjunktion liegen.

Mit meiner Simulation habe ich für alle Stellungen Erde-Sonne-Venus die beleuchtete Fläche (Phase) und den Abstand Erde-Venus berechnet und daraus deren Helligkeit abgeleitet.

Wie man anhand des Graphen erkennen kann, beträgt der Positionswinkel der Venus bei ihrer größten scheinbaren Helligkeit rund 160°. Die Erde befindet sich fix bei einem Positionswinkel von 180°. Dies deckt sich recht gut mit dem praktischen/echten Helligkeitsverlauf…