kleiner planetarer Simulator

Bernhard

Registriertes Mitglied
Mit den Formeln aus dem Wikipedia-Artikel zum Zweikörperproblem kann man daraus bereits die genaue Bahn berechnen:
E = -3.88e28 J, L = 2.82e34 Js, mu = 7.26e22 kg, epsilon = 0.094, p = 3.74e8 m, a = 3.77e8 m und b = 3.76e8 m
Hallo Kibo,

hier haben sich leider auch kleinere Fehler eingeschlichen. Deswegen noch die folgende Erklärung. Für die Berechnung der Integrationskonstanten und Bahnparameter berechnet man zuerst den Differenzvektor der Ortsvektoren und dessen Zeitableitung aus den Startbedingungen:

r = (KoordX1 - KoordX2, KoordY1 - KoordY2, KoordZ1 - KoordZ2)
v = (v_X1 - v_X2, v_Y1 - v_Y2, v_Z1 - v_Z2)

rb := Betrag von r
vb := Betrag von v

, sowie die reduzierte Masse mu = m1 * m2 / (m1 + m2)

Bei Verwendung von Schwerpunktskoordinaten berechnet sich die Gesamtenergie dann gemäß

E = 0.5 * mu * v^2 - G * m1 * m2 / rb

Der vektorielle Drehimpuls berechnet sich zu L = mu * r x v. Da nur der Betrag weiter verwendet wird, nutzt man die Eigenschaft, dass r und v aufeinander senkrecht stehen und berechnet den Betrag gemäß l = mu * rb * vb.

Es gilt dann

p := l^2 / (mu * G * m1 * m2)
epsilon = sqrt(1.0 + 2.0 * E * p / (G * m1 * m2) )
a = p / (1.0 - epsilon * epsilon)
b = p / sqrt(1.0 - epsilon * epsilon)

Für Tests der Programme mit zwei Massekörpern empfiehlt es sich also r und v so zu wählen, dass sie aufeinander senkrecht stehen. Man spart sich dann das Kreuzprodukt beim Impuls und kann relativ leicht die Bahnform berechnen. Die Formeln für epsilon, a und b sind nur für E < 0 gültig. Für E = 0 erhält man eine Parabel und für E > 0 eine Hyperbel.
MfG
 

Kibo

Registriertes Mitglied
Guten Abend,

Ich darf Fortschritte vermelden, denn Version 1.4 :

verwendet SI Einheiten (Meter, Kg, m/s²),
zeigt beim System Erde-Mond halbwegs realitätsnahe Ergebnisse (Besten Dank für die Daten),
gibt jetzt auch Simulierte Zeit, Schrittweite und Laufzeit aus,
enthält sowohl Runge-Kutta als auch Euler,
basiert auf überprüfbaren Daten anstatt rumprobieren. ;)


Bekannte Bugs:
Ausgegebene "Simulierte Zeit" stimmt nicht mehr bei im laufenden Betrieb veränderter Schrittweite
Runge-Kutta liefert bei großer Schrittweite immernoch eine hohe Bahninstabilität im Gegensatz zu Euler
Weiterhin fehlerhafte Texturierung der Körper

Download

Skydrive

Rapidshare

mfg

PS: Da die Planeten.txt nach Bernhards Daten in der Datei integriert ist, sollte sichergestellt werden, dass sich keine alte Planeten.txt im selben Ordner befindet.
 
Zuletzt bearbeitet:

Bernhard

Registriertes Mitglied
Hi Kibo,

ich denke, Du hast das Euler- und das Runge-Kutta-Verfahren jetzt korrekt implementiert :) . Gratuliere. Als "Belohnung" gibt es von mir eine weitere Planeten.txt.

Code:
KoordX	KoordY	Koordz	ImpulsX	ImpulsY	ImpulsZ	Masse		dichte	s	v	Name	s und v sind noch nicht implementiert
-7.48e10	0.0001	0.0001	0.0001	2.105e4	0.0001	1.99e30	0	0	Sonne1	1.408
7.48e10	0.0001	0.0001	0.0001	-2.105e4	0.0001	1.99e30	0	0	Sonne2	1.408

Das sind die Startbedingungen für ein Binärsystem aus zwei genau gleichen Sonnen im Abstand 1 AE. Jede einzelne Sonne hat die gleiche Masse und Dichte wie unsere Sonne (Sol). Beide Sonnen laufen auf einer Kreisbahn umeinander.

