PDA

Archiv verlassen und diese Seite im Standarddesign anzeigen : kleiner planetarer Simulator



Kibo
27.09.2012, 20:58
Wie im anderen Thread besprochen, eröffne ich einen neuen Thread um keinen zu stören.
Also es gibt eine neue Version (1.9.1)

Hier die Neuerungen:

Die Buttons wurden durch Kreise ersetzt was viele Vorteile hat.
Die Initialisierung mit dynamischen Objekten ist wesentlich einfacher, das findet sich im Code noch nicht wieder.
Neben dem ästhetischen Faktor wird über den Kugelradius auch die Masse der Planeten wiedergegeben. Das läuft über die nicht ganz realitätsgetreue Formel.:
r :=sqrt(m)+5

Das Programm läuft bis auf Weiteres noch in 2D und mit GDI+. Der Downloadlink ist immer noch der Gleiche über Skydrive
(https://skydrive.live.com/embed?cid=58397D4453EFC662&resid=58397D4453EFC662%21155&authkey=AMsx3LIPCBfxaMc)

mfg

Bernhard
28.09.2012, 08:52
Also es gibt eine neue Version (1.9.1)
Hi Kibo,

nettes Programm. Man sieht sogar eine Perihelbewegung des innersten Planeten. Die rote Linie gehört vermutlich also zu Merkur? Hast Du den BS-Integrator verwendet?


Das läuft über die nicht ganz realitätsgetreue Formel.: r :=sqrt(m)+5
Die Sonne würde ich etwas kleiner darstellen, weil die Merkurbahn schon fast die "Sonnenoberfläche" streift.
Gruß

Kibo
28.09.2012, 12:53
Hallo Bernhard,

Der Code für die physikalischen Berechnungen ist komplett selbst geschrieben, ich wollte einfach meine Kenntnisse auf dem Gebiet auffrischen und habe mir daher meine eigenen Gedanken gemacht weil man dabei mehr lernt.
Wäre ich auf größere Hürden gestoßen, hätte ich mir natürlich euren Quellcode genauer angeschaut. Wenn du möchtest, schicke ich dir mal meinen Code zu.
Was die Größenverhältnisse angeht, steht da noch nicht fest wie ich es am Ende haben möchte. Die Körper sollen halt alle gut sichtbar sein aber eben auch die Massenverhältnisse wiederspiegeln.

Das Beispielsystem soll nicht unser Sonnensystem darstellen, ich habe einfach ein bisschen mit Zahlen herumgespielt bis ein halbwegs stabiles System da stand.
Die Werte für die Planeten können in der Datei Planeten.txt nach belieben angepasst werden. Man kann auch noch mehr Objekte hinzufügen (wobei das noch nicht ganz nach meiner Zufriedenheit funktioniert).

mfg

Bernhard
28.09.2012, 14:29
Wenn du möchtest, schicke ich dir mal meinen Code zu.
Hi Kibo,

da ich momentan etwas im Zeitdruck bin, wäre es mir lieber, wenn Du kurz beschreiben würdest, welche Integrationsmethode Du verwendet hast.
Gruß

Kibo
28.09.2012, 17:22
Nun ich habe 2 Schleifen, dort nehme ich einen Planeten, berechne Entfernung, Masse und Richtung zu einem anderen Planeten und addiere diesen Vektor mit dem schon vorhandenen Impuls. Das mache ich so lange bis jeder Körper zu jedem anderen in Beziehung steht. Das ist so in etwa der Kern.

Kibo
28.09.2012, 20:39
Es gibt ein etwas größeres Update auf Version 1.10!

Es wurde ein kleines Formular zur direkten Manipulation der Himmelskörper integriert, ausserdem kann man jetzt mit den hoch und runter Tasten zoomen.

Die Änderungen gehen mit einem erheblichen Gewinn an Benutzerfreundlichkeit (und Spaß ) einher, von daher kann ich den Download nur wärmstens empfehlen!

Downloadlink

(https://skydrive.live.com/embed?cid=58397D4453EFC662&resid=58397D4453EFC662%21158&authkey=ALOT_odIptuPmsw)mfg

Kibo
12.10.2012, 22:41
Guten Abend,

Es gibt mal wieder ein Update, Ich bin mitlerweile bei Version 1.12 (https://skydrive.live.com/embed?cid=58397D4453EFC662&resid=58397D4453EFC662%21162&authkey=AE1hWr8SAHliMvw) angelangt.
Diese ist nun endlich in der Lage mit einer beliebigen Anzahl an Testkörpern umzugehen. Ich habe den Physikteil erheblich verbessert, sodass sich keine Ungenauigkeiten aufgrund von veränderten Planetenanzahlen ergeben sollten. Ich habe soweit Ich das sehe nun alle Bugs behoben, damit ist das jetzt der stabile Release von Version 1.
Das ist auch der letzte Release der GDI verwenden wird für Version 2 in 3d steige ich dann auf OpenGl um.

Falsl jemand Lust hat, würde Ich gerne den Integrator besprechen, Ich wäre interessiert mehr zu dem Thema zu erfahren, mein jetziger Integrator ist ja komplett von mir selbst zusammen geschustert ohne das ich dazu auch nur eine Seite Theorie zu gelesen habe, daher weiß ich nicht einmal auf welchem Niveau an Genauigkeit und Schnelligkeit der nun arbeitet.

Bernhard
12.10.2012, 23:54
Falsl jemand Lust hat, würde Ich gerne den Integrator besprechen, Ich wäre interessiert mehr zu dem Thema zu erfahren, mein jetziger Integrator ist ja komplett von mir selbst zusammen geschustert ohne das ich dazu auch nur eine Seite Theorie zu gelesen habe, daher weiß ich nicht einmal auf welchem Niveau an Genauigkeit und Schnelligkeit der nun arbeitet.
Hi Kibo,

ohne explizite Vorkenntnisse hast Du ziemlich sicher das Eulersche Polygonzugverfahren (http://de.wikipedia.org/wiki/Eulersches_Polygonzugverfahren) implementiert. Sobald man weiß, was die Ableitung einer Funktion einer Veränderlichen ist, kann man das mehr oder weniger intuitiv verstehen. Etwas ausgefeilter ist das Runge-Kutta-Verfahren (http://de.wikipedia.org/wiki/Runge-Kutta-Verfahren). Versuche also erst mal Dein aktuell eingesetztes Verfahren anhand dieser Links zu verbessern.

Kibo
13.10.2012, 08:50
Hallo Bernhard,

Ich weiß nicht so genau ob das Eulersche Polygonzugverfahren meinem entspricht. Ich fasse einfach mal mein Verfahren zusammen, sonst kann man ja schlecht drüber reden.

Also ich habe 2 Schleifen, die dienen zum bestimmen welcher Planet gerade mit der Berechnung dranne ist. Zum beispiel wird gerade der Planet a von b beinflusst.
Ich speichere zuerst den jetzigen Impuls von a (liegt in a.ix und iy) in die Hilfsvariablen a.hx und a.hy, wichtig später für die Trägheit. Dan kommt die Gravitationsberechnung die b auf a ausübt, die läuft über eine Methode die in der Planetenklasse steht. Ich bestimme den Vektor zwischen a und b mit
this.xd:=gegner.x-this.x
this.yd:=gegner.y-this.y

Daraus berechne ich den Abstand zwischen a und b.
abstand:=sqrt(this.xd*this.xd+this.yd*this.yd)

aus dem Vektor und dem Abstand kann na dann leicht den Einheitsvektor bestimmen.
this.xn:=this.xd/abstand
this.yn:=this.yd/abstand

Ich berechne dann die Stärke der Gravitation aus dem Abstand und der Masse der Planeten
this.s:=0.667384*gegner.m*this.m/(abstand*abstand)

Als Nächste schreibe ich die Gravitation noch in einen Vektor...
this.gx:=this.s*this.xn
this.gy:=this.s*this.yn

und addiere den neuen Impuls auf den alten auf.
%Planeten%.ix := %Planeten%.hx + %Planeten%.gx /(genauigkeit * %Planeten%.m)
%Planeten%.iy := %Planeten%.hy + %Planeten%.gy /(genauigkeit * %Planeten%.m)
Genauigkeit soll hier die Schrittweite bestimmen. %Planeten% ist der Dynamische Aufruf für den jetzigen Planeten also a (innerhalb der Klasse war das dann this für a und Gegner für b)

das war der inhalt der inneren Schleife, die wird nun für jeden anderen Planeten noch für a ausgeführt und am Schluss, wenn man alle Gravitationseinflüsse für a berechnet hat, werden die neuen Koordinaten für a bestimmt
%Planeten%.x:=%Planeten%.x+ %Planeten%.ix/genauigkeit
%Planeten%.y:=%Planeten%.y+ %Planeten%.iy/genauigkeit

das Ganze wird durch die äußere Schleife noch für alle anderen in die Planeten.txt eingetragenen Planeten wiederholt.

Wäre nett wenn du mal drüber schauen würdest Bernhard:)

Bernhard
13.10.2012, 11:18
aus dem Vektor und dem Abstand kann na dann leicht den Einheitsvektor bestimmen.
this.xn:=this.xd/abstand
this.yn:=this.yd/abstand
Hallo Kibo,

den Einheitsvektor braucht man eigentlich gar nicht. Besser man speichert die Differenz zwischen dem Ort des aktuellen Planeten und dem Ort des "Gegners" in einem Hilfsvektor h. Dann berechnet man über das Skalarprodukt dieses Vektors mit sich selbst das Quadrat des Abstandes sqr_abstand und zieht dann die Wurzel aus diesem Wert, um zuletzt den Abstand zu bekommen.


Ich berechne dann die Stärke der Gravitation aus dem Abstand und der Masse der Planeten
this.s:=0.667384*gegner.m*this.m/(abstand*abstand)
Mit den eben erwähnten drei Variablen sieht das dann folgendermaßen aus:

Beschleunigung(Planet) = 0.667384 * gegner.m * h / (sqr_abstand * abstand)

Man bekommt dann die Beschleunigung gleich als Vektor und vermeidet auch, dass man aus einer Größe zuerst die Wurzel zieht und diese dann wieder quadriert.

ZUsätzlich kann man anstelle von der Masse eines Planeten mit der kombinierten Variablen 0.667384 * gegner.m rechnen. Die Gravitationskonstante wird also einfach in die Masse fest rein multipliziert. Das spart haufenweise überflüssige Multiplikationen. Die Masse alleine wird im gesamten Programm nämlich gar nicht benötigt.



Als Nächste schreibe ich die Gravitation noch in einen Vektor...
this.gx:=this.s*this.xn
this.gy:=this.s*this.yn

Dieser Schritt entfällt dann.



und addiere den neuen Impuls auf den alten auf.
%Planeten%.ix := %Planeten%.hx + %Planeten%.gx /(genauigkeit * %Planeten%.m)
%Planeten%.iy := %Planeten%.hy + %Planeten%.gy /(genauigkeit * %Planeten%.m)
Genauigkeit soll hier die Schrittweite bestimmen. %Planeten% ist der Dynamische Aufruf für den jetzigen Planeten also a (innerhalb der Klasse war das dann this für a und Gegner für b)
Die Genauigkeit ist offensichtlich das Inverse der Schrittweite dt. Hier würde ich besser direkt mit der Schrittweite dt im Zähler rechnen. Das macht den Code übersichtlicher und vermeidet erneut unnötiges Ausrechnen. Wenn Du dem Anwender als Parameter lieber das Inverse der Schrittweite geben willst, musst Du eben beim Start des Programmes einmal den Kehrwert berechnen, um dt zu kennen.



das war der inhalt der inneren Schleife, die wird nun für jeden anderen Planeten noch für a ausgeführt und am Schluss, wenn man alle Gravitationseinflüsse für a berechnet hat, werden die neuen Koordinaten für a bestimmt
%Planeten%.x:=%Planeten%.x+ %Planeten%.ix/genauigkeit
%Planeten%.y:=%Planeten%.y+ %Planeten%.iy/genauigkeit

Das ist das genau das Euler-Verfahren, falls genauigkeit = 1.0 / dt.


das Ganze wird durch die äußere Schleife noch für alle anderen in die Planeten.txt eingetragenen Planeten wiederholt.
Hier fehlt noch das dritte newtonsche Gesetz (actio gegengleich reactio). Du berechnest in jedem Schritt sowohl actio, als auch reactio getrennt. Durch geschicktes Ordnen der Schleife berechnet man einmal actio (oder reactio) und verwendet diesen Wert dann sowohl mit positivem als auch mit negativem Vorzeichen. Man braucht dazu eigentlich nur für jeden Planeten eine vektorielle Variable für die Beschleunigung und läßt die äußere Schleife immer ab dem aktuellen Planeten starten. Die äußere Schleife wird dann mit fortschreitender innerer Schleife immer kleiner und man berechnet keine Zwischenergebnisse mehr doppelt.

Sämtliche Tips (mit unterschiedlichem Schwierigkeitsgrad) kann man in den genaueren Verfahren (Runge-Kutte, Bulirsch-Stoer, Symplektischer Integrator) beibehalten.
Gruß

Kibo
14.10.2012, 18:00
Hallo Bernhard,

Erst einmal danke für die Tipps. Die Masse außerhalb der Schleife mit der Gravitationskonstante zu verrechnen ist kein Problem, auch die kleine Anpassung mit der Schrittweite stellt mich vor keinem Hindernis, das mit der Beschleunigung direkt als Vektor ist allerdings neu für mich. Wie sieht denn so ein Beschleunigungsvektor aus, ist das ein Array?
Hab ich deinen letzten Tipp richtig verstanden? Du berechnest quasi gleichzeitig die Beeinflussung von Planet a durch b und b durch a ja? Das leuchtet mir ja ein, aber wozu ist da zwingend eine vektorielle Beschleunigungsvariable nötig.

Ich setze erstmal deine Tipps um, so wie ich sie verstanden habe, parallel dann zu OpenGl und Kollisionsdetektion.

mfg

Bernhard
14.10.2012, 18:42
Wie sieht denn so ein Beschleunigungsvektor aus, ist das ein Array?
Von der Mathematik her sind das nur drei reelle Zahlen für die drei Raumdimensionen. SRMeister hatte eine eigene Klasse PVektor dazu geschrieben. Wie Du das datentechnisch im Detail umsetzt, bleibt eigentlich Dir überlassen.


Hab ich deinen letzten Tip richtig verstanden? Du berechnest quasi gleichzeitig die Beeinflussung von Planet a durch b und b durch a ja?
Genau. Es geht allerdings hauptsächlich darum Zwischenergebnisse nicht doppelt zu berechnen, da sonst das Programm unnötig langsam wird.


Das leuchtet mir ja ein, aber wozu ist da zwingend eine vektorielleBeschleunigungsvariable nötig.
Der Beschleunigungsvektor dient nur der Übersichtlichkeit. Du musst das nicht so machen.
Gruß

SRMeister
14.10.2012, 22:19
Das Programm gefällt mir auch.
So wie du, hatte ich auch zu Beginn den intuitiven Ansatz gewählt und mich Schritt für Schritt weiter vertieft. Auch ich/wir haben nicht alles an einem Tag geschafft, also, es ist ganz normal, wenn man nicht alles sofort umsetzen kann. (Bernhard hat dir immerhin die komprimierte Schlussfolgerung aus mehreren Monaten gegeben)
Eine Anmerkung möchte ich noch zum Eulerverfahren machen: Prinzipbedingt summieren sich die Ungenauigkeiten "in eine Richtung" - das Bedeutet, die Fehler sind immer positiv und addieren sich ständig. Besser ist ein Verfahren welches zufällig Fehler verursacht, die mal positiv und mal negativ sind, damit sie sich in Summe eher aufheben.

So ein Verfahren ist das Midpoint-Verfahren und das sollte auch dein nächster Schritt sein. Hast du die Verallgemeinerung von Euler zu Midpoint verstanden/umgesetzt, sind die anderen Integratoren, wie Runge-Kutta und Bulirsch-Stoer nicht mehr weit entfernt.

Viel Spaß

Kibo
15.10.2012, 10:05
Danke SRMeister,

das wäre eh mein nächster Schritt gewesen, das Verfahren zu wechseln. Nun weiß ich dann auch welches Verfahren ich als Nächstes nehme.

mfg

Bernhard
15.10.2012, 13:53
Besser ist ein Verfahren welches zufällig Fehler verursacht, die mal positiv und mal negativ sind, damit sie sich in Summe eher aufheben.
Die Gesamtenergie erfährt beim Mittelpunktsverfahren angeblich ebenfalls eine Fehlerdrift, weswegen ein symplektischer Integrator zusätzliche Vorteile bringt. Diese Integratoren besitzen, so wie die anderen Verfahren auch, eine bestimmte Fehlerordnung, weswegen man auch beim symplektischen Integrator ein Extrapolationspolynom verwenden kann (s.a. Richardson-Extrapolation (http://en.wikipedia.org/wiki/Richardson_extrapolation)).

EDIT: Um "Berührungsängste" mit dem symplektischen Integrator zu vermeiden kann man sich vergegenwärtigen, dass es sich dabei auch nur um sehr geschickte Differenzenverfahren handelt.

Kibo
11.11.2012, 21:34
Guten Abend,

Ich habe mich dieses Wochenende mal wieder kurz ran setzen können.
Das (https://skydrive.live.com/embed?cid=58397D4453EFC662&resid=58397D4453EFC662%21163&authkey=AD7MpP7c6nxwvaQ) ist dabei raus gekommen. Wie man sieht, jetzt mit Open-Gl in 3 Dimensionen. Die Texturierung funktioniert noch nicht wie sie soll und bei der Platzierung der Sphären scheint es auch noch einen Fehler zu geben. Des weiteren gab es ein paar kleinere Probleme mit dem klassischen Runge-Kutta, weswegen in der Version das gute alte stabile Euler-verfahren werkelt, auf Grundlage von Bernhards Tips angepasst. Die Kamerasteuerung war für mich daher noch zweitrangig weswegen man mit WSAD nur zoomen und drehen kann (eher weniger als mehr bugfrei). Wie man sieht, hab ich im Moment nicht so viel Zeit, sobald mich die Muse küsst werden die aufgezählten Mängel demnächst behoben.

mfg

Bernhard
11.11.2012, 22:41
Das (https://skydrive.live.com/embed?cid=58397D4453EFC662&resid=58397D4453EFC662%21163&authkey=AD7MpP7c6nxwvaQ) ist dabei raus gekommen.
Nette Farben und Bahnen :) . Der Zentralkörper ist für meinen Geschmack noch immer (viel) zu groß, weil er größere Teile der Bahnen verdeckt. Hast Du bestimmte Exoplaneten im Sinn, die Du damit berechnen willst?

Kibo
27.11.2012, 23:10
Realitätstreue wäre programmiertechnisch kein Problem, blos dann erscheinen die Planeten im Vergleich zu den unendlichen Weiten des Weltalls doch wirklich sehr klein, nämlich um die Pixelgröße. Vielleicht bau ich noch einen Regler im GUI für ein. Wie gesagt, ich programmier das Teil nur so zum Spaß ohne ein bestimmtes Motiv dahinter.
Ich komm im Moment mit Quadric Texturen nicht wirklich klar, hab halt keine Ahnung von OpenGl ^^''
Habe im Moment relativ wenig Freizeit, so dass fürs Programmieren eher kaum was am Tag übrig bleibt. Ausserdem arbeite ich gerade noch an einer Client-Server Geschichte, die ist Momentan interessanter, von daher ruht der Simu jetzt etwas. Guter Hefeteig braucht ja auch viel Zeit zum aufgehen.:)

mfg

Kibo
08.02.2013, 15:55
Hallo,

Ich bin mal wieder dabei meinen Simulator ein bisschen aufzupolieren und rechne gerade am Beispielsystem Erde-Mond.

Bei mir hakt es im Moment an der Formel der Bahngeschwindigkeit.

wurzel((GM+m)(2/r-1/a))
G=Gravitationskonstante
M=Masse Erde
m= Masse Mond
r=a= große Halbachse

ich setze ein: G*(5,974*10^24+7,349*10^22)*(1/384400)= 32403

Richtige Lösung wäre meiner Meinung nach 1,02 km/s bzw 3672 km/h. Sieht wer meinen Fehler?

mfg

UMa
08.02.2013, 15:57
Hallo Kibo,

r und a in Metern statt km eingeben, dann passt es.

Grüße UMa

Kibo
08.02.2013, 17:19
Stimmt, danke

Kibo
09.02.2013, 15:04
Hmm erstaunlicherweise werden beim Runge Kutta die Bahnen ganz schön instabil im vergleich zum Euler.

Wer mal testen will hier die 2 verschiedenen Versionen

Mirrors

Rapidshare
Runge Kutta (http://rapidshare.com/files/2235517210/planeten%201.3%20RK.exe)
Euler (http://rapidshare.com/files/2334876395/planeten%201.3EU.exe)
Planeten Datei für Erde-Mond (http://rapidshare.com/files/2570060638/planeten.txt)


Skydrive
Alle 3 Dateien (https://skydrive.live.com/redir?resid=58397D4453EFC662!200&authkey=!ANztCmGQzaVHqzA)

Anleitung: Jeweilige exe und die Planeten.txt im selben Ordner herunterladen und exe starten. Wenn die Planeten.txt für Erde-Mond nicht bei der Exe Datei mit zu liegt, wird ein Standartsystem erstellt. Die Szenerie wird mit WASD QY und den Pfeiltasten gesteuert.

mfg

Bernhard
10.02.2013, 16:45
Hallo Kibo,


Hmm erstaunlicherweise werden beim Runge Kutta die Bahnen ganz schön instabil im vergleich zum Euler.
ich vermute da eher einen Programmierfehler. Poste hier, bei Bedarf, die entscheidenden Stellen deines Quelltextes. Dann kann man sich das genauer ansehen.



Runge Kutta (http://rapidshare.com/files/2235517210/planeten%201.3%20RK.exe)
Euler (http://rapidshare.com/files/2334876395/planeten%201.3EU.exe)
Planeten Datei für Erde-Mond (http://rapidshare.com/files/2570060638/planeten.txt)

Diese drei Links funktionieren bei mir nicht.
MfG

Kibo
10.02.2013, 18:31
Ok, Ist gefixt. Da ich die Links nicht mehr ändern kann muss ich sie halt noch mal neu rein stellen:

Rapidshare

Euler (http://rapidshare.com/files/2334876395/planeten%201.3EU.exe)

Runge Kutta (http://rapidshare.com/files/2235517210/planeten%201.3%20RK.exe)

Planeten.txt (http://rapidshare.com/files/2570060638/planeten.txt)

mfg

Bernhard
10.02.2013, 19:13
Ich komme mit der Bedienung des Programmes überhaupt nicht klar.
In welchen physikalischen Einheiten arbeitet das Programm?
Wie kann man die Startbedingungen vorgeben?
Was bedeuten Zeitangaben in Millisekunden (!)?
MfG

EDIT: Hier noch ein Link auf eine vorbildliche Lösung deiner Aufgabenstellung: http://phet.colorado.edu/sims/my-solar-system/my-solar-system_en.html .

Kibo
11.02.2013, 08:37
Hallo Bernhard,

Die Startbedingungen gibt man über die Datei Planeten.txt vor:

Koordinaten in x y z Impuls in x y z Masse s v Dichte des Planeten

Die Koordinaten sind in Kilometern anzugeben, der Impuls gibt quasi Startgeschwindigkeit und Richtung eines Planeten vor und bezieht sich im Moment noch auf das Kilometer/Berechnungsschritt. Die Masse sollte in Kilogramm angeben werden und ist im Moment noch der Übersichtlichkeit halber um 20 Nullen gekürzt, man kann auch solche Werte wie 5.7e11 eingeben, daher wird das in der nächsten Version angepasst. Die werte s v bitte auf 0 lassen, das funktioniert noch nicht so richtig. Die Dichte wird in g/cm³ angegeben.

Das Formular im Programm ist analog zu der Textdatei zu betrachten, und ist gedacht um die aktuell berechneten Werte abzulesen und anzupassen.

Das Feature Maus funktioniert nicht mehr, bitte ignorieren.
Der Slider Genauigkeit ermöglicht das Anpassen der Schrittweite, Der Slider geht von 1-100 und wird im Programm durch 5 gerechnet und dann als Kehrwert für die Schrittweite genommen. Das bedeutet Geneauigkeit 1 entspricht einer Schrittweite von 5 Einheiten (Kilometern?) und 100 entspricht der Schrittweite 1/ (100/5) 0.05.

Für die schlechte Beschriftung möchte ich mich entschuldigen, das Programm ist ja noch in Arbeit.

mfg

Bernhard
11.02.2013, 10:58
Hallo Kibo,


Die Startbedingungen gibt man über die Datei Planeten.txt vor:
die hatte mir gefehlt und damit wird jetzt alles gleich viel verständlicher :) . Nett wäre, wenn die Programme genau diese Datei selbst generieren, falls sie fehlt. Die aktuelle Default-Datei liefert ein Szenario, das eher an die Ergebnisse eines Teilchenbeschleunigers erinnert.


Die Masse sollte in Kilogramm angeben werden und ist im Moment noch der Übersichtlichkeit halber um 20 Nullen gekürzt, man kann auch solche Werte wie 5.7e11 eingeben, daher wird das in der nächsten Version angepasst.
Solche Spezialitäten solltest Du unbedingt vermeiden, da man die bei längeren Pausen sehr gerne vergisst. Es sei denn, Du machst Dir sehr genaue Notizen. Mir ist die letzten Tage selbst aufgefallen, wie wichtig Kommentare gerade in privaten Quelltexten sind.


Die werte s v bitte auf 0 lassen, das funktioniert noch nicht so richtig. Die Dichte wird in g/cm³ angegeben.

Das Formular im Programm ist analog zu der Textdatei zu betrachten, und ist gedacht um die aktuell berechneten Werte abzulesen und anzupassen.

Das Feature Maus funktioniert nicht mehr, bitte ignorieren.
Der Slider Genauigkeit ermöglicht das Anpassen der Schrittweite, Der Slider geht von 1-100 und wird im Programm durch 5 gerechnet und dann als Kehrwert für die Schrittweite genommen. Das bedeutet Geneauigkeit 1 entspricht einer Schrittweite von 5 Einheiten (Kilometern?) und 100 entspricht der Schrittweite 1/ (100/5) 0.05.

Für die schlechte Beschriftung möchte ich mich entschuldigen, das Programm ist ja noch in Arbeit.
OK, Rest ist erst mal Baustelle.

Jetzt zum Programm:
Bei dieser Auflösung (=visuell) und einem Umlauf sollte das Euler-Verfahren und das Runge-Kutta-Verfahren genau das gleiche Ergebnis berechnen. Du solltest den Quelltext des Runge-Kutta-Verfahrens nach kleinen aber wichtigen Fehlern durchsuchen. Wenn Du nicht weiterkommst, können wir den ersten Zeitschritt (ein Delta t) hier diskutieren, denn den kann man auch mit dem Taschenrechner nachrechnen.
MfG

Bernhard
11.02.2013, 14:41
Hallo Kibo,

um den Fehler in deinem Programm zu finden schreibe ich erst mal die DGL an (leider ohne LaTeX):



\dot{\vec{x}_1} = \vec{v}_1
\dot{\vec{v}_1} = G m_2 \frac{\vec{x_2}-\vec{x_1}}{\| \vec{x_2}-\vec{x_1} \|^3}

\dot{\vec{x}_2} = \vec{v}_2
\dot{\vec{v}_2} = G m_1 \frac{\vec{x_1}-\vec{x_2}}{\| \vec{x_1}-\vec{x_2} \|^3}

Es ist eine gewöhnliche DGL 1. Ordnung mit den 12 Variablen x_1, v_1, x_2 und v_2. Der erste Schritt berechnet den Wert dieser Variablen nach der Zeit dt. Bei einem ersten Test, um den Fehler zu finden, würde ich dt vergleichsweise groß wählen, z.B. 1 Tag = 86400 Sekunden.

Man kann dann die vier Hilfsvariablen k_1, k_2, k_3 und k_4 des Runge-Kutta-Verfahrens (http://de.wikipedia.org/wiki/Runge-Kutta-Verfahren) berechnen und testen, ob das Programm diese Variablen korrekt berechnet.

Beim ersten Überfliegen der Gleichungen sollte man bemerken, dass in allen vier Teilgleichungen der DGL die Gravitationskonstante immer als Produkt mit der Masse der Erde oder der Masse des Mondes vorkommt. Man kann die Gravitationskonstante also mit diesen zwei Massen schon vor der eigentlichen Berechnung multiplizieren und das Ergebnis als Masse verwenden. Die Gravitationskonstante taucht im Programm dann vorerst gar nicht mehr auf, was die Fehlersuche etwas vereinfacht.

Dann sollte man noch bemerken, dass man direkt mit Orts- und Geschwindigkeitsvektoren rechnen kann. Der Umweg über eine Impulsvariable ist momentan eher hinderlich. Die Werte aus der Datei planeten.txt kann ich zudem nicht nachvollziehen. Bei einer mittleren Orbitalgeschwindigkeit des Mondes von 1 km/s und einer Masse von 7.xeXX würde ich einen Wert von 7.xeXX oder etwas vergleichbares als Impuls erwarten.
MfG

Kibo
11.02.2013, 20:04
Hallo Bernhard,

Also deine Differentialgleichung kommt mir irgendwie spanisch vor :D Ich poste einfach mal den unterschiedlichen code. Ich hab ihn mal vonn Loops und Klassen befreit so dass man nicht viel umdenken muss, das ändert nichts am Ergebnis. Kommentare sind mit ; abgetrennt

Euler:


%Planeten%.hx:= %Planeten%.ix ; speichern des alten Impulses für trägheit
%Planeten%.hy:= %Planeten%.iy
%Planeten%.hz:= %Planeten%.iz

;steht in klasse this = planet um den es geht, gegner= beeinflussender planet

this.xd:=gegner.x-this.x ;Berechnung des Abstands zwischen den Planeten
this.yd:=gegner.y-this.y
this.zd:=gegner.z-this.z

abstand:=sqrt(this.xd*this.xd+this.yd*this.yd+this .zd*this.zd) ;bis hier

this.xn:=this.xd/abstand ; Einheitsvektor
this.yn:=this.yd/abstand
this.zn:=this.zd/abstand

this.s:=gegner.gm/(abstand*abstand) ; Berechnen der einwirkinden Gravitationskräfte

this.gx:=this.s*this.xn ; Verechnung der Grav mit der Richtung
this.gy:=this.s*this.yn
this.gz:=this.s*this.zn

;Ende Klassenteil

if (plzB=1)
%Planeten%.abstand:=abstand ; Speichert den Abstand zum Hauptkörper (idR Stern)


;aufaddieren des berechneten kraftvektors auf den vorhandenen impuls
%Planeten%.ix := %Planeten%.hx + %Planeten%.gx/genauigkeit
%Planeten%.iy := %Planeten%.hy + %Planeten%.gy/genauigkeit
%Planeten%.iz := %Planeten%.hz + %Planeten%.gz/genauigkeit

;Oben = Innerer schleifenteil, wird für jeden anderen Körper im System durchgeführt

%Planeten%.x:=%Planeten%.x+ %Planeten%.ix/genauigkeit ; Berechnung der neuen Koordinaten aus den summierten Kräften
%Planeten%.y:=%Planeten%.y+ %Planeten%.iy/genauigkeit
%Planeten%.z:=%Planeten%.z+ %Planeten%.iz/genauigkeit



Teil 2 folgt

Kibo
11.02.2013, 20:28
Nächste Ladung, klassisches Runge Kutta:




%Planeten%.hx:= %Planeten%.ix ; speichern des alten Impulses für trägheit
%Planeten%.hy:= %Planeten%.iy
%Planeten%.hz:= %Planeten%.iz

;Klassenteil ist der selbe daher siehe Teil 1

if (plzB=1)
%Planeten%.abstand:=abstand


%Planeten%.P1ix := %Planeten%.hx + %Planeten%.gx/(genauigkeit*2) ; %Planeten%.g/(genauigkeit*2 )=m1 prognose punkt 1 impuls (halbe schrittweite)
%Planeten%.P1iy := %Planeten%.hy + %Planeten%.gy/(genauigkeit*2)
%Planeten%.P1iz := %Planeten%.hz + %Planeten%.gz/(genauigkeit*2)

%Planeten%.P1x:=%Planeten%.x+ %Planeten%.P1ix/(genauigkeit*2) ;prognose punkt 1 koordinaten
%Planeten%.P1y:=%Planeten%.y+ %Planeten%.P1iy/(genauigkeit*2)
%Planeten%.P1z:=%Planeten%.z+ %Planeten%.P1iz/(genauigkeit*2)

%Planeten%.xd:=%PlanetenB%.x-%Planeten%.P1x ; abstand
%Planeten%.yd:=%PlanetenB%.y-%Planeten%.P1y
%Planeten%.zd:=%PlanetenB%.z-%Planeten%.P1z

abstand:=sqrt(%Planeten%.xd*%Planeten%.xd+%Planete n%.yd*%Planeten%.yd+%Planeten%.zd*%Planeten%.zd) ; immernoch abstand
%Planeten%.xn:=%Planeten%.xd/abstand
%Planeten%.yn:=%Planeten%.yd/abstand
%Planeten%.zn:=%Planeten%.zd/abstand

%Planeten%.s:=%PlanetenB%.gm/(abstand*abstand) ; gravitation

%Planeten%.g2x:=%Planeten%.s*%Planeten%.xn ;%Planeten%.g/(genauigkeit*2 )=m2 ; gravitationsvektor
%Planeten%.g2y:=%Planeten%.s*%Planeten%.yn
%Planeten%.g2z:=%Planeten%.s*%Planeten%.zn

%Planeten%.P2ix := %Planeten%.hx + %Planeten%.g2x/(genauigkeit*2) ; %Planeten%.g/(genauigkeit*2 )=m2 ;Prognosepunkt 2 impuls
%Planeten%.P2iy := %Planeten%.hy + %Planeten%.g2y/(genauigkeit*2)
%Planeten%.P2iz := %Planeten%.hz + %Planeten%.g2z/(genauigkeit*2)

%Planeten%.P2x:=%Planeten%.x+ %Planeten%.P2ix/(genauigkeit*2) ;teil 2 ; Prognosepunkt 2 koordninaten
%Planeten%.P2y:=%Planeten%.y+ %Planeten%.P2iy/(genauigkeit*2)
%Planeten%.P2z:=%Planeten%.z+ %Planeten%.P2iz/(genauigkeit*2)

%Planeten%.xd:=%PlanetenB%.x-%Planeten%.P2x ; siehe oben usw und so fort
%Planeten%.yd:=%PlanetenB%.y-%Planeten%.P2y
%Planeten%.zd:=%PlanetenB%.z-%Planeten%.P2z

abstand:=sqrt(%Planeten%.xd*%Planeten%.xd+%Planete n%.yd*%Planeten%.yd+%Planeten%.zd*%Planeten%.zd)
%Planeten%.xn:=%Planeten%.xd/abstand
%Planeten%.yn:=%Planeten%.yd/abstand
%Planeten%.zn:=%Planeten%.zd/abstand

%Planeten%.s:=%PlanetenB%.gm/(abstand*abstand)

%Planeten%.g3x:=%Planeten%.s*%Planeten%.xn ;%Planeten%.g/(genauigkeit*2 )=m3 teil 1
%Planeten%.g3y:=%Planeten%.s*%Planeten%.yn
%Planeten%.g3z:=%Planeten%.s*%Planeten%.zn

%Planeten%.P3ix := %Planeten%.hx + %Planeten%.g3x/genauigkeit ; %Planeten%.g/(genauigkeit*2 )=m3
%Planeten%.P3iy := %Planeten%.hy + %Planeten%.g3y/genauigkeit
%Planeten%.P3iz := %Planeten%.hz + %Planeten%.g3z/genauigkeit

%Planeten%.P3x:=%Planeten%.x+ %Planeten%.P3ix/genauigkeit ;teil 2
%Planeten%.P3y:=%Planeten%.y+ %Planeten%.P3iy/genauigkeit
%Planeten%.P3z:=%Planeten%.z+ %Planeten%.P3iz/genauigkeit

%Planeten%.xd:=%PlanetenB%.x-%Planeten%.P3x
%Planeten%.yd:=%PlanetenB%.y-%Planeten%.P3y
%Planeten%.zd:=%PlanetenB%.z-%Planeten%.P3z

abstand:=sqrt(%Planeten%.xd*%Planeten%.xd+%Planete n%.yd*%Planeten%.yd+%Planeten%.zd*%Planeten%.zd)
%Planeten%.xn:=%Planeten%.xd/abstand
%Planeten%.yn:=%Planeten%.yd/abstand
%Planeten%.zn:=%Planeten%.zd/abstand

%Planeten%.s:=%PlanetenB%.gm/(abstand*abstand)

%Planeten%.g4x:=%Planeten%.s*%Planeten%.xn ;%Planeten%.g/(genauigkeit*2 )=m4 teil 1
%Planeten%.g4y:=%Planeten%.s*%Planeten%.yn
%Planeten%.g4z:=%Planeten%.s*%Planeten%.zn

%Planeten%.P4ix := %Planeten%.hx + %Planeten%.g4x/genauigkeit ; %Planeten%.g/(genauigkeit*2 )=m4
%Planeten%.P4iy := %Planeten%.hy + %Planeten%.g4y/genauigkeit
%Planeten%.P4iz := %Planeten%.hz + %Planeten%.g4z/genauigkeit

%Planeten%.ix:=1/6*(%Planeten%.P1ix + 2 * %Planeten%.P2ix + 2*%Planeten%.P3ix + %Planeten%.P4ix) Verrechnung der 3 Prognosepunkte zum dann verwendeten Impuls
%Planeten%.iy:=1/6*(%Planeten%.P1iy + 2 * %Planeten%.P2iy + 2*%Planeten%.P3iy + %Planeten%.P4iy)
%Planeten%.iz:=1/6*(%Planeten%.P1iz + 2 * %Planeten%.P2iz + 2*%Planeten%.P3iz + %Planeten%.P4iz)

; ende innere schleife, ausgeführt für jeden gegnerischen Körper

%Planeten%.x:=%Planeten%.x+ %Planeten%.ix/genauigkeit ; Verrechnen des summierten neuen Impuls mit den Koordinaten
%Planeten%.y:=%Planeten%.y+ %Planeten%.iy/genauigkeit
%Planeten%.z:=%Planeten%.z+ %Planeten%.iz/genauigkeit



Ich hoffe es ist einigermaßen verständlich. Falls fragen sein sollten immer her damit.

Bernhard
12.02.2013, 00:25
Falls fragen sein sollten immer her damit.
Da ich im Quelltext keinen offensichtlichen Fehler finde, sollten wir uns besser die Startwerte genauer ansehen. Wie kommst Du auf den Wert -0.324 für die Variable ImpulsY,Mond?
MfG

Kibo
12.02.2013, 08:13
Ich muss gestehen, das war reines Ausprobieren, und zwar im Euler:o
Auch da geben die Werte ja keine perfekt kreisrunde Bahn, aber sie bleibt für viele Umläufe stabil.

Mir ist da noch eine kleine Idee gekommen. Beim RK berechne ich ja 3 Prognosepunkte von denen ich mir dann die Steigungen ausborge, Diese Prognosepunkte werden innerhalb des einen Rechenschritts nur für einen gegnerischen Planeten gerechnet, kann das das Ergebnis negativ beeinflussen? In dem benutzten Szenario, sollte es ja egal sein, da es ja nur 2 Körper im System gibt.

Ich werde mir mal ein ,paar Kontrollmechanismen überlegen, und mir meinen Impulsvektor genauer angucken, wird zeit das der mit reellen Zahlen rechnet. Empfiehlst du da sowas wie die Bahngeschwindigkeit in m/s oder Kraft in Newton?

mfg

mac
12.02.2013, 09:16
Hallo Kibo,

von hier: http://ssd.jpl.nasa.gov/horizons.cgi bekommst Du reale Daten.

Bei Ephemeris Type wählst Du VECTORS
Target Body die jeweiligen Objekte Sonne/Planeten/Monde/ Die mußt Du, soweit ich das bisher verstanden habe, leider alle einzeln abfragen. Und daß ich hier auch Sonne geschrieben habe, war kein Versehen.
Bei Coordinate Origin wählst Du Solar System Barycenter (SSB) (500@0)
Bei Time Span einen im Prinzip beliebigen Start und Stopwert. Grundsätzliche Fehler in Deiner Codierung werden sich aber bereits nach einer einzigen Iteration durch den Datenvergleich zeigen, also z.B. nach einer Sekunde oder einer Stunde.
Bei Table Settings wählst Du output units=KM-S m/s mußt Du daraus selber errechnen. Ich würde dabei, wo immer das möglich ist, die SI-Einheiten verwenden und wenn Du davon abweichen mußt, dann dokumentiere das im Programm incl. der Begründung, warum Du das tust.
Und dann noch Display/Output Da steht zur Auswahl HTML, plain text und download/save. Das mußt Du ausprobieren.

Du bekommst damit, soweit mir bekannt, die derzeit besten Daten für Position und Geschwindigkeit aller bekannten Objekte im Sonnensystem.

Als primären Test würde ich allerdings ein einziges schweres Zentralobjekt und einen einzigen masselosen Körper auf einer Kreisbahn wählen. Die Startwerte dafür kannst Du bekanntlich über die Formel zur 1. kosmischen Geschwindigkeit beliebig genau errechnen. http://de.wikipedia.org/wiki/Kosmische_Geschwindigkeiten

Als nächsten Test gibst Du dem Masselosen Körper eine Masse und bestimmst die Position beider Körper über http://de.wikipedia.org/wiki/Baryzentrum und ihre Geschwindigkeit über die Umlaufzeit um den so ermittelten Schwerpunkt des Systems (den Du dabei als Koordinatenursprung wählst) mit http://de.wikipedia.org/wiki/Satellitenbahn

Herzliche Grüße

MAC

Bernhard
12.02.2013, 11:02
Hallo Kibo,


Empfiehlst du da sowas wie die Bahngeschwindigkeit in m/s oder Kraft in Newton?
ich fände es prima, wenn Du alles auf SI-Einheiten umstellen würdest. Alle sechs Koordinaten also in m. Anstelle des Impulses sollte die Datei Planeten.txt die Startgeschwindigkeiten in m/s enthalten. Die Masse sollte man in kg angeben.

Die Datei Planeten.txt sähe dann beispielsweise so aus:


KoordX KoordY KoordZ v_X v_Y v_Z Masse s v Name dichte s und v sind noch nicht implementiert
-4.67e6 0.0 0.0 0.0 0.01 0.0 5.974e24 0 0 erde 5.515e3
3.794e8 0.0 0.0 0.0 -1.01e3 0.0 7.349e22 0 0 mond 3.341e3
Damit sollte das Programm Euler dann funktionieren. Falls es das nicht macht, müsste man erst mal dort nach Fehlern suchen. Man gibt damit ein System vor, das in etwa durch Schwerpunktkoordinaten vorgegeben wird. Bitte beachte dabei, dass ich in z-Richtung die Startkoordinaten und die Startgeschwindigkeiten gleich Null gesetzt habe. Das vereinfacht die Rechnung zusätzlich.

Auf die Daten der NASA-Horizons Seite können wir in diesem Stadium vorerst noch verzichten. Die Konstellation Erde-Mond ist für den Anfang ideal, weil sie sehr übersichtlich ist.
MfG

mac
12.02.2013, 11:56
Hallo Bernhard,


Auf die Daten der NASA-Horizons Seite können wir in diesem Stadium vorerst noch verzichten. Die Konstellation Erde-Mond ist für den Anfang ideal, weil sie sehr übersichtlich ist.
MfGDas Erde-Mond-System ist in meinen Augen nicht so gut geeignet, weil die überprüfbaren Daten nicht mit den Modelldaten übereinstimmen können. Besser wäre stattdessen der Test, den ich im letzten Absatz des Posts 33 vorgeschlagen hatte, der auch bei 2D noch genauso gut funktioniert und sehr leicht unabhängig überprüfbar ist.

Herzliche Grüße

MAC

Bernhard
12.02.2013, 12:30
Das Erde-Mond-System ist in meinen Augen nicht so gut geeignet, weil die überprüfbaren Daten nicht mit den Modelldaten übereinstimmen können.

Hallo MAC,

im Prinzip schlagen wir das gleiche vor, nämlich das Zweikörper-Problem. Mein Vorschlag für die Datei Planeten.txt setzt nur die mittlere Orbitalgeschwindigkeit des Mondes (s. Wikipedia) voraus und liefert über die Berechnung der Schwerpunktskoordinaten mit Hilfe der Distanz Erde-Mond (384.400 km) die vorgeschlagenen Daten.

Ich bin mir momentan auch nicht sicher, ob der Fehler anstelle im RK-Modul nicht im Euler-Modul liegt. Beide Module sollten aber im Display für einen Umlauf das gleiche Ergebnis liefern.
MfG

Bernhard
12.02.2013, 18:43
Die Datei Planeten.txt sähe dann beispielsweise so aus:


KoordX KoordY KoordZ v_X v_Y v_Z Masse s v Name dichte s und v sind noch nicht implementiert
-4.67e6 0.0 0.0 0.0 0.01 0.0 5.974e24 0 0 erde 5.515e3
3.794e8 0.0 0.0 0.0 -1.01e3 0.0 7.349e22 0 0 mond 3.341e3

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
Es ergibt sich also eine nahezu kreisförmige Bahn und das sollte sowohl die Berechnung nach Euler, als auch mit Runge-Kutta bestätigen.
MfG

Kibo
12.02.2013, 19:34
Vielen Dank für eure Vorschläge,

Meine Zeit für für das Programm ist Werktags arbeitstechnisch relativ knapp bemessen. Ich kann schon mal sagen, dass die Umstellung der Entfernungs- und Massenangaben überhaupt kein Problem darstellt, beim Impuls konnte ich aus Zeitmangel noch nicht viel machen. Werte von genau 0 sind in dem Programm schwierig umsetzbar, die führen aus mir nicht bekannten Gründen im kompilierten Programm zu reihenweise leeren Variablen, gleich nach dem ersten Start, daher vermeide Ich diese vorerst.

mfg

Bernhard
12.02.2013, 21:37
Werte von genau 0 sind in dem Programm schwierig umsetzbar, die führen aus mir nicht bekannten Gründen im kompilierten Programm zu reihenweise leeren Variablen, gleich nach dem ersten Start, daher vermeide Ich diese vorerst.
Das ist wirklich unschön. Gibt es da keinen Workaround? Kann man eventuell einen bestimmten String oder ein bestimmtes Sonderzeichen einlesen und im Programm die zugehörige Variable dann auf Null setzen?

Bernhard
13.02.2013, 00:35
Die Datei Planeten.txt sähe dann beispielsweise so aus:


KoordX KoordY KoordZ v_X v_Y v_Z Masse s v Name dichte s und v sind noch nicht implementiert
-4.67e6 0.0 0.0 0.0 0.01 0.0 5.974e24 0 0 erde 5.515e3
3.794e8 0.0 0.0 0.0 -1.01e3 0.0 7.349e22 0 0 mond 3.341e3

Hallo Kibo,

Bei v_Y_Erde habe ich mich verrechnet. Die neuen Zahlen für eine Kreisbahn lauten:


KoordX KoordY KoordZ v_X v_Y v_Z Masse s v Name dichte s und v sind noch nicht implementiert
-4.67e6 0.0 0.0 0.0 12.45 0.0 5.974e24 0 0 erde 5.515e3
3.797e8 0.0 0.0 0.0 -1.012e3 0.0 7.349e22 0 0 mond 3.341e3
Der Abstand zwischen Erde und Mond beträgt dann konstant 384.400 km. Die Ableitung des Differenzvektors der beiden Ortsvektoren (=v) habe ich mit sqrt( G * (m_1 + m_2) / r ) berechnet.

Bernhard
13.02.2013, 11:00
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
18.02.2013, 23:43
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 (http://sdrv.ms/Xl1UYc)

Rapidshare (http://rapidshare.com/files/2054815486/planeten%201.4.exe)

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.

Bernhard
19.02.2013, 21:03
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.



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
19.02.2013, 22:31
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
20.02.2013, 04:05
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

Bernhard
21.02.2013, 19:26
H_T = p1^2/(2*m1) + p2^2 / (2 * m2)
Der Kenner erkennt hier natürlich die Hamiltonfunktion zweier freier Massepunkte und löst das resultierende Gleichungssystem exakt, was man dann auch als Grund dafür hernehmen darf, warum der gesamte Integrator recht schnell arbeitet.

Kibo
24.02.2013, 19:26
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 (http://my.rapidshare.com/Kibo186/13327)

Skydrive (http://sdrv.ms/ZsLvTW)

mfg

Kibo
07.03.2013, 19:00
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

Bernhard
07.03.2013, 20:19
Hallo Kibo,

wie kann man das Erzeugen der Datei log.txt eigentlich verhindern?

Kibo
07.03.2013, 20:43
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?

Bernhard
07.03.2013, 21:04
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
07.03.2013, 21:13
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)

Bernhard
07.03.2013, 22:58
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:


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

Kibo
08.03.2013, 08:14
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
08.03.2013, 13:51
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

Bernhard
08.03.2013, 14:03
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 (http://iopscience.iop.org/1538-3881/116/4/2067/) 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
08.03.2013, 21:35
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
09.03.2013, 00:24
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%E2%80%93Kutta_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
10.03.2013, 23:08
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

Bernhard
11.03.2013, 08:27
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.

Bernhard
11.03.2013, 09:30
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)

Muss demnach zu


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


korrigiert werden.

Kibo
11.03.2013, 21:59
Wenn ich deinen Vektor richtig verstanden habe, dann ist bei dir

Die Position des Prognosepunkts P1 von Sonne1 (-7.48e10/ 1.0525e8/ 0.0)
und die "Steigung" m1 dann (2.9654e1, 2.105e4, 0.0)

Die von meinem Programm errechneten Werte sind ja
für P1 (-74799998516,438034.../10525000/0)
und für m1 (2,967124.../21050/0)

Es unterscheiden sich also der x-Wert von P1 und m1. Sehr wundert mich ersterer, da du selber schreibst


Die Beschleunigung als Vektor zeigt am Anfang also in Richtung der positiven x-Achse.


Also warum prognostiziert dein P1 diese Änderung nicht?

Wie wir übereinstimmend feststellten ist die durch die Gravitation verursachte Beschleunigung (rund) 5.93e-3 m/s^2, diese wirkt nur in Richtung x Achse. Da wir das für die halbe gewählte Schrittweite (500 s) betrachten kommen wir beide auf diese nahezu gleichen Werte, warum sie allerdings nicht 100%tig gleich sind, wüsste Ich schon gerne. Eine Kontrollrechnung in einer Exceltabelle, bestätigt die Werte meines Programms.

Nichtsdestotrotz muss mein Runge Kutta falsch sein, da bei einem Probelauf sich die Orbits der Sonnen schon nach kurzer Zeit immer weiter absenken sodass es zu Kollision kommt (Energieerhaltung???).
Euler hingegen, liefert selbst bei 800000 Sekunden Schrittweite noch nahezu perfekt kreisförmige Orbits, nach simulierten 1000 Jahren ist keine Destabilisierung erkennbar.

mfg

Bernhard
11.03.2013, 23:01
Also warum prognostiziert dein P1 diese Änderung nicht?
Die Gravitationskraft ist eine Beschleunigung und diese ändert ihrerseits im ersten Schritt nur die Geschwindkeit. Der Prognosepunkt enthält eben auch die Geschwindigkeit des Massepunktes 1 und die x-Komponente dieser Geschwindigkeit wird dann ja auch auf den Wert 29.7 m/s gesetzt.


warum sie allerdings nicht 100%tig gleich sind, wüsste Ich schon gerne.
Ich habe 10000 (zehntausend) Sekunden als Schrittweite verwendet. Wenn ich 1000 Sekunden verwende bekomme ich genau die gleichen Werte. Der Fehler liegt vermutlich wo anders. Deswegen rechne ich einfach mal weiter:

k2 = (2.9654e1, 2.105e4, 0.0, 5.9308e-3, -8.3452e-6, 0.0, -2.9654e1, -2.105e4, 0.0, -5.9308e-3, 8.3452e-6, 0.0)

sowie

y_n + 0.5 * h * k2 = (-7.48e10, 1.0525e8, 0.0, 2.9654e1, 2.105e4, 0.0, 7.48e10, -1.0525e8, 0.0, -2.9654e1, -2.105e4, 0.0)

und damit k3 gleich k2 in den hier aufgeschriebenen Dezimalstellen.

Kibo
12.03.2013, 07:24
Hmm, als ich den Runge Kutta implementiert habe, hab ich mich vornehmlich an diese Seite gehalten.

http://kohorst-lemgo.de/modell/runge.htm

Da die dortige Rechnung eine Änderung der Koordinaten impliziert, ist das dann dort falsch dargestellt?

mfg

Bernhard
12.03.2013, 09:26
Hallo Kibo,


Da die dortige Rechnung eine Änderung der Koordinaten impliziert, ist das dann dort falsch dargestellt?

das ist, so weit ich es verstehe, ein Mißverständnis. In der englischen Wikipedia wird anstelle von m1 bis m4 nur die Bezeichnung k1 bis k4 verwendet. Der Algorithmus ist ansonsten der selbe. Die Koordinaten sind immer euklidisch. Man muss nur beachten, dass man bei einem Zweikörperproblem die zwei Orte x_1 und x_2 und zusätzlich noch die zwei Geschwindigkeiten v_1 und v_2 praktisch als unabhängige Variablen betrachten muss. x_1, x_2, v_1 und v_2 sind dabei die bekannten, dreidimensionalen Vektoren im euklidischen Raum, die den Bewegungszustand der beiden Körper charakterisieren.

m1 bis m4 hat also jeweils zwölf Komponenten, bzw. Einträge. Bei den Prognosepunkten werden ebenfalls für alle vier Vektoren neue Werte berechnet und zwar wieder als zwölf reelle Zahlen (4 * 3 = 12).
MfG

Kibo
14.03.2013, 22:49
Nach einigem Kopfzerbrechen, gibt's ein Paar Fortschritte. Die Bahnen sind jetzt auch im RK Kreisrund wie sie sein sollten. Es existiert abe rnoch mindestens 1 weiterer Programmierfehler, da das System immernoch mit der Zeit seine Bahnenergie verliert (Todesspirale).

Bernhard
15.03.2013, 03:16
Ich vermute, dass immer noch ein Fehler in der Verwaltung, bzw. softwareseitigen Trennung der verschiedenen Variablen vorliegt. Bedenke, dass die Gravitationskraft nur die aktuellen Geschwindigkeiten ändert. Die Positionen (Orte) werden dagegen nur über die Geschwindigkeiten beeinflusst.

Kibo
16.03.2013, 16:44
Die 4 Sonnen, kreisen auch wieder wie gewohnt, es gab Inkompatibilitäten mit meiner Originalschleife und newton nr3. 1.7 rückt in greifbare nähe, jetzt muss ich noch den Spiralenfehler beseitigen.

mfg

Bernhard
17.03.2013, 14:56
Hallo Kibo,

anbei die letzten Zahlen für den ersten Schritt:

y + h * k_3 = (-7.48e10, 2.105e8, 0.0, 5.9308e1, 2.105e4, 0.0, 7.48e10, -2.105e8, 0.0, -5.9308e1, -2.105e4, 0.0)

k_4 = (5.9308e1, 2.105e4, 0.0, 5.9308e-3, -1.669e-5, 0.0, -5.9308e1, -2.105e4, 0.0, -5.9308e-3, 1.669e-5, 0.0)

Der neue Gesamtvektor (Ort + Geschwindigkeit für beide Sonnen) ergibt sich dann zu

y_1 = (-7.48e10, 2.105e4, 0.0, 5.9308e-3, 2.105e4, 0.0, 7.48e10, -2.105e4, 0.0, -5.9308e-3, -2.105e4, 0.0)

bei einer Schrittweite von 10.000 Sekunden bekommen die beiden y-Koordinaten eine Korrektur und die beiden Geschwindigkeitsvektoren bekommen bei der x-Komponente eine Korrektur. Die vier Korrekturen bestehen dabei aus zwei Werten, die mit beiden Vorzeichen eingehen.

Kibo
17.03.2013, 21:06
Also Bernhard, wenn ich das jetzt richtig verstanden habe, wird in einem Schritt bei dir folgendes gemacht:

Die Koordinaten der Körper werden aufgrund der Geschwindigkeit angepasst

DANACH

werden die Beschleunigungen ausgerechnet und die Geschwindigkeiten angepasst.

Sehe ich das so richtig? Wenn ja, warum nicht anders herum?

Bernhard
18.03.2013, 00:14
Sehe ich das so richtig?
So in etwa. Ich erkläre es aber gerne auch etwas genauer:

Die Gleichung F = m * a reicht auch zusammen mit dem newtonschen Gravitationsgesetz nicht aus, um Runge-Kutta zu implementieren, denn Runge-Kutta kann nur erste Ableitungen nähern. Da die Beschleunigung nun aber ungünstigerweise die zweite Ableitung des Ortes ist kann man nicht ausschließlich mit der Beschleunigung arbeiten.

Man verwendet vielmehr den Umweg über die Geschwindigkeit und das wie folgt: Die erste Ableitung des Ortes ist die Geschwindigkeit. Ferner gilt noch: Die erste Ableitung der Geschwindigkeit ist die Beschleunigung und die kann/soll/muss man aus dem newtonschen Gravitationsgesetz ausrechnen.

Dann denkt man sich den Ort und die Geschwindigkeit als eine große Variable (mit mehreren Komponenten) und die erste Ableitung dieser Variable ist jetzt bekannt, wie im zweiten Absatz eben gezeigt. Man muss dann nur noch zwischen der Ableitung und dem Wert der großen Variablen peinlich genau unterscheiden und ist vom Verständnis her fertig.

Zur Wiederholung:
a) In der zugrundeliegenden Differentialgleichung taucht die Geschwindigkeit einmal als Wert auf und gibt dort die erste Ableitung des Ortes an und
b) Die erste Ableitung der Geschwindigkeit wird gemäß F/m berechnet, wobei F die newtonsche Gravitationskraft auf den betreffenden Körper ist.
MfG

Kibo
18.03.2013, 20:56
Hmm,

also ich fass meine jetzige Lösung vielleicht am besten nochmal in Pseudocode zusammen:

Schleife bis alle Planeten durch sind
{
sichern des momentanen Geschwindigkeitsvektors in PlanetA.h & PlanetA.h

m1a[x,y,z] = gravitation(PlanetA koordinaten, PlanetB)
m1v = h + m1a*schrittweite Kommentar: Das entspricht nicht deiner Beschreibung, aber wie anders?
P1 = PlanetA koordinaten + m1v* Schrittweite*0,5
m2a[x,y,z] = gravitation(P1, PlanetB)
m2v = h + m2a*schrittweite
P2 = PlanetA koordinaten + m2v* Schrittweite*0,5
m3a[x,y,z] = gravitation(P2, PlanetB)
m3v = h + m3a*schrittweite
P3 = PlanetA koordinaten + m3v* Schrittweite
m3a[x,y,z] = gravitation(P3, PlanetB)
m3v = h + m3a*schrittweite

PlanetA.v:=( m1v + 2 * m2v + 2 * m3v + m4v)/6

PlanetA.a:= PlanetA.v - h Kommentar: Rückrechnen für Werte für PlanetB (Newton Nr3)
PlanetB.v:= -1 * PlanetA.a * PlanetA.m /PlanetB.m

}

gravitation(PlanetA koordinaten, PlanetB)
{
differenzvektor = koordninaten PlanetA - PlanetB
abstand = wurzel(Differenzvektor ²)
einheitsvektor für richtung n = differenzverktor/ abstand

Beschleunigung a = PlanetB.gm/abstand² Kommentar: PlanetB.gm = G*Masse bei Programmstart berechnet
grav Beschleunigungsvektor g = a * n
rückgabe g, abstand
}

Schleife für jeden Planeten
{
Planetkoordinaten:=Planetenkoordinaten + Planet.v*Schrittweite
}

Führt wie gesagt zu einer "Todesspirale"

MfG

Bernhard
18.03.2013, 23:56
sichern des momentanen Geschwindigkeitsvektors in PlanetA.h & PlanetA.h
Das ist ungünstig weil verwirrend. Warum nicht PlanetA.v für die Geschwindigkeit von PlanetA und PlanetA.x für den aktuellen Ort von PlanetA? Für die Geschwindigkeit nimmt man normalerweise immer das kleine v wegen englisch velocity.



m1a[x,y,z] = gravitation(PlanetA koordinaten, PlanetB)

muss heißen
m1a[x,y,z] = gravitation(PlanetB koordinaten, PlanetA)

Planet B zieht Planeten A in Richtung hin zu B. Die Kraft auf Planet A muss also von A nach B zeigen. Da haben wir also schon mal einen Vorzeichenfehler.



m1v = h + m1a*schrittweite Kommentar: Das entspricht nicht deiner Beschreibung, aber wie anders?

muss heißen:

m1v = aktuelle Geschwindigkeit des Planeten, also beispielsweise

m1v = PlanetA.v

Die Schrittweite h wird ferner nirgends irgendwo dazu addiert, weil die Zeit nirgends explizit vorkommt. h tritt dagegen überall nur multiplikativ auf. Jedes "h + ..." ist also falsch und muss korrigiert werden.

EDIT: Bei jedem Schleifendurchlauf berechnet man aus den aktuellen Positionen und Geschwindigkeiten der Planeten und der Schrittweite die neuen Positionen und Geschwindigkeiten der Planeten. Die Zeit kann man völlig seperat berechnen oder vorher festlegen. Man kann beispielsweise ein t0, ein t_stop und eine Anzahl von Schritten N vorgeben. Dann berechnet man h = (t_stop - t0) / N. Die Werte von t0 und t_stop gehen in die Rechnung dann gar nicht mit ein. Man muss nur die Schleife genau N mal durchlaufen und darin das eben berechnete h als Schrittweite verwenden.
MfG

Kibo
19.03.2013, 07:37
Guten Morgen Bernhard,

h ist die Hilfsvariable und speichert nur das alte v für die Trägheit. Die Schrittweite hab ich hier konkret mit schrittweite bezeichnet. Der letzte post von mir ist Pseudocode, und sollte eigentlich nur meine prinzipielle Vorgehensweise verständlich darstellen, ohne auf die Eigenheiten der verwendeten Programmiersprache einzugehen.
Da meine Planeten alle von einander angezogen und nicht abgestoßen werden, geh ich mal davon aus das der von dir gefundene Vorzeichenfehler auf meiner Unfähigkeit beruht meinen Code verständlich auszudrücken.:o

OK, du hast erst mal nur meine unverständlichen Variablennamen bemängelt, du siehst da also keinen logischen Fehler im Code an sich?

Ich lade mal das aktuell kompilierte Programm hoch (1.7alpha zu den Verbesserungen später mehr), damit du dir mal den Fehler ansehen kannst, bei 80000 Schrittweite sieht man es ganz gut. Sieht so ein bisschen nach Rundungsfehler aus, oder kann das hier ran liegen?:



Die Gleichung F = m * a reicht auch zusammen mit dem newtonschen Gravitationsgesetz nicht aus, um Runge-Kutta zu implementieren, denn Runge-Kutta kann nur erste Ableitungen nähern. Da die Beschleunigung nun aber ungünstigerweise die zweite Ableitung des Ortes ist kann man nicht ausschließlich mit der Beschleunigung arbeiten.



hier nochmal die Links:

Rapidshare (http://my.rapidshare.com/Kibo186/13327)

Skydrive (http://sdrv.ms/ZsLvTW)

mfg

Bernhard
19.03.2013, 09:10
Guten Morgen Kibo,

sinnvolle Bezeichnungen für Variablennamen sind extrem wichtig, damit man den Code auch nach längeren Pausen noch lesen kann. Ich schlage deswegen vor, wir suchen erst mal gemeinsam ein paar gute Bezeichnugen, um sinnvoll diskutieren zu können.

Das m1a und m1v gefällt mir schon mal gut. Es sollte für jeden Planeten exisitieren. Genauso m2a, m2v, m3a, m3v, m4a, m4v. Ferner benötigst Du für jeden Planeten einen Prognosepunkt und eine Prognosegeschwindigkeit.

Zusätzlich brauchst Du eine Funktion, welche die Differentialgleichung repräsentiert. SRMeister hatte diese Funktion in seinem Code "Evaluate" genannt. Diese Funktion erhält als Parameter die Schrittweite, den Prognosepunkt und die Prognosegeschwindigkeit und berechnet daraus die Werte für mia und miv, mit i = 1,2,3,4. Bei dieser Funktion helfe ich Dir dann gerne, wenn wir uns bei den benötigten Variablen geeinigt haben.
MfG

Kibo
19.03.2013, 10:28
Ok, ich fasse nochmal meine jetzigen Variablen zusammen:

%Planeten%.h(x,y oder z) hilfsvariable, enthält
%Planeten%.x, y oder z Koordinaten
%Planeten%.v(x,y oder z) Geschwindigkeit
%Planeten%.g(x,y oder z) Beschleunigung

Prognosevariablen
%Planeten%.m(1, 2, 3 oder 4)v(x,y oder z)
%Planeten%.m(1, 2, 3 oder 4)a(x,y oder z)
%Planeten%.P(1, 2 oder 3)(x,y oder z)

in der Funktion (auch alternativ als Klasse probiert) zur Gravitationsberechnung
%Planeten%.(x, y oder z)d
%Planeten%..(x, y oder z)n
%Planeten%.abstand
%Planeten%.a (ohne richtung)

Das müsste alles sein, bin gerade auf Arbeit, da ist das etwas schwierig:)

%Planeten% ist eine dynamische Variable und ändert sich innerhalb der schleife entsprechend in den Namen des Aktuellen Planeten. %PlanetenB% macht das selbe für den Gegnerplaneten. Wie erwähnt nach Newton Nr 3 die Beschleunigung für %PlanetenB% gleich mit gesetzt. Es werden in einem Durchlauf auch alle Planeten genau einmal berechnet, das hab ich kontrolliert um Fehler in der Schleife auszuschließen.

mfg

Bernhard
19.03.2013, 10:46
Da fehlt noch die Prognosevariable für die Geschwindigkeit der Planeten.

Kibo
19.03.2013, 11:11
Prognosevariablen
%Planeten%.m(1, 2, 3 oder 4)v(x,y oder z)
Ist das nicht die?

Oder schlägst du vor das RK umzubauen? zum Beispiel:

einmal RK zum integrieren Beschleunigung>Geschwindigkeit und nachträglich noch einmal RK für integrieren Geschwindigkeit>Koordinaten
Oder eins von beiden eulern?
oder ganz anders?

mfg

Bernhard
19.03.2013, 11:37
Ist das nicht die?
Ich würde die m-Variablen eher Steigungen nennen. Der Prognosepunkt und die Prognosegeschwindigkeit berechnen sich dann gemäß:

Px = x + Schrittweite / 2 * m1v
Pv = v + Schrittweite / 2 * m1a

usw., denn es gibt dann zu jedem m ein Prognosepaar.


Oder schlägst du vor das RK umzubauen?
Das Runge-Kutta-Verfahren ist fix. Es muss nur noch implementiert werden und zwar am besten fehlerfrei ;) .


Oder eins von beiden eulern?
Pfui.

Kibo
19.03.2013, 12:32
Gut,
näheres dazu dann nach 18 Uhr wenn ich zu Hause bin.

mfg

Kibo
19.03.2013, 21:42
Kann es sein, dass bei k1 - 4 sich Die Position der Planeten nicht ändert?
Das Internet gibt ja so einiges her zum RK4 blos das ist ja schon merkwürdig, wenn ich da jetzt einen Denkfehler (seit über einem Monat ) habe, korrigiere mich bitte:

Ich ging/gehe davon aus, dass das berechnen der Prognosekoordinaten auch irgendwas bringen muss, also rechne ich für P1
Gravitation aus = m1a
Daraus dann m1v = altes v +m1a*schrittweite
Daraus dann Prognosekoordinate P1 = altes x + m1v*schrittweite/2

dann nehme ich dieses P1 (koordinate) und bestimme die gravitation, die der Körper erfahren würde wenn er an dieser stelle steht.
Daraus dann wieder m2v = altes v + m2a*schrittweite
Daraus dann Prognosekoordinate P2 = altes x + m2v*schrittweite/2

Dann bestimme ich wieder die Gravitation an dieser neuen Koordinate und so weiter und so fort.
Da auf Grund dieser linearen Abhängigkeit ich die Koordinaten ja aus der Geschwindigkeit und diese aus der Beschleunigung berechnen muss kann ich die ja auch schlecht unabhängig von einander betrachten :confused:

Bernhard
20.03.2013, 01:29
Ich ging/gehe davon aus, dass das berechnen der Prognosekoordinaten auch irgendwas bringen muss, also rechne ich für P1
Gravitation aus = m1a
Daraus dann m1v = altes v +m1a*schrittweite
Daraus dann Prognosekoordinate P1 = altes x + m1a*schrittweite/2
Nein. Das wäre z.T. ja das Euler-Verfahren. Es gilt vielmehr:
m1v = altes v
und Prognosekoordinate
P1 = altes x + m1v*schrittweite/2.
Das ist der geschickt modifizierte Euler.
Und zusätzlich brauchst Du hier eben auch die Prognosegeschwindigkeit:
P1v = altes v + m1a * schrittweite/2.

Schaue Dir doch mal die Einheiten an! m1a ist eine Beschleunigung in m/s². Wenn man eine Beschleunigung mit der Schrittweite (Einheit Sekunden) multipliziert bekommt man m/s. Also eine Geschwindigkeit. Schon deswegen kann Deine Formel



P1 = altes x + m1a*schrittweite/2

nicht stimmen.



dann nehme ich dieses P1 (koordinate) und bestimme die gravitation, die der Körper erfahren würde wenn er an dieser stelle steht.

Prinzipiell richtig. Vergiss aber auch hier nicht die Prognosegeschwindigkeit.



Daraus dann wieder m2v = altes v + m2a*schrittweite
Daraus dann Prognosekoordinate P2 = altes x + m2a*schrittweite/2

Auch hier wieder der gleiche Fehler, den man schon aufgrund der physikalischen Einheiten aufspüren kann. Es gilt vielmehr:
m2v = P1v
m2a = Gravitation, bzw. Beschleunigung an der Stelle des Prognosepunktes P1.

Kibo
20.03.2013, 08:01
Ja, das war ein blöder Kopierfehler, aufgrund der späten Uhrzeit

Bernhard
20.03.2013, 08:48
Weißt Du jetzt, was die verbleibenden Fehler sind?

Kibo
20.03.2013, 09:00
Morgen Bernhard,

Ich habe deine Änderungen vorgenommen. Leider sind die Spiralen immernoch da, hat sich nicht verbessert.:( Hier mal der gesamte RK:



m1a:=Grav(%Planeten%.x,%Planeten%.y,%Planeten%.z,% PlanetenB%)

%Planeten%.P1v := %Planeten%.h + m1a *Schrittweite (ohne /2 das sonst Bahnen elliptisch werden)

%Planeten%.m1v := %Planeten%.h

%Planeten%.P1:=%Planeten%.x+ %Planeten%.m1v*Schrittweite/2

m2a:=Grav(%Planeten%.P1x,%Planeten%.P1y,%Planeten% .P1z,%PlanetenB%)

%Planeten%.P2v := %Planeten%.h + m2a *Schrittweite

%Planeten%.m2v:=%Planeten%.P1v

%Planeten%.P2:=%Planeten%.x+ %Planeten%.m2v *Schrittweite/2

m3a:=Grav(%Planeten%.P2x,%Planeten%.P2y,%Planeten% .P2z,%PlanetenB%)

%Planeten%.P3v := %Planeten%.x + m3a *Schrittweite

%Planeten%.m3v:=%Planeten%.P2v


%Planeten%.P3:=%Planeten%.x + %Planeten%.m3v *Schrittweite

m4a:=Grav(%Planeten%.P3x,%Planeten%.P3y,%Planeten% .P3z,%PlanetenB%)

%Planeten%.P4v := %Planeten%.h + m4a*Schrittweite

%Planeten%.v:=(%Planeten%.P1v + 2 * %Planeten%.P2v + 2*%Planeten%.P3v + %Planeten%.P4v)/6


mfg

Bernhard
20.03.2013, 09:25
OK. Das muss ich mir in Ruhe ansehen. Voraussichtlich erst heute abend.

UMa
20.03.2013, 09:39
Hallo Kibo,

so kann das nicht funktionieren. Wo wird z.B. m2a berechnet? Es sind viel zu viele Fehler drin.
Ich habe mal versucht, deine Bezeichnung ungefähr zu übernehmen. Ort steht für den Vektor der Positionen mit 3*n Komponenten. v für den Geschwindigkeitsvektor mit 3*n Komponenten, wobei n die Anzahl der Planeten ist. Wenn du die Planeten getrennt speicherst, musst du jedesmal wenn P1ort= ... P1v=... m1a=... usw. (also überall, außer bei dt=Schrittweite) steht eine Schleife über alle Planeten machen. Du kannst die Planeten nicht nacheinander berechnen.

Grüße UMa

Aus ort_alt und v_alt mache ort_neu und v_neu.
Grav() berechnet die Gravitationsbeschleunigung, abhängig von Ort und Masse.
Klassisches RKV:

dt=Schrittweite
P1ort=ort_alt
P1v =v_alt
m1a=Grav(P1ort)
P2ort=ort_alt + dt/2*P1v
P2v =v_alt + dt/2*m1a
m2a=Grav(P2ort)
P3ort=ort_alt + dt/2*P2v
P3v =v_alt + dt/2*m2a
m3a=Grav(P3ort)
P4ort=ort_alt + dt*P3v
P4v =v_alt + dt*m3a
m4a=Grav(P4ort)
ort_neu=ort_alt + dt/6*(P1v+2*P2v+2*P3v+P4v)
v_neu =v_alt + dt/6*(m1a+2*m2a+2*m3a+m4a)

Kibo
20.03.2013, 11:08
Hallo UMa,

ich hab die Zeile nachträglich eingefügt. Besten dank für deine Beteiligung


Du kannst die Planeten nicht nacheinander berechnen.


Du schlägst also vor nach der Berechnung von p2 gleich die Berechnung von p2 für alle anderen Planeten mit zu machen?

Auf deine weiteren Hinweise kann Ich erst heute Abend eingehen.



mfg

Bernhard
20.03.2013, 11:10
Du schlägst also vor nach der Berechnung von p1 gleich die Berechnung von p1 für alle anderen Planeten mit zu machen?

Ja. Genau so muss es gemacht werden :) .

Kibo
20.03.2013, 11:46
Das klingt einleuchtend, würde auch erklären warum der Euler funktioniert und Heun und RK4 nicht. Da setze ich mich heute Abend gleich ran, spannend:)

mgf

UMa
20.03.2013, 13:13
Hallo Kibo,

Du schlägst also vor nach der Berechnung von p2 gleich die Berechnung von p2 für alle anderen Planeten mit zu machen?
das ist für alle Verfahren, außer dem Eulerverfahren, zwingend erforderlich. Sonst berechnest du in der Zeile m2a=Grav(P2ort) die Gravitation nicht oder falsch. Du brauchst dafür die P2ort Positionen aller Planeten.

Grüße UMa

Bernhard
24.03.2013, 10:00
Hallo Kibo,

da die Arbeit an dem Programm scheinbar nur langsam vorangeht, hätte ich noch zwei Tips zum Lesen:

http://de.wikipedia.org/wiki/Raumkurve
mit
http://de.wikipedia.org/wiki/Frenetsche_Formeln

Ich weiß nicht, wie geläufig Dir diese Dinge sind, aber sie sollten als Hintergrundwissen eigentlich ganz hilfreich sein.
MfG

Kibo
24.03.2013, 10:53
Hallo Bernhard,

Danke schonmal. Ich schreib gerade große Teile um und probiere da auch andere Sachen aus. Das ändert aber nicht wirklich was an Genauigkeit oder Features sondern ist nur etwas sauberer und aufgeräumter programmiert. Da ich sowieso die ganze Schleife umstellen muss, bietet sich das ja an. Ich melde mich wenn ich fertig bin. Momentan sehe ich keine Baustelle die ich nicht selber schließen kann, sollte sich das ändern melde ich mich schon.;)

mfg

Kibo
26.03.2013, 22:19
Nabend, kleines Update aber keine guten Neuigkeiten.

Auch das Berechnen von P1v für alle Planeten vor P2v, sowie P2v vor P3v hat keine Verbesserung gebracht. Die Spirale ist nicht ein bisschen schmaler geworden. Sehr schade, ich habe wirklich gedacht das wäre der Fehler. Ich werde dann wohl mal ein paar Schritte komplett per Hand nachrechnen, vielleicht findet sich so eine Erklärung.:( Aber nicht mehr heute

gute Nacht

Kibo
29.03.2013, 16:53
Hmm,

Wenn ich UMa's Code anwende dann ist v um den Faktor 4 zu schwach.

Gleiche ich das aus mit:

ort_neu=ort_alt + dt/6*(P1v+2*P2v+2*P3v+P4v) * 4

Dann VERGRÖßERT sich der Abstand der Planeten zueinander nach jedem Orbit um einige %...


PS: Heun macht das selbe, nur ist der Effekt wesentlich schwächer

Bernhard
29.03.2013, 17:46
dt=Schrittweite
P1ort=ort_alt
P1v =v_alt
m1a=Grav(P1ort)
P2ort=ort_alt + dt/2*P1v
P2v =v_alt + dt/2*m1a
m2a=Grav(P2ort)
P3ort=ort_alt + dt/2*P2v
P3v =v_alt + dt/2*m2a
m3a=Grav(P3ort)
P4ort=ort_alt + dt*P3v
P4v =v_alt + dt*m3a
m4a=Grav(P4ort)
ort_neu=ort_alt + dt/6*(P1v+2*P2v+2*P3v+P4v)
v_neu =v_alt + dt/6*(m1a+2*m2a+2*m3a+m4a)
Hallo Kibo,

würdest Du bitte für jede Zeile den konkreten Wert in Deinem Programm bei dem ersten Schritt angeben. Also in der Art

dt = ...
P1ort = (..., ..., ...)
P1v = (..., ..., ...)
usw. bis
v_neu = (..., ..., ...)

dann kann ich nachprüfen, ob alles passt.
MfG

Kibo
29.03.2013, 19:14
dt= 1000
P1ort=(-7.48e10/0.0/0.0)
P1v =(-11.868006/-84200.040767/0.000000)
m1a=(-17.802000/0.003000/-0.000000)
P2ort=(74800000000.000000/-10525000.000000/0.000000)
P2v =(-17.802000/-105249.998000/0.000000)
m2a= (-0.029670/0.000006/-0.000000)
P3ort= (74799998516.500000/-10525000.000000/0.000000)
P3v = (-11.868000/-63149.998000/0.000000)
m3a= (-0.017802/0.000004/-0.000000)
P4ort= (74799997033.000000/-21049999.500000/0.000000)
P4v = (-5.934000/-21049.999000/0.000000)
m4a= (-0.005934/0.000002/-0.000000)
ort_neu= (74799988131.994064/-84200040.766666/0.000000)
v_neu = (-11.868006/-84200.040767/0.000000)

Nett dass du mal rüber schaust, danke

Bernhard
29.03.2013, 19:46
P1v =(-11.868006/-84200.040767/0.000000)
Sieht so aus, als wären einfach die Startbedingungen falsch:
P1v =(0.0/21060.0/0.0) für Sonne1

Bernhard
29.03.2013, 19:52
Der neue Gesamtvektor (Ort + Geschwindigkeit für beide Sonnen) ergibt sich dann zu

y_1 = (-7.48e10, 2.105e4, 0.0, 5.9308e-3, 2.105e4, 0.0, 7.48e10, -2.105e4, 0.0, -5.9308e-3, -2.105e4, 0.0)

bei einer Schrittweite von 10.000 Sekunden
Der Vollständigkeit halber und zur Kontrolle:

Ich habe das jetzt mal mit "meinem" Programm gerechnet und komme damit auf

y_1 = (-7.48e10, 2.106e8, 0.0, 5.94e1, 2.106e4, 0.0, 7.48e10, -2.106e8, 0.0, -5.94e1, -2.106e4, 0.0)

Bei einigen Komponenten habe ich die Multiplikation mit der Schrittweite übersehen.

Kibo
30.03.2013, 08:17
Guten Morgen,

Sorry Bernhard, eine der arithmetischen Vektorfunktionen hat anscheinend die v_alt überschrieben kurz bevor ich am ende des Schrittes alle Variablen abgefragt habe. Nach Korrektur des Fehler in der Funktion hat sich aber auch nichts an dem Verhaltend er Planeten gebessert. Dafür sollten die Ausgabewerte jetzt hoffentlich mal stimmen:


dt= 10000
P1ort=/-7.48e10/0.0/0.0
P1v =/0.0/-2.105e4/0.0
m1a=/-0.005934/-0.000000/-0.000000
P2ort=/74800000000.000000/-105250000.000000/0.000000
P2v =/-29.670000/-21050.000000/0.000000
m2a= /-0.005934/0.000008/-0.000000
P3ort= /74799851650.000000/-105250000.000000/0.000000
P3v = /-29.670000/-21049.960000/0.000000
m3a= /-0.005934/0.000008/-0.000000
P4ort= /74799703300.000000/-210499600.000000/0.000000
P4v = /-59.340000/-21049.920000/0.000000
m4a= /-0.005934/0.000017/-0.000000
ort_neu= 74798813199.406601/-841999354.332800/0.000000
v_neu = /-118.680059/-84199.935433/0.000000

Zu beachten ist, dass ich diesmal mit einer Schrittweite von 10000 gerechnet habe.

mfg

Bernhard
30.03.2013, 09:42
Guten Morgen Kibo,

die Hilfswerte stimmen schon mal alle :) . Man kann jetzt auch mit dem Taschenrechner nachprüfen, dass im letzten Schritt der y-Wert (zweiter Wert) falsch addiert wird und zwar sowohl bei ort_neu als auch v_neu. Bei v_neu ist zusätzlich der x-Wert falsch. Der x-Wert bei ort_neu ist dagegen korrekt.

Verwende bei diesem Schritt die korrekte Vektoraddition. Ganz offensichtlich wird in Deinem Programm im letzten Schritt die y-Komponente anders berechnet, als die x-Komponente.
MfG

Kibo
30.03.2013, 14:05
Das ist allerdings erstaunlich. Guck ich mir dann nach der Arbeit mal an. Danke

mfg

Kibo
01.04.2013, 17:23
Mal wieder gute Neuigkeiten:

Das RK läuft jetzt relativ rund. Grund für die Schwierigkeiten war v_neu =v_alt + dt/6*(m1a+2*m2a+2*m3a+m4a)
v_alt war einfach das alte v. Dieses wurde anscheinend irgendwo überschrieben, so dass die Gravitation immer ein ganzes Stück größer berechnet war als richtig gewesen wäre.

Was noch nicht funktioniert ist das Testsystem mit 4 Sonnen im RK, und Heun im generellen, außerdem will ich noch ein paar Features implementieren. Wenn alles gut läuft lade ich heute noch die fertige Version 1.7 hoch.

mfg

Bernhard
01.04.2013, 23:40
Was noch nicht funktioniert ist das Testsystem mit 4 Sonnen im RK
Man muss hier beachten, dass jede Sonne jeweils von den drei anderen Sonnen angezogen wird. Für die Berechnung der Hilfsgrößen m1a, m2a, m3a und m4a empfiehlt sich deswegen eine eigene Funktion, die alle Orte aller Körper kennt.

Bernhard
02.04.2013, 18:34
Hallo Kibo,

falls Du gleich weiter machen willst, hier die Ergebnisse für die vier Sonnen nach 10.000 Sekunden gemäß Runge-Kutta:

Sonne1
X = -82249763057.203293 Y = -963036708.834408 Z = 0.000000
v_x = 6043.466504 v_y = -96127.188943 v_z = 0.000000

Sonne2
X = -67349030295.003616 Y = 367277788.381377 Z = 0.000000
v_x = -5802.143831 v_y = 36551.513046 v_z = 0.000000

Sonne3
X = 67349030295.003616 Y = -367277788.381377 Z = 0.000000
v_x = 5802.143831 v_y = -36551.513046 v_z = 0.000000

Sonne4
X = 82249763057.203293 Y = 963036708.834408 Z = 0.000000
v_x = -6043.466504 v_y = 96127.188943 v_z = 0.000000

Schönen Gruß

Kibo
02.04.2013, 19:29
Danke, das wird sehr hilfreich sein.

Entweder ist es was Sprachspezifisches, dann muss ich einen Workaround schreiben, oder aber der Fehler liegt wirklich in der Verrechnung der einzelnen Kraftvektoren der einzelnen Gegnerplaneten innerhalb eines Prognoseschrittes. Ich hab das gestern einmal nachgerechnet und nichts Auffälliges gefunden, aber ich werde das noch mal mit deinen Zahlen abgleichen, das beschleunigt die Suche ungemein.
Ich gehe dennoch davon aus, das es nicht so schwierig sein kann das zu lösen. Prinzipiell funktioniert das RK, und meine Schleife liefert im Euler auch beim 4er System das geforderte Ergebnis. Wenn ich jetzt noch die notwendige Zeit zum programmieren finde.... (wohl das größte Hindernis momentan ;))

mfg

Bernhard
02.04.2013, 20:38
Ich gehe dennoch davon aus, das es nicht so schwierig sein kann das zu lösen.
Du musst vor allem beachten, dass man bei den "Runge-Kutta-Beschleunigungen" m1a bis m4a immer die zugehörigen, bzw. korrekten Prognosepositionen einsetzen muss. Prinzipiell haben wir das aber bereits erörtet (s.a. UMas kompakter Beitrag)- Viel Erfolg und lass Dir so viel Zeit, wie nötig.

Kibo
07.04.2013, 17:49
So, es ist mal wieder soweit. Ein neuer Release seit langem.

Die Version enthält einige Verbesserungen:
Bugfixes

- Fehler beim Öffnen einer anderen Planeten Datei beseitigt (Alte Körper wurden nicht sauber gelöscht)
- 2 logische Fehler im RK beseitigt
.............................................-> Falscher Ansatz zur gleichzeitigen Berechnung von a und v
.............................................-> Fehler bei der Verechnung der Kraftvektoren zwischen den Planeten
Kamerabeschleunigung sinnvoller gestaltet
Verbesserte Steuerung

Noch einmal danke an Bernhard und UMa.

Features
Verbessertes Logging: 1 Datei pro Planet. mit Tabulator separierte Zellen, in Excel und normalen Editor lesbar.
automatische Schrittweitenkontrolle
Links zum Download im Programm hinterlegt (Hilfe> Check_for_Updates)

Downloadbar wie immer über:

Rapidshare (http://my.rapidshare.com/Kibo186/13327)
Skydrive (http://sdrv.ms/ZsLvTW)

Das Heun-Verfahren läuft weiterhin nicht sauber, ist im Prinzip analog zum Runge Kutta aber nur mit zwei Schritten gemittelt, liefert trotzdem nicht das gewünschte Ergebnis.
Zum Runge Kutta möchte ich Bernhard nochmal bitten sich das kurz anzusehen. Die Bahnen in der Testdatei mit 4 Sonnen sehen am Anfang richtig aus, doch mit der Zeit werden sie sichtbar instabiler. Jetzt stellt sich mir die Frage ob da nun mein Verfahren falsch rechnet, oder das System prinzipiell nur im Euler stabil ist. Euler und RK werden ja sicher nicht immer zu den selben Ergebnissen kommen, dann hätte man ja auch gleich beim Euler bleiben können.

mfg

Bernhard
07.04.2013, 20:06
Kamerabeschleunigung sinnvoller gestaltet
Verbesserte Steuerung
Kamera und Koordinatensystem sind zweifelsfrei nette Features :) . Bei der Version 1.7.1 fällt mir aber auf, dass Oben <-> Unten anders funktioniert als Links <-> Rechts. Oben <-> Unten ist intuitiv verständlich. Links <-> Rechts ist irgendwie verdreht. Wenn man die Maus nach links bewegt, wird die Kamera nach rechts gesetzt.


Zum Runge Kutta möchte ich Bernhard nochmal bitten sich das kurz anzusehen. Die Bahnen in der Testdatei mit 4 Sonnen sehen am Anfang richtig aus, doch mit der Zeit werden sie sichtbar instabiler.
Da ist noch ein Fehler im Code. Deswegen folgender Vorschlag. Bitte liste hier wieder sämtliche Parameter auf, die zum ersten Schritt führen. Also ein Schritt mit beispielsweise 10000 Sekunden und dabei für jeden Planeten die Ausgabe aller Prognosepunkte und Prognosegeschwindigkeiten, sowie m1a bis m4a.


Euler und RK werden ja sicher nicht immer zu den selben Ergebnissen kommen, dann hätte man ja auch gleich beim Euler bleiben können.
Stimmt. RK ist genauer ;) .
MfG

Kibo
07.04.2013, 22:19
OK, ich hoffe das ist halbwegs verständlich:



Prognose 1:
Sonne1zuSonne2
P1/-82280000000.000000/0/0 zu /-67320000000.000000/0/0
P1v/0/-96392.0/0 zu /0/36816.0/0
P2v/2967.125000/-96392.000000/0.000000 zu /-2967.125000/36816.000000/0.000000
m1a/0.593425/0.000000/0.000000 zu /-0.593425/0.000000/0.000000
P2/-82280000000.000000/-481960000.000000/0.000000 zu /-67320000000.000000/184080000.000000/0.000000
Sonne1zuSonne3
P1/-82280000000.000000/0/0 zu /67320000000.000000/0/0
P1v/0/-96392.0/0 zu /0/-36816.0/0
P2v/29.670000/-96392.000000/0.000000 zu /-29.670000/-36816.000000/0.000000
m1a/0.599359/0.000000/0.000000 zu /-0.005934/0.000000/0.000000
P2/-82280000000.000000/-481960000.000000/0.000000 zu /67320000000.000000/-184080000.000000/0.000000
Sonne1zuSonne4
P1/-82280000000.000000/0/0 zu /82280000000.000000/0/0
P1v/0/-96392.0/0 zu /0/96392.0/0
P2v/24.520000/-96392.000000/0.000000 zu /-24.520000/96392.000000/0.000000
m1a/0.604263/0.000000/0.000000 zu /-0.004904/0.000000/0.000000
P2/-82280000000.000000/-481960000.000000/0.000000 zu /82280000000.000000/481960000.000000/0.000000
Sonne2zuSonne3
P1/-67320000000.000000/0/0 zu /67320000000.000000/0/0
P1v/0/36816.0/0 zu /0/-36816.0/0
P2v/36.630000/36816.000000/0.000000 zu /-36.630000/-36816.000000/0.000000
m1a/-0.586099/0.000000/0.000000 zu /-0.013260/0.000000/0.000000
P2/-67320000000.000000/184080000.000000/0.000000 zu /67320000000.000000/-184080000.000000/0.000000
Sonne2zuSonne4
P1/-67320000000.000000/0/0 zu /82280000000.000000/0/0
P1v/0/36816.0/0 zu /0/96392.0/0
P2v/29.670000/36816.000000/0.000000 zu /-29.670000/96392.000000/0.000000
m1a/-0.580165/0.000000/0.000000 zu /-0.010838/0.000000/0.000000
P2/-67320000000.000000/184080000.000000/0.000000 zu /82280000000.000000/481960000.000000/0.000000
Sonne3zuSonne4
P1/67320000000.000000/0/0 zu /82280000000.000000/0/0
P1v/0/-36816.0/0 zu /0/96392.0/0
P2v/2967.125000/-36816.000000/0.000000 zu /-2967.125000/96392.000000/0.000000
m1a/0.580165/0.000000/0.000000 zu /-0.604263/0.000000/0.000000
P2/67320000000.000000/-184080000.000000/0.000000 zu /82280000000.000000/481960000.000000/0.000000

Prognose 2:
Sonne1zuSonne2
P3v/2958.323358/-96260.292261/0.000000 zu /-2958.323358/36684.292261/0.000000
m2a/0.591665/0.026342/0.000000 zu /-0.591665/-0.026342/0.000000
P2/-82279877400.000000/-481960000.000000/0.000000 zu /-67319851650.000000/184080000.000000/0.000000
Sonne1zuSonne3
P3v/29.669941/-96391.940927/0.000000 zu /-29.669941/-36816.059073/0.000000
m2a/0.597599/0.026353/0.000000 zu /-0.005934/-0.000012/0.000000
P2/-82279877400.000000/-481960000.000000/0.000000 zu /67334835625.000000/-184080000.000000/0.000000
Sonne1zuSonne4
P3v/24.519583/-96391.856386/0.000000 zu /-24.519583/96391.856386/0.000000
m2a/0.602503/0.026382/0.000000 zu /-0.004904/-0.000029/0.000000
P2/-82279877400.000000/-481960000.000000/0.000000 zu /82265164375.000000/481960000.000000/0.000000
Sonne2zuSonne3
P3v/36.629853/36815.899854/0.000000 zu /-36.629853/-36815.899854/0.000000
m2a/-0.584339/-0.026362/0.000000 zu /-0.013260/0.000008/0.000000
P2/-67319851650.000000/184080000.000000/0.000000 zu /67334835625.000000/-184080000.000000/0.000000
Sonne2zuSonne4
P3v/29.669941/36816.059073/0.000000 zu /-29.669941/96391.940927/0.000000
m2a/-0.578405/-0.026350/0.000000 zu /-0.010838/-0.000041/0.000000
P2/-67319851650.000000/184080000.000000/0.000000 zu /82265164375.000000/481960000.000000/0.000000
Sonne3zuSonne4
P3v/2958.323358/-36684.292261/0.000000 zu /-2958.323358/96260.292261/0.000000
m2a/0.578405/0.026350/0.000000 zu /-0.602503/-0.026382/0.000000
P2/67334835625.000000/-184080000.000000/0.000000 zu /82265164375.000000/481960000.000000/0.000000

Prognose 3:
Sonne1zuSonne2
P34/5916.626735/-96128.585412/0.000000 zu /-5916.626735/36552.585412/0.000000
m3a/0.591663/0.026341/0.000000 zu /-0.591663/-0.026341/0.000000
P2/-82279754804.169998/-963918563.860000/0.000000 zu /-67319703300.589996/368160590.730000/0.000000
Sonne1zuSonne3
P34/59.329881/-96391.881874/0.000000 zu /-59.329881/-36816.118126/0.000000
m3a/0.597596/0.026353/0.000000 zu /-0.005933/-0.000012/0.000000
P2/-82279754804.169998/-963918563.860000/0.000000 zu /67349583233.580002/-366842922.610000/0.000000
Sonne1zuSonne4
P34/49.049166/-96391.712665/0.000000 zu /-49.049166/96391.712665/0.000000
m3a/0.602501/0.026382/0.000000 zu /-0.004905/-0.000029/0.000000
P2/-82279754804.169998/-963918563.860000/0.000000 zu /82250416766.419998/962602922.610000/0.000000
Sonne2zuSonne3
P34/73.249707/36815.799734/0.000000 zu /-73.249707/-36815.799734/0.000000
m3a/-0.584338/-0.026361/0.000000 zu /-0.013258/0.000008/0.000000
P2/-67319703300.589996/368160590.730000/0.000000 zu /67349583233.580002/-366842922.610000/0.000000
Sonne2zuSonne4
P34/59.349881/36816.118166/0.000000 zu /-59.349881/96391.881834/0.000000
m3a/-0.578403/-0.026350/0.000000 zu /-0.010840/-0.000041/0.000000
P2/-67319703300.589996/368160590.730000/0.000000 zu /82250416766.419998/962602922.610000/0.000000
Sonne3zuSonne4
P34/5940.109656/-36551.009673/0.000000 zu /-5940.109656/96127.009673/0.000000
m3a/0.580753/0.026507/0.000000 zu /-0.604851/-0.026540/0.000000
P2/67349583233.580002/-366842922.610000/0.000000 zu /82250416766.419998/962602922.610000/0.000000

Prognose 4:
Sonne1zuSonne2
m4a/0.586433/0.052218/0.000000 zu /-0.586433/-0.052218/0.000000
Sonne1zuSonne3
m4a/0.592365/0.052241/0.000000 zu /-0.005932/-0.000024/0.000000
Sonne1zuSonne4
m4a/0.597269/0.052299/0.000000 zu /-0.004905/-0.000057/0.000000
Sonne2zuSonne3
m4a/-0.579110/-0.052258/0.000000 zu /-0.013255/0.000016/0.000000
Sonne2zuSonne4
m4a/-0.573173/-0.052234/0.000000 zu /-0.010842/-0.000081/0.000000
Sonne3zuSonne4
m4a/0.577821/0.052752/0.000000 zu /-0.601917/-0.052817/0.000000


Die Beschleunigungen der Einzelplaneten werden in jedem Teilschritt einer Prognose aufaddiert, sodass m1a für Sonne1 bei Sonne1zuSonne4 der Gesamtvektor aus den vorherigen m1a dieses Schrittes ist

Kibo
07.04.2013, 22:26
Und hier noch einmal die Endergebnisse, damit es übersichtlich wird:


Sonne1

dt= 10000
P1ort= /-82280000000.000000/0/0
P1v = /0/-96392.0/0
m1a= /0.604263/0.000000/0.000000
P2ort= /-82280000000.000000/-481960000.000000/0.000000
P2v = /24.520000/-96392.000000/0.000000
m2a= /0.602503/0.026382/0.000000
P3ort= /-82279877400.000000/-481960000.000000/0.000000
P3v = /24.519583/-96391.856386/0.000000
m3a= /0.602501/0.026382/0.000000
P4ort= /-82279754804.169998/-963918563.860000/0.000000
P4v = /49.049166/-96391.712665/0.000000
m4a= /0.597269/0.052299/0.000000
ort_neu= -82279754785.622910/-963920970.233085/0.000000
v_neu = /6019.231080/-96128.955071/0.000000

Sonne2
dt= 10000
P1ort= /-67320000000.000000/0/0
P1v = /0/36816.0/0
m1a= /-0.580165/0.000000/0.000000
P2ort= /-67320000000.000000/184080000.000000/0.000000
P2v = /29.670000/36816.000000/0.000000
m2a= /-0.578405/-0.026350/0.000000
P3ort= /-67319851650.000000/184080000.000000/0.000000
P3v = /29.669941/36816.059073/0.000000
m3a= /-0.578403/-0.026350/0.000000
P4ort= /-67319703300.589996/368160590.730000/0.000000
P4v = /59.349881/36816.118166/0.000000
m4a= /-0.573173/-0.052234/0.000000
ort_neu= -67319703283.134903/368161130.174121/0.000000
v_neu = /-5778.254562/36553.278469/0.000000

Sonne3
dt= 10000
P1ort= /67320000000.000000/0/0
P1v = /0/-36816.0/0
m1a= /0.580165/0.000000/0.000000
P2ort= /67320000000.000000/-184080000.000000/0.000000
P2v = /2967.125000/-36816.000000/0.000000
m2a= /0.578405/0.026350/0.000000
P3ort= /67334835625.000000/-184080000.000000/0.000000
P3v = /2958.323358/-36684.292261/0.000000
m3a= /0.580753/0.026507/0.000000
P4ort= /67349583233.580002/-366842922.610000/0.000000
P4v = /5940.109656/-36551.009673/0.000000
m4a= /0.577821/0.052752/0.000000
ort_neu= 67349651736.590019/-367280058.216981/0.000000
v_neu = /5793.835026/-36551.889933/0.000000

Sonne4
dt= 10000
P1ort= /82280000000.000000/0/0
P1v = /0/96392.0/0
m1a= /-0.604263/0.000000/0.000000
P2ort= /82280000000.000000/481960000.000000/0.000000
P2v = /-2967.125000/96392.000000/0.000000
m2a= /-0.602503/-0.026382/0.000000
P3ort= /82265164375.000000/481960000.000000/0.000000
P3v = /-2958.323358/96260.292261/0.000000
m3a= /-0.604851/-0.026540/0.000000
P4ort= /82250416766.419998/962602922.610000/0.000000
P4v = /-5940.109656/96127.009673/0.000000
m4a= /-0.601917/-0.052817/0.000000
ort_neu= 82250348263.409973/963041249.736981/0.000000
v_neu = /-6034.811544/96127.566535/0.000000

Kibo
07.04.2013, 22:48
Hier noch schnell das Log nur für mXa aber mit Heraushebung des Verrechnungsvorgangs, wahrscheinlich liegt der Fehler ja da.
Ich kopier das noch schnell hier rein und dann geh ich aber auch ins Bett. Gute Nacht



Sonne1 zu Sonne2
m1a/0.593425/0.000000/0.000000=/0/0/0+/0.593425/0.000000/0.000000 zu /-0.593425/0.000000/0.000000=/0/0/0+/-0.593425/-0.000000/-0.000000
Sonne1 zu Sonne3
m1a/0.599359/0.000000/0.000000=/0.593425/0.000000/0.000000+/0.005934/0.000000/0.000000 zu /-0.005934/0.000000/0.000000=/0/0/0+/-0.005934/-0.000000/-0.000000
Sonne1 zu Sonne4
m1a/0.604263/0.000000/0.000000=/0.599359/0.000000/0.000000+/0.004904/0.000000/0.000000 zu /-0.004904/0.000000/0.000000=/0/0/0+/-0.004904/-0.000000/-0.000000
Sonne2 zu Sonne3
m1a/-0.586099/0.000000/0.000000=/-0.593425/0.000000/0.000000+/0.007326/0.000000/0.000000 zu /-0.013260/0.000000/0.000000=/-0.005934/0.000000/0.000000+/-0.007326/-0.000000/-0.000000
Sonne2 zu Sonne4
m1a/-0.580165/0.000000/0.000000=/-0.586099/0.000000/0.000000+/0.005934/0.000000/0.000000 zu /-0.010838/0.000000/0.000000=/-0.004904/0.000000/0.000000+/-0.005934/-0.000000/-0.000000
Sonne3 zu Sonne4
m1a/0.580165/0.000000/0.000000=/-0.013260/0.000000/0.000000+/0.593425/0.000000/0.000000 zu /-0.604263/0.000000/0.000000=/-0.010838/0.000000/0.000000+/-0.593425/-0.000000/-0.000000
Sonne1 zu Sonne2
m2a/0.591665/0.026342/0.000000=/0/0/0+/0.591665/0.026342/0.000000 zu /-0.591665/-0.026342/0.000000=/0/0/0+/-0.591665/-0.026342/-0.000000
Sonne1 zu Sonne3
m2a/0.597599/0.026353/0.000000=/0.591665/0.026342/0.000000+/0.005934/0.000012/0.000000 zu /-0.005934/-0.000012/0.000000=/0/0/0+/-0.005934/-0.000012/-0.000000
Sonne1 zu Sonne4
m2a/0.602503/0.026382/0.000000=/0.597599/0.026353/0.000000+/0.004904/0.000029/0.000000 zu /-0.004904/-0.000029/0.000000=/0/0/0+/-0.004904/-0.000029/-0.000000
Sonne2 zu Sonne3
m2a/-0.584339/-0.026362/0.000000=/-0.591665/-0.026342/0.000000+/0.007326/-0.000020/0.000000 zu /-0.013260/0.000008/0.000000=/-0.005934/-0.000012/0.000000+/-0.007326/0.000020/-0.000000
Sonne2 zu Sonne4
m2a/-0.578405/-0.026350/0.000000=/-0.584339/-0.026362/0.000000+/0.005934/0.000012/0.000000 zu /-0.010838/-0.000041/0.000000=/-0.004904/-0.000029/0.000000+/-0.005934/-0.000012/-0.000000
Sonne3 zu Sonne4
m2a/0.578405/0.026350/0.000000=/-0.013260/0.000008/0.000000+/0.591665/0.026342/0.000000 zu /-0.602503/-0.026382/0.000000=/-0.010838/-0.000041/0.000000+/-0.591665/-0.026342/-0.000000
Sonne1 zu Sonne2
m3a/0.591663/0.026341/0.000000=/0/0/0+/0.591663/0.026341/0.000000 zu /-0.591663/-0.026341/0.000000=/0/0/0+/-0.591663/-0.026341/-0.000000
Sonne1 zu Sonne3
m3a/0.597596/0.026353/0.000000=/0.591663/0.026341/0.000000+/0.005933/0.000012/0.000000 zu /-0.005933/-0.000012/0.000000=/0/0/0+/-0.005933/-0.000012/-0.000000
Sonne1 zu Sonne4
m3a/0.602501/0.026382/0.000000=/0.597596/0.026353/0.000000+/0.004905/0.000029/0.000000 zu /-0.004905/-0.000029/0.000000=/0/0/0+/-0.004905/-0.000029/-0.000000
Sonne2 zu Sonne3
m3a/-0.584338/-0.026361/0.000000=/-0.591663/-0.026341/0.000000+/0.007325/-0.000020/0.000000 zu /-0.013258/0.000008/0.000000=/-0.005933/-0.000012/0.000000+/-0.007325/0.000020/-0.000000
Sonne2 zu Sonne4
m3a/-0.578403/-0.026350/0.000000=/-0.584338/-0.026361/0.000000+/0.005935/0.000012/0.000000 zu /-0.010840/-0.000041/0.000000=/-0.004905/-0.000029/0.000000+/-0.005935/-0.000012/-0.000000
Sonne3 zu Sonne4
m3a/0.580753/0.026507/0.000000=/-0.013258/0.000008/0.000000+/0.594011/0.026499/0.000000 zu /-0.604851/-0.026540/0.000000=/-0.010840/-0.000041/0.000000+/-0.594011/-0.026499/-0.000000
Sonne1 zu Sonne2
m4a/0.586433/0.052218/0.000000=/0/0/0+/0.586433/0.052218/0.000000 zu /-0.586433/-0.052218/0.000000=/0/0/0+/-0.586433/-0.052218/-0.000000
Sonne1 zu Sonne3
m4a/0.592365/0.052241/0.000000=/0.586433/0.052218/0.000000+/0.005932/0.000024/0.000000 zu /-0.005932/-0.000024/0.000000=/0/0/0+/-0.005932/-0.000024/-0.000000
Sonne1 zu Sonne4
m4a/0.597269/0.052299/0.000000=/0.592365/0.052241/0.000000+/0.004905/0.000057/0.000000 zu /-0.004905/-0.000057/0.000000=/0/0/0+/-0.004905/-0.000057/-0.000000
Sonne2 zu Sonne3
m4a/-0.579110/-0.052258/0.000000=/-0.586433/-0.052218/0.000000+/0.007323/-0.000040/0.000000 zu /-0.013255/0.000016/0.000000=/-0.005932/-0.000024/0.000000+/-0.007323/0.000040/-0.000000
Sonne2 zu Sonne4
m4a/-0.573173/-0.052234/0.000000=/-0.579110/-0.052258/0.000000+/0.005937/0.000024/0.000000 zu /-0.010842/-0.000081/0.000000=/-0.004905/-0.000057/0.000000+/-0.005937/-0.000024/-0.000000
Sonne3 zu Sonne4
m4a/0.577821/0.052752/0.000000=/-0.013255/0.000016/0.000000+/0.591075/0.052736/0.000000 zu /-0.601917/-0.052817/0.000000=/-0.010842/-0.000081/0.000000+/-0.591075/-0.052736/-0.000000

Bernhard
08.04.2013, 02:28
Sonne1

dt= 10000
P1ort= /-82280000000.000000/0/0
P1v = /0/-96392.0/0
m1a= /0.604263/0.000000/0.000000
P2ort= /-82280000000.000000/-481960000.000000/0.000000
P2v = /24.520000/-96392.000000/0.000000
Hallo Kibo,

da sieht man sofort einen ziemlich groben Fehler, denn es soll ja gelten:
P2v = P1v + dt/2 * m1a

Schau Dir mal die x-Komponenten an. P1v_x = 0.0. Damit müsste gelten P2v_x = dt/2 * m1a_x = 5000 * m1a_x, mit m1a_x = 0.604263 ergibt das 3012.315 und nicht 24.52. Dieser Fehler zieht sich durch alle Prognosegeschwindigkeiten und sollte zuerst behoben werden.

Anschließend können wir uns bei Bedarf noch die mxa anschauen. Dort sind die Abweichungen zu meinem Programm aber relativ gering.
MfG

Kibo
08.04.2013, 22:17
So, das war endlich mal ein einfach zu beseitigender Fehler. Hab die Berechnung a->v aus der Berechnungsfunktion herausgenommen und stattdessen die eine Zeile nach der Berechnung des Gesamt-a der jeweiligen Prognose 1, 2, 3, 4 nochmal für alle Planeten eingefügt. Jetzt addiert er richtig und is sogar noch ein Stückchen schneller geworden.
Bernhard wenn du jetzt noch dein OK gibst, würde ich sagen Baustelle Runge-Kutta ist bis auf Weiteres geschlossen:D
Die Sache mit der Maus habe ich natürlich berichtigt, die neue Version heißt 1.7.2

mfg

Downloadlink ist wie immer:
Rapidshare (http://my.rapidshare.com/Kibo186/13327)
Skydrive (http://sdrv.ms/ZsLvTW)

Bernhard
08.04.2013, 23:07
Bernhard wenn du jetzt noch dein OK gibst, würde ich sagen Baustelle Runge-Kutta ist bis auf Weiteres geschlossen:D
Yep. Die zwei Beispiele werden jetzt korrekt dargestellt :) .

Bernhard
10.04.2013, 00:12
Bernhard wenn du jetzt noch dein OK gibst, würde ich sagen Baustelle Runge-Kutta ist bis auf Weiteres geschlossen:D
Hallo Kibo,

ich hoffe Du hast den Erfolg inzwischen ausreichend gefeiert, denn eine Idee hätte ich noch. Baue doch bitte die log-Datei erneut aus der Version 1.7.2 aus. Die Zwischenschritte der Berechnung, insbesondere die Werte von Hilfsvariablen wie m1a usw. sind nämlich nur für die Fehlersuche wichtig.

Abschließend wäre es noch interessant den RK-Integrator mit den Startbedingungen aus dem Thema "Sonnensystem: Position der Planeten" zu füttern, um letzte Mini-Fehler zu finden, bzw. um einen Vergleich zur Rechengenauigkeit von C++ zu haben. Dazu wäre es aber wichtig in einer Log-Datei zumindest die Endpositionen und Endgeschwindigkeiten aller Körper mit allen berechneten Stellen auszugeben. Wenn Du willst wäre auch ein komplettes Zeitlogging nett. Also Endpositionen und Endgeschwindigkeiten aller Körper nach jedem Zeitschritt mit Angabe der Zeit in Sekunden ab Start. Eine blockweise Ausgabe der Daten wäre dabei sinnvoll.
MfG

Kibo
10.04.2013, 09:50
Mach ich dann nach der Arbeit fertig :)

Kibo
10.04.2013, 20:58
So, deine Vorschläge sind in der 1.7.3 eingearbeitet (zu Finden unter Einstellungen>Featureeinstellungen>Logdatei v und k als Tabelle).
Sonst hat sich nichts verändert. Falls dir noch etwas falsch erscheinen sollte, sag bescheid.

mfg

Bernhard
11.04.2013, 20:02
Hallo Kibo,


Falls dir noch etwas falsch erscheinen sollte, sag bescheid.
schaue es Dir bitte selbst mal an:


KoordX KoordY Koordz ImpulsX ImpulsY ImpulsZ Masse dichte s v Name s und v sind noch nicht implementiert
-6.080091705285974E+08 2.454333708238292E+07 2.193321432139981E+06 2.809726412820165E+00 -1.036343441959919E+01 -4.204860196381065E-02 1.9897E+30 0 0 Sun 1.4
-5.331602529851193E+10 1.424323914756870E+10 6.000210156549231E+9 -2.277058581046753E+04 -4.495485791176564E+04 -1.582633849867808E+03 3.3032E23 0 0 Mercury 1.4
2.397967758924109E+10 -1.059717457755222E+11 -2.868972017005123E+09 3.388294521475002E+04 7.783094129350907E+03 -1.848596675496016E+03 4.8704E24 0 0 Venus 1.4
-1.474370019336122E+11 -2.789279589304088E+10 2.955826247744262E+06 5.079162483215040E+03 -2.939896910566963E+04 4.812093925554706E-01 5.976E24 0 0 Earth 1.4
-1.470514277349460E+11 -2.801509364690417E+10 3.730129039033502E+07 5.382156089456466E+03 -2.847733344864282E+04 2.235354455693894E+01 7.3505E22 0 0 Moon 1.4
2.035953497078640E+11 -3.456857769669849E+10 -5.736961786047079E+09 4.976094385716297E+03 2.594977520891352E+04 4.218104801061298E+02 6.421E23 0 0 Mars 1.4
7.115603707882032E+11 2.013825854433948E+11 -1.677091861778326E+10 -3.715952626588193E+03 1.319564076159422E+04 2.830947269851247E+01 1.8997E27 0 0 Jupiter 1.4
-1.396843425883871E+12 -3.381973875783179E+11 6.146826298318255E+10 1.753564063629373E+03 -9.410322369209542E+03 9.449519324100386E+01 5.6882E26 0 0 Saturn 1.4
3.003955040895567E+12 2.614531534364399E10 -3.882080718055668E+10 -1.089777236983593E2 6.492185531504123E+03 2.552972796165465E1 8.6875E25 0 0 Uranus 1.4
3.826621024919128E+12 -2.346728069113412E+12 -3.986079005702972E+10 2.805108135465661E+03 4.664771972182217E+03 -1.611737087776297E2 1.025E26 0 0 Neptune 1.4
4.599506779023092E+11 -4.750680975824932E+12 3.753056255249584E+11 5.510641913253943E+03 -5.401707214265152E-02 -1.536200853856072E+03 1.4717E22 0 0 Pluto 1.4

Anstelle schöner Planetenbahnen sieht man da seltsam gerade Bahnen :confused: .
MfG

Kibo
11.04.2013, 23:03
Ich habe die die Berechnung der Beschleunigung explizit auf Genauigkeit 10.15 für Floats gesetzt, jetzt sieht es relativ normal aus, soweit ich das beurteilen kann. Diese Änderung findet sich in 1.7.4

gute Nacht

Bernhard
12.04.2013, 21:56
Ich habe die die Berechnung der Beschleunigung explizit auf Genauigkeit 10.15 für Floats gesetzt, jetzt sieht es relativ normal aus
Hallo Kibo,

es ist etwas eigenartig, dass eine so kleine Änderung eine so große Wirkung zeigt. Wenn Du noch Lust hast, könntest Du deswegen den Masseparameter in der Eingabedatei durch das Produkt aus Gravitationskonstante mal Masse ersetzen. Horizons liefert für dieses Produkt sehr genaue Werte. Der Code wird dadurch auch etwas übersichtlicher, weil eben diese Multiplikation aus dem Code komplett entfernt werden kann. Die verwendeten Zahlen für Masse, Ort und Geschwindigkeit liegen dann im Exponenten auch nicht mehr so weit auseinander, wodurch Rundungsfehler vermieden werden.
MfG

Kibo
12.04.2013, 22:19
Die alleinige Masse wird im Schleifenteil so gar nicht mehr verwendet sondern stattdessen Variable GM die beim Start des Programmes einmalig berechnet wird. Der Knackpunkt ist an dieser Stelle Newton 3, denn:

Der Simulator berechnet zuerst die Beschleunigung von Körper 1 zu allen anderen. Da Körper 1 die Sonne in dem Falle ist, kommen da sehr geringe Beschleunigungen bei heraus wie 0.000000007014216 (in dem fall xWert für m1a für Sonne zu Merkur als Beispiel). Da wurde in der Programmiersprache für unsere Zwecke standartmäßig unvorteilhaft gerundet. Ich hab ihm dann einfach explizit gesagt das er möglichst genau rechnen soll und gut wars.

Edit: der Menge an Rechtschreibfehlern nach sollte ich ins Bett gehen
gute Nacht

Kibo
16.04.2013, 22:17
Version 1.7.5 ist fertig!

Bugfixes

Steuerung mit A, D hat nicht wie gewünscht funktioniert - behoben
Texturen wurden nicht ordnunggemäß angezeigt - behoben, es werden kleine Beispieltexturen erstellt

Features
Texturen lassen sich jetzt im img Ordner abändern, Namen kann man in settings.ini anpassen
Fadenkreuz hinzugefügt, im Settings.ini standartmäßig ein/ausschaltbar (1/0), im Menü Grafikeinstellungen im Laufenden Betrieb ein/aussschaltbar
Genauigkeit der Renderung von Kugelformen lässt sich über setting.ini und im Menü Grafikeinstellungen anpassen


Andere Anpassungen
In der Planetendatei wurden s und v entfernt, da obsolet. Planetendatein sollten an das neue Format angepasst werden


Beispiel Sonnesystem:
KoordX KoordY Koordz ImpulsX ImpulsY ImpulsZ Masse dichte s v Name s und v sind noch nicht implementiert
-6.080091705285974E+08 2.454333708238292E+07 2.193321432139981E+06 2.809726412820165E+00 -1.036343441959919E+01 -4.204860196381065E-02 1.9897E+30 Sun 1.4
-5.331602529851193E+10 1.424323914756870E+10 6.000210156549231E+9 -2.277058581046753E+04 -4.495485791176564E+04 -1.582633849867808E+03 3.3032E23 Mercury 1.4
2.397967758924109E+10 -1.059717457755222E+11 -2.868972017005123E+09 3.388294521475002E+04 7.783094129350907E+03 -1.848596675496016E+03 4.8704E24 Venus 1.4
-1.474370019336122E+11 -2.789279589304088E+10 2.955826247744262E+06 5.079162483215040E+03 -2.939896910566963E+04 4.812093925554706E-01 5.976E24 Earth 1.4
-1.470514277349460E+11 -2.801509364690417E+10 3.730129039033502E+07 5.382156089456466E+03 -2.847733344864282E+04 2.235354455693894E+01 7.3505E22 Moon 1.4
2.035953497078640E+11 -3.456857769669849E+10 -5.736961786047079E+09 4.976094385716297E+03 2.594977520891352E+04 4.218104801061298E+02 6.421E23 Mars 1.4
7.115603707882032E+11 2.013825854433948E+11 -1.677091861778326E+10 -3.715952626588193E+03 1.319564076159422E+04 2.830947269851247E+01 1.8997E27 Jupiter 1.4
-1.396843425883871E+12 -3.381973875783179E+11 6.146826298318255E+10 1.753564063629373E+03 -9.410322369209542E+03 9.449519324100386E+01 5.6882E26 Saturn 1.4
3.003955040895567E+12 2.614531534364399E10 -3.882080718055668E+10 -1.089777236983593E2 6.492185531504123E+03 2.552972796165465E1 8.6875E25 Uranus 1.4
3.826621024919128E+12 -2.346728069113412E+12 -3.986079005702972E+10 2.805108135465661E+03 4.664771972182217E+03 -1.611737087776297E2 1.025E26 Neptune 1.4
4.599506779023092E+11 -4.750680975824932E+12 3.753056255249584E+11 5.510641913253943E+03 -5.401707214265152E-02 -1.536200853856072E+03 1.4717E22 Pluto 1.4

Download wie immer:
Rapidshare (http://my.rapidshare.com/Kibo186/13327)
Skydrive (http://sdrv.ms/ZsLvTW)

Bernhard
17.04.2013, 21:27
Hallo Kibo,

es gibt bei mir immer wieder Probleme beim Start des Programmes. Kannst Du vielleicht noch die Lese-Funktion beim Programmstart so anpassen, dass die erste Zeile prinzipiell als Kommentar interpretiert wird. Für die Folgezeilen wäre ein Blank als Trennzeichen auch besser als ein Tabulator, weil gebräuchlicher.
MfG

Kibo
17.04.2013, 22:04
Hallo Bernhard,

Also prinzipiell wird die erste Zeile einer Planeten Datei bei mir immer als Kommentar behandelt, das war glaube ich sogar beid er ersten Version so. Jedenfalls ist in Version 1.8.1 das Trennzeichen frei in der settings.ini wählbar.
EDIT: Führe die frisch herunter geladene Exe-Datei am besten in einem komplett leeren Ordner aus, damit er alle Dateien einmal neu erstellt, die kannst du ja nach deinen Wünschen dann wieder anpassen.

Um Missverständnissen vorzubeugen erläutere Ich noch einmal kurz den Aufbau der der Datendatei, wie sie denn jetzt seit 1.7.5 nun ist:

Kommentarzeile, egal was steht
KoordX KoordY Koordz ImpulsX ImpulsY ImpulsZ Masse Name dichte
KoordX KoordY Koordz ImpulsX ImpulsY ImpulsZ Masse Name dichte
KoordX KoordY Koordz ImpulsX ImpulsY ImpulsZ Masse Name dichte
.
.
.


Rapidshare (http://my.rapidshare.com/Kibo186/13327)
Skydrive (http://sdrv.ms/ZsLvTW)

Ergänzung: wenn man nicht sichtbare Zeichen wie Tab oder space verwendet, sollten die in " " gesetzt werden um erkannt zu werden

Kibo
18.04.2013, 22:40
Schon wieder ein Update :)

Neben jedem Himmelskörper nun (Version 1.8.2) auch seine Bezeichnung.

mfg

Bernhard
19.04.2013, 00:38
Neben jedem Himmelskörper nun (Version 1.8.2) auch seine Bezeichnung.
Das sieht jetzt richtig fetzig aus :cool: .

Kibo
21.04.2013, 19:20
Mal wieder ein Update, mal wieder ein paar nette Features, die 1.9.1 mit sich bringt:

Das Seitenmenü wurde etwas aufgeräumt, alle Funktionen wurden in 3 Tabs untergebracht.

Neuer Verfolgungsmodus, man kann sich an einen Körper seiner Wahl ran hängen,verschafft einen, eine nette Aussicht auf den Mond von der erde aus :D

Mit dem Scrollrad wird jetzt die Bewegungsgeschwindigkeit der Kamera gesteuert, der ohnehin ziemlich sinnlose Zoom fällt, dafür weg.

Geschwindigkeit und Position werden jetzt im linken unteren Renderbereich mit angezeigt.

Die Maussteuerung wurde komplett neu gestaltet:
Mit klick in den Renderbereich verschwindet der Mauscursor und man kann sich jetzt per Mausbewegung umschauen ohne an irgendwelche Ränder zu stoßen. Mit einen weiteren Klick verlässt man diesen Modus wieder.

Unter Hilfe ist jetzt auch eine kleine Anleitung enthalten.

Downloadlinks wie immer:
Rapidshare (http://my.rapidshare.com/Kibo186/13327)
Skydrive (http://sdrv.ms/ZsLvTW)

Bernhard
21.04.2013, 21:34
Mal wieder ein Update, mal wieder ein paar nette Features, die 1.9.1 mit sich bringt:
Hallo Kibo,

ich habe auf sourceforge.net einen Account unter dem Nicknamen astronews334 (http://sourceforge.net/users/astronews334). Wenn Du das Programm inklusive Quelltext veröffentlichen willst, kann ich das unter diesem Account machen. Dann wird Dein Programm weltweit verteilt. Ohne Quelltext geht es allerdings nicht, da solche Projekte nach einer Meldung durch andere Teilnehmer entfernt werden.
Gruß

Kibo
23.04.2013, 21:04
Hi Bernhard,

du hast Post.


mfg

Bernhard
23.04.2013, 22:18
Hier ist der Link: http://sourceforge.net/projects/smallplanetarys/

Für die offenen Kommentarfelder, Beschreibungen usw. werden Vorschläge in Englisch gerne entgegen genommen.
MfG

Bernhard
24.04.2013, 12:02
Hallo Kibo,

wie kann man bei der neuen Version die Position der Kamera verändern? Das Scrollrad der Maus funktioniert leider nicht mehr und die Eingabefelder auf dem Karteireiter "Kamera" springen immer wieder auf den vorgegebenen Wert zurück.
MfG

Kibo
24.04.2013, 12:48
Hallo Bernhard,

Die Tasten für vorwärts, rückwärts, nach rechts und nach links sind mit W, S, A und D belegt. Das Scrollrad dient zum Einstellen der Geschwindigkeit in der sich die Kamera bewegt wenn man auf diese Tasten drückt.

Mit den Editfeldern scheint wirklich etwas nicht zu stimmen, die schaue ich mir noch mal genauer an. Aber erst nach der Arbeit.:)

mfg

Kibo
24.04.2013, 20:53
1.9.2 Steht am üblichen Link zum Download bereit, einzige Veränderung, ist, dass die Editfelder für die Kamera jetzt funktionieren.

mfg

Kibo
28.04.2013, 16:48
Wäre das Buch zu empfehlen, oder gibt es was Besseres für Symplectische Integration?
http://www.amazon.de/gp/product/3834804118/ref=s9_simh_gw_p14_d3_i1?pf_rd_m=A3JWKAKR8XB7XF&pf_rd_s=center-2&pf_rd_r=0G4XPZM0DXX6GNGGR36H&pf_rd_t=101&pf_rd_p=463375173&pf_rd_i=301128

mfg

Bernhard
28.04.2013, 18:17
Hallo Kibo,

besser wäre ein Lehr- und Übungsbuch zur hamiltonschen Mechanik, wie dieses hier (Amazon) (http://www.amazon.de/Mechanik-Lehrbuch-zur-Theoretischen-Physik/dp/3827421489/ref=sr_1_1?s=books&ie=UTF8&qid=1367165651&sr=1-1&keywords=torsten+flie%C3%9Fbach). Du kannst das Buch auch zum Selbststudium verwenden, weil es etliche Übungaufgaben enthält.
MfG

Bernhard
01.05.2013, 12:05
'Ich' hat im Quantenforum ein google-Buch verlinkt, was zu diesem Thema hier auch passt.

Sergei Kopeikin et al., "Relativistic Celestial Mechanics of the Solar System" (http://books.google.de/books?id=uN5_DQWSR14C&pg=PT860&lpg=PT860&dq=%22wave+vector%22+%22affine+parameter%22&source=bl&ots=bVV7nj-DuG&sig=AI6VCJ004k4QC7DzFTUmGFj2OHc&hl=de&sa=X&ei=cc5_UYaSG7PQ4QTWn4DYAw&sqi=2&redir_esc=y#v=onepage&q&f=false) (voreingestellte Seite hat für dieses Thema keine besondere Bedeutung).

Kibo
01.05.2013, 16:24
Sehr interessant, ich hab mir übrigends das von dir empfohlene Buch bestellt.
danke

Bernhard
01.05.2013, 19:23
ich hab mir übrigends das von dir empfohlene Buch bestellt.
OK. Damit kannst Du Dir im Idealfall zuerst die Lagrange-Mechanik erarbeiten. Der Übergang zur hamiltonschen Mechanik ist am Anfang dann entweder eine formale Angelegenheit oder Du schaust Dir zusätzlich auch die guten Wikipedia-Artikel zur Legendre-Transformation an. Diese Transformationen bewirken nämlich genau den Übergang von Lagrange nach Hamilton und viel mehr lernt man dann im Grundstudium Physik in der Mechanik-Vorlesung auch nicht.

Sollten Dir die Übungsaufgaben aus dem Fließbach nicht reichen oder zu leicht erscheinen, kannst Du zusätzlich/optional auch noch den Landau-Lifschitz, Band 1, "Mechanik" bearbeiten. Der verwendet aber dummerweise etwas andere Symbole als der Fließbach.
MfG

Kibo
09.05.2013, 20:59
Mal wieder ein kleines Update, (nein noch nichts mit symplektischen Integrator, das dauert noch^^).

Version 1.10 beinhaltet nun eine Kollisionsdetektion. Man kann sie per Häkchen im Simulations-tab aktivieren (standartmäßig deaktiviert).
Desweiteren kann man jetzt über " Hilfe> neueste Version " direkt die neueste Version von Sourceforge herunterladen.

Bugfixes:
Per Scrollrad übergebene Kamerageschwindigkeit wurde dauernd überschrieben.
Planetengröße wurde bei Echtzeit-Änderung der Masse nicht angepasst

Download (immer neueste Version):

Sourceforge (https://sourceforge.net/projects/smallplanetarys/files/latest/download?source=files)

Bernhard
10.05.2013, 01:31
Mal wieder ein kleines Update
Hallo Kibo,

damit die Dateien bei Sourceforge übersichtlich bleiben, wäre es gut, wenn Du auch bei den Quelltexten im Dateinamen immer eine Versionsnummer angeben würdest. Sonst weiß man irgendwann nicht mehr, welcher Quelltext zu welcher "exe" gehört. Sourceforge hätte sogar eine brauchbare Versionsverwaltung für den Quelltext, aber ich weiß leider nicht, wie man dort ein Repository anlegt.
MfG

Kibo
10.05.2013, 08:23
geändert

mfg

Bernhard
10.05.2013, 11:35
Prima. Sourceforge hat übrigens auch einen ganz guten Support. Soll ich mal nachfragen, ob ich das Projekt komplett an Deinen Account abgeben kann?

Boogi
10.05.2013, 13:17
Hi Kibo,

habe gerade Deinen Simulator ausprobiert. Finde ich klasse, dass Du Dich mit so was beschäftigst!
Allerdings bin ich nicht mit dem Programm zurecht gekommen (ich habe Version 1.10 oder 1.9.1 oder 1.1.9.4 verwendet. Wieso gibt's so viele Versionsnummern?)

- Nach dem Laden der Planeten habe ich 4 Sonnen, die aber 2x in der Liste auftauchen.
- Ich sehe aber nur eine Sonne.
- Nichts bewegt sich, obwohl die Simulation läuft.
- Die Steuerung der Blickrichtung ist äh... hakelig.
- Bewegungen mit WASD scheinen keinen Effekt zu haben.
- Das Mausrad scheint keinen Effekt zu haben.
- Was mache ich mit dem Kamera-Menü? Wieso kann ich dort keine Körper einstellen, auf die die Kamera dann zeigen soll?
- Wieso ging mein Browser auf, als ich Dein Programm laufen hatte?

Vielleicht bin ich auch nur zu doof dafür...

Gruß
Boogi

Kibo
10.05.2013, 14:40
Hallo Boogi,

Zunächst solltest du immer die neueste Version verwenden, die nächst Ältere liegt nur da falls die neue Macken hat. Dieses Programm befindet sich noch in der Entwicklung und ist dementsprechend auch noch nicht auf Anwenderfreundlichkeit getrimmt, das ändert sich natürlich nach und nach wenn die großen Baustellen geschlossen sind.
Die Kameraausrichtung kann man NUR mit der Maus verändern, um das zu tun muss man in den Renderbereich klicken. Die Simulation läuft relativ realitätsnah ab. Das bedeutet Körper die sehr weit weg sind, bewegen sich auf dem Bildschirm eben nur ein bisschen und sind auch sehr klein. Wir reden hier über Entfernungen von hunderttausenden von Kilometern. Das ist auch wahrscheinlich der Grund, warum du kaum eine Veränderung bemerkst wenn du versuchst dich zu bewegen. Man muss entsprechend der Entfernung auch eine relativ Große Geschwindigkeit wählen, sonst dauert eine Reise von Planet zu Planet, wie in der Realität auch mal Jahre. Am einfachsten verändert man seine Geschwindigkeit mit dem Scrollrad, du kannst deine Position und Geschwindigkeit im Renderbereich unten links ablesen.

Ich werde deine Kritik ernst nehmen und die Startwerte umgehend anpassen, sodass beim Erststart ein System erstellt wird, in dem man auch was sehen kann und in dem Man auch gleich eine Geschwindigkeit hat von der man auch was merkt:).

mfg

Boogi
10.05.2013, 15:19
Ok... hab's jetzt hinbekommen. Aus irgend einem Grund stand in der planeten.txt Müll. Deshalb waren auch die Sonnen doppelt vorhanden. Keine Ahnung warum.

Noch ein paar Ideen:
- Resetknopf für die Grafik: Die Kamera wird so eingestellt, dass alles im Blick ist.
- Wenn ich beispielsweise die Zeitschrittweite ändere wird die Kamera geändert. Warum?
- Die Maus läuft falsch rum, Wenn ich nach links fahre, dann will ich auch nach links schauen.
- Unser Planetensystem als default wäre cool.
Mehr fällt mir grad nicht auf...

Gruß
Boogi

Kibo
10.05.2013, 15:50
So, du kannst ja jetzt mal auf "neueste Version" klicken:)

Das mit der Maus müsste zumindest am Anfang in Ordnung sein, bei bestimmten Bewegungen kippt die Kamera einfach um. Die ausrichtung der Kamera ändert sich nach der Position der Maus zum Zeitpunkt des Klickens das ist nicht unbedingt gewollt aber jetzt auch nicht so schlimm und nicht so einfach zu beseitigen.

mfg

Boogi
10.05.2013, 15:58
Yay! Es geht!

Die Datei heißt jetzt aber "ERROR" statt planeten.txt

Kibo
10.05.2013, 16:27
Die Ladezeit der Settings.ini war zu kurz für deinen PC eingestellt, hab das ganze jetzt etwas konservativer eingestellt sodass es jetz gehen müsste.
Also einmal noch "neueste Version" und am besten alle vorherigen Dateien (also settings.ini, planeten.txt, alte exe) löschen damit es sauber neu erstellt wird.

Bernhard
11.05.2013, 18:13
(nein noch nichts mit symplektischen Integrator, das dauert noch^^).
Hallo Kibo,

ohne den Quelltext genauer angesehen zu haben, brauchst Du für den symplektischen Integrator als Vorbereitung eine Funktion, die als Eingabewerte alle Positionen und Geschwindigkeiten aller Körper und die Schrittweite dt (in Sekunden) bekommt. Diese Funktion soll dann die neuen Positionen und Geschwindigkeiten aller Körper nach dt Sekunden ausrechnen. Du kannst, falls noch nicht gemacht, den vorhandenen Runge-Kutta-Code entsprechend anpassen, dass es genau so eine Funktion gibt.

Für den symplektischen Integrator braucht man so eine Funktionsschnittstelle dann zweimal mit unterschiedlicher Funktionalität und unterschiedlichem Wert für die Schrittweite. EDIT: Die erste Funktion beschreibt dann die kinetische Energie aller Körper und kann sehr leicht programmiert werden. Die zweite Funktion beschreibt die potentielle Energie aller Körper und kann relativ leicht aus der im ersten Absatz beschriebenen Funktion des Runge-Kutta-Codes abgeleitet werden.
MfG

Kibo
20.05.2013, 12:36
Hallo Bernhard,


Soll ich mal nachfragen, ob ich das Projekt komplett an Deinen Account abgeben kann?
Kannst du machen, musst du aber nicht :)


Für den symplektischen Integrator braucht man so eine Funktionsschnittstelle dann zweimal mit unterschiedlicher Funktionalität und unterschiedlichem Wert für die Schrittweite.
Danke für den Tip, der SI kommt wenn ich mein Buch durchgearbeitet habe ;)

Kibo
20.05.2013, 13:36
Update auf 1.11

Features:

Man kann festlegen wie lange die Simulation laufen soll

beseitigte Bugs:

Fehlerhafte Ausgabe der simulierten Zeit
Geschwindigkeitsslider entfernt das zu großer Störfaktor
Stabileres Management der settings.ini


Sourceforge (https://sourceforge.net/projects/smallplanetarys/files/latest/download?source=files)

Kibo
23.05.2013, 14:00
Zu meinem großen Erstaunen und unter Zuhilfenahme einiger auf Youtube (http://www.youtube.com/watch?v=KelGNYzkarU) zur Verfügung gestellter Seminare eines gewissen Jörn Loviscach, Prof. für Ingenieurmathematik und technische Informatik, habe ich gerade festgestellt, dass das in meinem Programm implementierte Eulerverfahren bereits symplektisch ist.

Anscheinend hab ich das wohl unwissentlich durch Trial-and-Error so implementiert.
Eine kleine Umänderung der verwendeten Variablen hin zum normalen expliziten Euler führte erwartungsgemäß zu wesentlich instabileren Bahnverläufen.:)

Bernhard
25.05.2013, 13:38
Hallo Kibo,


Zu meinem großen Erstaunen und unter Zuhilfenahme einiger auf Youtube (http://www.youtube.com/watch?v=KelGNYzkarU) zur Verfügung gestellter Seminare eines gewissen Jörn Loviscach, Prof. für Ingenieurmathematik und technische Informatik, habe ich gerade festgestellt, dass das in meinem Programm implementierte Eulerverfahren bereits symplektisch ist.
Tolles Video. In der Wikipedia wird das gleiche hier: http://en.wikipedia.org/wiki/Symplectic_Euler_method beschrieben.

Ich habe die Dateien auf sourceforge.net neu sortiert und die aktuellen Parameterdateien neu hochgeladen. Den aktuellen Code habe ich so in eine zip-Datei gepackt, dass beim Entpacken automatisch ein neuer Ordner angelegt wird. Der Unterordner ist wichtig, falls man in einen Ordner entpackt, der bereits viele Dateien enthält.

Zusätzlich hätte ich noch eine Idee für eine neue Version. Beim Drücken auf den Start-Button sollte das Raster gleich so eingeblendet werden, dass man die Bahnen sofort sieht, ohne dass man die Kamera mit der Maus erst positionieren muss. Die Position der Kamera wird über die Datei settings.ini sowieso schon vorgegeben. Ein erneutes Positionieren mit der Maus stört deswegen.
MfG

Kibo
25.05.2013, 17:48
Hi Bernhard,

Dein Vorschlag ist gut, ich arbeite im Moment am Leapfrog. Wird noch etwas dauern bis das richtig läuft denk ich.
In der dann neuen Version wird's auch den nicht symplektischen Euler zum Vergleich geben.

mfg

Bernhard
25.05.2013, 18:48
In der dann neuen Version wird's auch den nicht symplektischen Euler zum Vergleich geben.
Interessant ist bei der Simulation der vier Sonnen ein Vergleich zwischen Runge-Kutta und Euler bei einer Schrittweite von 86400 Sekunden (=1Tag). Runge-Kutta wird da sehr schnell instabil, während das symplektische Eulerverfahren stabil bleibt.

Bernhard
26.05.2013, 08:36
ich arbeite im Moment am Leapfrog.
Der ist bei gleicher Genauigkeit im Vergleich zum Bulirsch-Stoer-Verfahren fast schon aberwitzig schnell. Ich teste das momentan an den Mondephemeriden und mache dabei Berechnungen über ein Jahr. Bei einer Schrittweite von 400 s braucht das Verlet-Verfahren (=Leapfrog) 8s Rechenzeit, während ich Bulirsch-Stoer nach über zwei Minuten abgebrochen habe. Bei einer Schrittweite von 86400s braucht BS ca. 4s, während man Verlet mit einer separaten Armbanduhr nicht mehr messen kann :) .

Nach welcher Anleitung gehst Du vor? Ich würde da nach wie vor die englische Wikipedia-Seite mit k=2 empfehlen. Es ist übrigens ganz interessant die dortige Anleitung auch für k=1 etwas genauer anzusehen. Man erhält das symplektische Euler-Verfahren aus dem von Dir verlinkten Clip nämlich nur mit einem (naheliegenden) kleinen Trick.
MfG

Kibo
26.05.2013, 09:07
Morgen Bernhard,

Schau mal hier:
http://www.youtube.com/watch?v=7VmYd9tJ68g
und weiter
http://www.youtube.com/watch?v=rbpCP89mgjI

dazu nehm ich noch Wikipedia

mfg

Bernhard
27.05.2013, 21:34
Bei einer Schrittweite von 86400s braucht BS ca. 4s, während man Verlet mit einer separaten Armbanduhr nicht mehr messen kann :) .
Hallo Kibo,

nachdem jetzt der lästige Fehler bei der Umrechnung in RA und DEC geklärt ist, habe ich nochmal das Verlet-Verfahren gegen BS antreten lassen und gesehen, dass Verlet bei großen Schrittweiten wie 86400s sehr ungenau wird. Bei Simulationen über ein Jahr muss man also wissen, welche Schrittweite am besten zu wählen ist. Dann sind beide Verfahren in etwa gleich schnell. Bei einer Schrittweite von 1000s bekommt man mit dem Leapfrog über ein Jahr sehr gute Ergebnisse.
MfG

Bernhard
01.06.2013, 15:17
Hallo Kibo,

Du wolltest weitere Literatur? Schau mal hier: http://lmgtfy.com/?q=symplektische+mechanik (Den Rahmen habe ich wegen einer persönlichen Vorliebe gewählt und soll keine Kritik bedeuten). Der dritte Link von L. Haug und A. Klaiber ist ganz interessant und hilfreich, wenn man das 1983er Paper von R. Ruth besser verstehen will: http://adsabs.harvard.edu/abs/1983ITNS...30.2669R
MfG

Kibo
02.06.2013, 08:28
Also eine mechanische Berechnung ist allgemein symplektisch wenn gilt:

E_pot + E_kin =konstant über den gesamten Ablauf der Simulation.

E_kin lässt sich wohl leicht ausrechnen über: E_kin=1/2m*v^2
Für E_pot habe ich beispielsweise die Formel gefunden: E_pot= -G*m*M/r

Ist das soweit erst einmal richtig?

mfg

Bernhard
02.06.2013, 10:07
Ist das soweit erst einmal richtig?
Hallo Kibo,

im Rahmen des aktuellen Themas ist das korrekt. Die vollständige Theorie ist aber allgemeiner aufgebaut. Wichtige Begriffe in diesem Zusammenhang sind beispielsweise:
http://de.wikipedia.org/wiki/Satz_von_Liouville_%28Physik%29
http://de.wikipedia.org/wiki/Phasenraum
MfG

EDIT: In dem von Dir verlinkten ersten Youtube-Clip wurde allerdings erwähnt, dass der symplektische Integrator die Gesamtenergie auch nur näherungsweise (und wesentlich besser als Euler) und nicht exakt konstant hält. Ein perfekter sympl. Integrator (Ordnung unendlich oder so ähnlich) sollte das aber prinzipiell machen.

Kibo
11.06.2013, 10:27
Hallo,

Ich möchte, dass mein Programm die Sternfarbe anhand seiner Masse anpasst (leichte rote Sterne, schwere blaue Sterne).

Zur Bestimmung der Sternenfarbe habe ich die sehr einfache Formel

Sternfarbe in Nanometern 3000000/temp in Kelvin

gefunden.

Ich konnte jedoch keine Formel zur Berechnung der Temperatur aus der Masse finden. Hat da jemand vielleicht etwas zur Hand?

mfg

Bernhard
11.06.2013, 21:47
Hallo Kibo,

das Hertzsprung-Russel-Diagramm verkoppelt zumindest die Farbe mit der absoluten Helligkeit. Vielleicht hilft das ja ein wenig weiter?
MfG

Kibo
12.06.2013, 08:10
Hallo Bernhard,

sowas hatte ich befürchtet, dann schätze ich eben eine Funktion anhand des HRD ab.

mfg und danke

joeydee
12.06.2013, 09:59
Erstmal zu deiner gefundenen Farbenformel: damit wäre die Sonne bei 3000000/5800=517nm grün?!
Es ist leider ein wenig komplexer: http://en.wikipedia.org/wiki/Planckian_locus#Approximation
Umrechnungen von LAB in RGB solltest du auch irgendwo finden.

Spektraltyp/Masse: Vielleicht reichen dir auch schon die Eckdaten aus dieser Tabelle: http://de.wikipedia.org/wiki/Spektralklasse#Einteilung - hier werden typische Massen für bestimmte Spektraltypen der Hauptreihensterne angegeben. Daraus kannst du eigentlich eine Wertetabelle zum Interpolieren basteln.

Kibo
12.06.2013, 19:10
Hallo Joey,

Ich finde das Ergebnis grün doch sehr treffend. Grün ist ja so ziemlich genau der Durchschnitt, des von der Sonne ausgestrahlten Lichts.
Sieht man hier (http://www.google.de/imgres?imgurl=http://www.helmholtz-berlin.de/media/media/oea/bildarchiv/forschung/solar/sonnenspektrum.jpg&imgrefurl=http://www.helmholtz-berlin.de/aktuell/pr/mediathek/bildarchiv1/forschung/solarenergie_de.html&h=529&w=800&sz=223&tbnid=fgXzeqdxj1SkWM:&tbnh=79&tbnw=120&zoom=1&usg=__wk6MT1lVHXs-iNPj0dAjEpRFskw=&docid=oKxZsFO3_1VU_M&sa=X&ei=D6q4UeGnMMmltAb9qYD4Bg&ved=0CDcQ9QEwAQ&dur=430) ganz gut. Jedenfalls, danke für die Links, leider helfen die mir aber auch nicht weiter wie ich von der Masse des Sterns denn nun zu seiner Oberflächentemperatur komme.

mfg

Bernhard
12.06.2013, 21:08
leider helfen die mir aber auch nicht weiter wie ich von der Masse des Sterns denn nun zu seiner Oberflächentemperatur komme.
Hallo Kibo,

in der ersten Tabelle von hier: http://de.wikipedia.org/wiki/Spektralklasse#Einteilung findet man in Spalte 5 die durchschnittliche Masse jedes Hauptreihensterntyps, in Spalte 3 die zugehörige Farbe, wie man sie mit dem Auge wahrnimmt und in Spalte 4 die durchschnittliche Temperatur. Damit sind 99% aller Sterne berücksichtigt.
MfG

joeydee
13.06.2013, 08:18
Ich finde das Ergebnis grün doch sehr treffend. Grün ist ja so ziemlich genau der Durchschnitt, des von der Sonne ausgestrahlten Lichts.
Wie kommst du damit auf die tatsächliche Sternfarbe?

Kibo
13.06.2013, 16:09
Wie kommst du damit auf die tatsächliche Sternfarbe?

Noch gar nicht, eines nach dem anderen.

Jedenfalls hab ich für die Beziehung Masse/Temperatur eine ganz gute Näherungsfunktion abgeschätzt:
f(m)=(m^0,51)*4950 ;m in Sonnenmassen

mfg

edit: hab ich, wie von Bernhard vorgeschlagen aus der Wikipedia-Tabelle erstellt.

Kibo
13.06.2013, 17:13
hab die funktion noch ein bisschen angepasst:
f(m)=(m-0,1)^0,5*5000+440
passt mehr oder weniger

Kibo
13.06.2013, 21:29
Färbung funktioniert wie von mir gewünscht. Ich arbeite noch am Leapfrog, dann gibt es das nächstes Release.

UMa
14.06.2013, 08:51
Hallo Kibo,

die in dieser Wikipediatabelle angegebenen Massen sind die maximalen in diese Spektralklasse, also z.B. für F0, G0, ... und passen zu der oberen Temperaturangabe, nicht der mittleren oder unteren. Die Sonne hat z.b. 5772 K bei m=1.0. Bei deiner Formel kommt viel weniger raus.

Grüße UMa

Kibo
14.06.2013, 20:22
Hallo Uma,

Wenn du jetzt schon meckerst sei froh dass du meinen Algorithmus Temperatur zu Farbe nicht gesehen hast ;)
Mir reicht die erste grobe Näherung, meine Sonne leuchtet schön gelb.

mfg

mac
14.06.2013, 22:54
Hallo Kibo,

im Laufe der Jahre sammelt sich was an:
http://adsbit.harvard.edu/cgi-bin/nph-iarticle_query?1992ApJS...81..865S
http://www.vendian.org/mncharity/dir3/starcolor/

Herzliche Grüße

MAC

Kibo
15.06.2013, 08:39
Wen vielleicht meine aktuelle Umsetzung interessiert:

Stemp:=(masse/1.989e30-0,1)**0.5 *5000 +440
spektrumdurchschnitt:=3000000/Stemp
rgbr:=-((spektrumdurchschnitt-650)/300)^2+1 ;blauer Anteil
rgbg:=-((spektrumdurchschnitt-525)/175)^2+1 ;grüner Anteil
rgbb:=-((spektrumdurchschnitt-475)/150)^2+1 ;roter Anteil

Kibo
25.06.2013, 22:55
Es gibt ein neues, im wahrsten Sinne des Wortes, schönes Update.

Neu sind die Farbberechnung des Sterns, OpenGl Schattenberechnung, Versionskontrolle der settings.txt, mehr oder weniger richtige Anfangsausrichtung der Kamera und wie gesagt explizites Euler.
Immernoch nicht Fertig ist der Leapfrog Algorithmus (gerade keine Lust mehr auf Troubleshooting).
Kann sein das ich die ein oder andere Änderung vergessen habe, letztes Update ist ja lange her.

mfg

Bernhard
27.06.2013, 17:42
Hallo Kibo,

hier eine Version des Leapfrog als Pseudocode a la UMa:

dt=Schrittweite

P1ort=ort_alt
P1v =v_alt
P1ort=ort_alt + dt/2*P1v
m1a=Grav(P1ort)

ort_neu=ort_alt + dt*v_alt + 1/2*dt*dt*m1a
v_neu =v_alt + dt*m1a

Kibo
30.06.2013, 14:51
Hallo Bernhard,



dt=Schrittweite

P1ort=ort_alt
P1v =v_alt
P1ort=ort_alt + dt/2*P1v
m1a=Grav(P1ort)

ort_neu=ort_alt + dt*v_alt + 1/2*dt*dt*m1a
v_neu =v_alt + dt*m1a

Bist du sicher das das so stimmt?

mfg

Bernhard
30.06.2013, 15:14
Bist du sicher das das so stimmt?
Hallo Kibo,

ich weiß, dass die Differenzenschemen in der Wikipedia etwas anders aussehen, aber ich bin mir sicher, dass die angeschriebene Variante zweiter Ordung ist. Es ist auch klar, dass es egal ist, ob die x-Variable den "Bock" bildet oder die v-Variable, denn es handelt sich insgesamt nur um einen zeitlichen Versatz von dt/2.

Kurz: Ich kann zwar nicht beweisen, dass dieses Schema symplektisch genannt werden darf, bin mir aber dessen ziemlich sicher, weil es eben aus dem englischen Artikel zum symplektischen Integrator stammt. Ich habe es in C++ bereits implementiert und damit in kürzester Rechenzeit sehr genaue Mond-Ephemeriden berechnen können. Zusammen mit der Richardson-Extrapolation ist dieses Verfahren bei langen Simulationszeiten schneller und genauer als Runge-Kutta und schneller als Bulirsch-Stoer. Nach dem gleichen Artikel habe ich mittlerweile auch das Verfahren dritter Ordnung implementiert und erfolgreich getestet.

In Deinem Programm sollte dieses Verfahren also eine deutliche Geschwindigkeitsverbesserung bei mindestens gleicher bis besserer Stabilität und Genauigkeit bringen.

Ich hätte auch noch einen Vorschlag für die Grafik des kleinen planetaren Simulators. Bei der Kameraposition würde ich bei der großen Zahl (z-Komponente?) statt einer 1 ganz vorne eine 5 verwenden, also fast die fünffache Entfernung. Für den Screenshot auf sourceforge habe ich genau diese Änderung verwendet. Man bekommt dann die Ansicht bis zur Marsbahn inklusive, was wesentlich netter aussieht als der "einsame" Merkur.
Gruß

Kibo
30.06.2013, 18:48
Dein Tip mal in die englische Wikipedia vom Leapfrog zu schauen hat's gebracht.

Hab die Variante:
x_i := x_(i-1) + v_(i-1/2) * delta_t
a_i := F(x_i)
v_(i+1/2) := v_(i-1/2) + a_i * delta_t

implementiert. Verglichen mit dem ewigen Rumprobieren mit dem was in der deutschen Wikipedia steht, ging das jetzt wirklich sehr schnell:)

Die 1.13 liegt dann wie immer auf sourceforge (link weiterhin in Hilfe zu finden)
Hab, wie von dir vorgeschlagen die Standart Kameraposition angepasst. Kann man wie man möchte weiterhin in der settings.txt nach eigenen Wünschen ändern.

mfg

Kibo
01.07.2013, 19:17
Noch ein kleines Update fürs Auge.

Ab Version 1.14 kann man wieder Zoomen. Und damit der Zoom nicht zu hässlich aussieht habe ich die Licht-Berechnung nochmal ein bisschen aufgehübscht. Hier mal ein Screenshot
https://a.fsdn.com/con/app/proj/smallplanetarys/screenshots/Screenshot%20sps3.png

Desweiteren gibt es jetzt eine funktionierende Ausgabe der Geschwindigkeit eines Körpers, und der Abstände der Körper zueinander.
Wie man sieht, habe ich hier eine andere Textur verwendet. Man kann die Texturen jederzeit ändern. Man muss nur die Texturen in den img-Ordner packen und entsprechend benennen.

mfg

Bernhard
01.07.2013, 21:00
Ab Version 1.14 kann man wieder Zoomen.
Das Scrollrad ist bei mir (mit WinXP) immer noch ohne Funktion :confused: .

Kibo
01.07.2013, 21:10
Das Scrollrad ist seit 1.10.2 dazu da, die Geschwindigkeit der Kamera zu regeln. Schau mal bitte ob sich die Geschwindigkeitsanzeige unten links im Renderfenster ändert wenn du scrollst.

Edit: Ich vergaß zu erwähnen das die Tastenbelegung zum Zoomen auf BildAuf/BildAb gelegt ist. Ich werd das mal mal gleich in die Hilfeseite einpflegen.

Bernhard
01.07.2013, 22:05
Edit: Ich vergaß zu erwähnen das die Tastenbelegung zum Zoomen auf BildAuf/BildAb gelegt ist. Ich werd das mal mal gleich in die Hilfeseite einpflegen.
Dann passt alles :) .
MfG

Kibo
04.07.2013, 21:48
Nabend Bernhard,

Kannst du die Richardson Interpolation kurz erklären? Ich kann mit dem Wikipedia-Artikel irgendwie nicht viel anfangen.

mfg

Bernhard
04.07.2013, 22:27
Kannst du die Richardson Interpolation kurz erklären?
Hallo Kibo,

kennt man die Ordnung eines Differenzenverfahrens so kennt man per Definition die Taylorentwicklung des Restglieds (=Fehler). Bei einem Verfahren zweiter Ordnung gilt beispielsweise y_exakt - y_numerisch = c_2 * dt^2 + c_3 * dt^3 + c_4 * dt^4 ....Vernachlässigt man die Koeffizienten c_3, c4, usw. erhält man eine Gleichung mit den zwei Unbekannten c_2 und y_exakt. Macht man den gleichen Rechenschritt mit zwei verschiedenen Schrittweiten bekommt zwei Gleichungen mit zwei Unbekannten und kann damit die zwei Unbekannten c_2 und y_exakt ausrechnen und erhält insgesamt den exakten Wert unter Vernachlässigung der Koeffizienten c_3, c_4 bis c_unendlich.

Grob erklärt ist das dann die Richardson-Interpolation.

Gegenfrage: Welche neuen Ziele verfolgst Du mit dem Programm, denn es geht mehr und mehr in Richtung Ephemeridenrechnung. Falls das Programm in dieser Richtung weiterentwickelt werden soll, könntest Du aus dem Nachbarprojekt http://sourceforge.net/projects/nbodysimulator/ die Umrechnung von den euklidischen Koordinaten in Rektaszension und Deklination verwenden und nach ahk übersetzen. Kurz, die grafische Darstellung zeigt beispielsweise zwischen dem Leapfrog und Runge-Kutta keinen Unterschied mehr. Der Unterschied wird allerdings dann wichtig, wenn man mit Teleskopen die genauen Positionen der berechneten Himmelskörper nachprüft. Man nennt das dann Astrometrie und kann das auch mit kleineren Teleskopen hobbymäßig betreiben.

Kibo
05.07.2013, 23:03
Hallo Bernhard, danke für die Antwort!

Zu deiner Frage: An sich entscheide ich das ganz spontan was ich noch einbaue (manche Sachen, verwerfe ich auch wieder beim Programmieren). Das Hauptziel ist seid Implementation der ersten wirklich genauen Algorithmen schon erreicht.
Ich plane gerade ein etwas größeres Update auf Planeten 2.0. Ich möchte die Modellierungsfunktion stark aufwerten.
Gerade fertig gestellt, habe ich einen Vorhersagemodus. Mit dem kann man sehen , wie sich Veränderungen an den Bahnbewegungen der Planeten auswirken.
Unter Umständen baue ich auch noch eine Planeten Drag&Drop Funktion ein. Dann kann man sich sein Sonnensystem recht schnell zurecht basteln wie man möchte und so simulieren. Ich habe mir auch vorgenommen die Einlesefunktion zu erweitern, sodass sie auch offiziellere Dateiformate liest, das ist sicher ganz praktisch wenn man direkt die aus dem Internet heruntergeladenen Dateien verwenden kann.
Deklination und Rektazension mit zu berechnen ist sicher keine schlechte Idee, ich denke das lässt sich recht fix realisieren.

mfg

Kibo
02.11.2013, 17:47
Das nächste Update steht zum Download bereit,

es ist allerdings nicht die 2.0 sondern nur 1.15. Das drag and drop ist noch nicht ganz ausgereift und das Äquatoriale Koordinatensystem ist noch nicht drinne, ihr könnt das ganze quasi als Beta betrachten.

Kurz zur Erklärung:
Die Planetendatei wurde leicht abgeändert, in der ersten Zeile stehen zusätzlich die Dateiversion und das Datum der Planetendaten.
Mit dem Datum kann man auch was anstellen, im Simulationsfenster kann man jetzt einstellen wie lange er simulieren soll, also beispielsweise 2000 Jahre oder bis zum 20.10.2019. Es gibt jetzt einen Vorhersagemodus, wenn man an bestimmten Planetendaten herumprobieren möchte (Beispielsweise spekulierte Asteroiden Flugbahnen), und wie gesagt gibt es jetzt auch eine Drag and Drop Funktion (noch sehr grob).

Downloadlink wie immer Sourceforge (https://sourceforge.net/projects/smallplanetarys/files/latest/download?source=files)

mfg

Bernhard
14.05.2014, 07:35
Ich möchte hier die interessanten Startparamater des Sitnikov-Problems festhalten:

KoordX KoordY Koordz v_X v_Y v_Z Masse Name dichte /1.15/20130721115300
-7.48e10 0.0 0.0 0.0 21068.5010849826 0.0 1.99e30 Sonne1 1.408
7.48e10 0.0 0.0 0.0 -21068.5010849826 0.0 1.99e30 Sonne2 1.408
0.0 0.0 1.0e4 0.0 0.0 0.0 1.0e4 Sitnikovia 1.408

1) Man lade sich von sourceforge die Datei planeten1.15.exe und starte diese einmal, um die Datei settings.ini zu erzeugen.
2) Nach dem Schließen des Programms ersetze man in der ebenfalls erzeugten Datei planeten.txt alle Eingaben durch die oben angegebenen Werte.
3) planeten1.15.exe erneut starten, Start-Button anklicken und unmittelbar danach mit der Maus gleich auf das kleine weiße Kreuz in der Bildmitte klicken.
4) Die Kreisbahnen der beiden Sonnen durch vorsichtiges Bewegen der Maus im Bild zentrieren und per erneutem Mausklick in der Position fixieren.

5) Beobachten, wie Sitnikovia nach mehreren Umläufen der Sonnen die Position in der Mitte des Bildes verlässt und später aus dem System katapultiert wird.

yetiholz
14.05.2014, 10:04
planeten1.15.exe
download:
https://sourceforge.net/projects/smallplanetarys/files/latest/download?source=files

(Ich hab einfach den Link von Kibo genommen.)
(Da brauchte ich nicht suchen.) :)

yetiholz
14.05.2014, 10:10
https://sourceforge.net/projects/smallplanetarys/files/latest/download?source=files

der editor von astronews hat da von selbst "..." ersetzt anstelle von:
... = "allplanetarys/files/latest/downloa"
(smallplanetarys/files/latest/download)
ich versuchs gleich nochmal den link von Kibo auszuschreiben

https://sourceforge.net/projects/smallplanetarys/files/latest/download?source=files

yetiholz
14.05.2014, 10:50
Hallo Bernhard, Hallo Kibo,
ja es hat gleich funktioniert mit der Anleitung von S. 19! :)
(super interessant)
(ich hab's mir gleich nochmal angesehen und bei der 2. Simulation ist Sitnikovia nach dem Katapultieren für einen kurzen Moment wieder in die Mitte und dann wieder raus)
(auf den Start button musste ich im Programm noch klicken, da hab ich eine Weile dafür gebraucht) :)
(bei der 3. Simulation ist Sitnikovia von der 1. Sonne zur 2. Sonne katapultiert worden) :)

Bernhard
14.05.2014, 11:06
(auf den Start button musste ich im Programm noch klicken, da hab ich eine Weile dafür gebraucht) :)
Habe den Punkt 3 entsprechend erweitert. Danke für den Hinweis.

Kibo
14.05.2014, 12:16
Hallo Yeti,

Die Instabilität dieses Systems lässt sich in dem Programm natürlich auf Rundungsfehler zurückführen. In der Realität wäre so ein System wohl stabiler, die Bahn des Planeten schwankt aber auch da unvorhersehbar. Auf Sitnikovia möchte ich nun wirklich nicht leben:
Jahrzehntelange eiskalte Winter, gefolgt von ebenso langen brütend heißen Sommern oder auch Jahreszeiten von nur wenigen Wochen, nahezu alles ist jederzeit möglich. Freut mich jedenfalls, dass dir mein Programm gefällt.

mfg

yetiholz
14.05.2014, 13:56
ja super :)
(jetzt bei der xten Simulation ist Sitnikovia von der 1. zur 2. Sonne und dann wieder zur 1. Sonne und dann raus)
(ich finde dass immer wieder total spannend) :)

Bernhard
14.05.2014, 14:44
(ich finde dass immer wieder total spannend) :)
Du kannst dann auch auch noch etwas mit den Startbedingungen in der Datei planeten.txt herumspielen. Welcher Wert was bedeutet ist ja in der Kommentarzeile ganz am Anfang beschrieben.

yetiholz
14.05.2014, 14:59
ich hab noch eine Sonne (3.) hinzugefügt :)
(großartig)
(leider hauts dann Sitnikovia raus, da mach ich noch was) :)

SRMeister
14.05.2014, 15:10
Ich schätze, mit 3 Sonnen wird das nix, aber mit 4 Sonnen könnte das wieder gehen :)

yetiholz
15.05.2014, 00:20
@ SRMeister
ja, mit 4 Sonnen geht's. :)
hier die Koordinaten:
KoordX KoordY Koordz v_X v_Y v_Z Masse Name dichte /1.15/20130721115300
-7.48e10 0.0 0.0 0.0 42137.0021699652 0.0 1.99e30 Sonne1 1.408
7.48e10 0.0 0.0 0.0 -42137.0021699652 0.0 1.99e30 Sonne2 1.408
0.0 -7.48e10 0.0 -42137.0021699652 0.0 0.0 1.99e30 Sonne3 1.408
0.0 7.48e10 0.0 42137.0021699652 0.0 0.0 1.99e30 Sonne4 1.408
0.0 0.0 1.0e4 0.0 0.0 0.0 1.0e4 Sitnikovia 1.408

(Wen es interessiert, einfach Koordinaten in Planeten.txt einfügen)
Programm von Kibo:
https://sourceforge.net/projects/smallplanetarys/files/latest/download?source=files
(super Anleitung hier auf S.19)

yetiholz
15.05.2014, 02:24
Hallo Bernhard, Hallo SRMeister
ich probierte noch ein bischen mit 3 Sonnen:
(ohne dass Sitnikovia gleich rausgeschleudert wird!) :)
(wobei ich Sitnikovia so schwer wie eine Zwergsonne machte)

KoordX KoordY Koordz v_X v_Y v_Z Masse Name dichte /1.15/20130721115300
-7.48e10 0.0 0.0 0.0 41068 0.0 0.99e28 Sonne1 1.408
7.48e10 0.0 0.0 0.0 -41068 0.0 0.99e28 Sonne2 1.408
0.0 0.0 0.0 0.0 0.0 0.0 1.99e30 Sonne3 1.408
2.0e10 0.0 0.0 0.0 68000 0.0 1.0e28 Sitnikovia 1.408

(Wen's interessiert, einfach die Parameter in Planeten.txt einfügen)
Programm von Kibo:
https://sourceforge.net/projects/smallplanetarys/files/latest/download?source=files
(super Anleitung hier auf S.19 von Bernhard)

Bernhard
15.05.2014, 09:28
ja, mit 4 Sonnen geht's.
hier die Koordinaten:
KoordX KoordY Koordz v_X v_Y v_Z Masse Name dichte /1.15/20130721115300
-7.48e10 0.0 0.0 0.0 42137.0021699652 0.0 1.99e30 Sonne1 1.408
....
Diese Startbedingungen führen auf relativ schlecht zentrierte Bahnen der Sonnen. Der Wert für die Startgeschwindigkeit kann aber gemäß:

v² = G * m / r * (1/sqrt(2) + 0.25) berechnet werden (Pythagoras).

Für r = 7.48e10, m=1.99e30 und G aus der Wikipedia ergibt das den Wert 41223.402162071053829891025355429
MfG

Bernhard
15.05.2014, 09:37
ich probierte noch ein bischen mit 3 Sonnen:
Hallo Yeti,

mit den angegebenen Startwerten steigt das Programm bei mir aus. Die Startbedingungen von Sonne3 sind zudem etwas auffällig. Geschwindigkeit Null und Start im Ursprung des Koordinatensystems macht das Problem sehr unsymmetrisch.

Zur Berechnung passender Werte benötigt man erneut etwas Trigonometrie und die Formel für die Zentrifugalkraft F_z = m * v² / r.
MfG

yetiholz
15.05.2014, 19:18
Hallo Bernhard,
ich hab gleich die von mir angegebenen Koordinaten nochmal in die Planeten.txt Datei eingefügt.
Das Programm läuft bei mir super. :)
(Vielleicht hatte Kibo eine verbesserte Version hochgeladen, die ich derzeit hab.)
Bisher hab ich nur probiert. Für die Berechnung: alle 3 gleichen Sonnen auf die selbe Kreisbahn zu bringen,
brauch ich etwas Zeit.
Vielen Dank für die Formeln. :)
(Damit hab ich's viel einfacher.)
v=41223... für 4 Sonnen probier ich gleich mal aus.
dankeschön :)

