TheiaSim

[SRMeister]

Zusammenfassung:
Es soll simuliert werden, wohin die Bruchstücke des Theiaeinschlags wandern, ob sie sich möglicherweise in bestimmten Regionen des Sonnensystems besonders häufen, oder ob sie auf Planeten niedergehen.
Das Programm dazu ist vorerst der von uns bereits programmierte Integrator nach Bulirsch Stoer. Dieser soll weiter an das Problem angepasst und verbessert werden. Wie genau, das ist noch offen. Möglicherweise werden adaptive Schrittweiten implementiert.
Die Objekte in der Simulation lassen sich in mehrere Klassen einteilen, je nachdem, welche gravitativen Kräfte auf sie wirken und von ihnen ausgehen.

So ich habe nun ein Sourceforge-Projekt mit SVN Repo angelegt.

Zum downloaden per SVN (ich empfehle TortoiseSVN) kurze Anleitung: Auf eigenem Rechner in den Explorer gehen und in einen Ordner gehen wohin ihr es runtergeladen haben wollt. Dann rechte Maustaste und "SVN Checkout..." wählen. Dann oben eingeben "https://theiasim.svn.sourceforge.net/svnroot/theiasim" und auf OK gehen. Das beinhaltet dann im Release-Unterordner auch eine kompilierte TheiaSim.exe, dass ist die Konsolenanwendung. Die Konsolenanwendung kann auch ohne SVN direkt von hier geladen werden.

Bernhard hat bereits Schreibrechte auf das SVN bekommen, wer sonst noch aktiv teilnehmen will bitte für Schreibrechte bei mir melden. Wir sind für jede Unterstützung dankbar.

Das Projekt ist so mit dem kostenlosen Visual C++ 2010 Express (Download) kompilierbar, sollte ohne Warnungen und Fehlermeldungen direkt funktionieren.

@Bernhard: Die vorgenommenen Anpassungen: RA,Dec entfernt, EIH entfernt, CNewtonSP entfernt, CNewton ist jetzt von der abstrakten Basisklasse PSystem abgeleitet. Des Weiteren wird die Zeitkonstante jetzt an DoStepBS übergeben, das ist notwendig für spätere DLL. Das Programm gibt die benötigte CPU Zeit aus (für Benchmark müsste man bei kurzen Integrationszeiten mehrere Werte mitteln, da es sonst stark schwankt)

Ich veröffentliche später das GUI.

Grüße

 

[Bernhard]

So ich habe nun ein Sourceforge-Projekt mit SVN Repo angelegt.

Hallo SRMeister,

genau so etwas in dieser Art schwebte mir auch vor. Das wäre also schon mal erledigt .

ABER:

RA,Dec entfernt

Hoffentlich war das nicht zu vorschnell. Diese paar kleinen Funktionen hätte ich persönlich drin gelassen. Das ist aber nur ein Detail und kann bei Bedarf wieder eingebaut werden (als Amateurastronom stehe ich natürlich auf RA und Dec).

EIH entfernt

OK. Aber auch diese Klasse stört nicht übermäßig, wenngleich es aus Performancegründen sicher nicht die erste Wahl für das neue Projekt TheiaSim ist.

CNewtonSP entfernt

Hm? Da steckt doch der Trick drin, die Position der Sonne über den konstanten Schwerpunkt des Sonnensystems zu berechnen. So etwas könnte für lange Simulationszeiten eventuell Stabilitätsvorteile bringen.

CNewton ist jetzt von der abstrakten Basisklasse PSystem abgeleitet. Des Weiteren wird die Zeitkonstante jetzt an DoStepBS übergeben, das ist notwendig für spätere DLL.

OK

Das Programm gibt die benötigte CPU Zeit aus (für Benchmark müsste man bei kurzen Integrationszeiten mehrere Werte mitteln, da es sonst stark schwankt)

Das gefällt mir.

Ich werde dann noch einen vergleichenden Blick in die mercury-Dokumentation werfen.
Gruß

 

[Bynaus]

Wow, das geht ja schnell. Ich bin echt gespannt, wie sich das entwickelt und was da rauskommt...

