PDA

Archiv verlassen und diese Seite im Standarddesign anzeigen : TheiaSim



Seiten : [1] 2

SRMeister
06.08.2012, 20:41
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 (http://theiasim.svn.sourceforge.net/) 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 (http://theiasim.svn.sourceforge.net/viewvc/theiasim/Release/) 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 (https://www.microsoft.com/germany/express/download/default.aspx#)) 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
06.08.2012, 21:54
So ich habe nun ein Sourceforge-Projekt (http://theiasim.svn.sourceforge.net/) 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
06.08.2012, 22:08
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
07.08.2012, 00:45
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
07.08.2012, 02:54
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
07.08.2012, 04:54
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
07.08.2012, 07:52
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
07.08.2012, 09:57
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
07.08.2012, 10:12
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 :D

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
07.08.2012, 10:35
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
07.08.2012, 11:00
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
07.08.2012, 11:35
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
07.08.2012, 11:59
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
07.08.2012, 12:09
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
07.08.2012, 12:33
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
07.08.2012, 12:37
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
07.08.2012, 13:32
Besser gleichverteilt über y [0...1] und dann theta=arccos(1-2y).

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

SRMeister
07.08.2012, 13:59
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
07.08.2012, 14:05
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 :confused: .

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.

ismion
07.08.2012, 14:15
Das würde aber doch die äquatoriale Richtung ein wenig bevorzugen und die Richtungen in der Nähe der beiden Koordinatenpole etwas ausdünnen :confused: .

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.

Kannst Du nicht erst einmal schlichtweg die Positionen mit Gnuplot darstellen? Schreibst Du für jedes "Partikel" ein einzelnes File raus?

Ich
07.08.2012, 14:25
Das würde aber doch die äquatoriale Richtung ein wenig bevorzugen und die Richtungen in der Nähe der beiden Koordinatenpole etwas ausdünnen
Nö, du willst ja konstante Dichte pro Raumwinkel (also pro Flächenelement auf der Kugeloberfläche). Die Metrik ist da in Umfangsrichtung ds=sin(theta)*dphi, von 0° bis 10° am Pol hat's also erheblich weniger Fläche als von 80° bis 90° am Äquator, das musst du berücksichtigen.

UMa
07.08.2012, 14:27
Hallo Bernhard,

Das würde aber doch die äquatoriale Richtung ein wenig bevorzugen und die Richtungen in der Nähe der beiden Koordinatenpole etwas ausdünnen :confused: .
die von "Ich" angegebene Formel ergibt eine exakte Gleichverteilung über die Sphäre. Eine Gleichverteilung von theta in 0..pi bevorzugt dagegen die Pole. Schaue dir z.B. einen Globus an. Zwischen 80-90 Grad Nord ist weniger Fläche, als zwischen 0-10 Grad Nord.

Grüße UMa

Bynaus
07.08.2012, 14:47
Ich sehe das auch so: Ich's Verteilung würde eine echte Gleichverteilung über die Sphäre bringen.

Allerdings hat auch SRMeister recht: Nach aller Wahrscheinlichkeit bewegten sich sowohl Theia als auch die Erde innerhalb der Ekliptik, entsprechend auch die Trümmer (hatte ich zuvor gar nicht bedacht). Also könnte man sich überlegen, ob man die Sphäre jeweils bei +-20° (oder so) beschneiden will (nur schon, um Rechenzeit zu sparen bzw. in der zur Verfügung stehenden Zeit mehr Fragmente zu simulieren).

Ich würde sagen, dass das durchaus eine vernünftige (verteidigbare) Annahme wäre. Was meint ihr?

ismion
07.08.2012, 14:57
Man könnte durchaus beide Fälle simulieren. Willst Du die Beschneidung nur in theta, oder auch in Phi richtung durchführen? In beiden Richtungen würde, meines Verständnisses nach, eher einem frühen Zeitpunkt nach einer Kollision entsprechen. Aber das alles ist ja nur eine Frage der gewählten Anfangsbedingungen.
Ich bin der Meinung, dass man solch eine Annahme durchaus gerechtfertigt machen darf!

Bynaus
07.08.2012, 15:25
Also wenn Theta der Winkel ist, der senkrecht auf der Ekliptik steht, dann würde ich nur diesen Winkel beschneiden. Den anderen (innerhalb der Ekliptik) würde ich aus den oben genannten Gründen bei 360° bzw. 2Pi belassen. Die Kollision hat wohl Trümmer v.a. in eine bevorzugte Richtung geschleudert - aber wir wissen nicht, welche das ist, also simulieren wir einfach alle und schauen, was rauskommt.

ismion
07.08.2012, 15:35
Also wenn Theta der Winkel ist, der senkrecht auf der Ekliptik steht, dann würde ich nur diesen Winkel beschneiden. Den anderen (innerhalb der Ekliptik) würde ich aus den oben genannten Gründen bei 360° bzw. 2Pi belassen. Die Kollision hat wohl Trümmer v.a. in eine bevorzugte Richtung geschleudert - aber wir wissen nicht, welche das ist, also simulieren wir einfach alle und schauen, was rauskommt.

Aber müsste dann nicht ein bestimmter Teil "schneller" sein? Du hast ja einen unterschiedlich starken Impulsübertrag des Impaktors, je nach Aufschlagswinkel. Oder ist dies schon berücksichtigt in der Annahme, dass die Fragmete unterschiedliche Anfangsgeschwindigkeiten erhalten?

Bernhard
07.08.2012, 15:39
Also wenn Theta der Winkel ist, der senkrecht auf der Ekliptik steht, dann würde ich nur diesen Winkel beschneiden.
Theta = 0 bedeutet senkrecht zur Ekliptik im Norden und Theta = 180° bedeutet senkrecht zur Ekliptik im Süden. In welcher Art soll also Theta behandelt werden? Ein bestimmter Bereich (70° ... 110° ) und zusätzlich Ichs arccos? Oder sind 40° Streuung zu wenig?

Bernhard
07.08.2012, 15:41
Oder ist dies schon berücksichtigt in der Annahme, dass die Fragmete unterschiedliche Anfangsgeschwindigkeiten erhalten?
Momentan ist die Anfangsgeschwindigkeit der Fragmente zwischen 1.0 * v_esc und 1.4 * v_esc gleichverteilt.

Bernhard
07.08.2012, 15:44
Schaue dir z.B. einen Globus an. Zwischen 80-90 Grad Nord ist weniger Fläche, als zwischen 0-10 Grad Nord.
Treffer pro Fläche. Damit kann ich was anfangen :) .

ismion
07.08.2012, 15:45
Okay, das klingt gut für den Anfang.

Bynaus, besaß die Urerde zu dem Zeitpunkt schon eine (evtl. dünne) Atmosphäre, weiß man da was? Könnte die eine Rolle spielen? Zumindest bei kleineren Fragmenten in Form von Reibung? (Ich glaube zwar selbst, dass das vernachlässigbar ist, aber man kann ja brainstormen).

Gibt es das Programm auch für Linux-User? Hab kein Windows :o

Bernhard
07.08.2012, 15:55
Gibt es das Programm auch für Linux-User? Hab kein Windows :o
Momentan Nein. Nachdem es aber noch keine grafischen Elemente enthält, bzw. als Ausgabe lediglich cout verwendet ist eine Migration nach Linux jetzt noch vergleichsweise leicht zu machen.

SRMeister
07.08.2012, 16:01
... ist eine Migration nach Linux jetzt noch vergleichsweise leicht zu machen.

wird auch so bleiben. die DLL Schnittstelle und die Konsolenausgabe können und sollen koexistieren. Vielleicht findet sich ja auch jemand, der zusätzlich irgendwie eine Anbindung an graphische Tools (gnuplot?) in Linux machen kann? Reicht ja ne passende Ausgabe in ne Datei.

Bernhard
07.08.2012, 16:04
Besser gleichverteilt über y [0...1] und dann theta=arccos(1-2y).
Habe es so eingebaut. Bei Bedarf kann der Bereich von y noch begrenzt werden.

Bernhard
07.08.2012, 16:10
Hab kein Windows :o
Wir können von mir aus auch mit Qt weitermachen. Das läuft auch unter Linux.

ismion
07.08.2012, 16:15
wird auch so bleiben. die DLL Schnittstelle und die Konsolenausgabe können und sollen koexistieren. Vielleicht findet sich ja auch jemand, der zusätzlich irgendwie eine Anbindung an graphische Tools (gnuplot?) in Linux machen kann? Reicht ja ne passende Ausgabe in ne Datei.

ich kenne mich mit dem programm nicht aus, aber eine art gnuplot/IDL/GDL interface zur darstellung der ergebnisse in linux könnte ich erstellen. jenachdem wie kompliziert das programm ist könnte ich den output-part da auch übernehmen oder ich bekomme einfach von euch eine test-datei, sodass ich den visualisierungsteil in linux machen kann.

Ich
07.08.2012, 16:26
Ich hab nicht so ganz mitgekriegt was ihr machen wollt, aber nach Drehimpulserhaltung lag die Verbindungslinie Erde-Theia kurz vor dem Einschlag nicht in der Ekliptik, sondern eher in der Ebene der Mondumlaufbahn. Dorthin sollte auch der Schutt.
Ich kenn' mich da nicht so wirklich aus, aber weil diese Ebene heute nicht in der Ekliptik liegt würde ich annehmen, dass sie das damals auch nicht tat, Präzession hin oder her. Was eine etwas größere Streuung als +- 20° zur Ekliptik nahelegen würde.

ismion
07.08.2012, 16:28
Ich hab nicht so ganz mitgekriegt was ihr machen wollt, aber nach Drehimpulserhaltung lag die Verbindungslinie Erde-Theia kurz vor dem Einschlag nicht in der Ekliptik, sondern eher in der Ebene der Mondumlaufbahn. Dorthin sollte auch der Schutt.
Ich kenn' mich da nicht so wirklich aus, aber weil diese Ebene heute nicht in der Ekliptik liegt würde ich annehmen, dass sie das damals auch nicht tat, Präzession hin oder her. Was eine etwas größere Streuung als +- 20° zur Ekliptik nahelegen würde.

Guter Punkt, den Du da angibst! Ich denke, sowas könnte auch untersucht werden, inwieweit die heutige Konfiguration reproduziert werden kann durch verschiedene AB

Bynaus
07.08.2012, 20:09
Ich denke, 20° Streuung um die Ekliptik sind realistisch (man kann das ja mal als freien Parameter gestalten und immer noch ändern, wenn nötig / EDIT: Bernhard hat das offenbar bereits gemacht).


Bynaus, besaß die Urerde zu dem Zeitpunkt schon eine (evtl. dünne) Atmosphäre, weiß man da was? Könnte die eine Rolle spielen? Zumindest bei kleineren Fragmenten in Form von Reibung? (Ich glaube zwar selbst, dass das vernachlässigbar ist, aber man kann ja brainstormen).

Brainstormen ist immer gut, und verrückte Fragen stellen auch, man weiss nie, was man übersieht (gilt natürlich für alle, die Mitlesen)! Nun, die Erde war kurz nach dem Giant Impact sowas wie eine superheisse Kugel aus Magma und verdampften Silikaten, von dem her... :) Ich denke, im Allgemeinen sind die Atmosphären wohl deutlich kleiner als der GF-Radius, daher darf man sie getrost vernachlässigen.


Ich hab nicht so ganz mitgekriegt was ihr machen wollt, aber nach Drehimpulserhaltung lag die Verbindungslinie Erde-Theia kurz vor dem Einschlag nicht in der Ekliptik, sondern eher in der Ebene der Mondumlaufbahn.

Würde man meinen, aber die Ebene der Mondumlaufbahn kann sich später auch durch Resonanzen ändern - tatsächlich hat man schon sehr viel zu der Frage vorgeschlagen und diskutiert, warum die Mondbahn weder in der Ekliptikebene noch in der Äquatorebene der Erde liegen könnte. Aber ich würde sagen, das ist (wenn auch interessant) "beyond the scope of this work".

Es geht bei diesem Projekt in ganz wenigen Worten darum, dass nach einem neuen Modell des Giant Impact eine Menge Trümmer (total bis zu einer Marsmasse) von Theia ins Sonnensystem geschleudert wurden, und wir simulieren möchten, was mit diesen geschieht, ob man z.B. erwarten darf, dass einige davon als Asteroiden im Gürtel gelandet sind (und wenn ja, wieviele).

Übrigens scheint mir, dass wir nicht vergessen sollten, eine Yarkovski-ähnliche Kraft einzubauen, sonst dürften wir den Anteil der Fragmente, der letztlich auf die Sonne fällt, deutlich überschätzen (es sei denn, sie sind ohnehin alle so schnell weg, dass es keine grosse Rolle spielt, was ich aber nicht glaube).

Bernhard
07.08.2012, 20:16
@SRMeister,

zwei Tips oder Bitten meinerseits damit der Code langfristig möglichst gut lesbar bleibt:

1) Eine Nummerierung, die auch für Anwender eine Rolle spielt (z.B. die Nummerierung der Testkörper) sollte immer bei 1 und nicht bei 0 beginnen.

2) Leerzeilen und Leerzeichen erhöhen die Übersichtlichkeit. Über diesen Punkt kann man zwar streiten, aber die Praxis bestätigt das eigentlich recht deutlich. Gerade wenn sehr unterschiedliche Programmierer an einem Projekt arbeiten sollte das Vorteile bringen.

Bernhard
07.08.2012, 20:19
Ich denke, 20° Streuung um die Ekliptik sind realistisch (man kann das ja mal als freien Parameter gestalten und immer noch ändern, wenn nötig / EDIT: Bernhard hat das offenbar bereits gemacht).
Hallo Bynaus,

ich habe vorerst nur den Tip von Ich eingebaut um die Kugelsymmetrie zu gewährleisten. Die Beschränkung auf 20° Streuung baue ich dann noch ein. Hier bietet sich ein define in der Datei Constants.h an.

SRMeister
07.08.2012, 22:02
ja Bernhard, werde ich ab jetzt so machen.
Wo wir gerade dabei sind: Nur Membervariable mit m_Xyz benennen, Parameter von Funktionen bitte als Xyz benennen.
Und wenn man merkt, dass in 2 klassen 90% gleich ist, dann am besten gleich beide von einer gemeinsamen Basisklasse ableiten. Das habe ich jetzt nochmal gemacht, also nicht wundern. Basisklasse CParticle abgeleitete Klassen CPlanet, CTestParticle.
Das musste ich jetzt auch machen, da die Schnittstelle nurnoch möglichst wenig geändert werden soll.
Ich will ja eigentlich in Deiner Baustelle nicht viel rumdoktorn, aber jetzt musste das nochmal sein :)

Das GUI ist bald fertig!

Bernhard
07.08.2012, 22:15
Wo wir gerade dabei sind: Nur Membervariable mit m_Xyz benennen, Parameter von Funktionen bitte als Xyz benennen.
Das kenne ich genauso und macht auch wirklich Sinn. Dummerweise haben wir vermutlich immer noch massenweise Membervariablen ohne das Prefix.



Basisklasse CParticle abgeleitete Klassen CPlanet, CTestParticle.

Wäre da nicht CBody besser?



Das GUI ist bald fertig!

Sehr gut.

SRMeister
07.08.2012, 22:36
okay, gib mir bitte bis morgen zeit, bevor du große Veränderungen einbaust, sonst kommen wir uns evtl irgendwo in die Quere. Ich mach jetzt alles ordentlich objektorientert (abstrakte Basisklassen) und vor allem versuche ich alle "doppelten Algorithmen" noch aufzulösen - also alles auf einen Ursprung zurückzuführen.

Bernhard
07.08.2012, 22:47
Gerne. Gib einfach Bescheid, wenn Du wieder pausieren willst, bzw. fertig bist.

SRMeister
08.08.2012, 00:27
so ich konnte jetzt alle redundanten Klassen und Funktionen aufräumen und das GUI erstmal damit lauffähig bekommen.

Link (https://sourceforge.net/projects/theiasim/files/)

Ich habe dabei auch die doppelten Integratorfunktionen wieder gelöscht. Deswegen ist die Komplexität momentan quadratisch steigend. Das kann man auch mit dem GUI beweisen, denn bei doppelt so vielen Körpern dauert die Berechnung pro Step viermal so lange.

Kannst du bitte versuchen in die CNewton::Evaluate() wieder Bedingungen einzubauen, damit die Berechnung für die Testpartikel wieder stimmt (und Komplexität wieder linear wird)?

Außerdem ist da auch ein Fehler, wodurch die Testpartikel wild durch die Kante springen, da ist wohl die Rechnung irgendwo falsch.

Außerdem werden die Testpartikel jetzt immer an der aktuellen Stelle der Erde (plus Entfernung) erstellt.
Es gibt noch 100 andere kleine Änderungen.
Ich bin soweit fertig mit dem Integrator und würde Dir das Feld wieder überlassen, aber worum ich Dich jetzt bitten würde:
Die Klasse PSystem benötigt vorerst keine Änderungen - das ist quasi die Schnittstelle für die DLL, das heißt, dass du bitte das Programm nur soweit verändern solltest, dass die dort genannten Funktionen ihre Gültigkeit behalten, also du kannst den Inhalt (Definition in CNewton) gern ändern, nur die Deklaration(Funktionsname in PSystem) nicht. Und bitte versuchen, neue Funktionalität immer "unterhalb" dieser Funktionen zu implementieren. Ich hoffe das ist jetzt verständlich erklärt, sonst nochmal bescheid sagen.


Wenn du jetzt Änderungen vornimmst und im GUI testen willst, brauchst du nur oben bei Projektkonfigurationen "Release - DLL" auswählen, dann erzeugt er die DLL, die musst du einfach in das Verzeichnis vom GUI kopieren.
Grüße

Bernhard
08.08.2012, 08:50
Ich habe dabei auch die doppelten Integratorfunktionen wieder gelöscht. Deswegen ist die Komplexität momentan quadratisch steigend.
Hallo SRMeister,

die Testkörper werden anders berechnet, als die Planeten, weil diese gravitativ nicht auf die Planeten einwirken. Du hast beim Aufräumen also eine ganze Menge Arbeit von mir zerstört!! Die Trennung in der Berechnung zwischen Planeten und Testkörper sollte erhalten bleiben, so wie es in Version 23 noch implementiert war. Ich hatte das hier im Forum übrigens auch so erwähnt. EDIT: Aber gut. Nachdem ich auch heute etwas Zeit über habe, baue ich die Funktionen halt wieder ein. Lass Sie zukünftig aber bitte drin. Bynaus hat diese Strukturierung ebenfalls gut geheißen.



Das kann man auch mit dem GUI beweisen, denn bei doppelt so vielen Körpern dauert die Berechnung pro Step viermal so lange.

Scheinbar geht die Laufzeit für die Berechnung der Planeten eben quadratisch. Hättest Du die alten Funktionen drin gelassen, würde die Laufzeit mit der Anzahl der Testkörper linear steigen



Außerdem ist da auch ein Fehler, wodurch die Testpartikel wild durch die Kante springen, da ist wohl die Rechnung irgendwo falsch.

Bei der Initialisierung der Geschwindigkeit der Testkörper habe ich die Geschwindigkeit der Erde vergessen. Das läßt sich relativ leicht beheben.



Außerdem werden die Testpartikel jetzt immer an der aktuellen Stelle der Erde (plus Entfernung) erstellt.

So soll es ja auch sein. Die Testpartikel sollen an der Erdoberfläche starten.



Und bitte versuchen, neue Funktionalität immer "unterhalb" dieser Funktionen zu implementieren.

So lange sich alles übersetzen läßt einverstanden.

Bernhard
08.08.2012, 10:07
Die Trennung in der Berechnung zwischen Planeten und Testkörper sollte erhalten bleiben, so wie es in Version 23 noch implementiert war.
Hallo SRMeister,

Du hast vergessen die Dateien der Implementierung von CBody hinzuzufügen -> rechte Maustaste TotoiseSVN/Add. Ich werde dann doch warten, bis das Projekt wieder übersetzbar ist. Vielleicht fügst Du dabei auch gleich die Funktionen EvaluateTP, DoStepMidpoint4TP, DoStepBS4TP wieder mit ein. Die unterscheiden sich nämlich von den Funktionen ohne TP.

Bernhard
08.08.2012, 10:52
Bedenke z.B. dass die Fragmente am Anfang eine Geschwindigkeit von ~<1.4 v_esc relativ zur Erde haben.
Hallo Bynaus,

erste Rechnungen mit einer lokalen Kopie des Programmes zeigen, dass sich bei einer relativen Startgeschwindigkeit von 1.0 * v_esc bis 1.4 * v_esc die Fragmente (bis auf statistische Ausnahmen) nach ein paar Jahren z.T. bei über 5e4 AE Abstand von der Sonne befinden, also wirklich sehr weit draußen. Kann das stimmen?

Bynaus
08.08.2012, 11:23
erste Rechnungen mit einer lokalen Kopie des Programmes zeigen, dass sich bei einer relativen Startgeschwindigkeit von 1.0 * v_esc bis 1.4 * v_esc die Fragmente (bis auf statistische Ausnahmen) nach ein paar Jahren z.T. bei über 5e4 AE Abstand von der Sonne befinden, also wirklich sehr weit draußen. Kann das stimmen?

Nun, schon möglich, je nach dem, in welche Richtung die Trümmer davonfliegen. Wenn die Erde eine Orbitalgeschwindigkeit von 30 km/s hat, dann ist die Fluchtgeschwindigkeit da sqrt(2)*30km/s, also etwa 42 km/s. Die Fluchtgeschwindigkeit der Erde beträgt 11 km/s, und 1.4 mal diese Geschwindigkeit ist 15 km/s. Das heisst, viele der Trümmer werden wohl mit lokaler Fluchtgeschwindigkeit (rel. zur Sonne) unterwegs sein, was durchaus erklären könnte, warum sie so schnell soweit weg sind (und wohl auch nie zurück kommen...).

Hm. Mit wievielen Teilchen hast du denn das bisher getestet? bzw., wieviele "statistische Ausnahmen" gibt es? Ich geh dem mal nach, 1 - 1.4 v_esc ist die Anfangsgeschwindigkeit von Theia relativ zur Erde im Modell, ich hab das jetzt einfach "ungebremst" auf die Trümmer übertragen, aber es kann natürlich sein, dass die Trümmer üblicherweise langsamer unterwegs sind. Anderseits können sie wohl kaum langsamer als mit 1 v_esc unterwegs sein, sonst wären sie wohl gebremst worden und zur Erde zurückgefallen.

Übrigens: hier könnte man sich auch noch inspirieren lassen: http://www.jstor.org/stable/2890263

EDIT: Aha, halt, natürlich (siehe Bemerkung 15 im oben verlinkten Paper). Um die Endgeschwindigkeit eines Teilchens, das ein Gravitationsfeld verlässt, zu berechnen, kann man natürlich nicht einfach addieren. Die Endgeschwindigkeit relativ zum Planeten wird v_inf = sqrt (v_0^2 - v_esc^2) betragen, wobei v_0 den oben genannten Wert von 1 - 1.4 v_esc haben soll. Ich denke, das sollte das Problem lösen.

Bernhard
08.08.2012, 11:30
In der Berechnung ist noch ein Fehler, weil die Ortsänderungen der Testkörper noch größer sind, als Startgeschwindigkeit * Simulationszeit. Ich versuche dann mal denn Fehler zu finden...

Bernhard
08.08.2012, 11:44
Hallo zusammen,

habe den Fehler schon gefunden. Bei der Berechnung der Fluchtgeschwindigkeit habe ich die Fluchtgeschwindigkeit einmal zu oft berücksichtigt und dadurch quadriert. Die neuen Zahlen über 10 Jahre sind sehr interessant. Die 10 Trümmer verteilen sich dabei auf ein Gebiet zwischen 0.66 und gut 6 AE!! :cool: :cool: .

Sobald SRMeister mit seinen Änderungen fertig ist, können wir nach Lust und Laune simulieren, was das Zeug hält... :)

SRMeister
08.08.2012, 12:02
okay ich hab es wieder eingesetzt, auch wenn es meiner Ansicht nach nicht schön ist, weil die Funktionen zu 99% identisch sind, wie ich sehe sind DoStepMidpoint und DoStepMidpoint4TP sogar komplett identisch, bis auf die Elemente für die sie ausgeführt werden. Schöner wäre es, wenn man jeweils nur eine Funktion hätte, und die Berechnung je nach (if Mass() == 0) unterschiedlich wäre. Aber das überlasse ich jetzt dir.

Übrigens habe ich einen Bug entdeckt, der möglicherweise die ganze Zeit da war, und der möglicherweise die letzten 60km Differenz erklären könnte, die wir damals noch zum JPL hatten.


void CNewton::Evaluate(double dt, int i)
{
...
for( m_source = m_sun; m_source != m_Last; m_source = m_source->next )
}

2 forSchleifen waren unvollständig, da sie bereits beim vorletzten Element der LinkedList aufhörten.
korrekt wäre (und war auch fast überall korrekt):

for( m_source = m_sun; m_source != NULL; m_source = m_source->next )

Die Testparticle fliegen immer noch wild umher. Kann man im GUI bei Detail View gut sehen. (Edit: oh okay dann bis gleich :) )

Grüße

SRMeister
08.08.2012, 12:34
ich überlasse den Fehler dir, bei mir fallen die TPs immer in die Sonne :)

Bernhard
08.08.2012, 12:41
Hallo SRMeister,

Nochmal zur generellen Struktur des Programmes:

a) Planeten und Testkörper sollten wirklich strikt voneinander getrennt werden. D.h. getrennte Listen und getrennte Funktionen damit die Performance nicht leidet. Wenn Du bei jeder Funktion die gemeinsame Liste erst filterst wird das Programm meiner Meinung nach zu langsam.

b) TP soll immer für Test Particle stehen und nicht für Planet. Sonst verliert man viel zu leicht den Überblick, was das Programm eigentlich berechnet.

c) Die Bedingung m=0 halte ich für eine Unterscheidung der zwei Listen deswegen für ungeeignet, weil ein Testkörper später vielleicht auch mal eine Masse haben kann/soll.

Die Fehler in der Bewegung der Testkörper (= Test particle) sind übrigens sehr schnell behoben, weswegen ich zuerst die allgemeine Struktur des Programmes klären möchte. Wir müssten auch noch klären, wer jetzt weiter machen soll/will.

Bernhard
08.08.2012, 12:42
2 forSchleifen waren unvollständig, da sie bereits beim vorletzten Element der LinkedList aufhörten.
korrekt wäre (und war auch fast überall korrekt):
Das ist mir heute auch schon aufgefallen.

SRMeister
08.08.2012, 13:15
Hallo SRMeister,

Nochmal zur generellen Struktur des Programmes:

a) Planeten und Testkörper sollten wirklich strikt voneinander getrennt werden. D.h. getrennte Listen und getrennte Funktionen damit die Performance nicht leidet. Wenn Du bei jeder Funktion die gemeinsame Liste erst filterst wird das Programm meiner Meinung nach zu langsam.
Das würde ich so nicht sagen. Getrennte Listen, vielleicht. Aber getrennte Funktionen? Schau mal, die sind fast identisch. Sowas ist keine ordentliche Programmierung. Die Performance leidet ganz sicher nicht, wenn man ganz zu Anfang der Fkt. einmal den Startpunkt oder Endpunkt festlegt (FirstTP = erster TP)
eine simple if abfrage wegen der Masse ist sogut wie kostenlos (die CPU schafft sowas in 1/2 Taktzyklus und der Compiler versteckt das wahrscheinlich sowieso weshalb es ganz kostenlos ist)


b) TP soll immer für Test Particle stehen und nicht für Planet. Hab ich es irgendwo falsch verwendet? Dachte eigentlich nicht.