Kibo
15.05.2014, 20:38
Hallo,

Ich habe jetzt die Testdaten mit 3 Sonnen auch probiert und hatte keine Probleme, blos von Sitnikovschen Verhalten habe ich leider nichts mehr sehen können. Nicht verwunderlich bei den Koordinaten :).

Kibo
15.05.2014, 21:03
Wenn dir so etwas gefällt, kannst du ja ein bisschen mit Mars, Phobos und Deimos (http://www.astronews.com/forum/showthread.php?7531-Wie-soll-Mars-seine-2-Monde-eingefangen-haben&p=102737#post102737) rumprobieren. Alle Planeten des Sonnensystems findest du in der Standartdatei.

mfg

Bernhard
15.05.2014, 21:48
Freut mich jedenfalls, dass dir mein Programm gefällt.
Das Display macht schon was her.


Ich habe jetzt die Testdaten mit 3 Sonnen auch probiert und hatte keine Probleme, blos von Sitnikovschen Verhalten habe ich leider nichts mehr sehen können. Nicht verwunderlich bei den Koordinaten :).
Habe es gerade noch mal ohne Probleme ausprobiert. Das Programm hatte sich beim letzten Mal wohl irgendwie "verschluckt" - ist egal. Mir ist jetzt auch klar, wie das mit den 3 Sonnen funktionieren soll. Sonne3 sitzt unbeweglich im Schwerpunkt des Systems. Sitnikovia stürzt im freien Fall auf die Sonne3. Um dieses Szenario sollen sich Sonne1 und Sonne2 auf einer Kreisbahn drehen.
MfG