Ich hoffe, ist ist okay, wenn ich mich auf die "wissenschaftliche" Seite des Projekts beschränke, als Inputs zu den Ausgangsbedingungen liefere und am Ende versuche, das ganze zu verschriftlichen (selbstverständlich unter angemessener Berücksichtigung aller Beiträge der Beteiligten). Meine Programmierkenntnisse sind zwar nicht gleich Null, aber ich habe das Gefühl, ich würde das Projekt eher verlangsamen, wenn ich meine Finger drinhätte... Ich werde aber natürlich hier im Thread sehr aktiv mitlesen und mich auch mal melden.

Ein möglicher Terminhorizont könnte übrigens Anfang Januar sein: dann ist die Deadline für Abstracts an der Lunar and Planetary Science Conference in Houston, TX. Nicht dass ich jetzt irgendjemanden stressen möchte!

 

[SRMeister]

Hoffentlich war das nicht zu vorschnell. Diese paar kleinen Funktionen hätte ich persönlich drin gelassen. Das ist aber nur ein Detail und kann bei Bedarf wieder eingebaut werden (als Amateurastronom stehe ich natürlich auf RA und Dec).

Würde auch sagen, lass uns das später wieder hinzufügen. Konzentrieren wir uns erstmal auf die eigentliche Aufgabe. Was mir in dem Zusammenhang noch einfällt: Aus demselben Grund habe ich auch RK-4 und Euler rausgenommen. Wir brauchen es ganz sicher vorerst nicht, und später kann man es zur Not wieder einfügen.

Hm? Da steckt doch der Trick drin, die Position der Sonne über den konstanten Schwerpunkt des Sonnensystems zu berechnen. So etwas könnte für lange Simulationszeiten eventuell Stabilitätsvorteile bringen.

Oh okay, dass hatte ich jetzt garnicht mehr im Kopf. Dann tuh ich die wieder mit rein erstmal, aber irgendwann sollten wir uns auf eine Version festlegen. In Bynaus Sinne, dass das Programm auch im wissenschaftlichen Kreis anerkannt werden kann, sollten wir denke ich soweit wie möglich auf eigene Experimente vorerst verzichten und auf bewährte und dokumentierte Methoden zurückgreifen. Dann wird es für ihn wahrscheinlich leichter die Software zu erklären.
Meine Idee eines Ansatzes wäre folgende: Wir suchen erstmal ein paar Papers zusammen, wo aktuelle Methoden vorgestellt werden, möglichst allgemein gehalten sind, so dass wir uns später daran orientieren können und Bynaus Quellen für die eingesetzten Methoden angeben kann.

Ich hoffe, ist ist okay, wenn ich mich auf die "wissenschaftliche" Seite des Projekts beschränke

Das ist natürlich völlig okay.Wo ich mich jetzt wieder, mit etwas zeitlichem Abstand der ganzen Sache nähere, fällt mir irgendwie auf dass wir uns möglicherweise beim letzen Mal etwas zu sehr in den Details der Programmierung verloren haben. Außerdem haben wir möglicherweise auch die Kompetenzen nicht so klar getrennt. Ich selbst würde mich also auch etwas auf das User interface spezialisieren, aber natürlich mit dem Integrator auch mithelfen wo ich kann, vielleicht aber nur theoretisch.

 

[SRMeister]

Hier noch ein N-Body-Integrator mit 4 verschiedenen Integratoren, der Erste ist dem in Mercury wohl sehr ähnlich. Programmiersprache C
http://www.astrobiology.ucla.edu/~varadi/NBI/NBI.html

Hier etwas über die Langzeitstabilität (Lyapunov time):
http://arxiv.org/abs/0708.2875 (Meine Schlussfolgerung: Bis ein paar Mio.Jahre ist wohl der einfache BS Integrator stabil/führt nicht zu Chaos)

We integrate the system of Jovian planets using only Newtonian gravity. The inner
planets are accounted for by adding their masses to the Sun and perturbing the Sun’s position
and linear momentum to equal that of the Sun+Mercucy+Venus+Earth+Mars system. This
ensures that the resonances between the outer planets is shifted by an amount that is second
order in this mass ratio, roughly 3 × 10^−11
(Murray and Holman 1999), which is far smaller
than the uncertainty in the outer-planet positions. We assume constant masses for all objects
and ignore many effects which are probably relevant over a 200My timescale (see for example
Laskar (1999)). We account for solar mass loss at a rate of ˙ m/m ≈ 10^−7
per million years
(Laskar 1999; Noerdlinger 2005), but note that we observe no noticable difference if we keep
the solar mass constant.