c) Die Bedingung m=0 halte ich für eine Unterscheidung der zwei Listen deswegen für ungeeignet, weil ein Testkörper später vielleicht auch mal eine Masse haben kann/soll.
Dann kann man es ja so machen, dass Mass() nur unter bestimmten Bedingungen eine Masse ausgibt. Schon jetzt kann man theoretisch dem TP eine Masse geben die aber nicht über Mass() ausgegeben wird. Das wäre die elegantere Variante und auch aus objektorienterter Sicht sinnvoll. Denn nur das Objekt selbst soll seine Eigenschaften bestimmen, nicht die Funktionen die es benutzen.



Die Fehler in der Bewegung der Testkörper (= Test particle) sind übrigens sehr schnell behoben, weswegen ich zuerst die allgemeine Struktur des Programmes klären möchte. Wir müssten auch noch klären, wer jetzt weiter machen soll/will.
Also mein Vorschlag wäre wirklich, dass wir einheitliche Funktionen für alle Klassen von Objekten benutzen. Es kann ja auch später noch andere Unterteilungen geben (zb mit Radius, ohne Radius usw) willst du dann wieder alle Integratorfunktionen verdoppeln?

a) Man könnte ja 2 getrennte LinkedLists machen und die beiden dann unabhängig durchlaufen.

b) Zu den TP-Fehler: Ich würde außerdem noch das Objekt Erde direkt an den TP-Konstruktor übergeben, dann kann er sich dort die Pos und Vel rausziehen. Ist doch besser als über 6 Argumente.

Wenn a+b okay ist, dann würde ich das noch schnell machen.

UMa
08.08.2012, 13:29
Hallo,

ich würde es so implementieren:
Keine Unterscheidung zwischen Planeten und Testköpern, aber alle massebehafteten Körper zuerst im Array. Dann braucht man nur die Beschleunigungsberechnung so zu ändern, dass dort die Schleife nur bis "nm" geht. Alle anderen Schleifen weiterhin bis "n". Wenn "n" die Gesamtzahl der Körper ist und "nm" die Anzahl der ersten "nm", die eine Masse haben und deren Gravitation berücksichtigt wird.

Grüße UMa

SRMeister
08.08.2012, 13:33
so wie du vorschlägst ist das Array(eig. Verkettete Liste) ja aufgebaut, aber es gibt trotzdem 2 Sätze von Integratorfunktionen, obwohl die fast identisch sind, bis auf eine.

Bernhard
08.08.2012, 13:39
Damit wir uns hier nicht endlosen Formfragen verlieren habe ich meine aktuelle Version eingespielt. Die Verteilung nach zehn Jahren geht doch etwas weiter raus, aber probiert es einfach selbst mal aus.

Bernhard
08.08.2012, 13:40
2 forSchleifen waren unvollständig, da sie bereits beim vorletzten Element der LinkedList aufhörten.
korrekt wäre (und war auch fast überall korrekt):
Das stimmt so, weil sich sonst Pluto selbst beschleunigen würde.

SRMeister
08.08.2012, 14:06
Die Fehler in der Bewegung der Testkörper (= Test particle) sind übrigens sehr schnell behoben, weswegen ich zuerst die allgemeine Struktur des Programmes klären möchte. Wir müssten auch noch klären, wer jetzt weiter machen soll/will.

Also das mit den klären hat sich dann wohl erledigt?

Ich hatte die Änderungen nicht ohne Grund vorgenommen. Die Schnittstelle und damit die Windowsanwendung kann mit den TPs nichts anfangen.
Deswegen hatte ich Planeten und TPs von Body abgeleitet.

Jetzt ist es so: 2 Klassen von fast identischen Objekten, 2 Sätze von fast identischen Integratorfunktionen, 2 Listen die nicht kompatibel sind.
Jetzt müsste ich, damit das Windows Programm klarkommt, auch fast die komplette Schnittstelle verdoppeln und das Windows Programm auch noch ewig anpassen.


Die letzte Version von mir (rev30) hatte die Integratorfunktionen doppelt und es hat auch alles funktioniert, also ich verstehe wirklich nicht warum du jetzt den Schritt zurück gemacht hast. Den ersten Fehler in den TPs hätte man auch dort beheben können.

Die Verteilung nach zehn Jahren geht doch etwas weiter raus, aber probiert es einfach selbst mal aus.

Die TPs sind immernoch nicht korrekt erstellt, dass wäre dir sicher aufgefallen, hättest du das im GUI in 1-Tages schritten mal verfolgt. Fehler: Die TPs müssen bei Erstellung auch die Geschwindigkeit der Erde(oder des Systems Erde-Mond) bekommen, denn dass ist ihr lokales Koordinatensystem. Sonst fallen sie einfach auf die Sonne. Edit ok doch das stimmt soweit, auch wenns etwas ungeschickt implementiert ist (siehe mein Hinweis oben zur Übergabe an den Konstruktor)

Bernhard
08.08.2012, 15:36
a) Man könnte ja 2 getrennte LinkedLists machen und die beiden dann unabhängig durchlaufen.

b) Zu den TP-Fehler: Ich würde außerdem noch das Objekt Erde direkt an den TP-Konstruktor übergeben, dann kann er sich dort die Pos und Vel rausziehen. Ist doch besser als über 6 Argumente.

Wenn a+b okay ist, dann würde ich das noch schnell machen.
zu a) Schau Dir mal bitte die aktuelle Version an und prüfe, ob das nicht vielleicht schon so gemacht ist :mad: .
zu b) Gute Idee.

SRMeister
08.08.2012, 15:58
zu a) Schau Dir mal bitte die aktuelle Version an und prüfe, ob das nicht vielleicht schon so gemacht ist

Ja das ist jetzt so, aber du hast alle Änderungen rückgängig gemacht, insbesondere die CBody klasse ist verschwunden.
Jetzt funktioniert nurnoch die Konsolenanwendung, der Windowsteil ist jetzt komplett unbrauchbar.
Wie schon vorhin beschrieben, waren die Änderungen notwendig und auch sinnvoll, um Redundanzen zu verringern, die Übersichtlichkeit zu erhöhen, den Code objektorientierter zu gestalten und vorallem um die Schnittstelle funktionsfähig zu halten. Die Schnittstelle war von Anfang an seit rev1 vorhanden nur durch das "verdoppeln" aller möglichen Funktionen/ Klassen hast du die leider komplett unbrauchbar gemacht.


Ich weis garnicht mehr was ich jetzt machen soll... ?

Vielleicht wäre es praktikabler wenn man eben doch 2 unabhängige Varianten erstellt... keine Ahnung??

Bernhard
08.08.2012, 16:35
Die Schnittstelle war von Anfang an seit rev1 vorhanden nur durch das "verdoppeln" aller möglichen Funktionen/ Klassen hast du die leider komplett unbrauchbar gemacht.
Hallo SRMeister,

ich werde dann mal versuchen ein eigenes svn-Projekt bei sourceforge.net aufzumachen. Langfristig sehe ich da momentan keine Alternative dazu.
Gruß

ralfkannenberg
08.08.2012, 16:54
Ich denke, 20° Streuung um die Ekliptik sind realistisch (man kann das ja mal als freien Parameter gestalten und immer noch ändern, wenn nötig / EDIT: Bernhard hat das offenbar bereits gemacht).
Hallo zusammen,

ich will mich ja nicht unnötig einmischen, aber ich würde diese Zahl robuster auslegen, auch wenn man nicht so genau weiss, warum. Bedenkt, dass sogar der grösste Zwergplanet eine Bahnneigung von fast 45° hat, was man eigentlich nicht erwartet hat.

Ist das denn sehr teuer (in Form von Rechenzeit), wenn man einen grösseren Wert verwendet ?


Freundliche Grüsse, Ralf

Bernhard
08.08.2012, 17:18
Ist das denn sehr teuer (in Form von Rechenzeit), wenn man einen grösseren Wert verwendet ?
Hi Ralf,

ich finde das ist ein guter Tip, denn bei meiner Programmversion (ich bin jetzt mal ein wenig gemein) hat der Streuwinkelbereich keinen Einfluss auf die Rechenzeit.

Was mich momentan ebenfalls sehr interessiert, ist die Anzahl an Abstürzen der Trümmer auf den Planeten. Eine triviale, geschwindigkeitsunabhängige Absturzbedingung läßt sich recht schnell implementieren und ich will jetzt einfach mal wissen, wieviele Abstürze die Simulation bei beispielsweise 1000 Jahren berechnet und wo diese stattfinden.
Gruß

SRMeister
08.08.2012, 18:43
ich versteh grad irgendwie nich, was eigentlich das Problem ist? Ich hatte selbst auch vorgeschlagen, wieder 2 Listen einzuführen. Ich habe auch angeboten das zu tun. Gut, du hast durch rückspielen einer alten Version dieses Problem gelöst aber die Zusammenarbeit beendet. Nur warum??? Da guck ich jetzt doof in die Röhre.

Okay, dann wird es eben keine Version kompatibel mit Konsole und GUI mehr geben...
Tolles Arbeiten...

ralfkannenberg
08.08.2012, 18:56
Okay, dann wird es eben keine Version kompatibel mit Konsole und GUI mehr geben...
Tolles Arbeiten...
Hallo SRMeister,

aus der Praxis (Industrie, Banken, etc.) weiss man, dass zwischen dem Erstellen einer ersten korrekt lauffähigen Software und ihrer Eingliederung in die vorhandene Systemlandschaft mindestens ein Faktor (!!) 10 an zusätzlichem Aufwand erforderlich ist.

Also: bitte nicht frustiert sein - bei einer gemeinsamen Arbeit, bei der die Leute nicht am selben Schreibtisch sitzen, müsst Ihr mit ähnlichen Faktoren rechnen.

Von meiner Seite einfach nur ein riesiges Kompliment, ich kann mich nicht erinnern, dass man sich hier im astronews in den vergangenen 7 Jahren die Mühe für ein solches Projekt gemacht hat !


Freundliche Grüsse, Ralf

Bernhard
08.08.2012, 22:16
Eine triviale, geschwindigkeitsunabhängige Absturzbedingung läßt sich recht schnell implementieren und ich will jetzt einfach mal wissen, wieviele Abstürze die Simulation bei beispielsweise 1000 Jahren berechnet und wo diese stattfinden.

Ganz so einfach ist das leider doch nicht und zwar aus dem folgenden Grund. Der Vorteil des Bulirsch-Stoer-Algorithmus ist gerade die große Schrittweite (bei uns 1 Tag). Da aber auf triviale Weise (Abstand) ein Absturz nur nach jedem Simulationsschritt abgefragt werden kann, verliert man hier leider viele interessante Informationen.

Mit der aktuellen Konsolenversion kann man in akzeptabler Zeit (und nach auskommentieren der TimeStamp-Ausgabe am Ende) 1000 Jahre mit 10 Testkörpern bereits simulieren. Bei mir ergab sich dabei folgendes Ergebnis: Ein Testkörper hat das Sonnensystem komplett verlassen. Mehrere haben sich im Erdorbit angesammelt. Der Rest verteilte sich auf das innere Sonnensystem bis etwa 7 AE. Interessant finde ich solche Ergebnisse allemal, allerdings lief das aus Zeitgründen mit einer für mich übersichtlicheren Version ohne verlinkte Listen, dafür aber mit Arrays aus Objekten.

Bynaus
09.08.2012, 07:41
Ich kann und will mich nicht in die Details eurer Diskussion bzw. eures Streits einmischen, aber ich glaube, es wäre absolut sinnlos, zwei verschiedene Simulationen zu dem einen Thema aufzubauen - dann machen wir die Arbeit doppelt, die Kontrolle bzw. das Review der Ergebnisse - und damit das zu erwartende Niveau - kann nur abnehmen. Ich kann euch nicht sagen, wie ihr dies lösen könnt, aber ich würde vorschlagen, mit einer klaren Rollenverteilung anzufangen.

Zurück zum eigentlichen Inhalt: Diese ersten Ergebnisse, Bernhard, klingen qualitativ schon mal ganz okay, dh, in etwa das, was man wohl erwarten würde. Die wichtigste Frage dazu muss lauten: können wir ihnen trauen? Dafür wird eine sehr sorgfältige Dokumentation nötig sein.

Ich würde mir das Programm (nicht unbedingt den Code, zu diesem Zeitpunkt) gerne mal auch selbst ansehen. Wie genau muss ich vorgehen?

ismion
09.08.2012, 09:54
Ich schließe mich Bynaus an und halte mich raus aus Eurer "Diskussion" :D Genauso sehe ich aber auch, dass Ihr das Kriegsbeil begraben müsst, weil das sonst nach hinten losgeht.

@Bernhard: Du sagtest, du hättest eine Konsolenversion. Ist die Linux-Geeignet oder ist das eher sowas wie DOS/Windows-Eingabeaufforderung? ;)

Wenn die Linux-basiert ist, kann ich gerne simulieren, weil mich das auch interessiert, was rauskommt. Ich habe (unseren Cluster ausgenommen) 6 CPU zur Verfügung. Läuft das Programm parallel via MPI?Wie lange läuft eine Simulation im Schnitt?

Beste Grüße,
Basti

Bernhard
09.08.2012, 10:09
@Bernhard: Du sagtest, du hättest eine Konsolenversion. Ist die Linux-Geeignet oder ist das eher sowas wie DOS/Windows-Eingabeaufforderung? ;)

Wenn die Linux-basiert ist, kann ich gerne simulieren, weil mich das auch interessiert, was rauskommt. Ich habe (unseren Cluster ausgenommen) 6 CPU zur Verfügung. Läuft das Programm parallel via MPI?Wie lange läuft eine Simulation im Schnitt?
Hallo Basti,

die Konsolenversion enthält nur noch marginale Windows-, resp. Microsoft-Technik und genau so wollte ich das auch haben, weil die Rechenzeiten mit den erweiterten Anforderungen natürlich sofort nach oben gehen.

Deswegen schlage ich folgendes vor:

@SRMeister: Darf ich meine aktuelle Version auf Deinen svn-Space aufspielen? Bynaus und Basti könnten sich dann den Quelltext ansehen und ausgiebig testen. Das würde das Projekt sehr voranbringen und auf das GUI können wir dabei momentan ohne Probleme verzichten. So etwas hält halt momentan eher auf. Idealerweise sollte die Simulation also komplett ohne grafischen Overhead auf einem Linux-Cluster laufen. So etwas wäre auch für mich interessant. Man müsste nur prüfen mit welchem Zusatzaufwand das zu machen wäre.
Gruß

Bernhard
09.08.2012, 10:40
Ich würde mir das Programm (nicht unbedingt den Code, zu diesem Zeitpunkt) gerne mal auch selbst ansehen. Wie genau muss ich vorgehen?
Hallo Bynaus,

schicke mir am besten eine PN mit Deiner eMail-Adresse. Ich schicke Dir den Code dann privat zu, weil der momentan nur auf meiner privaten Festplatte existiert.

SRMeister
10.08.2012, 14:02
Hallo SRMeister,

ich werde dann mal versuchen ein eigenes svn-Projekt bei sourceforge.net aufzumachen. Langfristig sehe ich da momentan keine Alternative dazu.
Gruß
okay

Dass die Windowsversion das Programm langsamer machen würde, oder Overhead erzeugen würde, ist totaler quatsch. Die DLL wird genausoschnell ausgeführt, wie die .exe Datei. Das konnte ich mit dem Intel Profiler und dem eingebauten Benchmark beweisen.


Ich kann und will mich nicht in die Details eurer Diskussion bzw. eures Streits einmischen, aber ich glaube, es wäre absolut sinnlos, zwei verschiedene Simulationen zu dem einen Thema aufzubauen
stimmt natürlich. Ich will/kann auch keine eigenen Integratoren erstellen, deswegen werde ich vorhandene, in der wissenschaftlichen Umgebung gut dokumentierte Software einsetzen.
Ich bin momentan noch auf der Suche nach einem passenden Ansatz, aber es kristallisiert sich bereits heraus, dass es ein symplektischer Integrator werden wird. das SWIFT Paket sieht da sehr vielversprechend aus. Darauf bauen viele wissenschaftlichen Arbeiten zum Thema Sonnensystem, Oortsche Wolken usw. auf.

Mit der letzten lauffähigen GUI-Version und den letzten von Bernhard veröffentlichten Änderungen habe ich momentan eine lauffähige Version. Damit habe ich ein Video erstellt, wo 2500 Fragmente simuliert werden. Habe es bei YouTube hochgeladen.

www.youtube.com/watch?v=RUd06bwsPMw (Am besten mit 480p Auflösung ansehen)

Grüße

ralfkannenberg
10.08.2012, 14:24
Dass die Windowsversion das Programm langsamer machen würde, oder Overhead erzeugen würde, ist totaler quatsch. Die DLL wird genausoschnell ausgeführt, wie die .exe Datei. Das konnte ich mit dem Intel Profiler und dem eingebauten Benchmark beweisen.
Hallo SRMeister,

nichts gegen schöne Benchmarks, aber hast Du es mal mit einem ganz konkreten Loadtest ausprobiert ? - Ich habe schon viel zu dem Thema gehört, aber wenn früher die Applikation, die ich programmieren sollte, zeitkritisch wurde, habe ich meist zuvor einen kleinen Loadtest gemacht, um zu schauen, welche Befehle mehr Zeit brauchen, und dann entsprechend programmiert. Mit sowas konnte man früher locker einen Faktor 5 - 10 herausholen.

Ok, ich habe zweimal ein Programm programmiert, dessen Laufzeit über 100000x länger als das Alter des Universums war - da helfen solche Tricks dann natürlich auch nicht mehr weiter.


Freundliche Grüsse, Ralf

SRMeister
10.08.2012, 15:00
das ist denke ich genau das, was ich mit dem Profiler gemacht habe.
Da siehst du genau, welche Funktion Zeit kostet, sogar bis hin zu einzelnen CPU-Instructions.
Es sind natürlich die mathematischen Operationen. Besonders die Division und noch mehr das Wurzelziehen.
Ich hatte mich damit ausgiebig beschäftigt. Es gibt von Intel genaue Anleitungen, wo die CPU Zyklen pro Befehl aufgeführt sind. Das was man dort findet, habe ich 1:1 im Profiler widergefunden.
Da konnte ich aber auch durch Assembler viel Zeit rausholen.

SRMeister
10.08.2012, 16:48
hier nochmal ein anderes Video, wo nur die ersten 2 Jahre zu sehen sind, dafür 5000 Partikel.
http://www.youtube.com/watch?v=O2h1FjA24uo

Bernhard
10.08.2012, 22:02
Dass die Windowsversion das Programm langsamer machen würde, oder Overhead erzeugen würde, ist totaler quatsch.
Das hat bei genauem Hinsehen auch niemand behauptet und mit Overhead wird gerne nur zusätzlicher Code bezeichnet.


www.youtube.com/watch?v=RUd06bwsPMw
Das sieht schonmal ganz interessant aus. Wieviel Jahre werden hier eigentlich simuliert?

SRMeister
10.08.2012, 22:08
für 350 Jahre(Schrittweite 8Tage) hat es etwa 6h gedauert, aber in das Video sind nur die ersten 65 Jahre eingeflossen, da sich dann wenig verändert. Man sieht nur dass die Exzentrizität runtergeht, durch Wechselwirkung mit den inneren Planeten.
Kollision habe ich ja nicht implementiert.

Was man aber schlussfolgern kann, ist, dass die Dichte der Objekte im Innern Sonnensystem viel größer ist und da dort auch mehr Planeten pro Volumen sind, wird wohl ein Großteil dort enden.

Viel mehr wird man wahrscheinlich mit dem BS nicht nachweisen können, vermute ich. Alles Andere spielt sich auf ganz anderen Zeitskalen ab.

Bynaus
13.08.2012, 08:29
hier nochmal ein anderes Video, wo nur die ersten 2 Jahre zu sehen sind, dafür 5000 Partikel.
http://www.youtube.com/watch?v=O2h1FjA24uo

Oh wow, das sieht schon mal sehr gut / interessant aus! Scheint so, als ob Merkur eine Menge Treffer abbekommt.... (zählst du die bereits? vermutlich nicht. EDIT: Nein, noch nicht.) Auch die Erde (natürlich), auch hier wäre es sehr interessant zu wissen, wieviel am Ende auf die Erde zurückfällt, weil das - wenn es sehr viel Masse ist - natürlich auch einen Einfluss auf isotopische Ähnlichkeit zwischen Erde und Mond haben wird (je mehr davon die Erde absorbiert, desto ähnlicher wird sie dem Mond in der Isotopenzusammensetzung, wobei der Effekt eher klein ist).

Ich glaube, die Kleinplaneten (insb jene jenseits des Neptuns) und Monde kann man bei der Berechnung weglassen (höchstens als mögliche "Ziele" für Trümmerteile drinbehalten). Das sollte die Simulation schneller machen und dürfte kaum einen Unterschied auf die Teilchenbahnen haben.


Viel mehr wird man wahrscheinlich mit dem BS nicht nachweisen können, vermute ich. Alles Andere spielt sich auf ganz anderen Zeitskalen ab.

Okay. Wäre eine mögliche Lösung dieses Problems, dass man am Anfang nur kleine Zeitschritte simuliert, die dann mit der Zeit (z.B. in umgekehrter Abhängigkeit der Anzahl Partikel, die noch da sind) immer grösser werden?

Für den Asteroidengürtel könnte man, damit man nicht hunderttausende von Asteroiden simulieren muss, vielleicht einfach davon ausgehen, dass wenn sich das Teilchen im Bererich zwischen a = 1.9 und a = 3.5 befindet (evtl noch Wahrscheinlichkeiten für verschieden dichte Bereiche des Gürtels, und auch Berücksichtigung der Bahnneigung des Teilchens), es eine bestimmte Wahrscheinlichkeit gibt, mit der es von einem "nur gedachten" Asteroiden absorbiert (und damit aus der Simulation entfernt) wird. Diese Wahrscheinlichkeit wäre dann aus der Asteroidendichte zu ermitteln (mein Job :) ). Wäre das eine Möglichkeit?

Weiter: wenn man an das Endprodukt denkt, dann wären sicher auch Exzentrizität vs. Inklination / Grosse Halbachse der einzelnen Teilchen interessant. Wäre z.B. toll, wenn man den zeitlichen Verlauf nicht nur in der xy-Ebene des Sonnensystems plotten könnte, sondern auch in diesen "Parameter-Ebenen".

Ich bin wirklich begeistert davon. Ich kann mir gut vorstellen, dass da etwas verwertbares (publizierbares) rauskommt...

ralfkannenberg
13.08.2012, 09:35
Hallo Bynaus,


Scheint so, als ob Merkur eine Menge Treffer abbekommt....
was ist mit der Venus ? Die hat ja mehr Masse als der Merkur, ist aber natürlich weiter von der Sonne entfernt, dafür zeitweise deutlich näher an der Erde/Mond-Trümmerwolke.


Ich glaube, die Kleinplaneten (insb jene jenseits des Neptuns) und Monde kann man bei der Berechnung weglassen (höchstens als mögliche "Ziele" für Trümmerteile drinbehalten). Das sollte die Simulation schneller machen und dürfte kaum einen Unterschied auf die Teilchenbahnen haben.
Ich denke auch, dass vor allem der Ziel-Aspekt interessant sein könnte, wobei ich nicht "glaube", dass die "Erde-Mond-Trümmerwolke" einen signifikanten Beitrag zu den mittel- und langfristigen Kometen im Kuipergürtel oder gar der Oort'schen Wolke leisten dürfte. Hier könnte man vielleicht auch anstelle einer Simulation einfach nur eine Wahrscheinlichkeit ausrechnen, ob es überhaupt grössere ("entdeckbare") Körper da draussen gibt, die genügend Anfangsenergie hatten, da überhaupt hinzukommen ?


Für den Asteroidengürtel könnte man, damit man nicht hunderttausende von Asteroiden simulieren muss, vielleicht einfach davon ausgehen, dass wenn sich das Teilchen im Bererich zwischen a = 1.9 und a = 3.5 befindet (evtl noch Wahrscheinlichkeiten für verschieden dichte Bereiche des Gürtels, und auch Berücksichtigung der Bahnneigung des Teilchens), es eine bestimmte Wahrscheinlichkeit gibt, mit der es von einem "nur gedachten" Asteroiden absorbiert (und damit aus der Simulation entfernt) wird. Diese Wahrscheinlichkeit wäre dann aus der Asteroidendichte zu ermitteln (mein Job :) ). Wäre das eine Möglichkeit?
Zumindest für den ersten Approach. Für den zweiten und dritten vermutlich auch völlig ausreichend.



Weiter: wenn man an das Endprodukt denkt, dann wären sicher auch Exzentrizität vs. Inklination / Grosse Halbachse der einzelnen Teilchen interessant. Wäre z.B. toll, wenn man den zeitlichen Verlauf nicht nur in der xy-Ebene des Sonnensystems plotten könnte, sondern auch in diesen "Parameter-Ebenen".
Das ist vermutlich ein erheblicher Mehraufwand aus Sicht des Rechenaufwandes, aber vermutlich auch unbedingt nötig. Hier sollte man m.E. nicht sparen.


Ich kann mir gut vorstellen, dass da etwas verwertbares (publizierbares) rauskommt...
Denke ich auch, wird aber sicherlich noch etwas Aufwand benötigen ...

Blöde Frage: Könnte es die Rechenzeit optimieren, wenn man den Mars einfach als grossen Planetoiden am inneren Ende eines erweiterten Planetoidengrütels simuliert und beispielsweise nur in der Nähe der Opposition als eigenständigen Planeten ?


Freundliche Grüsse, Ralf

Bernhard
13.08.2012, 09:59
Für den Asteroidengürtel könnte man, damit man nicht hunderttausende von Asteroiden simulieren muss, vielleicht einfach davon ausgehen, dass wenn sich das Teilchen im Bererich zwischen a = 1.9 und a = 3.5 befindet (evtl noch Wahrscheinlichkeiten für verschieden dichte Bereiche des Gürtels, und auch Berücksichtigung der Bahnneigung des Teilchens), es eine bestimmte Wahrscheinlichkeit gibt, mit der es von einem "nur gedachten" Asteroiden absorbiert (und damit aus der Simulation entfernt) wird. Diese Wahrscheinlichkeit wäre dann aus der Asteroidendichte zu ermitteln (mein Job :) ). Wäre das eine Möglichkeit?
Hallo Bynaus,

ja das wäre eine Möglichkeit, um den Code weiter zu entwickeln. Man bräuchte dabei voraussichtlich aber eine Wahrscheinlichkeit pro Weglänge. Je nachdem, wie lange (=Weglänge) sich das Trümmerteil in dem kritischen Bereich aufhält findet also ein Einfang statt oder nicht. Der Computer braucht hier eben eine konkrete Ja-Nein-Aussage, die sich aber vielleicht über die Weglänge realisieren ließe?
Gruß

Bynaus
13.08.2012, 10:08
was ist mit der Venus ?

Ja natürlich, auch die Venus bekommt einiges ab. Wegen der mühsamen Angwohnheit der Venus, ihre Oberfläche alle paar hundert Millionen Jahre mit Lavaströmen zu fluten, ist das für Planetenforscher allerdings nur am Rande interessant. :)


Könnte es die Rechenzeit optimieren, wenn man den Mars einfach als grossen Planetoiden am inneren Ende eines erweiterten Planetoidengrütels simuliert und beispielsweise nur in der Nähe der Opposition als eigenständigen Planeten ?

Für mich schwierig zu beurteilen - aber ich kann mir vorstellen, dass dem Mars eine wichtige Rolle beim "Einbremsen" der Trümmer in den Asteroidengürtel zukommt. Ich muss mal nachschauen, wie das Bottke et al. (2006) so gelöst haben.


Der Computer braucht hier eben eine konkrete Ja-Nein-Aussage, die sich aber vielleicht über die Weglänge realisieren ließe?

Ja, das klingt gut. Ich denke, für den Anfang dürfte eine uniforme Wahrscheinlichkeit für den ganzen Gürtel vollauf genügen. Wenn sich das als wichtig erweist, kann man immer noch mehr ins Detail gehen.