PS: Hier die Datei für 3 Sonnen und eine etwas zentrischere Bahn der Sonnen 1+2:

KoordX KoordY Koordz v_X v_Y v_Z Masse Name dichte /1.15/20130721115300
-7.48e10 0.0 0.0 0.0 47110.60061004 0.0 1.99e30 Sonne1 1.408
7.48e10 0.0 0.0 0.0 -47110.60061004 0.0 1.99e30 Sonne2 1.408
0.0 0.0 0.0 0.0 0.0 0.0 1.99e30 Sonne3 1.408
2.0e10 0.0 0.0 0.0 68000 0.0 1.0e4 Sitnikovia 1.408

Die Masse von Sitnikovia habe ich wieder auf Asteroidengröße reduziert. Läßt man das Programm länger laufen, verlässt Sonne3 den Mittelpunkt und kreist dann um Sonne1. Kurzfristig kehrt Sonne3 dann immer wieder mal zum Mittelpunkt zurück.
Interessante Simulation :) .

Bernhard
16.05.2014, 10:40
Läßt man das Programm länger laufen, verlässt Sonne3 den Mittelpunkt und kreist dann um Sonne1. Kurzfristig kehrt Sonne3 dann immer wieder mal zum Mittelpunkt zurück.
Interessante Simulation :) .
Reduziert man die vom Programm vorgegebene Schrittweite etwa um 50% zeigen sich weitere interessante Details dieser Konstellation. Sitnikovia kreist zuerst mit starker Periheldrehung in dem von Sonne3 gebildeten Potentialtopf. Nach vielen Umläufen von Sonne1+2 verlässt dann Sonne3 erneut den Koordinatenursprung und Sitnikovia wird wie bei einem Hammerwerfer von zwei Sonnen aus dem System herausgeschleudert.
MfG