Leider ist das Programm ein wenig langsam. An der Ausführungsgeschwindigkeit sollte man also noch arbeiten. Ich habe auch den Eindruck, dass die Nachleuchtspur nicht ganz passt und nicht exakt die Bahn nachzeichnet.

Es macht aber trotzdem schon Spaß das Programm zu testen.
MfG
 

Kibo

Registriertes Mitglied
Hallo Bernhard,

Erst einmal herzlichen Dank für's Testen.:) Mit der Geschwindigkeit gebe ich dir recht. Deine 2 Sonnen zeigen auch eines der noch nötigen Verbesserungen auf, die Möglichkeit, die Schrittweite komplett frei zu bestimmen (wenn nicht gar, je nach Situation automatisiert variabel zu gestalten), bei Entfernungen von 1 AE schläft man ja selbst bei niedrigster Genauigkeit ein. An meinem Rechner berechnet das Programm 3800 Schritte in 60 Sekunden, und das sowohl im RK als auch im Euler!
Das finde ich persönlich ziemlich erstaunlich, wenn man bedenkt wie viel größer der RK-Quelltext ist im Vergleich zum Euler. Eine Schrittweite von 1 entspricht übrigens genau einer Sekunde simulierter Zeit. Bedeutet bei niedrigster Genauigkeit, bei der Euler immerhin noch stabile wenn auch leicht exzentrische Bahnen bei Erde-Mond liefert (Abstand schwankt zwischen 382000 und 385000 km, die Bahn war ja entgegen der Realität kreisrund gedacht oder?), werden 2500 Sekunden pro Schritt, und damit ziemlich genau 110 Tage pro Minute berechnet. Mein PC hat einen Core2Quad 6700 mit 2,6 GHz getaktet. Performancetechnisch kann man sicher noch so Einiges aus dem Programm heraus holen. Es verwendet immer noch nicht das newtonsche Gesetz Nr 3 Actio= Reactio, das sollte die Leistung wohl noch mal verdoppeln (Gar nicht so einfach, wenn der Code selbst unabhängig von der vorgegeben Anzahl an simulierten Körpern sein soll). Von so Spielchen wie Multithreading, und Codeoptimierung gar nicht zu reden.

Auf jeden Fall gehört die von mir verwendete Sprache nicht zu den Schnellsten. Wäre vielleicht ganz interessant was man da noch an Geschwindigkeit herauskitzeln kann, hat aber erstmal noch wenig Priorität.
Die Grafik ist die größte Baustelle, Texturierung funktioniert noch nicht und du hast recht, die Bahnspur wird momentan nicht vom Mittelpunkt des Planeten gezeichnet. Des weiteren fehlen mir jegliche Hilfslinien zur Entfernungsbestimmung, und die Kamerasteuerung per Maus möchte ich auch gerne implementieren. Das sind alles Sachen die relativ schnell gehen sollten so denn ich denn Zeit und Muße habe.
Große Schwierigkeiten macht im Moment noch die Löschung von Testkörpern im laufenden Programm nach einer erkannten Kollision, das scheint ziemlich tricky, ist aber eine Aufgabe die noch mal gehörig Spaß für mich in das Projekt bringt sobald es da Fortschritte gibt.

Definitiv möchte ich noch mehr Integrationsalgorithmen ausprobieren, der RK liefert nicht die Ergebnisse, die Ich mir erhofft habe. Bernhard, wäre es möglich, dass du mal dieses symplektische Integrationsverfahren erklärst?

mfg
 

Bernhard

Registriertes Mitglied
Bernhard, wäre es möglich, dass du mal dieses symplektische Integrationsverfahren erklärst?
Hallo Kibo,

Du findest hier http://en.wikipedia.org/wiki/Symplectic_integrator eine sehr gute Einführung in dieses Thema und könntest einen einfachen symplektischen Integrator mit k=2 implementieren (Verlet-Methode). Dazu zerlegst Du bei dem besprochenen Zweikörperproblem die gesamte Hamiltonfunktion in die folgenden zwei Teile:

H_T = p1^2/(2*m1) + p2^2 / (2 * m2)

und

H_V = - G * m1 * m2 / r mit r := |r1 - r2|