Bernhard
13.08.2012, 10:15
Ja, das klingt gut. Ich denke, für den Anfang dürfte eine uniforme Wahrscheinlichkeit für den ganzen Gürtel vollauf genügen. Wenn sich das als wichtig erweist, kann man immer noch mehr ins Detail gehen.
Hallo Bynaus,

ich habe gerade gesehen, dass ich auf dem svn-repository keine Schreibrechte mehr habe. SRMeister hast mit mir scheinbar doch recht große Probleme. Ist natürlich die Frage, wie man unter solchen Voraussetzungen dann weiter machen soll?
Mit Gruß

Bernhard
13.08.2012, 10:18
für den Anfang dürfte eine uniforme Wahrscheinlichkeit für den ganzen Gürtel vollauf genügen.
Sobald das Trümmerteil in den Asteroidengürtel eintritt wird es also sofort eingefangen?

Bynaus
13.08.2012, 11:30
Sobald das Trümmerteil in den Asteroidengürtel eintritt wird es also sofort eingefangen?

Es kann eingefangen werden, ja. Sagen wir, der Weg durch die Zone des Asteroidengürtels ist 3 AU lang, und die Absorbtionswahrscheinlichkeit im Gürtel wäre z.B. 0.01/AU. Dann sollten nur 3 von 100 Teilchen, die den Asteroidengürtel auf einer solchen Bahn (3 AU Weglänge) durchqueren, absorbiert werden.

Was eure Differenzen angeht, da kann und will ich nicht urteilen. Ich kann höchstens anbieten, in einem Chat oder Videochat zu moderieren, um eine Lösung zu finden, die für beide Seiten stimmt.

ralfkannenberg
13.08.2012, 11:33
ich habe gerade gesehen, dass ich auf dem svn-repository keine Schreibrechte mehr habe. SRMeister hast mit mir scheinbar doch recht große Probleme. Ist natürlich die Frage, wie man unter solchen Voraussetzungen dann weiter machen soll?
Hallo Bernhard,

mal ein Blick in die vielgepriesene Praxis: in der Firma "verlieren" wir andauernd Schreibrechte, obgleich eigentlich jeder alles können sollte. Die Verantwortlichen haben noch nicht herausgefunden, woran das liegt, und bauen dann die Schreibrechte wieder mühsam ein.

Dann kommt irgendein Manager daher und kann 5 Euro pro Jahr sparen, wenn man auf einen anderen Server umzieht, worauf alle höheren Manager begeistert sind; natürlich werden die User nicht informiert und die Schreibrechte sind dann wieder weg ...

Ich würde deswegen erst einmal mit SRMeister per PN abklären, was da möglicherweise gelaufen sein könnte.


Freundliche Grüsse, Ralf

Bernhard
13.08.2012, 13:57
Ich würde deswegen erst einmal mit SRMeister per PN abklären, was da möglicherweise gelaufen sein könnte.
Hallo Ralf,

SRMeister antwortet eh nicht auf meine öffentlichen Fragen also was soll's. Es ist sein Thema und er will sowieso einen anderen Integrator einbauen.
Gruß

SRMeister
13.08.2012, 15:19
Ich glaube, die Kleinplaneten (insb jene jenseits des Neptuns) und Monde kann man bei der Berechnung weglassen (höchstens als mögliche "Ziele" für Trümmerteile drinbehalten). Das sollte die Simulation schneller machen und dürfte kaum einen Unterschied auf die Teilchenbahnen haben.
Das kann man machen. Man könnte sie möglicherweise auf eine Keplerbahn reduzieren, so dass der Rechenaufwand sehr klein wird und sie als Ziele nach wie vor vorhanden bleiben.



Okay. Wäre eine mögliche Lösung dieses Problems, dass man am Anfang nur kleine Zeitschritte simuliert, die dann mit der Zeit (z.B. in umgekehrter Abhängigkeit der Anzahl Partikel, die noch da sind) immer grösser werden? Die Integrationszeit braucht man gewöhnlich nicht auf diese Weise variabel zu machen. Durch die immer weniger werdenden Teile wird die Simulation sowieso immer schneller. Nachwievor müssen wir das Problem lösen, dass man einen Integrator braucht, der "close encounters" simulieren kann, also die Zeitschritte senkt oder die Genauigkeit erhöht, sobald sich 2 Körper nahe kommen. Sonst wird es schwer, die Kollision zuverlässig zu erkennen.
Allerdings glaube ich momentan nicht, dass wir das Problem mit dem einfachen BS-Integrator lösen können. Wir müssen zwangsläufig auf einen symplektischen Integrator gehen. Da ist wiederum auch nicht jeder geeignet, da gibt es welche die close encounters beherrschen und andere, die es nicht können. Dann ist das Problem, dass die meisten in Fortran-77 programmiert sind oder nicht open source vorliegen. Ich hatte mir einen sehr vielversprechenden ausgesucht, der genau unsere gewünschten Eigenschaften hat, sehr schnell ist, aktuell ist (2010) und der einigermaßen gut dokumentiert ist. Dieser nennt sich SCATR und ist eine Weiterentwicklung von SWIFT, der mMn. genauso verbreitet wie Mercury ist.
Link (http://arxiv.org/abs/1011.3830)
Die Autoren schreiben, man könne den Quellcode per EMail anfragen. Das habe ich vor einer Woche getan, aber leider keine Antwort erhalten. Möglicherweise fehlt mir dazu der wissenschaftliche Hintergrund (zb uni-email oder Mastertitel) um ernst genommen zu werden.
Vielleicht kannst du den Teil übernehmen, Bynaus, und versuchen, Kontakt aufzubauen.
Auch wenn der dann in Fortran programmiert ist, schätze ich ein, dass man ihn an ein eigenes Programm anschließen kann. Denn der ist so perfekt, dass wir wohl am Code selbst keine Änderungen mehr vornehmen müssten.


Weiter: wenn man an das Endprodukt denkt, dann wären sicher auch Exzentrizität vs. Inklination / Grosse Halbachse der einzelnen Teilchen interessant. Wäre z.B. toll, wenn man den zeitlichen Verlauf nicht nur in der xy-Ebene des Sonnensystems plotten könnte, sondern auch in diesen "Parameter-Ebenen".
Genau das möchte ich als nächstes tun. Am GUI habe ich schon entsprechende Funktionen vorbereitet. Die örtliche Verteilung ist ja letztendlich garnicht so interessant.


Das ist vermutlich ein erheblicher Mehraufwand aus Sicht des Rechenaufwandes, aber vermutlich auch unbedingt nötig. Hier sollte man m.E. nicht sparen.
Das denke ich nicht, das ist ein relativ kleiner Mehraufwand. Wir hatten auch schonmal eine Version, die die Exzentrizität berechnete. Da ist der Unterschied kaum spürbar, da die Werte wie Sonnenabstand sowieso während der Berechnung anfallen.



Blöde Frage: Könnte es die Rechenzeit optimieren, wenn man den Mars einfach als grossen Planetoiden am inneren Ende eines erweiterten Planetoidengrütels simuliert und beispielsweise nur in der Nähe der Opposition als eigenständigen Planeten ?
Ich würde eigentlich auf solche eigenen Varianten lieber verzichten, auch wenn es etwas helfen würde. Der Grund ist, wie Bynaus schonmal sagte, dass man sonst die komplette Software sozusagen wissenschaftlich verifizieren müsste, um den Einsatz in einer wissenschaftlichen Arbeit zu rechtfertigen.

Ich stehe momentan auf dem Punkt, wo ich gerne nur vorhandene, zum Problem passende Integratoren verwenden möchte.
Das hat wenig mit der Differenz zu Bernhard zu tun oder damit, dass ich selbst wohl keinen symplektischen Integrator programmieren könnte.



ich habe gerade gesehen, dass ich auf dem svn-repository keine Schreibrechte mehr habe.
Ja das ist richtig. Ich hatte deine Aussage, dass du eine eigene Variante weiterentwickeln möchtest, dementsprechend interpretiert. Ich wollte auf dem SVN space keinen Versionsstreit starten sondern einfach "meine Variante" weiter entwickeln.


SRMeister hast mit mir scheinbar doch recht große Probleme. Ist natürlich die Frage, wie man unter solchen Voraussetzungen dann weiter machen soll?
Mit Gruß
Genau die Frage hatte ich mir und hier öffentlich, bereits letzte Woche gestellt. Ich hatte dann aber den Eindruck gewonnen, du suchst keine Lösung die eine Zusammenarbeit beinhaltet. Wie ich geschrieben hatte, war ich da sehr verwirrt und wusste/weis bis jetzt eigentlich nicht so richtig, wie es nun weitergehen soll.
Ich habe die Schreibrechte nur deaktiviert, um zu vermeiden, dass du dort meine Arbeit wieder überschreibst.



SRMeister antwortet eh nicht auf meine öffentlichen Fragen also was soll's.
Ich habe auf deine Fragen doch immer geantwortet. Ich war am Wochenende nicht am PC, deswegen erst jetzt meine Rückmeldung.
Ich habe mehrere Vorschläge gemacht, wie man zusammen weiterarbeiten könnte, wie man die programmiertechnischen Probleme lösen könnte. Ich habe versucht, meine Änderungen am Programm zu rechtfertigen und angeboten, sie teilweise rückgängig zu machen, um eine gemeinsame Basis zu schaffen. Darauf hast du aber garnicht reagiert.

Gruß

ralfkannenberg
13.08.2012, 15:49
Ich habe auf deine Fragen doch immer geantwortet. Ich war am Wochenende nicht am PC, deswegen erst jetzt meine Rückmeldung.
Ich habe mehrere Vorschläge gemacht, wie man zusammen weiterarbeiten könnte, wie man die programmiertechnischen Probleme lösen könnte. Ich habe versucht, meine Änderungen am Programm zu rechtfertigen und angeboten, sie teilweise rückgängig zu machen, um eine gemeinsame Basis zu schaffen. Darauf hast du aber garnicht reagiert.
Ach herrje,

genau wie früher bei uns in der Firma.

Ihr wisst schon, dass man zuerst eine Spezifikation erstellt, diese dann reviewt und erst dann mit der Programmierung anfängt, nicht wahr ? - Wenn man schon erste Ergebnisse braucht, kann man ja - parallel dazu - einen Protoypen erstellen, der nur die Basisfunktionen enthält.

Und ja: den Prototypen löscht man dann wieder vor der richtigen Programmierung. Die Erfahrung zeigt, dass die endgültige Software dann bessere Qualität aufweist als wenn man den Prototypen solange "umbiegt", bis der auch ungefähr das kann was man will.


Freundliche Grüsse, Ralf

Bernhard
13.08.2012, 16:12
Ich habe die Schreibrechte nur deaktiviert, um zu vermeiden, dass du dort meine Arbeit wieder überschreibst.
Hallo SRMeister,

genau das hatte ich mir so auch gedacht und es ist auch völlig OK. Dazu folgenden Vorschlag. Ich habe mittlerweile eine Version ohne verlinkte Listen, dafür aber mit zwei Arrays jeweils für die Planeten (Objekte mit Masse) und die Testkörper (Objekte mit Masse). Die Rechnung selbst läuft unverändert, so dass es sich um rein formale Änderungen handelt. Weitere Features (adaptive Stepsize, Kollisionen) sähe ich gerne in dieser, meiner Meinung nach, deutlich übersichtlicheren Version.

Wir haben jetzt also drei Möglichkeiten:

1) Ich mache ein eigenes Projekt auf sourceforge.net auf und verzichte dabei auf jegliche Visualisierung (das ist mir einfach zu viel Arbeit)
2) Ich spiele diese Version in Dein Repository und wir versuchen gemeinsam eine Schnittstelle zu Deinem GUI (das mir sehr gut gefällt, wegen der Youtube-Clips) zu entwickeln.
3) Ich könnte Dir als Vorbereitung diese neue Version ebenfalls per eMail zuschicken. Es sind allerdings etliche Änderungen an den Methoden der Klasse CNewton dabei.

Ich persönlich würde Version 2 oder 3 bevorzugen.
Gruß

Bernhard
13.08.2012, 16:16
so dass es sich um rein formale Änderungen handelt
Die zwei Variablen m_source und m_target habe ich ebenfalls ausgebaut, weil es bei den Planeten selbst ja definitiv nicht um Streuprozesse geht. Das ständige "Buffern" der Bedeutungen war für mich einfach zu anstrengend.
Gruß

Bynaus
13.08.2012, 16:53
Vielleicht kannst du den Teil übernehmen, Bynaus, und versuchen, Kontakt aufzubauen.

Hab ich soeben getan. Ich glaube, die E-Mail-Adresse im Abstract war nicht mehr aktuell.

Ich würde vorschlagen, dass wir uns (ihr euch) als nächstes auf einen Kollisionsdetektions-Mechanismus einigt, der dann umgesetzt wird. Ich würde Zähler für jeden grösseren Körper vorschlagen, sowie einen speziellen Zähler (wie oben vorschlagen) für den Asteroidengürtel.

UMa
13.08.2012, 17:15
Hallo,

für den Kollisionsdetektions-Mechanismus würde ich folgendes vorschlagen:

Für jeden kleineren Körper nach jedem Zeitschritt prüfen, ob er in der Nähe eines größeren Körpers, z.B. näher als 10*Radius, ist.
Dann (nur für die wenigen Nahen!) die geringste Distanz unter Vernachlässigung anderer Körper analytisch ausrechnen. Wenn kleiner als Radius dann Kollision, sonst nicht.

Das umgeht das Problem, dass Körper während eines Zeitschrittes durch andere durchfliegen können und ermöglicht daher größere Zeitschritte.

Grüße UMa

Bernhard
13.08.2012, 17:50
Hallo UMa,


Für jeden kleineren Körper nach jedem Zeitschritt prüfen, ob er in der Nähe eines größeren Körpers, z.B. näher als 10*Radius, ist.
da ist der Testkörper dann aber eventuell schon durch den Planeten durchgeflogen.


Dann (nur für die wenigen Nahen!) die geringste Distanz unter Vernachlässigung anderer Körper analytisch ausrechnen.
Das ginge nur, wenn man innerhalb dieser Zeitspanne das äußere Gravitationsfeld konstant hält. Oder gibt es da noch bessere Möglichkeiten?

UMa
14.08.2012, 10:50
Hallo Bernhard

da ist der Testkörper dann aber eventuell schon durch den Planeten durchgeflogen.
Spielt keine Rolle. Ist das unten ausgerechnete Perizentrum q kleiner oder gleich dem Radius des großen Körpers, wird eine Kollision für den großen Körper gezählt und der kleine Körper wird aus der Rechnung entfernt.

Das ginge nur, wenn man innerhalb dieser Zeitspanne das äußere Gravitationsfeld konstant hält. Oder gibt es da noch bessere Möglichkeiten?
http://de.wikipedia.org/wiki/Zweik%C3%B6rperproblem#Bahnkurve

Ich habe aber sowohl den Drehimpuls als auch die Energie durch die Masse des Testkörpers geteilt, so dass im Vergleich zu den Wikipediaformeln µ wegzulassen bzw. =1 zu setzen ist.

Falls Abstand eines kleinen Körpers zu einem großen so klein ist, dass die anderen vernachlässigt werden können, dann folgendes:

gegeben:

W Geschwindigkeitsvektor des kleinen Körpers relativ zum nahen großen Körper
R Ortsvektor des kleinen Körpers relativ zum nahen großen Körper
D.h. Geschwindigkeit und Ort des großen Körpers von denen des kleinen Körpers abziehen, so dass großer Körper im Ursprung ruht.
GM Gravitationskonstante*Masse des großen Körpers

Rechenweg:

(spezif.) Drehimpuls l ausrechnen, Kreuzprodukt, dann Länge des Vektors
l=|WxR|

(spezif.) Energie E ausrechnen
E=w^2/2 - GM/r
wobei w=|W| und r=|R| die Beträge des Geschwindigkeits- und Ortsvektors sind.

Bahnparameter p ausrechnen
p=l^2/(GM)

Exentrizität e ausrechnen
e^2 = 1 + 2 E l^2 / (GM)^2

Perizentum q, d.h. Distanz des kleinsten Abstandes
q=p/(1+e)

q mit Radius des großen Körpers vergleichen, fertig.

Grüße UMa

PS: Das kleine "L" für den Drehimpuls ist leider nur schwer von Betragsstrichen oder großem "i" (welches nicht auftritt) zu unterscheiden.

Edit: Ganz am Anfang wird natürlich für alle eine Kollision mit der Erde registriert. Das muss irgendwie abgefangen werden.

Bernhard
14.08.2012, 12:06
Edit: Ganz am Anfang wird natürlich für alle eine Kollision mit der Erde registriert. Das muss irgendwie abgefangen werden.
Hallo UMa,

das ist im Programm bereits berücksicht, bzw. entsprechend angelegt. Ich werde mir im Laufe der Woche dann Deine weiteren Vorschläge und Formeln ansehen. Weitere Änderungen werde ich voraussichtlich in die Links aus dem Nachbarthema http://astronews.com/forum/showthread.php?6248-Entstehung-des-Mondes einbauen (müsssen). Währenddessen kann sich SRMeister dann in Fortran einarbeiten.
Schönen Gruß

Bernhard
14.08.2012, 14:27
Für den Asteroidengürtel könnte man, damit man nicht hunderttausende von Asteroiden simulieren muss, vielleicht einfach davon ausgehen, dass wenn sich das Teilchen im Bererich zwischen a = 1.9 und a = 3.5 befindet (evtl noch Wahrscheinlichkeiten für verschieden dichte Bereiche des Gürtels, und auch Berücksichtigung der Bahnneigung des Teilchens), es eine bestimmte Wahrscheinlichkeit gibt, mit der es von einem "nur gedachten" Asteroiden absorbiert (und damit aus der Simulation entfernt) wird. Diese Wahrscheinlichkeit wäre dann aus der Asteroidendichte zu ermitteln (mein Job :) ). Wäre das eine Möglichkeit?
Hallo Bynaus,

mir ist gerade aufgefallen, dass man mit Hilfe des internen Zufallsgenerators auch Wahrscheinlichkeiten auswerten/implementieren könnte. Diese dürfte auch geschwindigkeitsabhängig sein. Mir schwebt da ein IN- und OUT-Flag vor, die registriert, wann das Trümmerteil den Asteroidengürtel betritt, bzw. verlässt. Eine Ja/Nein-Einfangwahrscheinlichkeit, wie von mir oben gefordert, wäre hier zu unrealistisch.
Gruß

SRMeister
14.08.2012, 15:48
Hallo UMa,
deine Methode klingt sehr interessant. Da hast du ungefähr einen wichtigen Teil des Wisdom-Holman symplektischen Integrators beschrieben.

Man mus halt bedenken, dass die Keplerlösung nur eine Näherung ist. Insbesondere besteht die korrekte Bahn beim Wisdom Holman aus einem Keplerterm und einem Störungsterm.
Gut, wenn man dieses Verfahren nur im Nahbereich eines Körpers anwenden will, ergibt sich das Problem, dass die beiden Körper miteinander wahrscheinlich garkeine Keplerbeziehung haben, dass heist der Störungsterm dominiert.
Eine Keplerbeziehung kommt meiner Meinung nach nur in Frage, wenn sich der Testkörper innerhalb der Hillsphäre aufhält.
Man müsste also die Zeitschritte so klein machen, dass die Hillsphäre nicht übersprungen wird.

Man braucht wahrscheinlich also doch adaptive Zeitschritte.

Wozu die Methode aufjedenfall brauchbar wäre: Man könnte die Bahnparameter in Bezug auf die Sonne zu jedem Zeitpunkt errechnen. Es wäre also nicht nötig, einen ganzen Umlauf zu warten um die Exzentrizität zu errechnen.

Bernhard, ich würde natürlich deine aktuelle Version wieder integrieren, aber man müsste sich vorher diesmal auf eine Schnittstelle einigen, woran sich dann auch jeder hält. Wenn das Okay ist, dann würde ich hier eine Schnittstelle vorschlagen. Wenn du damit einverstanden bist, müsstest du dich aber auch darum bemühen, die Gültigkeit der Schnittstelle zu bewahren. Alles Andere kannst du dann so machen, wie du willst.

UMa
14.08.2012, 17:23
Hallo SRMeister,

Man müsste also die Zeitschritte so klein machen, dass die Hillsphäre nicht übersprungen wird.

Man braucht wahrscheinlich also doch adaptive Zeitschritte.
ich bin davon ausgegangen, dass adaptive Zeitschritte bereits implementiert sind. Ohne adaptive Zeitschritte geht es nicht. Sonst werden die Streuprozesse bei nahen Vorbeigängen nicht richtig berücksichtigt.
Der Fehlerschätzer für die adaptive Zeitschritte wird bei der Extrapolation quasi kostenlos mitgeliefert. Wert beim letzten Extrapolationsschritt minus Wert beim vorletzten.

Die Keplerlösung darf natürlich nur angewendet werden, wenn die Störungen der anderen Körper (Sonne) vernachlässigt werden können. Also nur weit innerhalb der Hillsphäre.

Grüße UMa

Bernhard
14.08.2012, 20:11
Bernhard, ich würde natürlich deine aktuelle Version wieder integrieren, aber man müsste sich vorher diesmal auf eine Schnittstelle einigen, woran sich dann auch jeder hält. Wenn das Okay ist, dann würde ich hier eine Schnittstelle vorschlagen. Wenn du damit einverstanden bist, müsstest du dich aber auch darum bemühen, die Gültigkeit der Schnittstelle zu bewahren. Alles Andere kannst du dann so machen, wie du willst.
Da Schnittstellen von Wünschen und Anforderungen abhängig sind, halte ich solche Forderungen prinzipiell für mehr oder weniger sinnlos.

Bernhard
14.08.2012, 20:21
Der Fehlerschätzer für die adaptive Zeitschritte wird bei der Extrapolation quasi kostenlos mitgeliefert. Wert beim letzten Extrapolationsschritt minus Wert beim vorletzten.
Kann man für den Fehlerschätzer einen maximalen, oberen Wert angeben, der gerade noch zulässig ist? Falls ja, wie groß ist dieser Wert dann typischerweise?

ralfkannenberg
14.08.2012, 20:42
Da Schnittstellen von Wünschen und Anforderungen abhängig sind, halte ich solche Forderungen prinzipiell für mehr oder weniger sinnlos.
Hallo Bernhard,

wie möchtest Du ein komplexes Programm, an dem mehrere Leute mitwirken, ohne einen solchen Konsens realisieren ?

Aber natürlich, zunächst einmal müssen diese Wünsche und Anforderungen vereinbart werden. Und äh, ja - dann muss man das auch mal testen, ob es irgendwo einen Schnittstellenfehler gibt. Solche Fehler sind lästig und sehr langweilig, die sollte man möglichst frühzeitig eliminieren, so lange das noch ohne grossen Aufwand geht.


Freundliche Grüsse, Ralf

Bernhard
14.08.2012, 20:58
wie möchtest Du ein komplexes Programm, an dem mehrere Leute mitwirken, ohne einen solchen Konsens realisieren ?
Hallo Ralf,

natürlich habe ich nichts gegen Vereinbarungen und halte es auch für selbstverständlich solche Vereinbarungen so weit es eben geht einzuhalten. Es bleibt dabei halt nur mal wieder die Frage, was getan wird, wenn sich die Realität anders entwickelt? Hat man dann beleidigte Leberwürste, die man nur noch mit der Zange anfassen darf oder kann man dann noch vernünftig miteinander reden? Genau das will ich vorher wissen, weil ich bereits viel Arbeit in das Projekt gesteckt habe und nur noch bedingt bereit bin hier den "Pausenkasper" zu spielen.

Mit Terminen halte ich es übrigens genau, wie in der Softwareentwicklung üblich. Es wird das gemacht, was geht. Und über das, was nicht geht, regt man sich dann weiter nicht auf. Sollte Euch diese Einstellung zu unverbindlich sein, kann die anfallende Tipparbeit auch gerne jemand anderes erledigen.

ralfkannenberg
14.08.2012, 21:37
Es bleibt dabei halt nur mal wieder die Frage, was getan wird, wenn sich die Realität anders entwickelt?
Hallo Bernhard,

das ist doch der Normalfall.

Gehen wir gar nicht so weit: seit vergangenem Dienstag führe ich hier in der Firma einen Text-Zyklus durch, die Zeit ist extrem knapp bemessen und ich bin kaum einen Tag vor 22 Uhr nach Hause. Nicht selten bin ich am morgen eine Stunde früher gekommen. Ich habe mir wirklich den Ar*** aufgerissen und alles Private hinten angestellt. Deadline ist heute Mitternacht und vor wenigen Minuten war meine letzte Auswertung abgeschlossen. Schön.

Nur: Heute nachmittag habe ich erfahren, dass wir eine veraltete Version getestet haben und die neue Funktionalität, für die dieser test-Zyklus extra durchgeführt wurde, mit der neuen Version testen sollen. Diese wird morgen geliefert.

Kurz und gut: die neue Version mag fehlerfrei sein, aber sie hat sehr schlechte Qualität, denn sie kommt einen Tag zu spät. Und auf das Testsystem des Kunden kommt eine Version, die wir gar nicht getestet haben.

Und dafür habe ich während einer Woche und einem Tag mein Privatleben "ausgeschaltet".


Was ich sagen will: auch wenn es Euch nicht bewusst ist, aber an sich habt Ihr ein Luxusproblem.


Freundliche Grüsse, Ralf

ralfkannenberg
14.08.2012, 21:38
doppelt, kann gelöscht werden

SRMeister
14.08.2012, 23:08
Jo das ist echt nen riesiger Witz langsam...
also ich habe nur vorgeschlagen dass wir zuallererst mal auf dem Papier bzw hier im Forum eine gemeinsame Grundlage schaffen, wenn du das nicht willst, also hier darüber zu reden wie wir zusammenarbeiten können, dann lassen wir es eben sein.

Bernhard
14.08.2012, 23:12
Was ich sagen will: auch wenn es Euch nicht bewusst ist, aber an sich habt Ihr ein Luxusproblem.
Hallo Ralf,

vielen Dank für den kleinen Einblick in Deine aktuelle Welt. So etwas ist natürlich unterhaltsam, bringt uns hier aber nur bedingt weiter, denn Freiwilligenprojekte laufen doch unter anderen Spielregeln, als kommerzielle Projekte. Ich werde jetzt erst mal abwarten in welche Richtung sich das Projekt weiterentwickelt. Mein Parallelprojekt hat seine Wirkung zum Glück nicht verfehlt. So viel Freiheit muss sein.

Dir persönlich würde ich noch raten auch auf Deine Gesundheit zu achten. In manchen Firmen gibt es deswegen eine maximale Tagesstundenzahl, was sicher keine schlechte Regelung ist. Ich weiß aber auch, dass Arbeitnehmer mit Familie in gewisser Weise erpressbar sind.
Mit Gruß

Bernhard
14.08.2012, 23:14
dann lassen wir es eben sein.
Bin schon gespannt, wie lange Du diese Aussage durchhältst.

SRMeister
14.08.2012, 23:18
Willst du jetzt mit aller Macht hier nen Streit provozieren?

Bernhard
15.08.2012, 00:22
Willst du jetzt mit aller Macht hier nen Streit provozieren?
Kann man sich mit Dir eigentlich auch normal unterhalten oder sind da Witze generell verboten.

FrankSpecht
15.08.2012, 01:16
@Bernhard,
@SRMeister,
das, was ihr hier gerade abzieht, ist Kinderkram!
:mad:

PS: Und völlig kontraproduktiv.
Vielleicht sollte ich euch meinen C-Code von '99 zur Verfügung stellen: N-Körper mit OpenGL, allerdings Windows-GUI?
Habt ihr da vielleicht einen gemeinsamen Nenner?