yetiholz
26.05.2014, 20:58
hallo allerseits,
hier die Parameter für 3 Sonnen auf einem Kreis:
(wen's interessiert, einfach die Parameter in Planeten.txt einfügen)
KoordX KoordY Koordz v_X v_Y v_Z Masse Name dichte /1.15/20130721115300
0.0 7.48e10 0.0 32116.4841 0.0 0.0 1.99e30 Sonne1 1.408
64778700203.0760107779 -3.74e10 0.0 -15844 -27690 0.0 1.99e30 Sonne2 1.408
-64778700203.0760107779 -3.74e10 0.0 -15844 27690 0.0 1.99e30 Sonne3 1.408
0.0 0.0 0.0 -600.0 0.0 0.0 1.0e4 Sitnikovia 1.408
(die Anfangspunkte müssten auf einem Kreis sein)
(leider sind die Vektoren nicht exakt, sondern nur über sukzessive iterative Approximation genähert)
(ich hab zwar einen konkreten Rechenweg für die Richtung der Vektoren, aber da muss ich die Wurzel von einer -Zahl ziehen)

hier die von Bernhard berichtigten Parameter für 3 Sonnen mit perfekter Kreisbahn:
(Schrittweite vor dem Start im Programm ~halbieren!)
KoordX KoordY Koordz v_X v_Y v_Z Masse Name dichte /1.15/20130721115300
-7.48e10 0.0 0.0 0.0 47110.60061004 0.0 1.99e30 Sonne1 1.408
7.48e10 0.0 0.0 0.0 -47110.60061004 0.0 1.99e30 Sonne2 1.408
0.0 0.0 0.0 0.0 0.0 0.0 1.99e30 Sonne3 1.408
2.0e10 0.0 0.0 0.0 68000 0.0 1.0e4 Sitnikovia 1.408

hier die Parameter für 4 Sonnen auf einem Kreis:
mit den richtigen Werten von Bernhard
(wen's interessiert, einfach die Parameter in Planeten.txt einfügen)
(lief bei mir ~800 Tage stabil)
KoordX KoordY Koordz v_X v_Y v_Z Masse Name dichte /1.15/20130721115300
-7.48e10 0.0 0.0 0.0 41223.4021620710 0.0 1.99e30 Sonne1 1.408
7.48e10 0.0 0.0 0.0 -41223.4021620710 0.0 1.99e30 Sonne2 1.408
0.0 -7.48e10 0.0 -41223.4021620710 0.0 0.0 1.99e30 Sonne3 1.408
0.0 7.48e10 0.0 41223.4021620710 0.0 0.0 1.99e30 Sonne4 1.408
0.0 0.0 1.0e4 0.0 0.0 0.0 1.0e4 Sitnikovia 1.408

(und ein kleines Horrorzenario :)
Schrittweite vor Start 2mal halbieren)
KoordX KoordY Koordz ImpulsX ImpulsY ImpulsZ Masse Name dichte /1.15/20130721115300
-6.080091705285974E+08 2.454333708238292E+07 2.193321432139981E+06 2.809726412820165E+00 -1.036343441959919E+01 -4.204860196381065E-02 1.9897E+30 Sun 1.4
-5.331602529851193E+10 1.424323914756870E+10 6.000210156549231E+9 -2.277058581046753E+04 -4.495485791176564E+04 -1.582633849867808E+03 3.3032E23 Mercury 5.4
2.397967758924109E+10 -1.059717457755222E+11 -2.868972017005123E+09 3.388294521475002E+04 7.783094129350907E+03 -1.848596675496016E+03 4.8704E24 Venus 5.2
-1.474370019336122E+11 -2.789279589304088E+10 2.955826247744262E+06 5.079162483215040E+03 -2.939896910566963E+04 4.812093925554706E-01 5.976E24 Earth 5.5
-1.470514277349460E+11 -2.801509364690417E+10 3.730129039033502E+07 5.382156089456466E+03 -2.847733344864282E+04 2.235354455693894E+01 7.3505E22 Moon 3.3
2.035953497078640E+11 -3.456857769669849E+10 -5.736961786047079E+09 4.976094385716297E+03 2.594977520891352E+04 4.218104801061298E+02 6.421E23 Mars 3.9
7.115603707882032E+11 2.013825854433948E+11 -1.677091861778326E+10 -3.715952626588193E+03 1.319564076159422E+04 2.830947269851247E+01 1.8997E27 Jupiter 1.4
-1.396843425883871E+12 -3.381973875783179E+11 6.146826298318255E+10 1.753564063629373E+03 -9.410322369209542E+03 9.449519324100386E+01 5.6882E26 Saturn 1.4
3.003955040895567E+12 2.614531534364399E10 -3.882080718055668E+10 -1.089777236983593E2 6.492185531504123E+03 2.552972796165465E1 8.6875E25 Uranus 1.4
3.826621024919128E+12 -2.346728069113412E+12 -3.986079005702972E+10 2.805108135465661E+03 4.664771972182217E+03 -1.611737087776297E2 1.025E26 Neptune 1.4
4.599506779023092E+11 -4.750680975824932E+12 3.753056255249584E+11 5.510641913253943E+03 -5.401707214265152E-02 -1.536200853856072E+03 1.4717E22 Pluto 1.8
-7.54080091705285974E+10 2.454333708238292E+07 2.193321432139981E+0 0.0 47068.5010849826 0.0 1.99e30 Sonne1 1.408
7.41919908294714026E+10 2.454333708238292E+07 2.193321432139981E+0 0.0 -47068.5010849826 0.0 1.99e30 Sonne2 1.408

(Wen's interessiert, einfach die Parameter in Planeten.txt einfügen)
brillantes Programm von Kibo:
https://sourceforge.net/projects/smallplanetarys/files/latest/download?source=files
(super Anleitung von Bernhard hier auf S.19)

Bernhard
26.05.2014, 22:31
und ein kleines Horrorzenario :)
Hallo Yeti,

um das Sonnensystem zu zerstören reicht auch eine Sonne, die einmal quer durch das System geschickt wird. Interessant ist auch, dass Jupiter bei Deinem Vorschlag auf einer festen Position mit einer praktisch verschwindenden Geschwindigkeit landet. War das bei Deinem Programmlauf auch so?
MfG

yetiholz
27.05.2014, 02:25
Hallo Bernhard,
dieses Phänomen hatte ich vermutlich noch nicht. (Wie sieht man dass?)
Bei mir kreisen nach 9000 Tagen nur noch Saturn, Neptun & Pluto und alle anderen sind fast sofort weg. (mit 10000er Schrittweite, Leapfrog)
Mit den Standardwerten (50000er Schrittweite & Leapfrog) sind nach 4000 Tagen Uranus, Saturn, Neptun & Pluto noch drin und nach ~7000 Tagen hauts die 3 Sonnen raus.
mfG

Bernhard
27.05.2014, 06:00
(Wie sieht man dass?)
@Yeti: Einfach die Schrittweite auf ca. 25000 stellen und dann lange laufen lassen (bis schätzungsweise 5000 Tage). Ich hatte den ruhenden Jupiter auch bei einem zweiten Lauf mit leicht veränderter Schrittweite, ebenfalls bei etwa 25000.


Mit den Standardwerten (50000er Schrittweite & Leapfrog) sind nach 4000 Tagen Uranus, Saturn, Neptun & Pluto noch drin
Die Planeten entfernen sich bei dieser Einstellung bei mir alle sehr schnell weit aus dem System. Ich vergrößere dann die z-Koordinate der Kamera um den Faktor 10. Dann sieht man die Ränder des Koordinatengitters und einige schwerere Planeten wie Uranus, Saturn und Neptun kreisen im Außenbereich langsam um den Ursprung (=Schwerpunkt des Systems).


und nach ~7000 Tagen hauts die 3 Sonnen raus.
Bei mir auch, so etwa nach 7400 Tagen.

Bernhard
27.05.2014, 06:01
(wie kann ich meinen letzten Beitrag löschen?)
Das kann nur der Webadministrator. Lass' es einfach so stehen.

[Und der hat den Beitrag inzwischen gelöscht]

yetiholz
27.05.2014, 13:54
Hallo Bernhard,
jetzt sah ich den "ruhenden" Jupiter auch. :)
Die 25000er Schrittweite hat's gebracht. (einfach 10000 oder 30000 wie ich wollte ging nicht)
mfG