Hier Mercury:
http://star.arm.ac.uk/~jec/home.html#publications

auch interessant ist der Taylor 1.4 Integrator (aber nicht für dieses Projekt geeignet)
http://www.maia.ub.es/~angel/taylor/

Finally, Taylor 1.4 allows the user to specify the machine arithmetic to use, including software arithmetics. Out-of-the-box, Taylor 1.4 supports the use of IEEE 754 double precision (64 bit representation with a 53 bit mantissa), Intel extended precision (80 bit representation with a 64-bit mantissa, giving a machine precision of about 10^−19, accessible as long double when using GCC on an Intel machine), the DoubleDouble datatype (Briggs 1996) which provides software quadruple precision in C++, and the GNU Multiple Precision Library, which allows arbitrary precision floating point numbers in C++. Most of our integations using Taylor 1.4 used Intel extended precision, which is almost as fast as double precision and gives about 19 decimal digits of accuracy. Over our 200-million-year integrations using Intel extended precision, Taylor typically had relative energy errors of less than 8 × 10^−14; the worst relative energy error observed in any of our integrations was 2 × 10^−13. Integrations began with the Solar System’s barycentre at the origin with zero velocity. After 200 million years the barycentre drifted a maximum of 3 × 10^−10 AU, while the z component of the angular momentum was always conserved to a relative accuracy better than 3×10^−14. We also performed a suite of quadruple
precision integrations, in which energy and angular momentum were each conserved to at least 26 significant digits over 200 million years.

 

[Bernhard]

Oh okay, dass hatte ich jetzt garnicht mehr im Kopf. Dann tuh ich die wieder mit rein erstmal, aber irgendwann sollten wir uns auf eine Version festlegen.

Wir können die Klasse CNewtonSP erst mal rauslassen und uns bei solchen Details auch an dem mercury-Code orientieren. Ich habe inzwischen die Klasse (oder Struktur) Testparticle hinzugefügt, das Projekt insgesamt etwas aufgeräumt und dabei auch etliche Kommentare und Bezeichnungen in's Englische übersetzt. Der Code sollte so übersichtlich und lesbar wie möglich bleiben, so dass man möglichst bald die Bewegung einiger Testkörper berechnen kann. An dem Programmablauf habe ich nichts geändert.

Die Funktion AddObject habe ich überladen. Die ursprüngliche Version fügt jetzt wie bisher immer Planeten hinzu. Die neue Version ohne den Masseparameter fügt immer Testkörper hinzu.

Der nächste Schritt sollte dann die Berechnung der Bahnen der Testkörper sein. Man kann bei Gelegenheit auch die Radien der Planeten einbauen (willst Du das machen?) und einen Zähler, der registriert, falls ein Testkörper das Innere des Planeten kreuzt (=Absturz des Testkörpers auf den Planeten).

 

[Bynaus]

Man kann bei Gelegenheit auch die Radien der Planeten einbauen (willst Du das machen?) und einen Zähler, der registriert, falls ein Testkörper das Innere des Planeten kreuzt (=Absturz des Testkörpers auf den Planeten).

Ich würde vorschlagen, hier nicht die physischen Radien zu nehmen, sondern den Radius der gravitativen Fokussierung (dh, der Bereich, aus dem ein Teilchen wegen der Gravitation mit dem Planeten kollidieren würde, auch wenn es ihn im gravitationsfreien Raum geometrisch nicht treffen würde). Leider ist dieser Radius von der Geschwindigkeit des Teilchens abhängig.

Eine Herleitung und Formel (Seite 28) findet sich z.B. hier: http://www.scribd.com/doc/52296339/29/Gravitational-focusing

 

[Bernhard]

Hallo SRMeister,