FrankSpecht
15.08.2012, 02:15
@Bernhard,
@SRMeister,
das, was ihr hier gerade abzieht, ist Kinderkram!
:mad:

PS: Und völlig kontraproduktiv.
Vielleicht sollte ich euch meinen C-Code von '99 zur Verfügung stellen: N-Körper mit OpenGL, allerdings Windows-GUI?
Habt ihr da vielleicht einen gemeinsamen Nenner?
http://www.frank-specht.de/astronomy/n-body.jpg

Bynaus
15.08.2012, 09:58
Ich habe von Nathan Kaib die SCATR-Dateien zugeschickt bekommen. Wer sie weitergeleitet haben will, soll mir eine PN mit E-Mail-Adresse schicken bzw. sich sonstwie melden.

Ich werd jetzt mal versuchen, eine grobe Spezifikation des Programms aus meiner Sicht festzuhalten, soweit wie sich die Diskussion bisher entwickelt hat. Mir ist bewusst, dass einiges (oder vieles) davon schon umgesetzt wurde, und ich werde mich nicht in technischen Details wie Arrays etc. verlieren - das müssen jene machen, die es direkt betrifft.

Wir wollen:
- Ein Programm, das eine grosse (definierbare) Anzahl masseloser Teilchen mit Anfangsgeschwindigkeiten 1 - 1.4 v_esc auf ihrem Weg durch das Sonnensystem (Gravitation der 8 Planeten und der Sonne) verfolgen kann, wobei man vorgeben kann, in welchem Raumwinkel die Teilchen von der Erde "abgestrahlt" werden.
- Weiter sollen auch einige massive Teilchen simuliert werden.
- Es soll ein Zähler implementiert werden, der Kollisionen zwischen allen Teilchen und den Planeten (inkl. der Erde) entdeckt (und die Teilchen aus der Simulation entfernt). Dieser soll den Effekt der gravitativen Fokussierung berücksichtigen.
- Für den Asteroidengürtel wird der Zähler indirekt implementiert, in dem ich eine Trefferwahrscheinlichkeit pro zurückgelegte Weglänge in einem bestimmten Bereich (gemessen in AU) vorgebe, und alle Teilchen, die diesen Bereich durchfliegen, gemäss dieser Wahrscheinlichkeit aus der Simulation entfernt werden, wobei jeweils der Zähler für den Asteroidengürtel um eins erhöht wird. Evtl. muss ein ähnlicher Prozess für die "Einbremsung" von Teilchen in den Asteroidengürtel entwickelt werden.
- Am Ende (Simulation für zumindest einige Millionen Jahre, wenn möglich) wollen wir wissen, wieviele Teilchen auf welchen Planeten (und im Asteroidengürtel) niedergegangen sind, und ob einige der Teilchen (insbesondere der massiven) in den Asteroidengürtel "eingebremst" wurden.
- Es sollte auch bestimmt werden, ob es eine Abhängigkeit zwischen dem anfänglichen Raumwinkel der "Abstrahlung" und der Trefferwahrscheinlichkeit bestimmter Planeten bzw. des Asteroidengürtels gibt.
- Das Programm sollte weiter in der Lage sein, Bilder bzw. Videos der Simulation auszugeben, und zwar (aus Anschaulichkeitsgründen) in der Ekliptikebene des Sonnensystems, sowie in den Parameterräumen Inklination vs. Exzentrizität vs. Grosse Halbachse.

Das schaffen wir, oder!? Ich auf jeden Fall fände es absolut grossartig. :)

Bernhard
15.08.2012, 11:46
Das schaffen wir, oder!? Ich auf jeden Fall fände es absolut grossartig. :)
Hallo Bynaus,

in welcher Sprache ist der Code von Nathan Kaib denn verfasst? Gibt es eine Dokumentation dazu?

Zweitens: Die neuen Aspekte von UMa und SRMeister sind prinzipieller Art. Inwieweit und in welcher Art Keplerbahnen in das Programm eingebaut werden können muss noch geklärt werden und das ist meiner Meinung nach keine triviale Angelegenheit.

Ob alle Punkte realisiert werden können, werden wir sehen. Ich bin da ein wenig skeptisch angesichts des Termins im Januar.

Ich mache also weiter, ABER es sollte hier bitte jeder berücksichtigen, dass momentan Ferienzeit ist und dass es auch an Wochenenden schon mal zu Verzögerungen kommen kann. Eigentlich Selbstverständlichkeiten, aber scheinbar muss man es doch immer wieder mal in Erinnerung rufen.
Gruß

Bynaus
15.08.2012, 13:03
in welcher Sprache ist der Code von Nathan Kaib denn verfasst? Gibt es eine Dokumentation dazu?

Im zip-File, das er mir geschickt hat, ist eine Dokumentation drin, aber ich habs mir bisher nicht angesehen. Soll ich es dir auch weiterleiten?


Ob alle Punkte realisiert werden können, werden wir sehen. Ich bin da ein wenig skeptisch angesichts des Termins im Januar.

Der Termin im Januar ist natürlich nicht sakrosankt. Es gibt immer wieder solche Termine bzw. Konferenzen. Ich wollte euch nur darüber informieren, ohne Anspruch, ihn zu erreichen. Wir sind fertig, wenn wir fertig sind...

Mir ist natürlich voll und ganz bewusst, dass das ein Freizeit-Projekt ist, das nur von eurem Enthusiasmus dafür lebt. Ich kann zur Zeit leider nicht besonders viel beitragen, aber dafür wird mein Beitrag später, wenn es darum geht, dies alles zu wissenschaftlichem Papier zu bringen, umso wichtiger werden.

SRMeister
15.08.2012, 13:20
Das Projekt ist in Fortran-77 geschrieben.
Eigentlich ist es für den kostenlosen ifort (Intel Fortran Compiler) ausgelegt, den es nur in Linux gibt. Aber ich habe es bereits in Visual Studio als Konsolenanwendung kompilieren können.
Da man außerdem eine Bibliothek erhält, bin ich zuversichtlich, dass ich SCATR mittelfristig an ein C++ oder C# Projekt anschließen kann.(zb an das GUI)

Bernhard
15.08.2012, 19:18
Soll ich es dir auch weiterleiten?
Ja. Bitte schicke mir das Paket. Falls Du meine eMail-Adresse nicht mehr hast, bitte eine kurze Mitteilung per PN.

Bernhard
15.08.2012, 19:26
Da man außerdem eine Bibliothek erhält, bin ich zuversichtlich, dass ich SCATR mittelfristig an ein C++ oder C# Projekt anschließen kann.(zb an das GUI)
Das klingt interessant. Unter Linux gibt es auch noch Übersetzungsprogramme von Fortran nach C. Wenn wir Anpassungen an SCATR machen wollen, ist so was eventuell ganz hilfreich.

@SRMeister: Ich denke es wäre gut, wenn Du vielleicht doch mal eine Schnittstelle zu Deinem GUI definieren würdest. Irgendwie kommen wir da schon zusammen. Sollte sich die Schnittstelle aufgrund erweiterter Anforderungen dann als zu schmal darstellen, müsstest Du natürlich nachbessern.

ralfkannenberg
15.08.2012, 22:40
Sollte sich die Schnittstelle aufgrund erweiterter Anforderungen dann als zu schmal darstellen, müsstest Du natürlich nachbessern.
Hallo Ihr beide,

sorry dass ich wieder aus der Praxis komme ...

Es ist ja zu erwarten, dass die Schnittstelle im ersten Entwurf unvollständig sein wird, allein schon deswegen weil sich aufgrund Eurer Erfahrungen auch die Anforderungen ändern werden.

Deswegen solltet Ihr von vornherein einplanen, dass die Schnittstelle erweitert werden wird. Mit dem Ansatz braucht man dann die Schnittstelle nur zu ergänzen und nicht "nachzubessern".


Freundliche Grüsse, Ralf

Bernhard
16.08.2012, 12:28
Das Projekt ist in Fortran-77 geschrieben.
Eigentlich ist es für den kostenlosen ifort (Intel Fortran Compiler) ausgelegt, den es nur in Linux gibt.
Ich kann SCATR mittlerweile unter OpenSuSE 12.1 übersetzen und kann den Code damit zumindest prinzipiell an unsere Bedürfnisse anpassen und die Bibliothek libscatr.a erzeugen. Zusätzlich sollte man natürlich auch nachsehen/herausfinden, welche Funktionen diese Bibliothek standardmäßig exportiert und wie man diese verwenden kann.

Bernhard
16.08.2012, 18:04
Aber ich habe es bereits in Visual Studio als Konsolenanwendung kompilieren können.
Hallo SRMeister,

kannst Du bitte beschreiben, wie Du das genau gemacht hast. Kann man dabei einfach die bestehenden Fortran-Dateien in eine bestehende C++-Konsolenanwendung einfügen und übersetzen? Ich würde hier zumindest ein paar spezielle Compiler-Einstellungen erwarten.

Ich habe mir unter Linux inzwischen das Beispiel mit der Datei mvs.in mit 5 Planeten angesehen. Scheint alles ziemlich gut zu funktionieren. Der nächste größere Schritt wäre dann, den Testkörpern wieder beizubringen mit der vorgegebenen Zufallsverteilung bei der Erde zu starten. Da ich Fortran bisher aber auch nur am Rande mitbekommen habe wird das schon eine Weile dauern.

SRMeister
16.08.2012, 23:06
Also ich habe es mit Intel Composer XE 2011 gemacht. Das fügt sich perfekt in Visual Studio 2010 (Professional oder höher) ein. Ist natürlich beides dann nichtmehr Freeware.
Dann hat man aber trotzdem 2 Projekte, also pro Programmiersprache ein Projekt.
Das SCATR kompiliere ich als Fortran-DLL und den C++ Teil als Konsolen-exe. Später soll der C++ Teil aber wieder als Schnittstelle zwischen der SCATR-DLL und dem C# Programm laufen. Somit hat man dann wieder 3 Programmiersprachen und 2 DLL-Bibliotheken. Aber zumindest für eine Konsolenanwendung benötigst du nur die Fortran-DLL und eine zugehörige von mir zu erstellende c++ include Datei. Ich werde meine Fortschritte bezüglich Kombination Fortran-DLL und C++ Konsolenanwendung verfügbar machen.

Ich möchte, dass der C++ Teil genau den Teil des Programms enthält, der momentan in scatr.f, also im Hauptprogramm zu finden ist. Damit mache ich momentan schon ganz gute Fortschritte. Ich muss dazu DLL-Interfaces für etwa 10 Funktionen schreiben, damit diese dann von C++ aufgerufen werden können. Leider hat Fortran einige nervige Eigenheiten und die SCATR Funktionen haben auch jeweils ca. 20-24 Parameter daher ist das alles erstmal nicht ganz so einfach. Aber ich denke, am Ende hat man von der C++ Seite genug Möglichkeiten, Einfluss zu nehmen, insbesondere auf die Datenquelle(Startbedingungen), sowie auf alle Parameter und Statistiken der Planets und TPs nach jedem Zeitschritt. Das Programm zählt übrigens genau wie von Bynaus gefordert, welcher TP wo drauffällt oder endet. Außerdem ist das Verfahren der Close Encounters dem von UMa beschriebenen sehr ähnlich, wenn nicht gleich.
Ich denke am Integrator selbst muss man keine Änderungen mehr vornehmen.

PS: Ich bin wirklich gespannt wieviel schneller dieses Programm sein wird.

Bernhard
17.08.2012, 00:32
Ich möchte, dass der C++ Teil genau den Teil des Programms enthält, der momentan in scatr.f, also im Hauptprogramm zu finden ist. Damit mache ich momentan schon ganz gute Fortschritte. Ich muss dazu DLL-Interfaces für etwa 10 Funktionen schreiben, damit diese dann von C++ aufgerufen werden können. Leider hat Fortran einige nervige Eigenheiten und die SCATR Funktionen haben auch jeweils ca. 20-24 Parameter daher ist das alles erstmal nicht ganz so einfach. Aber ich denke, am Ende hat man von der C++ Seite genug Möglichkeiten, Einfluss zu nehmen, insbesondere auf die Datenquelle(Startbedingungen), sowie auf alle Parameter und Statistiken der Planets und TPs nach jedem Zeitschritt. Das Programm zählt übrigens genau wie von Bynaus gefordert, welcher TP wo drauffällt oder endet. Außerdem ist das Verfahren der Close Encounters dem von UMa beschriebenen sehr ähnlich, wenn nicht gleich.
Ich denke am Integrator selbst muss man keine Änderungen mehr vornehmen.
Viel Erfolg bei der Umsetzung, denn damit wäre Bynaus' Anforderung weitgehend bis komplett umgesetzt. Ein Visual Studio Professional habe ich leider nicht zur Verfügung, weswegen ich mich dann erst mal zurücknehmen werde.

SRMeister
17.08.2012, 00:32
An Ralf,
danke für deine immer nützlichen Beiträge. So habe ich mir jetzt gedacht, bevor man garnicht weis, was man überhaupt für ein Programm will, ist es besser erst garkeine Anforderungen festzulegen. Alles was jetzt geschieht, ist zwanglos und ohne irgendeinen Anspruch. Alles was man daraus lernt, kann später bei einer konkret formulierten Umsetzung nützlich sein.

An Frank,
auch wenn ich dein Programm jetzt nicht als Grundlage nutzen möchte, würde ich es mir doch gerne einmal anschauen. Rein interessehalber. Ist das MFC? Ich hab früher auch mit MFC und OpenGL programmiert, ist aber schon ne Weile her. Heute würde ich eine optisch ansprechende Darstellung von Planeten/ Weltall nur mit einem Raytracer machen. Der ist für sowas perfekt(aufgrund mehrere Aspekte) und kann nahezu realistische Bilder mit wenig Aufwand erzeugen. Das wäre vielleicht mal etwas für die Zukunft :)

SRMeister
17.08.2012, 00:41
Viel Erfolg bei der Umsetzung, denn damit wäre Bynaus' Anforderung weitgehend bis komplett umgesetzt. Ein Visual Studio Professional habe ich leider nicht zur Verfügung, weswegen ich mich dann erst mal zurücknehmen werde.

Naja, aber die DLL mit include Datei kann man ja dann auch in Express benutzen. Da ist ja dann kein Fortrancompiler mehr nötig. Und falls SCATR später geupdatet werden sollte, bräuchte man wohl nur die DLL Datei auszutauschen und die darauf basierende Anwendung bleibt funktionsfähig.

Bernhard
17.08.2012, 01:00
Naja, aber die DLL mit include Datei kann man ja dann auch in Express benutzen.
Diese dll nur zu anzuwenden finde ich ehrlich gesagt wenig reizvoll. Da sehe ich mir lieber die Veröffentlichungen ausgehend von der Arbeit von Wisdom&Holman bis Nathan Kaib an.

Bynaus
17.08.2012, 11:00
@SRMeister: Das klingt ja alles sehr vielversprechend!

Ich hab mich mal auf die Suche nach Bottke's Simulationen der "Einbremsung" von einigen Planetesimalen in den Asteroidengürtel gemacht. Es geht um dieses Nature-Paper (http://www.nature.com/nature/journal/v439/n7078/full/nature04536.html), das soweit ich gesehen habe, leider nicht frei verfügbar ist (ausser man sitzt an einer Uni ;) ). Hier ist allerdings ein (frei verfügbares) Abstract der LPSC-Konferenz in Houston, das dem Paper vorausging: http://www.lpi.usra.edu/meetings/lpsc2006/pdf/1388.pdf Er geht hier offenbar davon aus, dass es im Asteroidengürtel damals viele Embryos (mit jeweils 0.04 M_E) gab, die bei der Einbremsung halfen. Ich bin nicht sicher, ob wir das zur Zeit der Mond-Bildung immer noch annehmen dürfen - wohl eher nicht mehr. Das macht es natürlich schwieriger. Ich glaube auch nicht, dass das Einbeziehen einiger grosser Asteroidengürtelobjekte (wie Ceres, Vesta, Pallas etc.) da gross helfen würde.

Was den Yarkovski-Effekt angeht, so hängt dieser ja stark von der Grösse des Fragments sowie seiner Beschaffenheit ab. Ausserdem haben kleine Fragmente eine sogenannte "collisional lifetime" (clt), eine typische Zeit, die vergeht, bis sie durch Kollisionen zerstört werden (genauer: nach der clt sind nur noch 1/e der ursprünglichen Fragmente vorhanden), die sich etwa zu clt = 16.8 * sqrt(R) Millionen Jahre berechnet, mit R in Metern. Bottke verweist übrigens auch auf die Arbeit von Vokrouhlicky und Farinella (2000) (http://adsabs.harvard.edu/abs/2000Natur.407..606V). Dieser Artikel ist hier frei verfügbar (http://astro.mff.cuni.cz/davok/papers/vf_nat00.pdf). Die beiden modellierten, wie der Yarkovski-Effekt sich auf die Fragmente von Asteroidenkollisionen auswirkt. Sie interessierten sich zwar für den Asteroidengürtel, aber ihre Schlüsse könnten auch für uns wichtig sein (z.B. wollen wir ja nicht den Anteil der Fragmente, die wieder auf die Erde fallen, überschätzen, nur weil wir den Yarkovski-Effekt vernachlässigen).

Das muss natürlich nicht alles sofort ins Programm. Das sind nur denkbare Erweiterungen und Einwände, die noch kommen könnten.

Was noch fehlt, ist die Absorbtionswahrscheinlichkeit im Gürtel. Ich geh dieser Frage als nächstes nach.

Bynaus
17.08.2012, 17:22
Okay, hier meine Hausaufgaben: ich lade alle gerne ein, das zu überprüfen und mich zu korrigieren, wenn ich hier völligen Mist rechnen sollte...

Für die Asteroiden mit 50 km Durchmesser und mehr (nur bei diesen kann man offenbar davon ausgehen, dass sie auch schon bei der Bildung des Sonnensystems zugegen waren - die kleineren sind Kollisionsfragmente) - es sind deren 700, die könnte man vielleicht sogar aktiv simulieren, wenn man ihre gegenseitigen gravitativen Beeinflussungen vernachlässigt... (? was meint ihr? Man könnte dann sogar die richtigen Bahnelemente verwenden) - habe ich eine kumulierte Querschnittsfläche von etwa 10+-5 Mio Quadratkilometern bestimmt. Das klingt nach viel, aber die Querschnittsfläche des Asteroidengürtels bei 3 AU (+-10° unter und über der Ekliptik) beträgt etwa 4.5e15 Quadratkilometer. Verdichtet man alle Masse auf diese Fläche bei 3 AU, käme eine Kollisionswahrscheinlichkeit bei Durchquerung dieser Fläche von ca. ~2e-11 raus. Nur schon daran sieht man, dass die Trefferwahrscheinlichkeit winzig ist.

Alternativ könnte man natürlich die 700 Objekte homogen über den Bereich 2 bis 4 AU (+-10° unter und über der Ekliptik) verteilen. Das gibt dann ein Volumen von 8e25 Kubikkilometern, oder 1e23 Kubikkilometern pro Asteroid. Wenn ich das jetzt mal mit der mittleren freien Weglänge L annähere, dann gäbe das:

L = 1/(n*W)

Mit n der Anzahl Asteroiden pro Volumeneinheit und W dem Wirkungsquerschnitt pro Asteroid. Wenn wir einen Kubikkilometer als Grundlage nehmen, haben wir also n = 1/1e23 und W = 1e7 / 700, womit die freie Weglänge 7e18 km wäre. Oder, mit anderen Worten, pro km beträgt die Absorptionswahrscheinlichkeit etwa 1.4e-19 oder etwa 2e-11 pro AU (das ist ziemlich genau der Wert, den wir oben schon erhalten hatten).

Das würde bedeuten, dass ein Fragment, das im Asteroidengürtel auf einer kreisrunden Bahn mit a=3 AU kreist, dementsprechend einen Bahnumfang von 19 AU und eine Periode von 5 Jahren hat, wobei es auf 100% dieser Bahn absorbiert werden könnte, ca. 13 Milliarden Jahre lang kreisen müsste, bis es absorbiert wird (im Mittel). Okay... Vielleicht kann man das wirklich vernachlässigen.

Bernhard
17.08.2012, 19:05
Für die Asteroiden mit 50 km Durchmesser und mehr (nur bei diesen kann man offenbar davon ausgehen, dass sie auch schon bei der Bildung des Sonnensystems zugegen waren - die kleineren sind Kollisionsfragmente) - es sind deren 700, die könnte man vielleicht sogar aktiv simulieren, wenn man ihre gegenseitigen gravitativen Beeinflussungen vernachlässigt... (? was meint ihr? Man könnte dann sogar die richtigen Bahnelemente verwenden)
Wenn man auch noch die gravitative Wirkung der Asteroiden auf die Planeten vernachlässigt, könnte man die Asteroiden einfach zu den Testkörpern zählen.


- habe ich eine kumulierte Querschnittsfläche von etwa 10+-5 Mio Quadratkilometern bestimmt. Das klingt nach viel, aber die Querschnittsfläche des Asteroidengürtels bei 3 AU (+-10° unter und über der Ekliptik) beträgt etwa 4.5e15 Quadratkilometer.
Um die Zahlen korrekt zu haben: Das sind 4.4e17 Quadratkilometer.


Verdichtet man alle Masse auf diese Fläche bei 3 AU, käme eine Kollisionswahrscheinlichkeit bei Durchquerung dieser Fläche von ca. ~2e-11 raus.
Die Wahrscheinlichkeit passt dann wieder :) .


Vielleicht kann man das wirklich vernachlässigen.
Von mir aus momentan natürlich ein ja dazu.

BTW: Ich habe mich heute ein wenig in den SCATR-Code eingelesen, denn ich würde gerne Vergleiche mit unserem BS-Code rechnen. Man muss aber Fortran schon recht gut lesen können, um die Input-Dateien so modifizieren zu können, dass "unsere" Startbedingungen (=01.04.2011 oder so ähnlich) für alle Planeten gerechnet werden. Die Anzahl der Testkörper kann man dabei erst mal auf Null setzen. Mich würde interessieren, welche Ergebnisse SCATR beispielsweise für ein Jahr berechnet und wie groß die Abweichungen zu den Ergebnissen mit dem BS-Code sind.

Bernhard
18.08.2012, 13:11
@Bernhard,
@SRMeister,
das, was ihr hier gerade abzieht, ist Kinderkram!
:mad:
Hallo Frank,

da hast Du Recht und ich danke Dir genaugenommen für diesen "Tritt in den A...", zeigt es doch, dass an dem Thema ein gewisses Interesse besteht. Das laufende Projekt wurde durch den SCATR-Code für mich wieder deutlich interessanter. Wenn Du den Fortran-Code ebenfalls haben willst, sende ich ihn Dir gerne zu. Der ist nämlich wirklich gut gemacht und sehr lehrreich.
Schönen Gruß

Bernhard
20.08.2012, 08:53
Weiß jemand, wie eine solare Masse von 2.959e-4 zu verstehen ist? Aufgrund der anderen Angaben sollte dieser Wert proportional zur SI-Masse sein. Geometrische Einheiten sind es nicht.

SRMeister
20.08.2012, 18:05
ich bin erstmal im Urlaub.
Bernhard, vielleicht in Sonnenmassen?

Bernhard
20.08.2012, 19:00
Weiß jemand, wie eine solare Masse von 2.959e-4 zu verstehen ist?
Das ist missverständlich ausgedrückt. Die Masse der Sonne wird im Programm mit dem Wert 2.959e-4 verwendet.

Bernhard
20.08.2012, 19:10
Der nächste größere Schritt wäre dann, den Testkörpern wieder beizubringen mit der vorgegebenen Zufallsverteilung bei der Erde zu starten.
Einige Zwischenergebnisse:

a) Die Startbedingungen der Testkörper werden über die Datei tp.in gesteuert. Diese Datei enthält allerdings etliche Parameter, die mir noch nicht ganz klar sind.

b) Die Datei param.in steuert den Programmablauf und gibt aktuell vermutlich eine Simulation über 1 Millionen Jahre vor. Ich habe dazu einige Beispiele mit Laufzeiten über 1 Jahr und 100 Jahre gerechnet.

c) Bei den Startbedingungen der Planeten werden in der Datei pl.in die heliozentrischen Koordinaten in AU und die Geschwindigkeiten in AU/Tag angegeben.

Bynaus
20.08.2012, 19:28
Die Masse der Sonne wird im Programm mit dem Wert 2.959e-4 verwendet.

Darauf kann ich mir keinen Reim machen.

Im Allgemeinen würde ich vorschlagen, konsequent SI-Einheiten zu benutzen, selbst wenn sie sperrig gross werden. Denn da macht man bestimmt keine Umrechnungsfehler.

Bernhard
20.08.2012, 20:24
Darauf kann ich mir keinen Reim machen.
Schade.


Im Allgemeinen würde ich vorschlagen, konsequent SI-Einheiten zu benutzen, selbst wenn sie sperrig gross werden.
In manchen Fällen muss man Zahlen in einem bestimmten Bereich halten, um Rechen- und Rundungsfehler zu vermeiden. Bevor wir oder ich den Code nicht besser verstehe(n), kann man das also erst mal so lassen. Die Einheit für die Masse sollte allerdings geklärt werden. Man könnte dazu, um Nathan Kaib nicht unnötig zu belästigen, eventuell auch einen SWIFT-Autor kontaktieren.

BTW: Ich denke momentan zusätzlich auch über Möglichkeiten nach den gesamten Code in C/C++ zu übersetzen. Das Linux-Tool f2c scheitert zwar an einigen Passagen, könnte dabei aber Tipparbeit abnehmen. Hilfreich wäre so eine Übersetzung auf jeden Fall.
Gruß

SRMeister
21.08.2012, 10:14
Vielleicht hilft das?

The code requires units in which the grav. constant G is unity.
Any combo of lengths, masses and times that keeps that true is O.K.
For example, the pl.in file given in the examples dir. uses lengths
in AU and time in days (where one year is exactly 365.25 days). This
forces the Solar mass to be about 2.96e-4

ralfkannenberg
21.08.2012, 10:20
In manchen Fällen muss man Zahlen in einem bestimmten Bereich halten, um Rechen- und Rundungsfehler zu vermeiden.
Hallo Bernhard,

selbstverständlich, aber dieses Problem kann man auch lösen, indem man die SI-Einheiten "normalisiert", d.h. in Mantisse und Zehnerpotenz aufspaltet; die Mantissen sind im "bestimmten Bereich" und die Zehnerpotenzen kann man seperat behandeln, konkret durch triviale Addition und Subtraktion ihrer Exponenten. So kann man Rundungsfehler ebenfalls ganz einfach vermeiden und arbeitet trotzdem mit SI-Einheiten.


Freundliche Grüsse, Ralf

Bernhard
21.08.2012, 10:59
Vielleicht hilft das?
Das erklärt es dann.

Bynaus
24.08.2012, 11:43
Hier ist übrigens ein Abstract von der LPSC von letztem Jahr - eingereicht von einem Masterstudenten. Es ist eine wirklich nette, einfache Arbeit, und nicht allzu weit weg von dem, was wir hier versuchen. Sollte unser Projekt gelingen, würde mir ein solcher Abstract als erster Schritt zu einer Publikation vorschweben (das ist eben diese Konferenz, deren Deadline im frühen Januar ist).

http://www.lpi.usra.edu/meetings/lpsc2012/pdf/2467.pdf

Bernhard
24.08.2012, 14:07
Hallo Bynaus,

gerade für große Zeiträume bietet SCATR in der Tat interessante Ansätze. Einem BS-Integrator würde ich als Referee im Bereich von Millionen Jahren selbst nicht trauen.

Ich arbeite mich also weiter in den Code ein und habe im Prinzip jetzt auch alles in C. Allerdings kommt der automatische Übersetzer mit den Ein- und Ausgabefunktionen nur sehr schlecht zurecht. Diese Teile muss man selber übersetzen. Dabei lernt man aber auch etwas über die Bedienung und die verschiedenen Optionen des Programmes. Ich melde mich spätestens, wenn ich eine lauffähige Version in C habe. Das interessiert mich unabhängig von SRMeisters Ansätzen für ein mehrsprachiges Projekt.
Gruß