@Webadministrator
vielen Dank :)

yetiholz
09.06.2014, 15:16
Hallo,
Daten (Schrittweite runter auf 10000 vor Start) für 8 Sonnen & Sitnikovia: :)
KoordX KoordY Koordz v_X v_Y v_Z Masse Name dichte /1.15/20130721115300
-7.48e10 0.0 0.0 0.0 70570 0.0 1.99e30 Sonne1 1.408
-52891587232.753754825183158285443 52891587232.753754825183158285443 0.0 49900.525548334658796967586633739 49900.525548334658796967586633739 0.0 1.99e30 Sonne2 1.408
0.0 7.48e10 0.0 70570 0.0 0.0 1.99e30 Sonne3 1.408
52891587232.753754825183158285443 52891587232.753754825183158285443 0.0 49900.525548334658796967586633739 -49900.525548334658796967586633739 0.0 1.99e30 Sonne4 1.408
7.48e10 0.0 0.0 0.0 -70570 0.0 1.99e30 Sonne5 1.408
52891587232.753754825183158285443 -52891587232.753754825183158285443 0.0 -49900.525548334658796967586633739 -49900.525548334658796967586633739 0.0 1.99e30 Sonne6 1.408
0.0 -7.48e10 0.0 -70570 0.0 0.0 1.99e30 Sonne7 1.408
-52891587232.753754825183158285443 -52891587232.753754825183158285443 0.0 -49900.525548334658796967586633739 49900.525548334658796967586633739 0.0 1.99e30 Sonne8 1.408
0.0 0.0 0.0 0.0 0.0 0.0 1.0e4 Sitnikovia 1.408