ich fürchte ich habe gerade einen bösen Fehler gemacht. Ich dachte die Dateien TheiaSim.suo und TheiaSim.user sind anwenderspezifische Dateien. Damit die nicht bei jedem Commit von uns gegenseitig geändert werden, habe ich beide Dateien gelöscht. Seitdem hat Visual Studio ein größeres Problem mit dem vorkompilierten Header und bringt ziemlich viele Warnings und auch Fehler.

Kannst Du diese beiden Dateien bitte aus Deinen lokalen Quellen wieder einspielen, oder den vorhanden Fehler beheben?

 

[ismion]

Leider ist dieser Radius von der Geschwindigkeit des Teilchens abhängig.

Ich kenn diesen Ausdruck als "Runaway Growth". Hab den zuletzt in meiner "Geophysik des Sonnensystems" - Vorlesung gesehen

Der Wirkungsquerschnitt ist von der Geschwindigkeit "im Unendlichen" (Sigma) abhängig und von der Fluchtgeschwindigkeit. Die Fluchtgeschwindigkeit ist proportional der Wurzel der Masse und proportional zu Wurzel aus 1/Radius des Planeten. Für Akkretion der Körper muss v_esc immer wieder geupdated werden.

Aber, worauf ich eigentlich hinaus möchte. Für kleine Körper ist v_esc relativ niedrig. Bei der Erde ist das schon "nur" 11 km/s. Kann man nicht die extrem vereinfachende Annahme einführen, dass v_esc << sigma => (v_esc/sigma)^2 ~ 0? Somit wäre die Geschwindigkeitsabhängig dahin.

 

[Bynaus]

Für kleine Körper ist v_esc relativ niedrig.

Nur dass wir uns richtig verstehen: ich will ja nicht die Fragmente (untereinander) akkretieren. Aber wir wollen wissen, wie oft die Fragmente von (den bestehenden oder allenfalls zusätzlichen, hypothetisch eingefügten) Planeten akkretiert werden. Die Fluchtgeschwindigkeit von Planeten ist jedoch durchaus in der Grössenordnung der Bahngeschwindigkeiten (also v_esc ~ sigma), deshalb kann man das nicht vernachlässigen. Bedenke z.B. dass die Fragmente am Anfang eine Geschwindigkeit von ~<1.4 v_esc relativ zur Erde haben.

Bei Asteroiden im Gürtel mag es jedoch sein, dass v_esc << sigma ist (z.B. Vesta, v_esc = 360 m/s - Bahngeschwindigkeit von Vesta = 19 km/s, typische Kollisionsgeschwindigkeiten im Gürtel ~5 km/s => (v_esc/sigma)^2 = ~<1%), das heisst, hier könnte man das wohl vernachlässigen und deren phsische Radien nehmen, um abzuschätzen, wieviele der Fragmente von den Asteroiden akkretiert werden.

 

[ismion]

Ach so, okay. Gut zu wissen das mit den Geschwindigkeiten. Danke für die Aufklärung!

Wenn wir also den geschwindigkeitsabhängigen Term nehmen, dann wird das Problem leider aufwändiger. Du musst dann in jedem Schritt prüfen, ob ein Fragment innerhalb der Einflussregion liegt und ob es zu ner Kollision kommt. Wenn eine Kollision zustande gekommen ist, muss die Masse erhöht werden, was dann wiederum den Einflussradius für die Berechnung der weiteren Fragmente erhöht. Oder man umgeht das Problem indem man nur eine Kollision pro "akkretierendem Objekt" pro Zeitschritt erlaubt.

 

[Bynaus]

Danke für die Aufklärung!

Danke dir für den Hinweis, dass man diesen Effekt im Asteroidengürtel wohl vernachlässigen kann! Damit muss man das wirklich nur für die acht grossen Planeten bestimmen.

ob ein Fragment innerhalb der Einflussregion liegt und ob es zu ner Kollision kommt

Jetzt kommts drauf an, was du mit Einflussregion meinst. Es gibt eine Region, bei der die Bahn des Teilchens abgelenkt wird, ohne dass es mit dem Planeten kollidiert.