Bei der Verlet-Methode integriert man dann (beispielsweise mit der Runge-Kutta-Methode) zuerst das System mit H_T um eine halbe Schrittweite (also mit 0.5 * dt), dann H_V um eine ganze Schrittweite (also dt) und zuletzt nochmal H_T um eine halbe Schrittweite. Bei der Integration von H_V nimmt man bei den Ortsvariablen diejenigen Orte, die man zuvor mit der ersten Integration über 0,5 * dt berechnet hat. Das ist der eigentliche Trick an dieser Methode.

Da H_T nicht vom Ort abhängt und H_V nicht von den Impulsen, vereinfachen sich die hamiltonschen Gleichungen so stark, dass deren Integration (z.B. mit Runge-Kutta) entsprechend schnell vom Rechner durchgeführt werden kann.

Bevor man diese Methode implementiert muss/sollte man sich also die hamiltonschen Gleichungen für H_T und für H_V getrennt notieren und für beide Systeme dann einen normalen Integrator (z.B. Runge-Kutta) mit beliebiger Schrittweite implementieren.

Das n-Körper-Problem geht dann ganz analog. H_T und H_V enthalten dann aber entsprechend mehr Terme.
MfG
 

Kibo

Registriertes Mitglied
Hallo,

Es gibt mal wieder ein kleines Update. Ich habe Newton Nr 3 implementiert, dies hat erstaunlicherweise zu einer veränderten Bahnberechnung geführt, die Bahnen bei RK sind stabiler geworden. Ich weiß nicht so recht ob mich das freuen soll (da unvorhergesehen). Des weiteren gibt es ein paar zusätzliche Features:

Schrittweite ist jetzt alternativ auch übe Edit-feld unbegrenzt wählbar
Hotkeys werden nur abgefragt wen das Fenster auch aktiv ist, so behindert das Programm einem nicht mehr bei anderen Arbeiten

Es gibt eine Menüleiste mit verschiedenen Funktionen:
-Man kann Planetendurchläufe speichern und auch Dateien neu laden
-Man kann die Anzahl der Frame-Aktualisierungen pro Sekunde anpassen
-Man kann die Länge der Bahnspur selbst bestimmen

Ich werde als Download jetzt diese öffentlichen Ordner bereit stellen, sodass sich der Downloadlink jetzt nicht mehr ständig ändert.

Rapidshare

Skydrive


mfg
 

Kibo

Registriertes Mitglied
Version 1.6 ist jetzt im Ordner verfügbar.

- Fenstergröße ist beliebig veränderbar
- Hilfsgitter zur Orientieung (An/Aus über Grafikeinstellungen)
- Man kann einen Hilfsorbit (Kreis) um einen Körper mit beliebigen Radius zeichnen (An/Aus über Grafikeinstellungen)
- Steuerung wurde stark verändert. Blickrichtung der Kamera ist jetzt mit der Maus steuerbar, WSAD steuerung hängt von der Blickrichtrung ab (Firstperson Steuerung, noch buggy)
- Kamera beschleunigt um so mehr, je öfter/länger man eine WSAD-Taste gedrückt hält, beschleunigung wird durch drücken einer anderen Taste zurückgesetzt
- Man kann den Radius der Himmelskörper um einne beliebigen Faktor anpassen (über Grafikeinstellungen)
mfg
 
Zuletzt bearbeitet:

Kibo

Registriertes Mitglied
Nabend Bernhard,

Wie du dir sicher denken kannst, war das nicht so geplant, sorry. Die Logdatei ist jetzt (Version 1.6.2) standardmäßig aus und kann über den Punkt Feature-Einstellungen bei Bedarf aktiviert werden.

mfg

btw hast du eigentlich noch ein paar hübsche Systeme kreiert?
 
Zuletzt bearbeitet:

Bernhard

Registriertes Mitglied
btw hast du eigentlich noch ein paar hübsche Systeme kreiert?
Hallo Kibo,