yetiholz
11.06.2014, 14:45
Hallo, :)
hier Sitnikovia, mit 4 Sonnen von Bernhard aus Beitrag 53:
(nach 120 Tagen, sobald sich Sitnikovia anfängt zu bewegen muss man die Schrittweite runter auf 10000 stellen, sonst sieht man nichts)
KoordX KoordY Koordz v_X v_Y v_Z Masse Name dichte /1.15/20130721115300
-8.228e10 0.0 0.0 0.0 -96392.0 0.0 1.99e30 Sonne1 1.408
-6.732e10 0.0 0.0 0.0 36816.0 0.0 1.99e30 Sonne2 1.408
6.732e10 0.0 0.0 0.0 -36816.0 0.0 1.99e30 Sonne3 1.408
8.228e10 0.0 0.0 0.0 96392.0 0.0 1.99e30 Sonne4 1.408
0.0 0.0 0.0 0.0 0.0 0.0 1.0e4 Sitnikovia 1.408

yetiholz
11.06.2014, 15:04
(Wen's interessiert, einfach die Parameter in Planeten.txt einfügen)
brillantes Programm von Kibo:
https://sourceforge.net/projects/smallplanetarys/files/latest/download?source=files
(super Anleitung von Bernhard hier auf S.19)

http://sourceforge.net/projects/smallplanetarys/
http://sourceforge.net/projects/smallplanetarys/files/latest/download
https://sourceforge.net/projects/smallplanetarys/files/latest/download?source=files

AlphaCentauri
12.06.2014, 11:51
(Wen's interessiert, einfach die Parameter in Planeten.txt einfügen)
brillantes Programm von Kibo:
https://sourceforge.net/projects/smallplanetarys/files/latest/download?source=files
(super Anleitung von Bernhard hier auf S.19)

http://sourceforge.net/projects/smallplanetarys/
http://sourceforge.net/projects/smallplanetarys/files/latest/download
https://sourceforge.net/projects/smallplanetarys/files/latest/download?source=files

Danke für die Links! :)

yetiholz
22.06.2014, 01:30
@AlphaCentauri, gerne :)