Bernhard
24.08.2012, 16:25
INFO: Ich habe SCATR testweise eben mit den angepassten Startbedingungen aus dem alten BS-Programm (Sonne + Planeten + 4 Asteroiden) und Schrittweite ein Tag 365 Tage vom 01.04.2011 aus rechnen lassen. Die Position von Merkur weicht von den Ergebnissen vom NASA-Horizons-Interface um ca. 0.01 AU ab. Beim Mond bekommt eine ähnliche Abweichung von 0.017 AU. Übermäßig präzise ist das nicht, aber für den Anfang auch nicht ganz verkehrt.

EDIT: Hier noch der Link auf Referenz [5] des obigen pdfs: http://adsabs.harvard.edu/abs/1965ApJ...141.1536A

SRMeister
25.08.2012, 10:29
Hallo Bernhard,
SCATR verlangt heliozentrische Koordinaten, unsere waren glaube ich baryzentrisch. Aber aus Horizons bekommt man auch heliozentrische

Bernhard
25.08.2012, 15:07
Hallo Bernhard,
SCATR verlangt heliozentrische Koordinaten, unsere waren glaube ich baryzentrisch.
Da bei "bar" (also bei unseren Startbedingungen) die Position und Geschwindigkeit der Sonne bekannt ist, ist die Umrechnung von bar -> hel lediglich je eine Vektoraddition für Position und Geschwindigkeit und genau so habe ich es berücksichtigt. Das ergibt z.B. die folgende Datei für pl.in:


15 15
2.959122082855912e-004
0.000000000000000e+000 0.000000000000000e+000 0.000000000000000e+000
0.000000000000000e+000 0.000000000000000e+000 0.000000000000000e+000
4.912547342993006e-011
-3.523313258706317e-001 9.504611091594724e-002 4.009426609758516e-002
-1.315273650151517e-002 -2.595761761117310e-002 -9.140232477760470e-004
7.243452490032174e-010
1.643585342906148e-001 -7.085414292529864e-001 -1.919255484837590e-002
1.956741558333183e-002 4.501098380608712e-003 -1.067629632867861e-003
8.887692392411485e-010
-9.814911942588033e-001 -1.866158863169087e-001 5.097029871362697e-006
2.931839043977104e-003 -1.697334005001157e-002 3.022067795326298e-007
1.093189659159437e-011
-9.789137899355652e-001 -1.874333963075148e-001 2.346822772010697e-004
3.106832494474352e-003 -1.644105091782470e-002 1.293453737002694e-005
9.549535048237397e-011
1.365015143164586e+000 -2.312407313953970e-001 -3.836388232646479e-002
2.872312236658535e-003 1.499323465249643e-002 2.436402223642883e-004
1.404726574197985e-013
2.310378879836105e+000 -1.818355029949789e+000 -4.824841588729188e-001
5.907611012671520e-003 7.483884537718916e-003 -8.551195036579472e-004
4.196673818013239e-015
-2.815351112852315e+000 2.412702045753047e-001 5.892739737641274e-002
-3.216956785527108e-003 -9.132876949748969e-003 2.202761475540308e-003
3.188506350957331e-014
8.656984163765726e-001 -2.699411709384842e+000 1.793408702321638e+000
8.215311727150669e-003 7.979516547040460e-004 -1.235707923749165e-003
3.968910003289545e-014
1.214545584869641e-001 -2.158557655801899e+000 5.016793101419336e-002
1.208600933373512e-002 3.155876262806481e-004 -1.477945481111473e-003
2.825345909640263e-007
4.760551581845321e+000 1.345995375310020e+000 -1.121213280759918e-001
-2.147764976969137e-003 7.627105635028521e-003 1.637437370629982e-005
8.459715186095996e-008
-9.333257286778627e+000 -2.260873964001870e+000 4.108752977421108e-001
1.011145239145509e-003 -5.428930560404771e-003 5.459982590329613e-005
1.292024916765944e-008
2.008426347372372e+001 1.746065761892764e-001 -2.595157292190287e-001
-6.456265483588176e-005 3.755536279097650e-003 1.476893678293213e-005
1.524358900811198e-008
2.558344591678675e+001 -1.568707229328016e+001 -2.664675853628982e-001
1.618462759020618e-003 2.700116634445506e-003 -9.306132082543796e-005
2.188700015808439e-012
3.078644668841170e+000 -3.175650493699054e+001 2.508748490001770e+000
3.181039267103034e-003 5.954186745479077e-006 -8.872059485933566e-004


Der Umrechnungsfaktor von unserem GM auf den Masseparameter ist 86400^2/(AU^3). Den Wert für eine AU habe ich aus Horizons genommen.

SRMeister
25.08.2012, 17:40
ja gut die Werte sehen soweit ganz gut aus, ich finde Differenzen zu Horizons im Bereich ab der 9. Stelle hinterm Komma, aber ich habe nur ein paar Werte geprüft.
Was mir noch einfällt, man kann irgendwo die Genauigkeit einstellen bei SCATR. Entweder in den input-files, oder bei den Preprozessor-Direktiven im Quellcode.
Man soll die Schrittweite so wählen, dass sie etwa 1/100 der Periode des Sonnennächsten Planeten ist, also sollte es daran wohl nicht liegen.
Wenn ich ausm Urlaub zurück bin, kann ich auch wieder etwas intensiver mitmachen.

Bernhard
30.08.2012, 11:03
Allerdings kommt der automatische Übersetzer mit den Ein- und Ausgabefunktionen nur sehr schlecht zurecht.
Leider sind die automatisch übersetzten Dateien bis auf die Funktionsschnittstellen kaum zu gebrauchen. Die Übersetzung wird also noch etwas mehr Zeit in Anspruch nehmen. Momentan habe ich in etwa die Hälfte der Dateien geschafft. Die Überprüfung des Programmes wird dann auch noch mal etwas dauern, weil die Standardindizes bei Fortran bei 1 und bei C bei 0 starten. Zum Glück findet man im Netz relativ viele Quellen zu Fortran77, so dass man alles ohne Bücher machen kann.

Bernhard
03.09.2012, 13:23
Hallo Bynaus, SRMeister,

die Übersetzung des SCATR-Codes ist fertig, aber noch nicht fehlerfrei. Ich habe derzeit ein Visual Studio 6.0 C++-Projekt, weil die Express-Version von 2010 mit der Ordnerstruktur irgendwie nicht klar kommt. In dem 6.0er-Projekt EDIT: befindet sich aktuell noch ein syntaktischer Fehler.
Gruß

Bynaus
03.09.2012, 14:42
die Übersetzung des SCATR-Codes ist fertig, aber noch nicht fehlerfrei.

Wow, super, danke. Das war bestimmt eine grosse Arbeit (wie wohl auch das debugging, das noch anzustehen scheint). Lasst mich einfach wissen, wenn ihr einen Input von mir braucht. Zurzeit kann ich ja, wie es scheint, nicht viel beitragen.

Bernhard
03.09.2012, 17:08
@Bynaus,


Das war bestimmt eine grosse Arbeit
ich schätze mal so 40-50 Stunden. Es hat mich aber einfach interessiert Fortran mal ein wenig besser kennen zu lernen :) . Bei der Verwaltung von Matrizen funktioniert Fortran scheinbar etwas anders, als C, was an zwei, drei Stellen noch einige unschöne Zusatzvariablen nötig gemacht hat, die man langfristig eventuell durch eine geschicktere Programmierung wieder herausnehmen kann.


wie wohl auch das debugging, das noch anzustehen scheint.
Das Programm sollte in C natürlich die gleichen Ergebnisse berechnen, wie mit Fortran. Das kann nochmal einige Arbeit bedeuten, um kleine, versteckte Tippfehler zu finden. Aber prinzipiell sollte ich das auch alleine hinbekommen.

Den fertigen Code würde ich dann vorerst vor allem an SRMeister weitergeben wollen, damit er eine Visualisierung (Youtube-Clip oder etwas in der Art) damit programmieren kann. Sollte daran anschließend wirklich ein Paper entstehen, würde ich den Code voraussichtlich auch an die Originalautoren weitergeben.
Gruß

Bernhard
03.09.2012, 19:18
Das Programm sollte in C natürlich die gleichen Ergebnisse berechnen, wie mit Fortran. Das kann nochmal einige Arbeit bedeuten, um kleine, versteckte Tippfehler zu finden. Aber prinzipiell sollte ich das auch alleine hinbekommen.
Visual Studio 6.0 ist wegen der Größe des Projektes ausgestiegen, aber ich habe jetzt eine Visual Studio Express Version die man auch debuggen kann. Die maximale Anzahl der Testkörper musste ich dazu erst von 30001 auf 3001 und dann noch auf 301 reduzieren. Das Limit im Debug-Modus liegt also irgendwo zwischen 301 und 3001. Der Release-Modus läßt hoffentlich noch mehr Testpartikel zu.
Gruß

SRMeister
03.09.2012, 23:36
das klingt doch schonmal ganz gut!
ich würde mir den Code gerne anschauen, vllt. kannst du mir den per EMail schicken?

Grüße

Bernhard
04.09.2012, 00:48
ich würde mir den Code gerne anschauen, vllt. kannst du mir den per EMail schicken?
Hallo SRMeister,

eigentlich würde ich den Code gerne noch so weit debuggen, dass er die gleichen Ergebnisse liefert, wie das Original. Wenn Dir das zu lange dauert, kannst Du mir per PN aber auch Deine eMail-Adresse zuschicken.

Ich habe beispielsweise gerade bemerkt, dass man CR+LF bei den drei Beispieldateien *.in möglichst in die DOS-Konvention umwandeln sollte, weil sonst die ganzen io-Funktionen praktisch nur noch "Schrott" liefern.

Für unser Projekt kann man später die drei Eingabedateien auch durch Code ersetzen. Aktuell möchte ich aber zuerst eine möglichst genaue Übersetzung machen, damit man diese Version dann bei Bedarf an die Autoren des Fortran-Codes schicken kann, sozusagen als Dank für das Original.
Gruß

Bernhard
04.09.2012, 18:19
Hallo SRMeister,

das Programm läuft im Debug-Modus jetzt durch und erzeugt dabei zumindest korrekt formatierte Ergebnisse. Die berechneten Positionen sind nicht komplett verkehrt aber leider nicht identisch mit den Originalzahlen.

Da ich die nächsten Tage gerne ein wenig Abstand von dem Programm haben will, sende ich Dir den Code jetzt einfach mal zu. Im Release-Modus liefert das Programm recht eigenartige Ausgaben in der DOS-Box, die im Debug-Modus nicht erscheinen :confused: .

Du kannst den Code jetzt also querlesen und mit dem Programm einfach mal ein bisschen herumspielen.
Gruß

EDIT: Wenn wir uns ausführlich austauschen und über jede Änderung am Code unterrichten, können wir auch gemeinsam den Code debuggen und dann zu neuen Versionen zusammen mergen. Du besitzt für das Debugging die idealen Entwicklungswerkzeuge und es wäre schön, wenn Du diese Möglichkeiten nutzen könntest. Ich werde in der Zwischenzeit mal nachsehen, ob meine ifort-Installation unter Linux eine grafische Oberfläche mit Debugger besitzt. Vor einem zeilenweisen Debugging mit write-Ausgaben graust es mir schon ein wenig.

Bernhard
05.09.2012, 18:22
Im Release-Modus liefert das Programm recht eigenartige Ausgaben in der DOS-Box, die im Debug-Modus nicht erscheinen :confused: .

Dieser Punkt hat sich gerade geklärt. Ich hatte einfach vergessen die drei Eingabedateien in das Arbeitsverzeichnis zu kopieren ("selbst hau"). Hier sollte/kann man ein paar Meldungen einbauen, um so triviale Fehler gleich abzufangen.

Kibo
06.09.2012, 09:30
Aufgrund eures Projektes wurde bei mir das Interesse an solchen Simulationen geweckt. Mit meinen bescheidenen Kenntnissen ist das hier (https://skydrive.live.com/?cid=58397D4453EFC662&id=58397D4453EFC662%21132) dabei heraus gekommen.

mfg

Bernhard
06.09.2012, 09:55
Hallo Kibo,

schön, dass Du Dich für dieses Thema hier interessierst. Dein Pong scheint ganz gut zu funktionieren :) . Das Programm planeten oop1.exe muss wohl noch etwas weiter entwickelt werden :confused: Wie wäre es, wenn Du einfach "unseren" BS-Astronews-Integrator (http://sourceforge.net/projects/originofthemoon/?source=directory) einbauen würdest? Da laufen dann keine Planeten mehr aus dem Display und man kann auch zusätzlich Testkörper simulieren. Dieser Integrator wurde von SRMeister bereits erfolgreich verwendet.

SRMeister
06.09.2012, 12:11
Hallo Kibo, auch ich finde das sehr interressant. Auto Hotkey war mir auch völlig unbekannt. Scheint ja sehr leicht damit, Windows Progs zu schreiben.

Irgendwie ist aber dein Integrator noch etwas fehlerhaft: Denn wenn ich die Maus so lege, dass "a" einen auf dem Bildschirm sichtbaren Orbit um "b" einnimmt, dann fällt auf, dass die Bahn keine Keplerbahn ist. Die Perihel- und Aphel Punkte sind nicht bei jedem Umlauf an der selben Stelle sondern schwanken sehr stark. Ich vermute dass das Abstands-Kraftgesetz nicht richtig forumliert ist.
Ansonsten aber sehr schön! Die Idee mit dem Mauszeiger gefällt mir gut, so bekommt man irgendwie ein Gefühl für die wirkenden Kräfte.

Grüße

Kibo
06.09.2012, 12:53
Danke für die Antworten,

Die .ahk Dateien lassen mit jedem Texteditor öffnen. Ich hab versucht objektorientierte Elemente mit einfließen zu lassen, einfach um das mal auszuprobieren. Normal programmiere ich lieber strukturiert.

Gut möglich, dass das physikalisch nicht ganz richtig umgesetzt ist, das Ganze ist halt noch in Arbeit. Mir war es erstmal wichtig das sich das Ganze halbwegs richtig bewegt :D.
Fehlen tut da noch eine ganze Menge z.B. richtige Masse, mehr als 2 Körper. Die ganze Impulsberechnung muss dafür noch umgestellt werden. Langfristig möchte ich noch in die dritte Dimension und dafür OpenGL in Autohotkey nutzen (das ist durchaus machbar). Da ist naoch ganz viel Platz nach oben offen.

Den Integrator schau ich mir auf jeden Fall noch genauer an, da lässt sich sicher viel daraus lernen. Ansonsten geht es mir hierbei nur ums Verständnis und da hilft selber Programmieren am meisten:)

edit: Das Pong ist auch noch alles andere als Fertig. Mir ging es bei dem Projekt um die Echtzeit-Datenübertragung zwischen Client und Server, welche wie angedeutet noch nicht so ganz funktioniert. Je nach Lust und Laune wird die Grafik auch noch etwas aufgepeppt

Bernhard
06.09.2012, 21:59
Hallo SRMeister,

ich habe meine Zwischenpause wieder beendet und mache jetzt mit dem Debugging des Cpp-Codes weiter. Ich bin gerade einem Fehler in der Berechnung der Testkörper auf der Spur.

SRMeister
07.09.2012, 00:12
Das klingt echt super.
Ich hab leider momentan grad etwas beruflichen Stress, aber bin bald wieder mehr dabei, solange bleibe ich mal passiv.
Ein Problem auf das ich gestoßen bin als ich ein Interface bauen wollte, waren die Strings. In Fortran werden die immer mit verstecker Länge übergeben, dass heist die Funktionen haben in C entsprechend mehr Parameter als in Fortran angegeben, pro char* noch ein int mit der Länge. In C sind die ja mit \0 terminiert, also ist die Längenangabe überflüssig solande der String nich geändert wird. Bist du da auf Probleme gestoßen?

Die ganze Sache mit den globalen Arrays für alle Zustände widerspricht zwar einerseits meiner objektorientierten Einstellung, andererseits finde ich das in dem Fall garnicht schlecht, so können wir schon recht einfach Einfluss auf die Parameter nehmen und unsere eigenen Importfunktionen schreiben.

Bernhard
07.09.2012, 09:15
Bist du da auf Probleme gestoßen?
Da haben wir nichts zu befürchten, weil es hier nur um drei Dateinamen geht, die ich für unsere Zwecke sowieso schon fest in den Code eingebaut habe. Bei der Version, die ich Dir geschickt habe habe ich als Dateinamen param.in, pl.in und tp.in verwendet. Die Eingabe der Startparameter über Dateien kann und sollte für unsere Zwecke eh komplett ausgebaut werden. Dieses Feature bringt für uns nichts, weil wir mit festen Startbedinungen recht gut auskommen und zusätzliche Änderungen per eMail, bzw. Internet vereinbaren können. Bei den restlichen Stringpararmetern kann man ohne wesentliche Performanceverluste auch in den Funktionen nach dem Stringendezeichen 0 suchen.


Die ganze Sache mit den globalen Arrays für alle Zustände widerspricht zwar einerseits meiner objektorientierten Einstellung, andererseits finde ich das in dem Fall garnicht schlecht, so können wir schon recht einfach Einfluss auf die Parameter nehmen und unsere eigenen Importfunktionen schreiben.
Es gibt im Fortran-Code sowieso nur ein globales Array (r2hill in der Datei scatr.f) und das wird nirgends verwendet. Im Cpp-Code habe ich diese Variable deswegen gar nicht mehr angelegt.

Für das Debugging brauche ich noch einen grafischen Debugger, weil sich das Programm im Detail doch sehr stark verästelt. Ich habe mir gestern mal FTN95 angesehen. Das ist in der Personal-Edition aber nicht wirklich ernst zu nehmen. Heute möchte ich mir den Intel Composer ansehen. Wenn es da eine 30-Tage-Testversion gibt wäre das schon mal sehr hilfreich. Die Linux-Version habe ich ja schon, aber ich habe da aktuell nur den Kommandozeilen-Compiler ifort gefunden und der reicht für unsere Zwecke einfach nicht aus.

EDIT: Nachdem der Intel Composer nur mit Visual Studio Professional läuft, wird sich das Debugging von meiner Seite aus wohl wieder etwas verzögern.

Bernhard
10.09.2012, 19:28
EDIT:

Ich habe jetzt Visual Studio und den Intel Compiler als Testversion installiert und habe das Fortran-Projekt zum Laufen gebracht. Prinzipiell habe ich damit jetzt rund 27 Tage zum debuggen.
Gruß

Bernhard
15.09.2012, 22:01
Hallo Bynaus und SRMeister,

ich bin mit dem Cpp-Code wieder ein Stückchen weiter gekommen. Bei einer Simulationszeit von einer Millionen Jahre gab es bisher immer einen "Planeten", der dabei instabil wurde. Im Debugger habe ich jetzt gesehen, dass dieser "Planet" der Mond ist. Entfernt man den Mond aus den Startbedingungen läuft das Programm ohne Fehler durch.

Ich werde jetzt das Programm also mit Testkörpern laufen lassen. Man kann in der jetzigen Version bereits über eine Variable festlegen, wieviel Trümmerteile von der Erdoberfläche aus starten sollen, allerdings bekommt man damit wieder Fehlermeldungen. Das Debugging kann/muss also weiter fortgesetzt werden.
Gruß

EDIT: Der Cpp-Code läuft jetzt mit 20 Trümmerteilen schon mal für T_Sim = 1 Mio. Jahre durch, allerdings gefallen mir die Ergebnisse für die Trümmerteile noch nicht. Bei dem Fortran-Code verlassen in diesem Fall bereits etliche Trümmerteile das Sonnensystem und das sollte das übersetzte Programm meiner Meinung nach ebenfalls zeigen.

Kibo
18.09.2012, 08:34
Guten Morgen,

mir ist gerade eine Idee gekommen was man noch berechnen könnte.

Szenario:

In 10 Jahren sind alle für Sateliten lohnenswerten Erdorbits mit Milliarden von Trümmerteilchen zugekleistert. Da die Teilchenmenge zu hoch ist um sie per Laser zu entfernen, beschließen die Regierungen auf der Erde seinen Orbits eine Radikalkur zu verpassen :)
Ein Asteroid mit ausreichender Größe muss her und an der Erde haarscharf vorbei ziehen. Seine Eigengravitation soll die Bahnen der Teilchen destabilisieren.
Ein Riskanter Plan, aber lohnenswert?

----------------------------------------------------------------------------------------

Mein Programm ist übrigends auch ein Stück voran geschritten (link hat sich nicht verändert). Ich habe gerade etwas mehr Zeit, vielleicht kommt heute noch das große Update bei dem dann alle Objekte sich gegenseitig gravitativ beeinflussen, atm kreisen alle um eines.

mfg

Bynaus
18.09.2012, 08:55
@Bernhard und SRMeister: Ich freue mich, dass ihr Fortschritte macht.
Ich hab mir die Sache mit den Simulationen von massereichen Partikeln nochmals durch den Kopf gehen lassen und auch die eigentlichen Kollisionssimulationen nochmals angeschaut. Ich würde sagen, dass wir zwei Dinge simulieren könnten:
1) Eine Wolke masseloser Partikel (so viel wie die Rechenzeit eben zulässt), um zu sehen, wohin es Trümmer allgemein verschlagen kann
2) Ein paar wenige, massive Trümmer, die sich im Bereich von einigen Mondmassen befinden. Solche entstehen bei den Kollisionen nämlich recht häufig.

Wäre interessant zu sehen, wo die Trümmer aus 2) am Ende hinkommen, ob sie z.B. alle mit der Erde kollidieren (dann müsste man das bei diesem Szenario künftig sicher berücksichtigen), oder ob sie das System verlassen, mit anderen Planeten zusammenstossen... etc.

@Kibo: Ich bezweifle, dass der Vorbeizug eines Asteroiden auch nur die geringste Auswirkung auf solche Teilchen hätte (ausser jene, die sich in extremer Nähe befinden). Die Hill-Sphäre eines 10 km-Asteroiden (Dichte 3 g/ccm) auf einem parabolischen Orbit in 36000 km Entfernung von der Erdoberfläche beträgt gerade mal r = a * Wurzel (m_asteroid/(M_erde * 3)) oder etwa 400 m (noch weniger wenn er noch näher an der Erde vorbeizieht), kleiner als sein Radius also. Das heisst, dass an jedem Punkt ausserhalb der Asteroidenoberfläche die Gravitation der Erde über jene des Asteroiden dominiert.

Kibo
18.09.2012, 09:30
Ich gehe davon aus, das eine Dominanz auch nicht nötig sein muss. Ich hoffe eher, dass durch den Vorbeizug, die Orbits der Teilchen elliptischer werden, und dann in die Erdatmospghäre stürzen. Ich werd jedenfalls versuchen das zu simulieren

Bernhard
18.09.2012, 10:16
STATUS WinSCATR:

1.) Die Ergebnisse des Cpp-Codes nähern sich immer weiter an die Ergebnisse des Fortran77-Codes an.
2.) Die Rechengenauigkeit des Cpp-Codes ist bei Standardfunktionen, wie der Quadratwurzel, bei Cpp deutlich höher, als bei Fortran77.
3.) Einige offensichtliche Fehler beim Datenlogging und den Ausgabefunktionen sind noch im Cpp-Code enthalten. Ich werde natürlich versuchen diese zu beheben, solange meine Testversionen noch benutzbar sind.
4.) Die Funktion der globalen Variable r2hill hatte ich übersehen und ist mittlerweile eingebaut. Damit können Testkörper von Planeten eingefangen werden.
5.) In den ersten hunderttausend Jahren der Simulation werden bei dem Cpp-Code noch deutlich mehr (15 von 20) Testkörper eingefangen, als bei dem Fortran-Code (2 von 20). Auch hier können sich neben der unterschiedlichen Rechengenauigkeit noch kleinere Programmierfehler befinden, die aber über die Ausgabefunktionen gefunden werden können.
6.) Die Schnittstelle der Step-Funktion rmvs3_step (Schiebt das Gesamtsystem in der Zeit um dt nach vorne) musste ich um einige Cpp-spezifische Parameter erweitern, damit die Laufzeit nicht unnötig verlängert wird. Die Ausgabefunktionen und damit die Schnittstellen zur Visualisierung sind davon aber nicht betroffen. Die Laufzeit ist damit bei beiden Programmversionen in etwa gleich groß.
7.) Die Rechenergebnisse sämtlicher bisher verschickten Versionen (per PN) sind nicht wirklich ernst zu nehmen :) .
Gruß

Kibo
18.09.2012, 10:51
Zitat von Microsoft


Sie können Ihre Testphase auf bis zu 90 Tage verlängern, wenn Sie sich registrieren und einen kostenlosen Product Key für die Testversion erhalten.

So wie ich das sehe, auf Hilfe und dann Testversion registrieren gehen, dann kriegt man insgesamt 90 Tage :)

Bernhard
18.09.2012, 11:25
So wie ich das sehe, auf Hilfe und dann Testversion registrieren gehen, dann kriegt man insgesamt 90 Tage :)
Danke für den Tip. Für Visual Studio habe ich mittlerweile auch eine ziemlich günstige Lizenz bei eBay bestellt. Bliebe noch der Intel Compiler.

Bynaus
18.09.2012, 11:36
@Bernhard:


1.) Die Ergebnisse des Cpp-Codes nähern sich immer weiter an die Ergebnisse des Fortran77-Codes an.
2.) Die Rechengenauigkeit des Cpp-Codes ist bei Standardfunktionen, wie der Quadratwurzel, bei Cpp deutlich höher, als bei Fortran77.

Schön! Beim Vergleichen der Ergebnisse sollten wir die Rechengenauigkeit nicht vergessen: nicht, dass wir versuchen, die Ergebnisse des Fortran77-Codes zu replizieren, aber an dessen geringerer Rechengenauigkeit scheitern!


Die Funktion der globalen Variable r2hill hatte ich übersehen und ist mittlerweile eingebaut. Damit können Testkörper von Planeten eingefangen werden.
5.) In den ersten hunderttausend Jahren der Simulation werden bei dem Cpp-Code noch deutlich mehr (15 von 20) Testkörper eingefangen, als bei dem Fortran-Code (2 von 20). Auch hier können sich neben der unterschiedlichen Rechengenauigkeit noch kleinere Programmierfehler befinden, die aber über die Ausgabefunktionen gefunden werden können.

Was heisst, "eingefangen"? Der permanente Einfang eines masselosen Teilchens sollte eigentlich (aus Symetriegründen) nicht möglich sein. Oder meinst du eine Kollision? Implementiert der Fortran-Code eigentlich auch Kollisionen?

@Kibo: du kannst das gerne simulieren, aber ich sage voraus, dass der Vorbeiflug keinen nennenswerten Einfluss auf die Teilchenbahnen haben wird... :)

Bernhard
18.09.2012, 11:41
Was heisst, "eingefangen"? Der permanente Einfang eines masselosen Teilchens sollte eigentlich (aus Symetriegründen) nicht möglich sein. Oder meinst du eine Kollision? Implementiert der Fortran-Code eigentlich auch Kollisionen?
Das sind gute Fragen und man wird sehen müssen, wieweit man bei der Klärung hier alleine mit dem vorhandenen Quelltext kommt. Das Programm spuckt auch Meldungen aus, wenn sich ein Testkörper über eine bestimmte Entfernung von der Sonne hinausbewegt. Wo die Schranke dabei liegt, kann man im Quelltext nachlesen. Ich werde auch versuchen, die anderen Fragen zu klären. Es kann nur sein, dass man für ein 100%iges Verständnis auch die Hilfe der Entwickler benötigt.
Gruß