vorerst Nein. Die Version 1.6 zeigt übrigens keine Orbits mehr, wenn man obige Planeten.txt mit den zwei gleich schweren Massepunkten verwendet :( . Das solltest Du unbedingt fixen. Die neuen Linien kann man ohne Anleitung leider kaum zuordnen, bzw. nutzen. Das Programm ist in dieser Hinsicht also nicht selbsterklärend.

Leider noch ein Kritikpunkt: Solange das Programm keine Startwerte mit Null einlesen kann, ist meine Motivation neue Planeten.txt-Dateien anzulegen eher gleich Null, weil man so praktisch gar nicht mehr einschätzen kann, wo welche Fehler wirken.
MfG
 

Kibo

Registriertes Mitglied
Hmm, also bei mir zeigt er die Bahnen einwandfrei an, sowohl bei Erde-Mond als auch bei den 2 Sonnen von dir.

Zu den Linien, das sind Hilfslinien um die Abstände einschätzen zu können. Sie liegen jeweils genau 1 AE auseinander.

EDIT: Der 0er Bug ist soweit ich das sehe jetzt gefixt (1.6.3)

Ich habe übrigens noch ein Feature unterschlagen: Man kann die Koordinaten der Kamera im linken Formular anzeigen lassen (obere Felder) und verändern (die unteren)
 
Zuletzt bearbeitet:

Bernhard

Registriertes Mitglied
EDIT:

Hallo Kibo,

ich habe eben die Positionierung der Kamera bei der Version 1.6.3 ausprobiert. Wenn man in die ersten drei Kästchen die Zahlen 0, 0, und 150000000000 einsetzt und dann zusätzlich noch den Zeiger der Maus korrekt positioniert, sieht man wieder die berechneten Bahnen. Dabei habe ich leider festgestellt, dass die Bahnen mit dem Runge-Kutta-Verfahren noch ziemlich falsch sind. Bei dem Beispiel mit den zwei Sonnen entstehen keine kreisrunden Bahnen. Bei dem Euler-Verfahren passt es. Dort ergeben sich die kreisrunden Bahnen.

Als weiteren Test habe ich geeignete Startbedingungen für vier Sonnen berechnet:
Code:
KoordX	KoordY	Koordz	ImpulsX	ImpulsY	ImpulsZ	Masse		dichte	s	v	Name	s und v sind noch nicht implementiert
-8.228e10	0.0	0.0	0.0	-96392.0	0.0	1.99e30	0	0	Sonne1	1.408
-6.732e10	0.0	0.0	0.0	36816.0	0.0	1.99e30	0	0	Sonne2	1.408
6.732e10	0.0	0.0	0.0	-36816.0	0.0	1.99e30	0	0	Sonne3	1.408
8.228e10	0.0	0.0	0.0	96392.0	0.0	1.99e30	0	0	Sonne4	1.408

Wenn man das mit dem Euler-Verfahren betrachtet, sieht man, dass die Symmetrie des Systems nicht erhalten bleibt. Da befindet sich also auch noch ein Programmierfehler.

Es wäre auch schön, wenn man mit dem Slider noch größere Schrittweiten vorgeben könnte. 2500s ist doch etwas klein für astronomische Systeme. Ich würde das Maximum auf 86400s setzen. Runge-Kutta kommt mit einer Schrittweite von einem Tag eigentlich sehr gut zurecht.
MfG
 
Zuletzt bearbeitet:

Kibo

Registriertes Mitglied
Guten Morgen,

Alternativ kannst du in dem Editfeld unter dem Slider eine Schrittweite festlegen. Ich werde mal schauen ob ich die Bahnen nicht kreisrund kriege.

Mfg
 

Bernhard

Registriertes Mitglied
Guten Nachmittag.

Alternativ kannst du in dem Editfeld unter dem Slider eine Schrittweite festlegen.
Bei einer Schrittweite von 86400 Sekunden spiralen die beiden Sonnen nach einigen wenigen Umläufen in Richtung gemeinsamer Schwerpunkt. D.h. da ist noch irgendwo ein Fehler im Runge-Kutta-Code.

Ich werde mal schauen ob ich die Bahnen nicht kreisrund kriege.
Bei den vier Sonnen, die sich paarweise gegen den Uhrzeigersinn drehen, ist beim Euler-Verfahren die Sonne Nr. 4 gravitativ zu schwach an den Partner (Sonne 3) gekoppelt und bewegt sich tangential zu gerade weiter. Sonne 1 wird hier richtig berechnet und kann deswegen beim Debuggen als "Vorbild" für Sonne4 verwendet werden. Es müssen nur die Vorzeichen von Pos_Y und Vel_Y vertauscht werden. Der Runge-Kutta-Code steigt bei den vier Sonnen leider so ziemlich komplett aus (grins).
MfG
 
Zuletzt bearbeitet:

Bernhard

Registriertes Mitglied
Bernhard, wäre es möglich, dass du mal dieses symplektische Integrationsverfahren erklärst?
Hallo Kibo,

ich habe mir eben das Swifter-Paper von Duncan, Levison und Lee nochmal angesehen. Gleichung 10 ist hilfreich, wenn man den Leapfrog-Integrator implementieren will. Man muss bei Deinem Programm allerdings beachten, dass Gleichung 10 hierbei für unterschiedlich viele Körper auch unterschiedliche Gleichungen impliziert.
MfG
 

Kibo

Registriertes Mitglied
OK, danke für die Hinweise Bernhard. Ohne dich wäre ich wahrscheinlich noch beim Rumeulern von undefinierten Murmeln:)
Ich denke es kann noch ein bisschen dauern, bis Version 1.7 fertig ist. An die Formeln muss Ich mich erst mal heranarbeiten.