hier noch eine Sitnikovia mit 3 Sonnen, abgeleitet von Bernhards 4 Sonnen von Beitrag 53:
(die Zahlen schrieb ich nach Gefühl, das Rechnen gab ich auf, wie man schon an den vorherigen Beispielen sah)

KoordX KoordY Koordz v_X v_Y v_Z Masse Name dichte /1.15/20130721115300
-8.228e10 0.0 0.0 0.0 -96392.0 0.0 1.99e30 Sonne1 1.408
-6.732e10 0.0 0.0 0.0 36816.0 0.0 1.99e30 Sonne2 1.408
7.48e10 0.0 0.0 0.0 29775.0 0.0 3.98e30 Sonne3 1.408
0.0 0.0 0.0 0.0 0.0 0.0 1.0e4 Sitnikovia 1.408

Bernhard
22.06.2014, 09:05
das Rechnen gab ich auf
Hallo Yeti,

dazu braucht man etwas Vektorrechnung. Ich skizziere mal kurz, wie man das Problem von n kreisenden Sonnen zerlegen kann.

1) Das Problem muss symmetrisch sein. D.h. die Sonnen müssen alle den gleichen Abstand voneinander haben. D.h. Der Winkel zwischen Sonne1, Ursprung, Sonne2 ist gleich 360°/n. Ferner sollen alle Sonnen die gleiche Masse haben.