EDIT: Das Programm gibt während des Laufes immer wieder Meldungen aus, wieviele der anfangs definierten Testkörper noch "aktiv" (was das genau heißt, weiß ich noch nicht) sind. Sinkt diese Zahl auf Null bricht das Programm aktuell die weitere Berechnung ab.

Bernhard
18.09.2012, 12:45
Was heisst, "eingefangen"? Der permanente Einfang eines masselosen Teilchens sollte eigentlich (aus Symetriegründen) nicht möglich sein.
Stimmt. "Eingefangen" ist erst mal Unsinn. Man müsste aber klären, wie das Programm die Hill-Sphären (s.a. Wikipedia) berücksichtigt und ob das dann Kollisionen sind.

Man müsste auch versuchen zu klären wie das Programm mit Swing-By-Effekten umgeht. Per Swing-By können in der Realität Testkörper ja deutlich beschleunigt und auch abgebremst werden. Für den Anfang ist die Funktion der Hill-Sphären aber natürlich erst mal wichtiger.

SRMeister
18.09.2012, 12:53
die Testkörper können durch verschiedene Szenarien verschwinden. Ihre "Geschichte" wird in 2 Arrays festgehalten: istat und rstat. istat ist 2-dimensional und die erste Dimension ist der Index des Testpartikels TP.
Beispiel: ist istat[TP][0] == 0 dann bedeutet dass, das der TP noch aktiv ist. Dann kann man istat[TP][1] auswerten: ist der Wert 0, dann ist der TP momentan nicht im Nahbereich eines anderen massebehafteten Körpers (Planets). Ist der Wert +n, dann ist der TP im äußeren Bereich des Planeten mit der Nr. n. Ist der Wert -n, dann ist der TP im inneren Bereich des Planeten n.
istat[TP][2] enthält die Nr des Planeten mit dem es zuletzt ein "encounter" gab. Encounter = innerhalb Hillsphäre.
Man kann da noch weiter ins Detail gehen, wie zb, wie oft war der TP bei welchem Planeten usw... aber ich denke das ist erstmal nebensächlich. Interessant ist noch folgendes:
Ist istat[TP][0] == 1 dann bedeutet dass, das der TP inaktiv ist. Dann steht die Zahl in istat[TP][1] für das, was mit ihm passiert ist:

==> if istat[TP][1] = -1 Danby did not converge
==> if istat[TP][1] = -2 a<0 & r_helio>rmax_a
==> if istat[TP][1] = -3 r_helio>rmax
==> if istat[TP][1] = -4 q<qmin
==> if istat[TP][1] = 1 r_helio<rmin
==> if istat[TP][1] = n too close to planet n
Bitte nicht mich fragen was Danby ist.
SCATR hat aber noch andere Werte in diesen Status arrays, nachzulesen in README.stat

rstat enthält im Wesentlichen zu den Informationen in istat die Abstände zu den entsprechenden Planeten

Bynaus
18.09.2012, 13:25
Man müsste aber klären, wie das Programm die Hill-Sphären (s.a. Wikipedia) berücksichtigt und ob das dann Kollisionen sind.

Der Name des Parameters (r2hill) legt offenbar nahe, dass er irgendwas mit der Hillsphäre zu tun hat. Nicht jeder Körper, der in die Hillsphäre eindringt, sollte zu einer Kollision führen (aber natürlich sind Swing-Bys möglich und sogar wichtig). Ob es zu einer Kollison kommt, sollte vom Radius der Gravitativen Fokussierung abhängen.

@SRMeister: r ist vermutlich die grosse Halbachse, q der Perihel der Bahn - irgendwo kann man dafür wohl Maximal/Minimalwerte festlegen. Danby wird irgend eine Näherungsfunktion sein. Interessant wäre auch zu sehen, was "too close to planet n" meint, dh, ab wann ein Teilchen als "too close" betrachtet wird.

Bernhard
18.09.2012, 13:33
Bitte nicht mich fragen was Danby ist.
Danke für die Erklärungen zu den Status-Flags.

Bernhard
18.09.2012, 13:36
Der Name des Parameters (r2hill) legt offenbar nahe, dass er irgendwas mit der Hillsphäre zu tun hat.
r2hill wird mit der Formel für den Hill-Radius (http://de.wikipedia.org/wiki/Hill-Sph%C3%A4re) berechnet. Es ist ein programmglobales Array und wird für jeden Planeten bei der Initialisierung einmal berechnet.

Bernhard
18.09.2012, 19:24
Danby wird irgend eine Näherungsfunktion sein.
Danby ist der Nachname des Autors, bzw. des Erfinders des zugehörigen Algorithmus:

/* ALGORITHM: Based on pp. 70-72 of Fitzpatrick's book "Principles of */
/* Cel. Mech. ". Quartic convergence from Danby's book. */

SRMeister
18.09.2012, 21:40
Hallo Bernhard, ich vermute die verbliebenen Probleme haben ihre Ursache hier:
Link (http://www.yolinux.com/TUTORIALS/LinuxTutorialMixingFortranAndC.html)

Order of multi dimensional arrays in C/C++ is the opposite of FORTRAN.

Bernhard
19.09.2012, 07:07
Hallo Bernhard, ich vermute die verbliebenen Probleme haben ihre Ursache hier:
Link (http://www.yolinux.com/TUTORIALS/LinuxTutorialMixingFortranAndC.html)
Hallo SRMeister,

das hatte ich beim Debuggen schon bemerkt und entsprechend angepasst. Es sind sehr wahrscheinlich noch Fehler in der Verwaltung der Status-Flags, denn die Ausgaben sehen an manchen Stellen noch ziemlich verdächtig aus. Ich werde/muss die nächsten paar Tage wieder eine "kreative Pause" einlegen, um den Kopf wieder frei zu bekommen. Das Programm ist zwar sehr interessant aber insgesamt auch ziemlich groß.
Gruß

SRMeister
19.09.2012, 10:22
alles klar, ich bin gerade dabei die hybride Version zu machen. Hänge aber auch grad etwas fest.

SRMeister
19.09.2012, 14:11
Hallo Bernhard, ich habe zwar nicht deine neueste Version, aber vllt hast du diesen Fehler ja noch nicht gefunden:


int rmvs3_hel(long &nbod, long &npl, long &ntp,
double &time, double *mass, double &j2rp2, double &
j4rp4, double *xh, double *yh, double *zh, double *
vxh, double *vyh, double *vzh, double *xht, double *
yht, double *zht, double *vxht, double *vyht, double *
vzht, long istat[NTPMAX][NSTAT], double &dt, double &tinc, long *
iclose, long vnclose, double &rcrit);

long vnclose sollte sein long &nclose

SRMeister
19.09.2012, 14:30
Man darf in C/C++ die großen Arrays nicht auf dem Stack erstellen(das ist nur ein kleiner Speicherbereich direkt auf der CPU - dafür sehr schnell). Deshalb hast du vermutlich die Grenze von Testpartikel auf 301 runtergesetzt. Mit dem Schlüsselwort "static" werden die Arrays im RAM angelegt, dann können die beliebig groß sein. Beispiel:

static int istat[NSTAT][NTPMAX] = {0};
Mit ={0} wird das Array mit Nullen gefüllt (ist wahrscheinlich zwar überflüssig aber beugt Fehlern vor)

Die Hybridversion ist nun fertig und lauffähig. Sie erzeugt bei mir die gleichen Ausgaben wie die reine Fortranversion. Ich habe mich gegen eine dynamische und für eine statische Bibliothek entschieden. Das bedeutet man muss lediglich die .lib und die .h zu dem C++Projekt hinzufügen. Man erhält am Ende also eine einzige .exe Datei.(bei mir wirds aber wahrscheinlich eine .DLL die ich dann in das Windows-Prog einbinde)

Bernhard
19.09.2012, 14:38
Mit dem Schlüsselwort "static" werden die Arrays im RAM angelegt, dann können die beliebig groß sein
Das ist interessant und ich sollte mal ausprobieren, ob ich das nutzen kann. In meiner aktuellen Version wird der Speicher im Hauptprogramm mit new angelegt und die Zeiger an die entsprechenden Funktionen übergeben. Das funktioniert zwar, ist aber ein wenig unübersichtlich, weil dann mehrere Zeiger über verschiedene Funktionen weiter gegeben werden müssen.

Bernhard
19.09.2012, 18:42
long vnclose sollte sein long &nclose
Danke für den Tip, aber diese Variable wird in der Funktion nicht weiter verwendet.

Bernhard
19.09.2012, 21:14
Die Hybridversion ist nun fertig und lauffähig. Sie erzeugt bei mir die gleichen Ausgaben wie die reine Fortranversion. Ich habe mich gegen eine dynamische und für eine statische Bibliothek entschieden. Das bedeutet man muss lediglich die .lib und die .h zu dem C++Projekt hinzufügen. Man erhält am Ende also eine einzige .exe Datei.(bei mir wirds aber wahrscheinlich eine .DLL die ich dann in das Windows-Prog einbinde)
Prima. Schaffst Du zusätzlich auch den Einbau der Startbedingungen? Falls Nein, könnte ich mir das auch anschauen. Wird aber etwas dauern, weil ich die nächsten Tage auch noch mit der Cpp-Version experimentieren/debuggen möchte.

SRMeister
19.09.2012, 21:30
Das würde bedeuten, wir müssten Funktionen haben, um einzelne Körper hinzuzufügen oder? Also so wie die alten AddPlanet() Funktionen? Und wenn wir das haben, müssen wir nur aus der alten Version den vorhandenen Algorithmus für die TPs kopieren, richtig?
Was mir da kopfzerbrechen bereitet sind vor allem die Einheiten. Am besten wäre es doch wenn wir alles über das metrische System eingeben könnten und es dann intern umgewandelt wird ? Also dass AddPlanet() mit Meter und Meter/Sekunde arbeitet?
Ich habe aber noch nicht nachgedacht, wie eine etwaige Umwandlung aussieht.

Wenn du dein Okay gibst, würde ich meine Version dir zuschicken und eine Anleitung für Visual C++ Express Edition schreiben, aber vllt. verwendest du ja auch lieber eine andere Version (Professional), also gib mir nochmal bescheid.

Bernhard
19.09.2012, 22:46
Das würde bedeuten, wir müssten Funktionen haben, um einzelne Körper hinzuzufügen oder?
Es würde sich hier sehr stark anbieten, den Fortran-Code einfach etwas zu erweitern. Dort gibt es eine Variable npl (=number of plantes) für die Anzahl der Planeten und eine andere Variable ntp (=number of test particles) für die Anzahl der Testkörper. Die Dynamik der Planeten und Testkörper wird dann über Arrays für die dreidimensionale Position und Geschwindigkeit modelliert (=beschrieben+berechnet). ntp gibt also an, wieviele Werte der Arrays für Position und Geschwindigkeit zu berücksichtigen sind. Diese müssen am Anfang mit den korrekten Werten belegt werden und werden dann vom Programm für jeden Zeitschritt neu berechnet, ganz so wie bei unserem Bulirsch-Stoer-Code und eigentlich sehr komfortabel. Die maximale Größe aller Arrays wird über die Defines in den Header-Dateien festgelegt und darf bei der Initialisierung nicht überschritten werden.


Also so wie die alten AddPlanet() Funktionen? Und wenn wir das haben, müssen wir nur aus der alten Version den vorhandenen Algorithmus für die TPs kopieren, richtig?
Aus der alten Version brauchen wir nur die Berechnung der Startbedingungen. Alles andere haben die SCATR-Autoren sehr wahrscheinlich deutlich "besser", bzw. fortschrittlicher gemacht :o .


Was mir da kopfzerbrechen bereitet sind vor allem die Einheiten. Am besten wäre es doch wenn wir alles über das metrische System eingeben könnten und es dann intern umgewandelt wird ? Also dass AddPlanet() mit Meter und Meter/Sekunde arbeitet?
Umgekehrt wäre es halt deutlich weniger Arbeit. Die Startbedingungen in SCATR-Einheiten auszurechnen, habe ich ja bereits gemacht. Das ist eigentlich ziemlich trivial und das Umwandeln der Ergebnisse in SI-Einheiten ist ebenfalls trivial.


Wenn du dein Okay gibst, würde ich meine Version dir zuschicken und eine Anleitung für Visual C++ Express Edition schreiben, aber vllt. verwendest du ja auch lieber eine andere Version (Professional), also gib mir nochmal bescheid.
Warum nicht? Schicke mir bitte die Version für Visual C++ Express zu.
Gruß

Kibo
19.09.2012, 23:49
Nabend,

Wer möchte kann sich mal meinen kleinen Simulator anschauen

hier zu finden

(https://skydrive.live.com/embed?cid=58397D4453EFC662&resid=58397D4453EFC662%21154&authkey=AIevkdvIeTreL_g)Das Programm erzeugt eine Planeten.txt mit Beispielwerten die angepasst werden kann.

mfg

Bernhard
20.09.2012, 10:11
"STATUSBERICHT:"

Hallo SRMeister, (CC Bynaus),

habe heute im Cpp-Code noch einen Fehler bei der Verwaltung der Array-Indizes gefunden und behoben. Die Ergebnisse nähern sich damit wieder ein gutes Stückchen an die Fortran-Ergebnisse. Ganz perfekt ist der Code aber vermutlich immer noch nicht, was unser gemeinsames Projekt aber nicht unbedingt aufhalten muss, weil ich voraussichtlich spätestens nächste Woche auch den Fortran-Code für unsere Bedürfnisse anpassen kann. Unsere speziellen Startbedingungen machen Änderungen lediglich an zwei Dateien notwendig. Die Schnittstellen zur Visualisierung werden sich dadurch auch nicht mehr ändern ;) . Du könntest dann einfach ein Update auf die Hybridversion machen und diese dann neu compilieren. Wir hätten dann die gewünschten Berechnungsfunktionen und könnten dann erste Modelltests rechnen und mit Bynaus diskutieren.

BTW: Vielen Dank noch für den static Tip von oben. Genau so etwas hatte ich gesucht :o .

SRMeister
20.09.2012, 12:55
Download Link (Sourceforge) (http://sourceforge.net/projects/theiasim/files/scatr-cpp.zip/download)
Hier die Vorgehensweise für Visual C++ 2010 Express:
Neues Projekt erstellen: Win32 Konsolenanwendung
In den Projektordner die Zipdatei kopieren und entpacken.
Zum Projekt die Dateien scatr.cpp und scatr-lib.h hinzufügen und die originale .cpp Datei löschen, sonst gibt es 2 main-funktionen.
Dann die .lib Datei zum Projekt hinzufügen: Projekt-> ***-Projekteigenschaften -> Linker -> Eingabe -> Zusätzliche Abhängigkeiten -> dort die Datei Scatr-StaticLib.lib eintragen
Dann die .lib und .h Dateien der Intelbibliotheken hinzufügen:
Projekt-> ***-Projekteigenschaften -> VC++-Verzeichnisse -> Includeverzeichnisse : Includepfad hinzufügen
Projekt-> ***-Projekteigenschaften -> VC++-Verzeichnisse -> Bibliothekenverzeichnisse : .Libpfad hinzufügen
Includepfad (bei mir): C:\Program Files %28x86%29\Intel\Composer XE 2013\compiler\include\ia32;$(IncludePath)
Libpfad (bei mir): C:\Program Files %28x86%29\Intel\Composer XE 2013\compiler\lib\ia32;$(LibraryPath)

-Es ist auch möglich, statt den beim Intel Composer mitgelieferten Libraries, die kostenlos frei verfügbaren "Intel redistributable Libraries" zu verwenden. Link (http://software.intel.com/en-us/articles/redistributable-libraries-of-the-intel-c-and-fortran-compiler-for-windows/)
-Alle Änderungen in den Projekteigenschaften müssen getrennt für Debug und Release vorgenommen werden.
-Die .Lib habe ich mit maximaler Gleitkommagenauigkeit erstellt.
-Es sind nicht alle SCATR-Funktionen exportiert, nur die bis jetzt benötigten, aber andere hinzuzufügen ist einfach.
-Es liegen Bernhard seine Theia-Startbedingungen bei und ist voreingestellt auf 1Mio Jahre mit 200TPs.(Dauer bei mir etwa 5min?)

SRMeister
21.09.2012, 09:30
Bernhard hat mich darauf hingewiesen, dass ich die Parameterdateien für die Startbedingungen vergessen hatte. Ich habe die .zip Datei aktualisiert. Die Dateien müssen in den selben Ordner entpackt werden, in dem auch die Quellcodedateien sind. Dann sollte das Programm lauffähig sein.

Ich würde mich wirklich über Rückmeldungen freuen.

Stefan

Bynaus
21.09.2012, 09:56
Hallo SRMeister/Stefan, vielen Dank. Ich werde mir das Programm morgen gerne anschauen, da hab ich (hoffentlich) Zeit.

Bernhard
21.09.2012, 13:51
Ich werde mir das Programm morgen gerne anschauen, da hab ich (hoffentlich) Zeit.
Hallo Bynaus,

wenn Dich der Quelltext nicht näher interessiert, wäre eine bereits compilierte (Release-)Version vermutlich besser. Das wären dann nur vier Dateien. Die ausführbare Datei WinSCATR.exe und die drei Parameter- und Startwerte-Dateien. Du könntest damit unmittelbar testen. Gib mir einfach per PN bescheid, wenn ich Dir diese vier Dateien separat zuschicken soll.
Gruß

Bernhard
22.09.2012, 10:22
Hallo SRMeister,

kennst Du eventuell brauchbare Fortran77-Funktionen zum Erzeugen von Zufallszahlen. Ich habe etliche aus dem Internet (google) ausprobiert, aber die kennt der Intel-Compiler scheinbar alle nicht. Aktuell verwende ich srand() zum Initialisieren und rand() zum Erzeugen der Zufallszahl. Dabei bekommt man dann die folgende Fehlermeldung

error #6404: This name does not have a type, and must have an explicit type. [RAND]
Definiert man den Typ der Funktion als real*8 kann man zwar alles compilieren, aber der Debugger bricht dann bei der Funktion rand() kommentarlos ab. Mich wundert dabei, dass die Funktion srand akzeptiert wird :confused:
Gruß

SRMeister
22.09.2012, 14:08
der wichtige Teil im folgenden Code ist der Befehl "USE IFPORT"

program test_rand
USE IFPORT
implicit none
write (*,*) rand()
end

Bernhard
22.09.2012, 20:45
der wichtige Teil im folgenden Code ist der Befehl "USE IFPORT"
Dazu gehören laut Intel Homepage auch einige Compiler-Einstellungen. Dummerweise musste ich heute meine Version von Visual Studio 2010 Professional wieder deinstallieren, weil im Hintergrund soviele Prozesse gestartet wurden, dass man kaum noch die Maus bewegen konnte.

SRMeister
23.09.2012, 08:01
Bernhard, ich bin mir nicht sicher ob die ganzen Hintergrunddienste wirklich vom Studio kamen. Aber du hast schon recht, das Ding installiert allerhand Unsinn. Vor allem irgendwelche SQL Server usw..
Ich konnte aber auf meinem alten Laptop das Teil auch zum Laufen kriegen, ohne ruckeleien, wenn man den entsprechenden Unrat entsorgt.
Was hast du denn eigentlich für ein System? Meins: Intel Q6600 auf 3Ghz, 6 GB RAM und Win7x64 (die Hardware ist von 2007)

Bernhard
23.09.2012, 08:33
Was hast du denn eigentlich für ein System? Meins: Intel Q6600 auf 3Ghz, 6 GB RAM und Win7x64 (die Hardware ist von 2007)
Ich verwende ein WindowsXP-Home (32-Bit) auf einem Athlon-LE-1640-Rechner (64-Bit) mit 1,5GB RAM.

SRMeister
24.09.2012, 16:10
Das ist ungefähr das selbe wie mein alter Laptop und dort habe ich auch starke Probleme mit VS2010.
Vielleicht können wir ja eine gute kostenlose Entwicklungsumgebung finden? Die haben ja zumeist wenig features und sind dadurch auch schneller. In der Schule haben wir mal mit Dev C++ (http://de.wikipedia.org/wiki/Dev-C%2B%2B) gearbeitet. Was hälst du davon?

Bernhard
24.09.2012, 18:02
Was hälst du davon?
Hallo SRMeister,

das ist ein interessanter Vorschlag. Der Gnu-Cpp-Compiler ist langfristig eventuell wirklich die bessere Wahl. Allerdings läuft inzwischen Visual Studio auf meinem Rechner auch ganz gut. Die erneute Installation hat scheinbar irgendetwas verändert. Eine weitere, wichtige Frage wäre auch, ob Dev-C++ mit Fortran klar kommt. BTW: Hast Du meine eMail erhalten? Ich habe Dir zwei Fortran-Dateien auf Deine eMail-Adresse geschickt.

Ich habe heute übrigens noch das Buch von Danby, "Fundamentals of Celestial Mechanics" bestellt, um den Einstieg in die Thematik etwas zu erleichtern. Da das Buch aus USA verschickt wird, dauert der Versand aber voraussichtlich bis Mitte Oktober.
Gruß

SRMeister
24.09.2012, 20:42
ob Dev-C++ mit Fortran klar kommt.
Nein ich denke nicht, aber irgendeine IDE, die Fortran und C++ kann, muss es doch geben?? (Außer Visual Studio mit Intel Composer)

BTW: Hast Du meine eMail erhalten? Ja die habe ich erhalten, ich könnte im Prinzip die Funktion mit in die Bibliothek kompilieren, aber eigentlich wäre es doch gut, wenn man die Funktion in C++ hätte? Sagmal, bist du jetzt auf Fortran umgestiegen??? :) :)

Zur Problematik x64 und x86: ich kompiliere meine Programme generell für x86, die laufen auf Win7x64 ja auch. Bynaus und allen anderen stillen Mitlesern könnte man ja eine x86 .exe zur Verfügung stellen? ich glaube auch nicht dass das Prog mit 64 Bit merkbar schneller wird.
Gruß

@Kibo: Das Programm gefällt mir richtig gut, ich bin gespannt wie das mit OpenGL aussehen wird! Sehe ich das richtig, dass du bis jetzt nur in 2 Dimensionen rechnest?

Kibo
24.09.2012, 22:00
Ja das siehst du richtig, ist ja im Moment noch ein Prototyp, sollte kein großes Problem sein da noch eine extra Dimension mit einzufügen.
Die Linien sind gdi+, im Prinzip ist der Umstieg auf OpenGl kein Problem und für eine ansprechende 3D-Ausgabe wohl auch nötig. Im Moment habe ich noch ein andere Projekt am laufen, daher läuft der Simu eher auf Sparflamme.
Ich plane für das Programm noch ein GUI an der linken Seite zum manipulieren der Planeten in Echtzeit, sowie mehrere Einstellungsmöglichkeiten für Ansicht, Geschwindigkeit, Zentrierung, Zoom und natürlich wären beschriftete Gitternetzlinien auch nicht schlecht. Man soll ja die Umlaufbahnen auch mit der Realität irgendwie vergleichen können.

Bernhard
24.09.2012, 22:39
Nein ich denke nicht, aber irgendeine IDE, die Fortran und C++ kann, muss es doch geben?? (Außer Visual Studio mit Intel Composer)
Das wäre wichtig, weil wir so schnell auf das Fortran-Programm nicht verzichten können. Das Debuggen des Cpp-Codes ist mittlerweile so kompliziert geworden (alle offensichtlichen Fehler sind behoben), dass ich mich erst mal mit dem Programm selbst beschäftigen werde, angefangen bei dem Buch von Danby


Zur Problematik x64 und x86: ich kompiliere meine Programme generell für x86, die laufen auf Win7x64 ja auch.
Das wäre gut. Ich habe da allerdings auch schon Gegenteiliges erlebt.


Bynaus und allen anderen stillen Mitlesern könnte man ja eine x86 .exe zur Verfügung stellen?
Bynaus habe ich bereits eine x86.exe geschickt und ich hoffe, dass er damit bald erste Tests machen kann.


ich glaube auch nicht dass das Prog mit 64 Bit merkbar schneller wird.
Um Geschwindigkeit geht es gar nicht, sondern darum möglichst bald erste wissenschaftlich verwertbare Ergebnisse zu bekommen.

@Bynaus: Lass mich bitte kurz zusammenfassen, was wir bisher erreicht haben:

Wir haben ein Programm (SCATR), welches das Sonnensystem für mindestens eine Millionen Jahre simulieren kann und das maximal 30 Tausend Testkörper mitberechnen kann. Über die Startbedingungen können wir die Anzahl der Testkörper festlegen. Wir können den Eintrittswinkel von Theia relativ zur Ekliptik festlegen und wir können den maximalen Streuwinkel der Trümmer festlegen.

Welche Ziele wollen wir uns ausgehend davon setzen?
a) Soll die Simulationszeit vergrößert werden? -> In den Veröffentlichungen seit etwa 1980 wurden Zeiten bis 3 Milliarden Jahre berechnet.
b) Wie tiefgreifend müssen wir uns mit der Stabilität des Sonnensystems über diese Zeiträume beschäftigen? -> Etliche Veröffentlichungen zeigen bereits sehr detaillierte Ergebnisse zu diesem Thema, wie z.B. die Stabilität der Plutobahn, die Stabilität der inneren Planeten usw.

Bernhard
24.09.2012, 22:44
@Kibo: Das Programm gefällt mir richtig gut, ich bin gespannt wie das mit OpenGL aussehen wird! Sehe ich das richtig, dass du bis jetzt nur in 2 Dimensionen rechnest?
Hallo,

könntet Ihr dieses Parallelthema bitte in einem eigenen, bzw. neuen Thema verfolgen? Das Engagement von Kibo ist zwar sehr zu begrüßen, aber ich empfinde es ehrlich gesagt als ziemlich störend, wenn hier immer wieder themenfremde Einschübe kommen.
Gruß

Bernhard
24.09.2012, 23:12
Zur Problematik x64 und x86: ich kompiliere meine Programme generell für x86, die laufen auf Win7x64 ja auch. Bynaus und allen anderen stillen Mitlesern könnte man ja eine x86 .exe zur Verfügung stellen?
Hallo SRMeister,

es wäre auch für unser Projekt hier sehr hilfreich auf der sourceforge-Seite sowohl eine 32-Bit, als auch eine 64-Bit Variante des Programmes zu veröffentlichen, weil meine 32-Bit-Version laut Bynaus nicht auf seinem 64-Bit Rechner läuft. Es wäre also hilfreich, wenn Du mit Deiner 64-Bit-Software ebenfalls eine Release-exe erzeugen könntest. Ich schicke Dir dann meine 32-Bit-Release-exe zu und Du könntest dann beide Dateien auf Deiner sourceforge-Seite veröffentlichen.

Würdest Du das machen?

Und wie gesagt. Solche "Problemchen" zwischen 32-Bit und 64-Bit sind alles andere als ungewöhnlich.

SRMeister
24.09.2012, 23:51
Würdest Du das machen?
klar, dann schick mal her das Teil :)
Ich mach dann mal eine 64 bit Version mit den Dateien die du mir zugeschickt hast.

SRMeister
25.09.2012, 00:09
ich konnte das Problem mit der rand() Funtion lösen. Hier mal alle Zeilen die was damit zu tun haben:


....
USE IFPORT
....
integer iseed
real*8 rran
....
iseed = 35791246
call srand(iseed)
....
rran = rand()
y = rran * (y_max - y_min) + y_min
....

Ich musste keine Compilereinstellungen ändern.

Bernhard
25.09.2012, 07:31
Hallo SRMeister,

die genaue Position der USE IFPORT-Direktive war für meinen Compiler sehr wichtig. Könntest Du deswegen bitte die gesamte Datei als Forenbeitrag mit den Code-Tags veröffentlichen?
Gruß

Bynaus
25.09.2012, 08:06
@Bynaus: Lass mich bitte kurz zusammenfassen, was wir bisher erreicht haben:

Wir haben ein Programm (SCATR), welches das Sonnensystem für mindestens eine Millionen Jahre simulieren kann und das maximal 30 Tausend Testkörper mitberechnen kann. Über die Startbedingungen können wir die Anzahl der Testkörper festlegen. Wir können den Eintrittswinkel von Theia relativ zur Ekliptik festlegen und wir können den maximalen Streuwinkel der Trümmer festlegen.

Welche Ziele wollen wir uns ausgehend davon setzen?
a) Soll die Simulationszeit vergrößert werden? -> In den Veröffentlichungen seit etwa 1980 wurden Zeiten bis 3 Milliarden Jahre berechnet.
b) Wie tiefgreifend müssen wir uns mit der Stabilität des Sonnensystems über diese Zeiträume beschäftigen? -> Etliche Veröffentlichungen zeigen bereits sehr detaillierte Ergebnisse zu diesem Thema, wie z.B. die Stabilität der Plutobahn, die Stabilität der inneren Planeten usw.