Das graviational focussing (im Folgenden GF) hat damit aber nicht so viel zu tun - es erhöht bloss den effektiven Radius des Planeten. Man muss also nur für jedes Teilchen bei gegebener Geschwindigkeit im Unendlichen (aus den Bahnelementen) prüfen, ob es weniger als einen (um einen GF-Faktor erweiterten) Radius des Planeten an dessen Schwerezentrum vorbeizieht. Wenn ja, wird es akkretiert, aus der Simulation entfernt, und der Teilchenzähler des Planeten um eins erhöht. Wenn nicht, wird es nur abgelenkt.

Wenn eine Kollision zustande gekommen ist, muss die Masse erhöht werden, was dann wiederum den Einflussradius für die Berechnung der weiteren Fragmente erhöht.

Das ist, denke ich, unnötig. Zunächst einmal sind die meisten Fragmente masselos, und die paar wenigen, die es nicht sind, haben sehr geringe Massen relativ zu den Planeten, so dass die zusätzliche Masse des Trümmerstücks kaum ins Gewicht fällt. Die Planeten sind zur Zeit des Giant Impact (ca. 30 - 130 Mio Jahre nach dem Anfang des Sonnensystems) auch ohnehin schon praktisch fertig gebildet.

Es interessiert ja nur, wo die Trümmer von Theia gelandet sein könnten und ob wir einen Teil davon im Asteroidengürtel erwarten können.

 

[ismion]

Mich würde noch interessieren, ob sich präferierte "Positionen" für die Trümmer einstellen, trotz radialsymmetrischer Fragment-Anfangspositionen. Aber das sind alles kleine Spielereien.

 

[Bynaus]

Das ist sicher auch eine interessante Frage, auch aus meiner Sicht, und zwar aus folgendem Grund:

Man kann die radiale Symetrie auch so verstehen, dass wir nicht wissen, in welche Richtung die Trümmer von Theia davongeschleudert wurden. Durch Abdecken aller Richtungen bekommen wir trotzdem einen Eindruck davon, wo die Trümmer gelandet sein könnten. Ergeben sich starke Asymetrien, könnte das z.B. das Fehlen (oder Vorhandensein) von Theia-Trümmern im Asteroidengürtel erklären, weil es dann sein könnte, dass Theia einfach "im falschen (richtigen) Winkel" mit der proto-Erde kollidiert ist.

 

[Bernhard]

Kannst Du diese beiden Dateien bitte aus Deinen lokalen Quellen wieder einspielen, oder den vorhanden Fehler beheben?

Das hat sich erledigt. Bei den betreffenden Dateien musste ledigliche ein #include "stdafx.h" eingefügt werden. Die Konsolenanwendung läßt sich damit wieder fehlerfrei übersetzen.

 

[Bernhard]

Man kann die radiale Symetrie auch so verstehen, dass wir nicht wissen, in welche Richtung die Trümmer von Theia davongeschleudert wurden.

Der mittlerweile eingebaute Zufallsgenerator erzeugt aktuell einfach gleichverteilte Zufallswinkel theta und phi im Bereich 0..pi, bzw. 0..2*pi. Diese Bereiche lassen sich sehr leicht anpassen, so dass ein bevorzugter Streubereich für die Trümmer simuliert werden könnte.

 

[Ich]

Besser gleichverteilt über y [0...1] und dann theta=arccos(1-2y).

 

[ismion]

Man könnte die Trümmer auch einem gegebenem Dichteprofil gehorchend verteilen, wenn wir schon dabei sind, die Verteilung zu diskutieren.

 

[SRMeister]

ich denke man kann den Großteil der Trümmer in die Ebene setzen, der Einschlag hat ja sicher auch in etwa in der Ebene stattgefunden. Also vielleicht theta -20° bis +20° oder so ?

 

[Bernhard]

Besser gleichverteilt über y [0...1] und dann theta=arccos(1-2y).

Das würde aber doch die äquatoriale Richtung ein wenig bevorzugen und die Richtungen in der Nähe der beiden Koordinatenpole etwas ausdünnen .

BTW: Der erste Prototyp mit 10 Testpartikeln ist bei mir gerade ohne Fehler durchgelaufen . Ein Visualisierungstool wäre für die zu erwartenden Bahnen im Asteroidengürtel für den Anfang sehr hilfreich.