2) Da das Problem symmetrisch ist, genügt es das Kräftegleichgewicht für eine Sonne auszurechnen

3) Die Zentrifugalkraft auf die herausgegriffene Sonne ist gleich mv²/r. r ist dabei der Abstand der Sonne zum Ursprung. Diese Formel liefert am Ende die gesuchte Verknüpfung zwischen den Startparametern r_x, v_x und m.

4) Die Gravitationskraft auf die herausgegriffene Sonne wird von den anderen n-1 Sonnen und der Geometrie des Problems festgelegt und kann per Vektoraddition berechnet werden. Die Abstände der herausgegriffene Sonne zu den restlichen Sonnen wird trigonometrisch berechnet, also mit Pythagoras und/oder den trigonometrischen Funktionen, wie sin, cos, tan usw.

5) Kennt man den einzelnen Abstand d zwischen der herausgegriffenen Sonne und der gerade betrachteten Sonne aus dem n-1-Pool, so kennt man über Gm²/d² auch die Stärke der einzelnen Gravitationskraft. Zusammen mit der Richtung dieser einzelnen Kraft kann man für jede der n-1 Sonnen den zugehörigen Vektor der Gravitationskraft ausrechnen.

6) Zuletzt müssen alle diese Vektoren addiert werden. Aufgrund der Symmetrie des Problems muss dieser Vektor genau in Richtung des Koordinatenursprung zeigen.

7) Jetzt muss der Betrag des resultierenden Vektors ausgerechnet werden und mit mv²/r gleichgesetzt werden. Fertig.

(8) Für n=2, 3 und 4 ist das Problem relativ übersichtlich und relativ leicht auszurechnen, wenn man sich genau an die angegebenen Schritte hält. Zusammen mit einer geeigneten Skizze sollte diese Aufgabe bereits ein interessierter Abiturient bewältigen können.

(9) Sitnikovia ist ein Testkörper, der in der hier beschriebenen Rechnung besser nichts zu suchen hat.
MfG

Bernhard
23.06.2014, 18:55
Hallo Yeti,

zu meiner Entlastung gebe ich hier das Ergebnis für drei kreisende Sonnen an:

KoordX KoordY Koordz v_X v_Y v_Z Masse Name dichte /1.15/20130721115300
7.48e10 0.0 0.0 0.0 32017.1979351182 0.0 1.99e30 Sonne1 1.408
-3.74e10 6.477870020308e10 0.0 -27727.7067698070344 -16008.5989675591 0.0 1.99e30 Sonne2 1.408
-3.74e10 -6.477870020308e10 0.0 27727.7067698070344 -16008.5989675591 0.0 1.99e30 Sonne3 1.408

Der Betrag der Geschwindgkeit jeder Sonne berechnet sich mit:

v = G * m / ( r * sqrt(3) )

Man muss bei dieser Aufgabe schon etwas rechnen, weil man immer wieder den 30° Winkel in der Startbedingung berücksichtigen muss.

Die Simulation ist aber auch nett anzusehen. Schrittweite als Empfehlung so 10000 bis 20000.
MfG

yetiholz
25.06.2014, 11:35
Hallo Bernhard, ich könnte mir dass stundenlang anschauen. :)

(Wen's interessiert, einfach die Parameter in Planeten.txt einfügen)
brillantes Programm von Kibo:
https://sourceforge.net/projects/smallplanetarys/files/latest/download?source=files
(super Anleitung von Bernhard hier auf S.19[Beitrag#190{Programm}] & S.22[Beitrag#220{Berechnung}])
http://sourceforge.net/projects/smallplanetarys/
http://sourceforge.net/projects/smallplanetarys/files/latest/download
https://sourceforge.net/projects/smallplanetarys/files/latest/download?source=files

Bernhard
25.06.2014, 14:55
Hallo Yeti,


Hallo Bernhard, ich könnte mir dass stundenlang anschauen. :)
ja, diese ursprünglich von SRMeister aufgeworfene Anregung ist wirklich ein netter Spezialfall des n-Körper-Problems. Ich finde es auch interessant, wie nach etlichen Umrundungen die numerischen Störungen das System verzerren. Die Sonnen vollführen dann recht ansprechende, elliptische Pirouetten.

Wer sich für die Mathematik dieses Spezialfalles interessiert und gerne rechnet, kann auch noch eine zusätzliche Sonne mit neuer Masse in den Ursprung setzen, oder (wie oben bereits erwähnt) einen Sitnikov-Planeten, der dann senkrecht zu der Ebene der Sonnen hin und her pendelt.
MfG

Bernhard
01.07.2014, 13:44
Hallo zusammen,

hier die Formel für 8 Sonnen auf einer Kreisbahn:

v² = G*m/r * ( 1/4 + 1 / sqrt(2) + sqrt( 2 + sqrt(2) ) )

Sie ist (logischerweise) eine Erweiterung der Formel für 4, bzw. 2 Sonnen.

Wer es sich in Kibos Programm ansehen will, braucht nur die folgenden Zeilen in die Datei Planeten.txt zu kopieren:



KoordX KoordY Koordz v_X v_Y v_Z Masse Name dichte
-7.48e10 0.0 0.0 0.0 70569.929193167 0.0 1.99e30 Sonne1 1.408
-52891587232.75375482518 52891587232.75375482518 0.0 49900.4754803429 49900.4754803429 0.0 1.99e30 Sonne2 1.408
0.0 7.48e10 0.0 70569.929193167 0.0 0.0 1.99e30 Sonne3 1.408
52891587232.75375482518 52891587232.75375482518 0.0 49900.4754803429 -49900.4754803429 0.0 1.99e30 Sonne4 1.408
7.48e10 0.0 0.0 0.0 -70569.929193167 0.0 1.99e30 Sonne5 1.408
52891587232.75375482518 -52891587232.75375482518 0.0 -49900.4754803429 -49900.4754803429 0.0 1.99e30 Sonne6 1.408
0.0 -7.48e10 0.0 -70569.929193167 0.0 0.0 1.99e30 Sonne7 1.408
-52891587232.75375482518 -52891587232.75375482518 0.0 -49900.4754803429 49900.4754803429 0.0 1.99e30 Sonne8 1.408
0.0 0.0 1.0e4 0.0 0.0 0.0 1.0e4 Sitnikov 1.408


BTW: Bei der Formel in Beitrag #221, also ganz oben, fehlt ein Quadrat. Es muss dort also v² = ... heißen.
MfG

yetiholz
03.07.2014, 18:18
hier eine kleine Spielerei mit 12 Sonnen:
(das Programm ist eigentlich viel zu gut, um damit zu spielen :) )
(falls es nicht funktioniert, einfach Sitnikovia weglassen)
(Schrittweite runter auf 10000, vor dem Start)

KoordX KoordY Koordz v_X v_Y v_Z Masse Name dichte
-14.96e10 0.0 0.0 0.0 80000 0.0 1.99e30 Sonne1 1.408
-105783174465.50750965036631657089 105783174465.50750965036631657089 0.0 56568.542494923801952067548968388 56568.542494923801952067548968388 0.0 1.99e30 Sonne2 1.408
0.0 14.96e10 0.0 80000 0.0 0.0 1.99e30 Sonne3 1.408
105783174465.50750965036631657089 105783174465.50750965036631657089 0.0 56568.542494923801952067548968388 -56568.542494923801952067548968388 0.0 1.99e30 Sonne4 1.408
14.96e10 0.0 0.0 0.0 -80000 0.0 1.99e30 Sonne5 1.408
105783174465.50750965036631657089 -105783174465.50750965036631657089 0.0 -56568.542494923801952067548968388 -56568.542494923801952067548968388 0.0 1.99e30 Sonne6 1.408
0.0 -14.96e10 0.0 -80000 0.0 0.0 1.99e30 Sonne7 1.408
-105783174465.50750965036631657089 -105783174465.50750965036631657089 0.0 -56568.542494923801952067548968388 56568.542494923801952067548968388 0.0 1.99e30 Sonne8 1.408
-4.48e10 0.0 0.0 0.0 -50000 0.0 1.99e30 Sonne01 1.408
4.48e10 0.0 0.0 0.0 50000 0.0 1.99e30 Sonne02 1.408
0.0 -4.48e10 0.0 50000 0.0 0.0 1.99e30 Sonne03 1.408
0.0 4.48e10 0.0 -50000 0.0 0.0 1.99e30 Sonne04 1.408
0.0 0.0 0.0 0.0 0.0 0.0 1.0e4 Sitnikovia 1.408

Bernhard
03.07.2014, 18:42
hier eine kleine Spielerei mit 12 Sonnen:
Läßt man die Anzahl der Sonnen gegen unendlich gehen, sollte man eine eindimensionale Massendichte, bzw. die Gesamtmasse des eindimensionalen Massenringes ableiten können :cool: .

yetiholz
03.07.2014, 20:55
Hallo Bernhard,
~2h brauchte ich, bis ich dass kapierte. :)
(durchgehenden, kontinuierlichen, ununterbrochenen Ring)
mit freundlichen Grüßen
yetiholz

Kibo
25.11.2014, 16:52
Hat jemand Koordinaten/Geschwindigkeitsvektoren für das Jupiter-Mond-System? bitte?

Bernhard
25.11.2014, 17:22
Hallo Kibo,

da bekommt man vermutlich die Bahnelemente der galileischen Monde etwas leichter und bekommt damit dann auch stabile Keplerbahnen. Das "Gefummel" mit den Horizons-Vektoren tue ich mir auf jeden Fall nicht an, und zwar auch deshalb, weil ich da neulich für Ralfs Thema ganz eigenartige Werte bekommen habe.
MfG

Kibo
25.11.2014, 17:32
Hallo Bernhard,

ja, auf der Horizons Seite war ich heute morgen. Nach 5 Minuten habe ich diese Möglichkeit dann auch ersteinmal verworfen. Wenn es hart auf hart kommt und keiner Bahndaten hat, werde ich wohl was programmieren müssen. Das wäre eigentlich ziemlich geil, ein Programm, das aus Ephemeriden stabile Bahndaten ableitet. (Oder besser "aufleitet" XD)

Bernhard
25.11.2014, 19:33
das aus Ephemeriden stabile Bahndaten ableitet. (Oder besser "aufleitet" XD)
Hallo Kibo,

die Ephemeriden sind die Zahlen, die man unmittelbar in eine Sternkarte eintragen kann, also Rektaszension und Deklination. Die wirklichen Eingabedaten sind eigentlich das Datum, sowie das Äquinoktium der Sternkarte, eventuell noch der Beobachtungsort und Zeitzone des Beobachters, je nachdem wieviel Arbeit man sich machen will.

Am Anfang ist man für gewöhnlich schon mal froh, wenn man aus den Bahnelementen, die Draufsicht auf Jupiter in Abhängigkeit von der Zeit korrekt ausrechnet und das ist auch nicht sooo schwer:
Beschreibung des Berechnungsverfahrens (http://de.wikibooks.org/wiki/Astronomische_Berechnungen_f%C3%BCr_Amateure/_Druckversion/_Himmelsmechanik)
Siehe auch: http://de.wikipedia.org/wiki/Bahnelement
sowie die Bahnelemente der vier Monde:
http://de.wikipedia.org/wiki/Io_%28Mond%29
http://de.wikipedia.org/wiki/Europa_%28Mond%29
http://de.wikipedia.org/wiki/Ganymed_%28Mond%29
http://de.wikipedia.org/wiki/Kallisto_%28Mond%29

Erste Aufgabe:
Berechne die wahre Anomalie aus der mittleren Anomalie per Iteration
Zweite Aufgabe:
Übertrage den Ort des Mondes in ein geeignetes Koordinatensystem (deutlich aufwendiger mit offenem Ende)

EDIT: Will man diese ganze Rechnung deutlich abkürzen, kann man die Bahnen vorerst auch als Kreisbahnen annehmen und die berechneten Positionen mit Cartes Du Ciel o.ä. vergleichen. Für kleinere Zeiträume sollte das gut funktionieren, da alle Exzentrizitäten der Bahnen relativ klein sind.

Kibo
25.11.2014, 22:27
OK bei der ersten aufgabe (EDIT: richtig??):o:

Berechnung von der exzentrischen Anomalie E(t)

iterative Nullstellenberechnung nach Newtonverfahren

t=1; t0=0
e=0,041 (für Io)
U=1,769 tage


f(E) = E - e *sin(E) - 2*pi*(t-t0)/U
f'(E) = 1-0,041*cos(E)

Newtonverfahren:

x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}

startwert:E=1

3,553063682
3,55394685
3,553942082
3,553942108
3,553942108
3,553942108
3,553942108
3,553942108
3,553942108
3,553942108
3,553942108
3,553942108

EDIT: sieht die Nullstelle da richtig aus?

Kibo
25.11.2014, 22:50
für T:
tan (T/2) = sqrt((1+e)/(1-e)) *tan(E/2)

tan (T/2) = sqrt((1+0,041)/(1-0,041)) *tan(3,553942108/2)

tan (T/2)= -5,190166064 |arctan

T/2 = -1,380456687
T= -0,690228344


In der Zwischenzeit habe ich vergessen was ich da eigentlich rechne. Stimmt es denn wenigstens?

Kibo
25.11.2014, 23:11
Oh, deinen Link "Beschreibung des Berechnungsverfahrens" habe ich übersehen. Gucke ich mir dann morgen an.

gute Nacht

Bernhard
26.11.2014, 10:23
iterative Nullstellenberechnung nach Newtonverfahren
Hallo Kibo,

die Formeln hast Du korrekt angegeben, aber mit den Zahlenwerten für E geht schon beim ersten Wert etwas schief. Das korrekte Ergebnis für E lautet 3.53607147955... und gibt an, wo sich Io gerade auf der Ellipse befindet. Bei E = Pi befindet sich Io gerade im Aphel, bei E = 0 im Perihel der Ellipse. Werte dazwischen geben den Winkel bezüglich der großen Achse an. Man sollte das aber auch in den verlinkten Seiten nachlesen können.
MfG

Kibo
13.05.2015, 23:38
Hi,

Die Ephemerien habe ich derzeit wieder vergessen, aber dafür gibt es ein Update, dass die Performance verbessert in dem es einen Teil der Berechnung in eine c dll auslagert. Ausserdem kann man die CPU Auslastung in den settings unter Performance einstellen.

mfg