Yay, ich kann auch wieder was beitragen! :)

Zuerst eine Frage: Was meinst du genau mit dem "Eintrittswinkel" von Theia bzw. dem Streuwinkel der Trümmer? Wie genau gehen sie in die Simulation ein? Verstehe ich das richtig: das Programm simuliert nun also nicht mehr eine sphärische Verteilung aller Trümmer von der Position der Erde aus? Die Sache ist natürlich die, dass wir keine Ahnung haben (können) in welchem Winkel zur Ekliptik Theia die Erde getroffen hat, bzw. in welchem Streuwinkel die Trümmer davonflogen. Je nach Szenario wird da was anderes naheliegen, aber belegen kann man nichts davon. Deshalb würde ich vorschlagen, grundsätzlich eine sphärische Verteilung der Trümmer anzunehmen. Nicht, weil das realistisch wäre, sondern weil uns dies erlaubt, alle möglichen Fälle abzudecken (sonst müssten wir ja alle möglichen Eintrittswinkel und Streuwinkel separat durchprobieren). Dann sehen wir, ob diese Winkel überhaupt irgend eine Konsequenz haben (z.B.: "Die meisten Trümmer, die den Mars treffen, kommen von der sonnenzugewandten Seite der Trümmerwolke" vs. "Alle Trümmer haben etwa die gleiche Chance, Mars zu treffen").

Dann zu deinen Fragen.

a) Nein, die meisten Trümmer werden ja eher klein sein, selbst wenn sie in Asteroidenform vorliegen. Das heisst, die Sonne wird ihre Bahnen in wenigen Millionen Jahren auf einen Kollisionskurs mit sich selbst zwingen (Yarkovsky, YORP Effekte). Die Simulation sollte deshalb an sich nicht über ein paar Millionen Jahre hinaus gehen, andernfalls müssen wir diese Kräfte einbeziehen, und dann wirds wirklich kompliziert.
b) Da es ja um die Frage geht, was mit den Trümmern passiert ist, muss uns die Stabilität des Planetensystems nicht weiter gross beschäftigen (es sei denn, wir könnten zeigen, dass massive Trümmer einer bestimmten Grösse dieses destabilisieren - das wäre natürlich ein wichtiges Ergebnis).

Ich könnte mir einen Fall vorstellen, wo die Simulation über längere Zeit nützlich sein könnte: Wir sollten, um das obere Ende des Trümmerspektrums abzudecken, ein paar Simulationen mit massiven Trümmern machen (einzeln, oder bis zu maximal 10 davon): vielleicht von der Masse des Mondes, dh, 1% der Erdmasse, und schauen, welchen Einfluss diese auf das Planetensystem haben bzw mit welchen Himmelskörpern sie wann ungefähr kollidieren. Eine solche Simulation müsste dann wohl über deutlich längere Zeiträume laufen, evtl. 4.5 Mrd Jahre oder zumindest so lange, bis alle massiven Trümmer das System verlassen haben oder mit der sonne kollidiert sind.

Bernhard
25.09.2012, 09:19
Was meinst du genau mit dem "Eintrittswinkel" von Theia bzw. dem Streuwinkel der Trümmer? Wie genau gehen sie in die Simulation ein? Verstehe ich das richtig: das Programm simuliert nun also nicht mehr eine sphärische Verteilung aller Trümmer von der Position der Erde aus? Die Sache ist natürlich die, dass wir keine Ahnung haben (können) in welchem Winkel zur Ekliptik Theia die Erde getroffen hat, bzw. in welchem Streuwinkel die Trümmer davonflogen. Je nach Szenario wird da was anderes naheliegen, aber belegen kann man nichts davon. Deshalb würde ich vorschlagen, grundsätzlich eine sphärische Verteilung der Trümmer anzunehmen. Nicht, weil das realistisch wäre, sondern weil uns dies erlaubt, alle möglichen Fälle abzudecken (sonst müssten wir ja alle möglichen Eintrittswinkel und Streuwinkel separat durchprobieren). Dann sehen wir, ob diese Winkel überhaupt irgend eine Konsequenz haben (z.B.: "Die meisten Trümmer, die den Mars treffen, kommen von der sonnenzugewandten Seite der Trümmerwolke" vs. "Alle Trümmer haben etwa die gleiche Chance, Mars zu treffen").
OK. Das ist eine gute Idee. Momentan werden die Trümmer symmetrisch zur Ekliptik in einem Bereich von +/- 20° gestreut. Aufgrund von Ralfs Einwand, dass manche Bahnen von heutigen Kleinplaneten relativ stark von der Ekliptik abweichen, sollte man diesen Bereich also auf +/- 40° vergrößern.


Ich könnte mir einen Fall vorstellen, wo die Simulation über längere Zeit nützlich sein könnte: Wir sollten, um das obere Ende des Trümmerspektrums abzudecken, ein paar Simulationen mit massiven Trümmern machen (einzeln, oder bis zu maximal 10 davon): vielleicht von der Masse des Mondes, dh, 1% der Erdmasse, und schauen, welchen Einfluss diese auf das Planetensystem haben bzw mit welchen Himmelskörpern sie wann ungefähr kollidieren. Eine solche Simulation müsste dann wohl über deutlich längere Zeiträume laufen, evtl. 4.5 Mrd Jahre oder zumindest so lange, bis alle massiven Trümmer das System verlassen haben oder mit der sonne kollidiert sind.
Ich werde dann nochmal testen, wie gut der SCATR-Code mit dem Mond zurecht kommt. Falls solche gebundenen Bahnen ebenfalls stabil berechnet werden, kann ich x Trümmer bei der Liste der Planeten einbauen.

Bynaus
25.09.2012, 11:08
Momentan werden die Trümmer symmetrisch zur Ekliptik in einem Bereich von +/- 20° gestreut. Aufgrund von Ralfs Einwand, dass manche Bahnen von heutigen Kleinplaneten relativ stark von der Ekliptik abweichen, sollte man diesen Bereich also auf +/- 40° vergrößern.

Ja, das kann man, wenn genügend Rechenzeit zur Verfügung steht, nur schon um zu sehen, ob das einen Unterschied im Ergebnis macht. Ich erinnere mich, dass Planetenembryos - zumindest in einschlägigen Simulationen - gelegentlich auch relativ hohe Inklinationen haben können.


Ich werde dann nochmal testen, wie gut der SCATR-Code mit dem Mond zurecht kommt. Falls solche gebundenen Bahnen ebenfalls stabil berechnet werden, kann ich x Trümmer bei der Liste der Planeten einbauen.

Okay. Ich würde bei den massiven Trümmern vorschlagen, dass wir deren Anfangsparameter zufällig auswählen.

Vielleicht könnte man drei Serien von Simulationen machen:

1) Nur masselose Trümmer
2) Nur massive Trümmer
3) Sowohl masselose als auch massive Trümmer

Wenn wir eine grobe Anhnung davon haben, was in allen Fällen passiert, kann man entscheiden, wie wir weiter vorgehen wollen, ob wir uns etwa auf einen Aspekt fokussieren wollen, etc.

Bernhard
25.09.2012, 11:52
Okay. Ich würde bei den massiven Trümmern vorschlagen, dass wir deren Anfangsparameter zufällig auswählen.
D.h. Verteilung der Anfangsorte und Anfangsgeschwindigkeiten wie bei den masselosen Trümmern?
Wieviel Trümmer mit Masse sollen es sein?
Welche genaue Massenverteilung sollen die Trümmer mit Masse dann haben?

SRMeister
25.09.2012, 12:24
Dein 32bit Programm läuft bei mir fehlerfrei, ohne dass es nötig wäre, spezielle Kompatibilitätsoptionen zu setzen.
Ich kann mir im Moment nicht so recht vorstellen, woran es liegen könnte, dass es auf einem anderen x64 System nicht läuft.
Benutzt das Programm die gleiche theia_tp_200.in, die ich auch benutze? Weil bei meinem "Kompilat" bleibt eine andere Anzahl Körper übrig.
Ich werde jetzt deine 32bit und meine 64bit Version veröffentlichen, wobei die 64 Bit Version momentan die TPs über deine Funktion selbst erzeugt, die Datei theia_tp.in ist also erstmal überflüssig.
Als nächsten Teilschritt habe ich mir überlegt, möchte ich die Planetenkollision aktivieren. Die ist ja leider noch ausgeschaltet und r(Planet) sind auch noch nicht vorgegeben. Das wird wohl auch der Grund sein warum so viele TPs das System verlassen: sie tunneln vorher durch Planeten und bekommen dabei nen ordentlichen Schubs.

Zur Initfunktion:
Sonst habe ich noch die Verwendung von rand() verändert, so dass bei jeder Benutzung eine neue Zufallsvariable bereit steht.


c Last revision: 9/25/12 S.R.

subroutine io_init_tp_theia(ntp,npl,x_earth,y_earth,z_earth,
& vx_earth,vy_earth,vz_earth,
& xht,yht,zht,vxht,vyht,vzht,istat,rstat,frame)
USE IFPORT
include '../scatr.inc'
include 'io.inc'

c... Input
integer ntp,npl
real*8 x_earth,y_earth,z_earth
real*8 vx_earth,vy_earth,vz_earth
character*80 frame

c... Output
real*8 xht(NTPMAX),yht(NTPMAX),zht(NTPMAX)
real*8 vxht(NTPMAX),vyht(NTPMAX),vzht(NTPMAX)
real*8 rstat(NTPMAX,NSTATR)
integer istat(NTPMAX,NSTAT)

c... Internal
integer iseed
real*8 rran
integer i,j,ierr,ns,ftest
real*8 v_esc,r_earth
real*8 y,phi,vn,y_min,y_max
real*8 nx,ny,nz
integer SCATTERING_ANGLE

c-----
c... Executable code

c Initialize
c 1 AU=149597870691 m

c Initialize random generator
iseed = 35791246
call srand(iseed)

SCATTERING_ANGLE = 20
v_esc = 6.4604555902823027300785018765024e-3 ! AU/day = 11186.0 m/s
r_earth = 4.2563891254523551325458216979349e-5; ! mean radius of earth in AU (=6367467.5 m)

y_min = 0.5 * (1.0 - cos((90.0 - SCATTERING_ANGLE)/PI*180.0))
y_max = 0.5 * (1.0 - cos((90.0 + SCATTERING_ANGLE)/PI*180.0))

if(ntp.gt.NTPMAX) then
write(*,*) ' SWIFT ERROR: in io_init_tp_theia: '
write(*,*) ' The number of test bodies,',ntp,','
write(*,*) ' is too large, it must be less than',NTPMAX
call util_exit(1)
endif

write(*,*) ' '
write(*,*) 'ntp : ',ntp

if(ntp.eq.0) then
write(*,*) ' '
return
endif ! <===== NOTE

c... Determine the number of istat and rstat variables. In what follows,
c... we assume that they are the same.
ns = npl + 2

if(ns.ne.NSTAT) then
write(*,*) 'Warning: The size of istat and rstat arrays is '
write(*,*) ' not NSTAT=',NSTAT,', but is ',ns
endif

c Read in the x's and v's and istat(*,*)
write(*,*) ' '
do i=1,ntp
rran = rand()
y = rran * (y_max - y_min) + y_min ! add random number instead of 0.5d0
rran = rand()
phi = rran * TWOPI ! add random number instead of 0.5d0
rran = rand()
vn = rran * 0.4d0 + 1.0d0 ! add random number instead of 0.5d0

nx = 2.0d0 * sqrt(y - y * y) * cos(phi)
ny = 2.0d0 * sqrt(y - y * y) * sin(phi)
nz = 1.0d0 - 2.0d0 * y

xht(i) = -9.814911942588033e-001 + nx * r_earth
yht(i) = -1.866158863169087e-001 + ny * r_earth
zht(i) = 5.097029871362697e-006 + nz * r_earth

vxht(i) = 2.931839043977104e-003 + vn * v_esc * nx
vyht(i) = -1.697334005001157e-002 + vn * v_esc * ny
vzht(i) = 3.022067795326298e-007 + vn * v_esc * nz

do j=1,ns
istat(i,j) = 0
enddo

do j=ns+1,NSTAT
istat(i,j) = 0
if (j.eq.NSTAT-2.and.
& (frame(1:3).eq.'hel'.or.frame(1:3).eq.'HEL')) then
istat(i,j) = 1
endif
enddo

do j=1,NSTATR
rstat(i,j) = 0.0d0
enddo
enddo

return
end ! io_init_tp_theia.f
c-----------------------------------------------------------------

SRMeister
25.09.2012, 13:00
So die beiden Varianten sind jetzt bei Sourceforge (https://sourceforge.net/projects/theiasim/files/).

Bernhard
25.09.2012, 13:31
Dein 32bit Programm läuft bei mir fehlerfrei, ohne dass es nötig wäre, spezielle Kompatibilitätsoptionen zu setzen.

@Bynaus: Dann fehlen bei Dir sehr wahrscheinlich noch weitere Bibliotheken oder Du hast das falsche Redistributable-Paket erwischt. Wir können vorerst aber auch so weitermachen.

Die Ausgaben des Programmes sind bereits so aussagekräftig, dass man jetzt auch mal nachsehen sollte, wo die masselosen Trümmer im Laufe einer Million Jahre "herunterkommen". Ansonsten erscheint mir ein großes Bruchstück mit der Masse 1% der Erde als sehr wichtig. Es wäre sehr interessant, wenn dieser "Mondkeim" in der Simulation ebenfalls Trümmer einfangen würde. Die Startgeschwindigkeit dieses "Keimes" müsste man in etwa so groß wählen, wie die Geschwindigkeit von Theia. Ich habe da aber wenige Ideen, wie groß diese Geschwindigkeit sein sollte.

BTW: Das Thema nimmt so langsam richtig Fahrt auf, WOW.

Bynaus
25.09.2012, 14:35
D.h. Verteilung der Anfangsorte und Anfangsgeschwindigkeiten wie bei den masselosen Trümmern?
Wieviel Trümmer mit Masse sollen es sein?
Welche genaue Massenverteilung sollen die Trümmer mit Masse dann haben?

1) Ja, so dass wir verschiedene Richtungen mal durchprobieren und schauen, was rauskommt. Evtl nichts, aber das sollten wir zuerst überprüfen. Aber natürlich sollten wir pro Simulationsrun nicht hunderte von massiven Fragmenten haben, sondern nur soviele wie einigermassen realistisch.
2) & 3) Die totale Masse der massiven Trümmer sollte 1 Marsmasse = 10 Mondmassen = 0.1 Erdmassen nicht überschreiten. Wie gross sie allerdings genau sein sollten, ist schwierig zu sagen. Fragmentation bei einem solchen Impakt ist ein schwieriges Thema, das zudem kaum erforscht ist. Ich hätte jetzt - vom Blick auf die Hit-and-run-Impakt-Videos - gesagt, dass die grössten Fragmente stets etwas kleiner als 1 Mondmasse sind. Also nehmen wir doch 1 Mondmasse für das grösste Fragment. Und dann vielleicht noch vier mit 0.1 Mondmassen dazu. Ich würde darüber hinaus die Bahnen (Abstrahlwinkel & Geschwindigkeit) dieser Fragmente einigermassen ähnlich setzen, innerhalb eines Dec/Rect Winkels von nicht mehr als 15°. Das heisst, in jeder Simulation könnte zuerst eine Hauptrichtung festgelegt werden, um die dann die tatsächlichen Fragmentbahnen geringfügig variiert werden. Die Geschwindigkeit der Trümmer sollte sich im Bereich zwischen 1 und 1.4 * v_esc bewegen.

EDIT:


@Bynaus: Dann fehlen bei Dir sehr wahrscheinlich noch weitere Bibliotheken oder Du hast das falsche Redistributable-Paket erwischt.

Ja, mir scheint "MSVCR100.dll" zu fehlen (Meldung nach Installation der 64bit Version auf sourceforge).

SRMeister
25.09.2012, 14:42
ups, das ist Teil des Visual C++ Redistributable Package. Kann man von Microsoft (http://www.microsoft.com/de-de/download/details.aspx?id=14632)hier laden.

Bernhard
25.09.2012, 21:53
2) & 3) Die totale Masse der massiven Trümmer sollte 1 Marsmasse = 10 Mondmassen = 0.1 Erdmassen nicht überschreiten.
Hallo Bynaus,

ich habe gerade den SCATR-Code mit den Startbedingungen für Merkur-Venus-Erde-Mond-Mars-Jupiter-Saturn-Uranus-Neptun vom 01.04.2011, sowie Null Testkörper gestartet und 1 Millionen Jahre in die Zukunft rechnen lassen. Ergebnis: Mond und Mars verlassen das Sonnensystem und vagabundieren irgendwo bei 1.0e5 bis 1.0e6 AU durch das freie All :confused: . Genau so etwas hatte ich ein wenig befürchtet, weil der Ansatz des Wisdom/Holman-Algorithmus davon ausgeht, dass die Bewegung der Planeten in erster Linie von der Sonne bestimmt wird und der Rest lediglich Bahnstörungen sind. Beim Mond trifft das aber gar nicht zu. Dessen Bewegung wird natürlich maßgeblich von der Erde mitbestimmt, so dass man für eine korrekte Bahnberechnung vermutlich den SCATR-Code entsprechend anpassen müsste und das ist keine leichte Aufgabe. Dass auch der Planet Mars das Sonnensystem verläßt finde ich etwas eigenartig, aber vielleicht fehlt hier ja der Asteroidengürtel als Stabilisator?
Gruß

SRMeister
25.09.2012, 22:35
Bernhard, das Problem momentan ist folgendes: In der Literatur findet man die Angabe, dass für eine über Mrd. Jahre reproduzierbare Integration die Zeitschrittweite so eingestellt werden muss, dass der innere Planet eine Umrundung der Sonne in ca. 10 - 20 Zeitschritten macht. Das würde bedeuten, wir müssten etwa 10 Tage einstellen.
Momentan ist aber 1 Jahr eingestellt. Daher kommen die extremen Fehler.
Diese Grundeinstellungen haben die SCATR Programmierer genommen, da sie die Oortsche Wolke plus die 4 Gasriesen simuliert haben, wo Jupiter der innerste ist und 1/10 der Periode eben 1 Jahr ist.
Um die Periode auf 10 Tage einzustellen, musst du den letzten Wert in der erste Zeile in theia_param.in, der momentan 10 ist, auf 360 anheben.
Ich habe das probiert und die Ergebnisse sehen vielversprechend aus. Es werden über die 1 Mio Jahre nurnoch ca. 5% der TPs aus dem System geschleudert.

Zur Problematik mit dem Mond: soweit ich den algorithmus richtig verstanden habe, sucht er für jeden Körper selbstständig den dominierenden Gegenspieler.

Bernhard
25.09.2012, 23:39
Um die Periode auf 10 Tage einzustellen, musst du den letzten Wert in der erste Zeile in theia_param.in, der momentan 10 ist, auf 360 anheben.

Das ist ein korrekter Hinweis, aber es sind falsche Zahlen ;) . Parameter drei in Zeile eins ist die Schrittweite dt in Tagen. Parameter 4 ( = tinc = t_increment) gibt die Änderung der Schrittweite an, falls die Schrittweite angepasst wird (adaptive stepsize). Ich habe Parameter drei also auf 10 (Tage) und Parameter 4 auf 0.1 gesetzt. Damit bleiben Mond und Mars dann auch nach einer Millionen Jahre auf vernünftigen Bahnen.

Bernhard
25.09.2012, 23:46
Also nehmen wir doch 1 Mondmasse für das grösste Fragment. Und dann vielleicht noch vier mit 0.1 Mondmassen dazu.
Frage zum besseren Verständnis: Wie kann aus so einer Anfangsverteilung der heutige Mond mit einer Mondmasse entstehen? Es müsste bei dem Entstehungsprozess des Mondes eine ganze Menge an Material (=0.4 Mondmassen) den Mond wieder verlassen. Kann man diesen Masseverlust über Meteoriteneinschläge erklären?

Bynaus
26.09.2012, 07:47
@Bernhard: Mit unserem Mond haben die Fragmente nichts zu tun. Die Fragmente sind Trümmer von Theia, die zu schnell unterwegs sind, als dass sie im Erde-Mond-System verbleiben könnten. Sie segeln also weiter und umkreisen die Sonne. Die wichtigste Frage hinter diesem Projekt ist, was denn eigentlich mit diesen Trümmern passiert.

In der Zwischenzeit bildet sich aus denjenigen Trümmern, die das Erde-Mond-System (ich nenne es jetzt so, obwohl sich der Mond zu diesem Zeitpunkt erst noch bilden muss) nicht verlassen haben, eine Scheibe rund um die Erde. Aus dieser Scheibe verklumpt sich dann, nach ein paar 100 bis 1000 Jahren vielleicht (? man weiss es nicht so genau) der Mond. Das brauchen wir alles nicht direkt zu simulieren, weil das andere Modelle schon mit viel höherer Auflösung tun. Uns interessiert nur, was mit der ~Marsmasse an Trümmern geschieht, die das Erde-Mond-System verlassen und die Sonne umkreisen. Es geht um folgende Fragen:

1) Wieviel kehrt zur Erde zurück? Je nach dem müsste man das in isotopischen Mischrechnungen zwischen Erde und Theia nämlich berücksichtigen (solche Rechnungen macht man, weil man die Ähnlichkeit der Isotopie des Mondes und der Erde erklären will)
2) Wieviel landet auf anderen Planeten oder allenfalls im Asteroidengürtel? Das wird uns eine Ahnung geben, ob wir in unserer Meteoritensammlung auch Trümmer von Theia erwarten sollten.
3) Welche Auswirkungen haben die grössten Trümmer auf die Bahnen der Benachbarten Planeten? Wenn sich z.B. herausstellen würde, dass die Bahn der Venus dadurch stets viel exzentrischer wird als beobachtet, könnte man dieses Szenario ausschliessen bzw. die Masse des grössten Fragments nach oben begrenzen. Vielleicht stellt sich ja auch heraus, dass das grösste Fragment immer mit dem Mars kollidiert - dann ist es vielleicht für das Borealis-Becken verantwortlich? (ich spinne jetzt nur rum, um zu zeigen, wohin es führen könnte - so etwas zu belegen wäre natürlich eine ganz andere Aufgabe).

SRMeister
26.09.2012, 12:37
Das ist ein korrekter Hinweis, aber es sind falsche Zahlen ;)

Parameter 3 (dt) ist die Schrittweiter außerhalb von 300 AU (Siehe "rcrit" in param.in). Da haben die Autoren berechtigterweise 10 Jahre vorgegeben (3600 stand da).
Außerhalb von "rcrit" wird der Wisdom Holman Integrator verwendet, der sich auf das Barycenter bezieht.
Innerhalb von "rcrit" also unter 300AU wird ein Mixed-Variable-Symplectic Integrator namens RMVS3 verwendet. Die Zeitschritte in diesem Bereich ergeben sich aus "Wert 3" geteilt durch "Wert 4". Deswegen muss Wert 4 auch ein Integerwert sein, da sonst keine runde Teilung möglich ist.

dt is the timestep used when particles are far from the Sun
tinc is the fraction of dt used as an integration step when particles
are near the Sun (must be integer)
Ich versuche das mal zu übersetzen:

dt ist der Zeitschritt, der benutzt wird, wenn die Partikel weit von der Sonne entfernt sind.
tinc ist der Bruchteil von dt, der als Integrationsschritt verwendet wird, wenn die Partikel der Sonne nah sind. (Muss ein Integer sein)

Ich habe mir so meine Gedanken gemacht, ob wir vllt. davon profitieren können. Da wir keinen Kuipergürtel haben, könnten wir den Bereich hinter Neptun, sagen wir mal, ab 35 AU, als den "äußeren Bereich" festlegen. Dort wären dann Zeitschritte gemäß dt = Periode/10 von ca. 15 Jahren sinnvoll. Dort sind ja sogut wie nie irgendwelche Teilchen außer Pluto. Und im inneren Sonnensystem nehmen wir dann 10 Tage. Was hälst du davon?

Bernhard
26.09.2012, 13:34
Da wir keinen Kuipergürtel haben, könnten wir den Bereich hinter Neptun, sagen wir mal, ab 35 AU, als den "äußeren Bereich" festlegen. Dort wären dann Zeitschritte gemäß dt = Periode/10 von ca. 15 Jahren sinnvoll. Dort sind ja sogut wie nie irgendwelche Teilchen außer Pluto. Und im inneren Sonnensystem nehmen wir dann 10 Tage. Was hälst du davon?
Das passt schon. Ich hatte meine Informationen direkt aus dem Quelltext abgeleitet und dabei scheinbar einiges missverstanden. Als nächsten Schritt sollte man jetzt noch die Vorgaben von Byanus mit den fünf massiven Teilkörpern berücksichtigen. Wir können dabei auch gerne parallel arbeiten, weil das relativ schnell gemacht ist. Ein paar Extra-Simulationen mit unterschiedlichen Startwerten können auch nicht schaden.

@Bynaus: Falls Du mit der Software nicht weiter kommst und trotzdem mittesten willst, empfehle ich die Installation von Visual Studio 2010 und den Intel Composer. Du hättest dann die gleiche Software wie wir und mindestens 30 Tage Zeit (wegen des Intels Compilers), um damit zu experimentieren. Wenn das dann auch nicht funktioniert, müsstest Du wohl oder übel das Betriebssystem neu installieren :rolleyes: .

ralfkannenberg
26.09.2012, 15:46
Ergebnis: Mond und Mars verlassen das Sonnensystem und vagabundieren irgendwo bei 1.0e5 bis 1.0e6 AU durch das freie All :confused: .
Hallo zusammen,

mein Kompliment - inzwischen seid Ihr so weit, dass Ihr bereits mit den ersten Testläufen Fehler in den Anfangsbedingungen eruieren könnt. Klasse !!

Ich glaub', das wird bald spannend :)


Freundliche Grüsse, Ralf

Bernhard
26.09.2012, 18:20
Um die Periode auf 10 Tage einzustellen, musst du den letzten Wert in der erste Zeile in theia_param.in, der momentan 10 ist, auf 360 anheben.

Hallo SRMeister,

rechne das vielleicht auch selber noch mal nach, aber bei mir "diffundiert" der Mond bei diesen Einstellungen aus dem Sonnensystem.

Bernhard
26.09.2012, 18:23
mein Kompliment...Klasse !!
Danke.


Ich glaub', das wird bald spannend :)
Ich finde es bereits jetzt recht spannend, aber die Ergebnisse müssen natürlich noch aussagekräftiger werden.
Mit Gruß

Bernhard
27.09.2012, 09:26
Also nehmen wir doch 1 Mondmasse für das grösste Fragment. Und dann vielleicht noch vier mit 0.1 Mondmassen dazu. Ich würde darüber hinaus die Bahnen (Abstrahlwinkel & Geschwindigkeit) dieser Fragmente einigermassen ähnlich setzen, innerhalb eines Dec/Rect Winkels von nicht mehr als 15°. Das heisst, in jeder Simulation könnte zuerst eine Hauptrichtung festgelegt werden, um die dann die tatsächlichen Fragmentbahnen geringfügig variiert werden. Die Geschwindigkeit der Trümmer sollte sich im Bereich zwischen 1 und 1.4 * v_esc bewegen.