mfg
 

Bernhard

Registriertes Mitglied
An die Formeln muss Ich mich erst mal heranarbeiten.
Hallo Kibo,

ich möchte dazu einen Teil des ersten Schrittes des Runge-Kutta-Verfahrens vorrechnen und verwende dazu die Bezeichnungen von hier: http://en.wikipedia.org/wiki/Runge–Kutta_methods#Common_fourth-order_Runge.E2.80.93Kutta_method, sowie das Beispiel mit den zwei Sonnen von oben, wobei ich anstelle 0.001 0.0 verwende, damit die Zahlen übersichtlich bleiben.

y ist ein zwölfdimensionaler Vektor, den man sinnvollerweise wie folgt definiert: y = ( x_1, v_1, x_2, v_2 ). Jeder der vier Einträge ist selbst ein dreidimensionaler Vektor.

Für k1 ergibt sich dann der folgende Vektor:
k1 = ( 0.0, 2.105e4, 0.0, 5.9308e3, 0.0, 0.0, 0.0, -2.105e4, 0.0, -5.9308e3, 0.0, 0.0 )

Bei einer Schrittweite von 10000 Sekunden ergibt sich:
y_n + 0.5 * h * k1 = (-7.48e10, 1.0525e8, 0.0, 2.9654e7, 2.105e4, 0.0, 7.48e10, -1.0525e8, 0.0, -2.9654e7, -2.105e4, 0.0)

Diese Werte solltest Du auch bei dem ersten Schritt innerhalb der Implementierung finden. Mit dem letzten 12er-Zahlentupel kann man dann k2 berechnen und dann analog den ganzen ersten Schritt berechnen.
MfG
 

Kibo

Registriertes Mitglied
Nabend Bernhard,

Also für K1 ergeben sich bei mir folgende Werte
Steigung m1 für Sonne1: (2,967124.../21050/0)

Daraus berechneter Prognosepunkt P1: (-74799998516,438034.../10525000/0)

Das sieht beides für mich richtig aus, obwohl es anscheinend nicht mit deinen Werten übereinzustimmen scheint.
x-wert von m1 ist: Startgeschwindigkeit 0 + berechnete Gravitation von 0,0059342478...m/s² * Schrittweite 1000 /2
 
Zuletzt bearbeitet:

Bernhard

Registriertes Mitglied
Nabend Bernhard,

Also für K1 ergeben sich bei mir folgende Werte
Steigung m1 für Sonne1: (2,967124.../21050/0)
Guten Morgen Kibo,

dann lass uns mal den Wert 2,967124... ansehen. Die Ableitung der Geschwindigkeit von Sonne 1 ist gleich der Kraft durch die Masse von Sonne 1 wegen F_1 = m_1 * a_1. Die Kraft berechnet sich gemäß dem newtonschen Gravitationsgesetz. Es gilt also insgesamt:

a_1 = F_1 / m_1 = G * m_2 / r_12^3 * (r_2 - r_1). Die Beschleunigung als Vektor zeigt am Anfang also in Richtung der positiven x-Achse.

Fehlt noch der Wert in SI-Einheiten:

a_1 = 6.67e-11 * 1.99e30 / (149.6e9)^2 = 5.93e-3 m/s^2

Damit habe ich zwar einen Fehler im Exponenten, aber die Zahl 5.93 stimmt und sollte so auch im Programmablauf bei K1 zu finden sein.
MfG

EDIT: Mein erster Prognosepunkt ist damit leider falsch. Die Korrektur muss ich mir später ansehen.
 
Zuletzt bearbeitet:
Oben