Guten Morgen,

ich habe das versuchsweise mal mit einer sehr kleinen Schrittweite (par3 = 10, par4 = 0.1) und 200 masselosen Trümmern gerechnet. Die Ergebnisse sehen vom Überfliegen her so aus, als würden sich sämtliche Trümmer über die Erdbahn verteilen (mit Masse auf jeden Fall und ohne Masse mit Fragezeichen, weil ich nicht alle durchsuchen konnte). Inwieweit die Trümmer zusammenklumpen oder sich gleichmäßig verteilen kann ich im momentan nicht sagen. Keines der Trümmer hat das Sonnensystem verlassen.

Ich zweifle momentan auch ein wenig an manchen Informationen aus der Readme.first-Datei, weil diese Dateien nicht unbedingt immer aktuell sein müssen.

Bynaus
27.09.2012, 09:47
Hallo Bernhard,
Ich hab das Programm noch nicht zum Laufen gekriegt, aber ich arbeite daran (wenn auch langsam...). Aber das klingt schon mal interessant - wie verteilen sie sich denn über die Erdbahn? Es müsste ja vermutlich auch solche geben, deren Bahnen diejenigen der anderen Planeten kreuzen? Kannst du mal ein Video davon (falls möglich?) auf YouTube laden (bzw. mir per Mail zukommen lassen)?

Die "Selbstverklumpung" der Trümmer würde ich jetzt mal ignorieren. Es geht ja erst mal darum, wo sie am Ende alle hinkommen. Wichtig wäre in diesem Zusammenhang die Implementierung von Kollisionen mit den Planeten und allenfalls dem Asteroidengürtel.

Bernhard
27.09.2012, 10:19
Kannst du mal ein Video davon (falls möglich?) auf YouTube laden (bzw. mir per Mail zukommen lassen)?
Da muss ich aktuell leider passen, obwohl mich das selber sehr interessieren würde. Ich hatte ein wenig gehofft, dass das SRMeister machen würde, denn die Cpp-Schnittstellen sind ja jetzt vorhanden.


Wichtig wäre in diesem Zusammenhang die Implementierung von Kollisionen mit den Planeten und allenfalls dem Asteroidengürtel.
Ich muss mir da den Quelltext nochmal genau ansehen, ob so etwas nicht bereits implementiert ist. Falls man Kollisionen mit Asteroiden berechnen will, muss man diese bei den Startbedingungen definieren und das macht das Programm dann mehr und mehr rechenintensiv.

@SRMeister: Wie sieht es denn mit einem weiteren YouTube-Clip aus? Könntest Du so etwas mit dem vorhandenen Material machen? Die drei Parameter-Dateien kann ich Dir gerne zuschicken.

Bynaus
27.09.2012, 11:01
Falls man Kollisionen mit Asteroiden berechnen will, muss man diese bei den Startbedingungen definieren und das macht das Programm dann mehr und mehr rechenintensiv.

Man könnte die Asteroiden ja als masselos annehmen und nur als "Zielkörper" (die mit ihrer Querschnittsfläche einfangen) in die Simulation einbauen. Aber ich würde sagen, das hat vorerst eine geringere Priorität.

Bernhard
28.09.2012, 09:40
Kannst du mal ein Video davon (falls möglich?) auf YouTube laden (bzw. mir per Mail zukommen lassen)?
Hallo Bynaus,

ich werde in der verbleibenden Zeit, wo der Intel Composer bei mir noch läuft versuchen eine Visualisierung der SCATR-Ergebnisse zu programmieren. Die Library von SRMeister enthält dazu bereits recht brauchbare Funktionen.

@SRMeister: Könntest Du mir mal den gesamten Quelltext zur Erstellung der Library schicken? Besser wäre natürlich der gesamte Code Deiner Visualisierungstools. Ich würde mir dann alles mal durchsehen, damit wir hier bei der Interpretation der bereits vorhandenen Ergebnisse weiter kommen.
Gruß

Bernhard
29.09.2012, 11:36
ich habe das versuchsweise mal mit einer sehr kleinen Schrittweite (par3 = 10, par4 = 0.1) und 200 masselosen Trümmern gerechnet. Die Ergebnisse sehen vom Überfliegen her so aus, als würden sich sämtliche Trümmer über die Erdbahn verteilen (mit Masse auf jeden Fall und ohne Masse mit Fragezeichen, weil ich nicht alle durchsuchen konnte). Inwieweit die Trümmer zusammenklumpen oder sich gleichmäßig verteilen kann ich im momentan nicht sagen. Keines der Trümmer hat das Sonnensystem verlassen.

Ich zweifle momentan auch ein wenig an manchen Informationen aus der Readme.first-Datei, weil diese Dateien nicht unbedingt immer aktuell sein müssen.
Der Debugger zeigt, dass der Wert par4=0.1 tatsächlich sehr schlecht gewählt ist. Die Position beispielsweise von Merkur friert damit einfach ein.

EDIT: Mit par3 = 10.0 und par4 = 1.0 bekommt man halbwegs vernünftige Werte. Weitere Tests mit dem Debugger zeigen, dass par3 die Schrittweite in Tagen sein sollte. Anders lassen sich die Ergebnisse nach 1, 2, 3, .... Schritten eigentlich nicht erklären.

Bernhard
30.09.2012, 09:12
Weitere Tests mit dem Debugger zeigen, dass par3 die Schrittweite in Tagen sein sollte. Anders lassen sich die Ergebnisse nach 1, 2, 3, .... Schritten eigentlich nicht erklären.
Mir ist dabei noch nicht klar, ob es besser ist die Schrittweite und den Teiler tinc generell klein zu wählen (z.B. 10 Tage und 1) oder beide groß. Da bei unseren Rechnungen der Teiler praktisch immer verwendet (alle relevanten Distanzen sind kleiner als 300 AU) wird, sollte diese Frage schon geklärt werden. Bei kleinen Schrittweiten sollte der Rechenfehler eigentlich kleiner sein, als bei großen Schrittweiten.

Viele und zum Teil wichtige Dokumenationen zum SCATR-Code und einige Papers zur Entstehung des Mondes findet man übrigens auf der akademischen Seite von Jack Wisdom:

http://groups.csail.mit.edu/mac/users/wisdom/

EDIT: Habe gerade gesehen, dass bei einer Simulation der wichtigsten Planeten mit Erdmond mit dt,tinc = (3600.0, 360.0) der Erdmond nicht von der Erde gebunden wird. Bei (10.0, 1.0) entfernt sich der Erdmond auf etwa 2,7 AU.

Bernhard
02.10.2012, 00:14
Hallo Bynaus und SRMeister,

ich habe heute unseren alten Code (s. Position der Planeten) so erweitert, dass der dortige Integrator zu einem symplektischen Integrator wird. Es waren dazu nur ein paar kleinere Anpassungen nötig. Als Anleitung habe ich diesen Wikipedia-Artikel verwendet: http://en.wikipedia.org/wiki/Symplectic_integrator . Die Mittelpunktsregel nennt sich dann Verlet-Methode. Der Programmablauf wird dadurch mindestens um einen Faktor 1.8 schneller!! Wir könnten vorerst also auch mit diesem Code weitermachen. Der ist im Vergleich zum SCATR-Code zwar noch immer eine "lahme Krücke", aber vielleicht bleibt der Mond bei diesem Programm ja zur Abwechslung mal auf seiner Bahn um die Erde :) . Mein Rechner braucht für eine Millionen Jahre knapp 27 Stunden Rechenzeit. Eine Berechnung über 100.000 Jahre wäre mit knapp 3 Stunden problemlos machbar. Ich denke, ich werde so eine Simulation in den nächsten Tagen versuchsweise mal laufen lassen.
Gruß

Bynaus
02.10.2012, 15:23
Okay. Könnte es nicht sein, dass der Mond seine Bahn verlässt, weil ein Umlauf eben "nur" ein Monat dauert, die Zeitschritte bei der Berechnung also einen grossen Bruchteil seiner Umlaufszeit ausmachen? Der Mond ist in dieser Simulation in dem Sinn nicht so besonders wichtig, er hatte sich ja anfangs noch gar nicht gebildet und später wird er (wegen dem geringen Querschnitt) nur sehr wenige der Trümmer auffangen können. Wenn die Ergebnisse ohne Mond gut aussehen, dh zumindest die Planetenbahnen vernünftig berechnet werden, dann sehe ich keinen Grund, die Pferde umzusatteln.

Bynaus
07.10.2012, 12:39
Ich war nun einige Tage auf Exkursion und während der langen Autofahrten hab ich u.a. über dieses Projekt ein bisschen nachdenken können. Ich weiss, dass meisten "Probleme" im Moment eher technischer Natur sind, aber ich versuche, hier eine gewisse Langzeitperspektive zu entwickeln. Für ein Programm, das die Bahnen von vorwiegend masselosen Teilchen durch das Sonnensystem simulieren kann, gibt es unterschiedliche Verwendungsarten, und ich hab grad ein paar Ideen für Projekte, die man dereinst damit realisieren könnte (die über die "TheiaSim" hinaus gehen, auf die ich zur Zeit aber noch nicht im Detail eingehen kann).

Ein möglicherweise wichtiger Punkt wäre, dass man für jedes (?) "fixe" Objekt im Sonnensystem (Planeten, Monde, grosse Asteroiden) wählen kann, ob man ihre Bahnen gemäss gravitativen Interaktionen rechnen will, oder ob man sie einfach als gravitativ auf die masselosen Teilchen agierende, aber nicht untereinander interagierende, "Zielkörper" auf festgelegten Bahnen bewegen will. Wenn die Rechenschritte sich dann zu sehr an die Umlaufzeiten annähern (Faktor 10-100?), sollte automatisch auf festgelegte Bahnen gewechselt werden. So kann man z.B. den Mond in der Simulation behalten: man lässt ihn einfach auf seiner bekannten Bahn um das Schwerezentrum des Erde-Mond-Systems kreisen.

Weiter wäre es wichtig, dass man für jedes "fixe" Objekt (gem. oben) einzeln festlegen kann, bei welcher Distanz ein masseloses Teilchen aus der Simulation entfernt werden soll. Bei der Erde ist z.B. eine Distanz von etwa 1000 km zur Oberfläche ein guter Wert, weil man da annehmen kann, dass die Bahn von masesarmen Teilchen durch die Interaktion mit der Exosphäre stark verändert wird (stärker, als man von der Gravitation alleine vermuten würde). Für einige der Projekte, die ich angedacht habe, wäre es auch wichtig, wenn man Geschwindigkeitsvektor und -richtung zum Zeitpunkt der Annäherung an die Erdoberfläche auf 1000 km (also dann, wenn das Teilchen aus der Simulation entfernt wird) speichern könnte.

Der letzte Punkt, der mir wichtig wäre, ist wohl bereits weitgehend realisiert: Ein Input-File, in dem die Positionen und anfänglichen Geschwindigkeiten (und Vektoren) aller zu simulierenden Teilchen gespeichert sind.

Bernhard
07.10.2012, 13:33
Ich weiss, dass meisten "Probleme" im Moment eher technischer Natur sind, aber ich versuche, hier eine gewisse Langzeitperspektive zu entwickeln.
Hallo Bynaus,

ich habe mittlerweile mal probiert, das SCATR-Projekt mit Eclipse (Photran 7.0 mit Juno-Release) zu übersetzen, aber das hat leider nicht funktioniert. Eclipse liefert hier einige ziemlich kryptische Fehlermeldungen in den generierten Make-Dateien, was für mich auf irgendwelche internen Probleme von Eclipse hindeutet.

Der erweiterte Bulirsh-Stoer-Algorithmus liefert mittlerweile auch über eine Million Jahre ganz interessante Ergebnisse. Bei einer Schrittweite von 10 Tagen braucht das Programm nur einige Minuten für die Berechnung.
Gruß

Bernhard
07.10.2012, 14:05
Ein möglicherweise wichtiger Punkt wäre, dass man für jedes (?) "fixe" Objekt im Sonnensystem (Planeten, Monde, grosse Asteroiden) wählen kann, ob man ihre Bahnen gemäss gravitativen Interaktionen rechnen will, oder ob man sie einfach als gravitativ auf die masselosen Teilchen agierende, aber nicht untereinander interagierende, "Zielkörper" auf festgelegten Bahnen bewegen will.
Nur so bleibt die Fehlersuche auch meiner Meinung nach im Programm übersichtlich, denn es macht bei den Formeln schon einen recht großen Unterschied, ob der Körper ein gravitatives Feld erzeugt oder nicht. Im SCATR-Code ist das (so weit ich es überblicken kann) sehr deutlich getrennt und in den BS-Code habe ich ja unter "Protest" von SRMeister ebenfalls so eingebaut.


Weiter wäre es wichtig, dass man für jedes "fixe" Objekt (gem. oben) einzeln festlegen kann, bei welcher Distanz ein masseloses Teilchen aus der Simulation entfernt werden soll. Bei der Erde ist z.B. eine Distanz von etwa 1000 km zur Oberfläche ein guter Wert, weil man da annehmen kann, dass die Bahn von masesarmen Teilchen durch die Interaktion mit der Exosphäre stark verändert wird (stärker, als man von der Gravitation alleine vermuten würde). Für einige der Projekte, die ich angedacht habe, wäre es auch wichtig, wenn man Geschwindigkeitsvektor und -richtung zum Zeitpunkt der Annäherung an die Erdoberfläche auf 1000 km (also dann, wenn das Teilchen aus der Simulation entfernt wird) speichern könnte.
In den BS-Code habe ich das so eingebaut. Man kann dort einen kritischen Radius für jeden massebehafteten Körper angeben, ab dem eine Kollision stattfindet.


Der letzte Punkt, der mir wichtig wäre, ist wohl bereits weitgehend realisiert: Ein Input-File, in dem die Positionen und anfänglichen Geschwindigkeiten (und Vektoren) aller zu simulierenden Teilchen gespeichert sind.
Bei dem SCATR-Code ist das genau so gemacht. Allerdings kann man bei den masselosen Testkörpern auch gut darauf verzichten, weil diese Startbedingungen in unserem Fall sehr stark von der Startbedingung der Erde abhängen. Die zugehörigen Listen werden deswegen bei sehr vielen Testkörpern (>100) sehr wenig aussagekräftig und wären ledigleich als Detailinformation für die Dokumentation wichtig. Ich persönlich fände hier einen Eingabeparameter für die Anzahl der masselosen Trümmerteile als ausreichend. Um zusätzlich den Streuwinkel als Eingabeparameter zu bekommmen müsste man den SCATR-Code weiter anpassen und das macht dann wieder eine Fortran-Compiler notwendig, womit wir wieder bei den technischen Problemen wären. Man könnte eventuell auch mal versuchen hier Crowdfunding anzuwenden, um an eine Intel-Lizenz zu kommen. Mit 1900 Dollar wären wir aus dem Schneider.

Bynaus
07.10.2012, 16:15
Der erweiterte Bulirsh-Stoer-Algorithmus liefert mittlerweile auch über eine Million Jahre ganz interessante Ergebnisse. Bei einer Schrittweite von 10 Tagen braucht das Programm nur einige Minuten für die Berechnung.

Das klingt doch schon ganz gut.


Nur so bleibt die Fehlersuche auch meiner Meinung nach im Programm übersichtlich, denn es macht bei den Formeln schon einen recht großen Unterschied, ob der Körper ein gravitatives Feld erzeugt oder nicht.

Nicht dass hier Missverständnisse entstehen: ein gravitatives Feld sollte der Körper natürlich immer erzeugen (dh, er sollte immer eine Wirkung auf die masselosen Teilchen und die Planeten mit langen Umlaufszeiten haben). Aber seine Bahn sollte, sobald die Schrittweite sich seiner Umlaufzeit annähert (innerhalb etwa einer Grössenordnung: also bei 10 Tagen ist es für Merkur und den Mond bereits recht gefährlich, während die Venus wohl nicht so sehr unter diesem Problem "leidet"), nicht mehr durch Integration berechnet, sondern vorgegeben werden: sonst passiert eben das, dass Merkur und Mond aus dem System geworfen werden. Wir müssen so verfahren, weil wir ja weiterhin berechnen wollen, was mit den Teilchen passiert, die in der Nähe von Merkur rumdümpeln - deshalb können wir Merkur nicht einfach aus der Simulation entfernen.


Man kann dort einen kritischen Radius für jeden massebehafteten Körper angeben, ab dem eine Kollision stattfindet.

Sehr schön. Wie gesagt wäre es für spätere Projekte gut, wenn man auch wüsste, an welcher Stelle der Planet bzw. der gewählte Abstand getroffen wird und in welche Richtung das Teilchen sich bewegt.


Allerdings kann man bei den masselosen Testkörpern auch gut darauf verzichten, weil diese Startbedingungen in unserem Fall sehr stark von der Startbedingung der Erde abhängen.

Ja, in dem Fall. Aber in anderen denkbaren Simulationen, die man mit einem solchen Programm machen könnte, sieht es vielleicht anders aus. Ausserdem: wenn man für jedes Teilchen einzeln eine Position und Geschwindigkeit festlegt, kann man dieses Teilchen später (angenommen, es zeigt ein interessantes Ergebnis) "klonen" und variieren, um eine grössere Auflösung für dieses Ergebnis zu bekommen, ohne dass man dabei Abermillionen von uninteressanten Teilchen simulieren muss. Also sagen wir, ein Teilchen trifft den Mars, dann können wir in einem anschliessenden Durchlauf dieses Teilchen sagen wir hunderttausend Mal klonen und seine Anfangsparameter nur ein klein wenig variieren, um eine bessere Beschreibung der Wahrscheinlichkeit eines Mars-Treffers zu erhalten.


Man könnte eventuell auch mal versuchen hier Crowdfunding anzuwenden, um an eine Intel-Lizenz zu kommen. Mit 1900 Dollar wären wir aus dem Schneider.

Hm, ich wäre erstaunt, wenn das Erfolg hätte... Was sind denn, gesamthaft gesehen, die Optionen und Probleme? Ich habe mittlerweile etwas die Übersicht verloren...

Bernhard
07.10.2012, 18:17
Sehr schön. Wie gesagt wäre es für spätere Projekte gut, wenn man auch wüsste, an welcher Stelle der Planet bzw. der gewählte Abstand getroffen wird und in welche Richtung das Teilchen sich bewegt.
Ich werde versuchen, da mal etwas Brauchbares in den BSSI-Code einzubauen. (BSSI = Bulirsh-Stoer-Symplectic-Integrator)


Ja, in dem Fall. Aber in anderen denkbaren Simulationen, die man mit einem solchen Programm machen könnte, sieht es vielleicht anders aus. Ausserdem: wenn man für jedes Teilchen einzeln eine Position und Geschwindigkeit festlegt, kann man dieses Teilchen später (angenommen, es zeigt ein interessantes Ergebnis) "klonen" und variieren, um eine grössere Auflösung für dieses Ergebnis zu bekommen, ohne dass man dabei Abermillionen von uninteressanten Teilchen simulieren muss. Also sagen wir, ein Teilchen trifft den Mars, dann können wir in einem anschliessenden Durchlauf dieses Teilchen sagen wir hunderttausend Mal klonen und seine Anfangsparameter nur ein klein wenig variieren, um eine bessere Beschreibung der Wahrscheinlichkeit eines Mars-Treffers zu erhalten.
OK

Was sind denn, gesamthaft gesehen, die Optionen und Probleme? Ich habe mittlerweile etwas die Übersicht verloren...
Für 1900 $ bekäme man eine uneingeschränkt nutzbare Intel-Lizenz für Privatanwender.

Bernhard
07.10.2012, 19:33
Ich werde versuchen, da mal etwas Brauchbares in den BSSI-Code einzubauen. (BSSI = Bulirsh-Stoer-Symplectic-Integrator)
Ich habe eben eine Funktion eingebaut, welche bei jedem Schritt die Mindestdistanz der Trümmer zu den Planeten ausrechnet. Bei einer Rechnung über 10k Jahre und 20 Testkörpern fliegt beispielsweise ein Testkörper mit einer Entfernung von nur 16 Jupiterradien an Jupiter vorbei. 2 der 20 Testkörper verlassen das Sonnensystem und die 5 massiven Trümmer sammeln sich im Bereich 0.4 - 1.5 AE. Ein anderer Testkörper fliegt in einem Abstand von nur 40 Sonnenradien an der Sonne vorbei usw.

Ich könnte als nächstes die Schrittweite anpassen. Wenn sich z.B. ein Testkörper einem Planeten kleiner als einem bestimmten kritischen Wert nähert, könnte ich die Schrittweite von 10 Tagen auf einen Tag senken, um die genaue Mindestdistanz zum Planeten zu bestimmen. Während dieser Phasen mit erhöhter Genauigkeit könnte der Vorbeiflug bzw. die Kollision mit erhöhter Genauigkeit analysiert und die zugehörigen Daten (Nummer des Testkörpers, Position, Geschwindigkeit, Winkel zwischen den Bahnen usw.) abgespeichert werden.

EDIT: Habe mich leider bei der Schrittweite etwas vertan, so dass die Simulationszeit nur 10k Jahre anstelle von 100k Jahre beträgt.

SRMeister
08.10.2012, 16:04
Bernhard, wie hast du das Problem mit den Testpartikeln bei dem BSSI gelöst? Soweit ich das verstehe lassen sich masselose Partikel nicht trivial in eine symplektische Form bringen, da es zu Unendlichkeiten kommt.

Bynaus, eigentlich sollte soetwas nicht passieren:

Nicht dass hier Missverständnisse entstehen: ein gravitatives Feld sollte der Körper natürlich immer erzeugen (dh, er sollte immer eine Wirkung auf die masselosen Teilchen und die Planeten mit langen Umlaufszeiten haben). Aber seine Bahn sollte, sobald die Schrittweite sich seiner Umlaufzeit annähert (innerhalb etwa einer Grössenordnung: also bei 10 Tagen ist es für Merkur und den Mond bereits recht gefährlich, während die Venus wohl nicht so sehr unter diesem Problem "leidet"), nicht mehr durch Integration berechnet, sondern vorgegeben werden: sonst passiert eben das, dass Merkur und Mond aus dem System geworfen werden
und meinst du mit "vorgegeben", dass man die Planetenbahnen durch reine Keplerbahnen abschätzt?


Für 1900 $ bekäme man eine uneingeschränkt nutzbare Intel-Lizenz für Privatanwender.
Es gäbe die Option: Der Intel compiler ist für Linux kostenlos. Alternativ wäre ich immernoch an einer kostenlosen Windows Lösung interessiert, vielleicht mit einem GNU Fortran Compiler.

Bernhard
08.10.2012, 17:42
Bernhard, wie hast du das Problem mit den Testpartikeln bei dem BSSI gelöst? Soweit ich das verstehe lassen sich masselose Partikel nicht trivial in eine symplektische Form bringen, da es zu Unendlichkeiten kommt.
Ich sehe momentan kein Problem darin die normalen Hamilton-Gleichungen zu verwenden. Die Massen der Testkörper kürzen sich dort heraus. Man kann also den gleichen Algorithmus für die Testkörper verwenden, sofern man berücksichtigt, dass diese Testkörper kein gravitatives Feld erzeugen.

BTW: Die variable Schrittweite erscheint mir momentan recht interessant. Aktuell läuft bei mir im Hintergrund eine Version, bei der jeder Planet (inklusive Sonne) eine Art kritische Hülle von 70 Jupiterradien besitzt. Sobald diese Hülle von einem Trümmerteil durchstoßen wird, wird die Schrittweite von einem Tag auf 8640s und ab 7 Jupiterradien auf 864s gesetzt. Die Ergebnisse einer Simulation über 1000 Jahre mit 120 Trümmerteilen zeigen allerdings, dass die Zeitschritte vermutlich immer noch zu groß gewählt sind.

Bynaus
09.10.2012, 08:37
@SRMeister:


Bynaus, eigentlich sollte soetwas nicht passieren:

Du meinst die Destabilisierung der Bahnen? Warum nicht? Bei so grossen Schrittweiten relativ zur Bahn können sich (würde ich jetzt vermuten) kleine Rundungsfehler zu beträchtlichen "Kräften" auswachsen, die die Bahn verändern. Dass man Bahnen mit Umlaufzeiten in der Nähe der Schrittweite von der Simulation ausklammert, ist gang und gäbe (sieht man auch in anderen Arbeiten regelmässig). Aber kann auch sein, dass ich da etwas verwechsle.


und meinst du mit "vorgegeben", dass man die Planetenbahnen durch reine Keplerbahnen abschätzt?

Ja. Das ist zumindest nicht grob falsch.

@Bernhard:


Aktuell läuft bei mir im Hintergrund eine Version, bei der jeder Planet (inklusive Sonne) eine Art kritische Hülle von 70 Jupiterradien besitzt.

Wäre es da nicht besser, die kritische Hülle von der Masse & Entfernung des Planeten abhängig zu machen? Mit anderen Worten: die Hill-Sphäre zu verwenden?

Betreffend variable Schrittweite: Wird denn da das ganze Programm auf diese langsame Schrittweite heruntergefahren? Oder stoppt man das Programm und klammert das Trümmerteil aus, übergibt es an eine andere Instanz, die mit geringer Schrittweite den Vorbeiflug berechnet, die Endposition zurückgibt und dann mit normaler Geschwindigkeit die Simulation fortsetzt?

Bernhard
09.10.2012, 12:51
Wäre es da nicht besser, die kritische Hülle von der Masse & Entfernung des Planeten abhängig zu machen? Mit anderen Worten: die Hill-Sphäre zu verwenden?
Die Hill-Sphäre ist bei masselosen Testkörpern gleich Null und damit eigentlich ungeeignet. Wichtiger erscheint mir vielmehr eine Berücksichtigung der Geschwindigkeit des Testkörpers zu sein, um lineare Extrapolationen bei der Bahn machen zu können.


Betreffend variable Schrittweite: Wird denn da das ganze Programm auf diese langsame Schrittweite heruntergefahren? Oder stoppt man das Programm und klammert das Trümmerteil aus, übergibt es an eine andere Instanz, die mit geringer Schrittweite den Vorbeiflug berechnet, die Endposition zurückgibt und dann mit normaler Geschwindigkeit die Simulation fortsetzt?
An eine isolierte Betrachtung des lokalen Zweikörperproblems habe ich auch schon gedacht. Mir ist dabei nur nicht klar in welchem Bezugssystem man diese Rechnung dann durchführt.

Bin momentan ziemlich erkältet, was zu Verzögerungen führen wird. Das Buch von Danby ist heute aber eingetroffen.

SRMeister
17.10.2012, 01:58
Also wie ich das sehe, müssten wir den SCATR so modifizieren, dass Mond und Merkur durch reine Keplerbahnen abgeschätzt werden. Dann kann man 1/10tel der Periodendauer von Venus als Schrittweite nehmen und wir sind auf der sicheren Seite.
Korrekt?

Bernhard
17.10.2012, 10:22
Also wie ich das sehe, müssten wir den SCATR so modifizieren, dass Mond und Merkur durch reine Keplerbahnen abgeschätzt werden. Dann kann man 1/10tel der Periodendauer von Venus als Schrittweite nehmen und wir sind auf der sicheren Seite.
Korrekt?
@SRMeister: Ich denke, Du legst hier den Finger genau in die wunden Stellen des SCATR-Projektes, denn meiner Meinung nach ist das SCATR-Projekt viel zu schlecht dokumentiert, um damit sinnvoll arbeiten zu können. Sogar mit dem Buch von Danby wird praktisch überhaupt nicht klar, welche Ideen die Entwickler genau verfolgt haben, weswegen ich dieses Projekt vorerst nicht weiter unterstützen will. Auch wenn es mir selber leid tut, aber ohne weitere Dokumentation ist das Thema für mich erst mal erledigt...