Sonnensystem: Position der Planeten?

Runzelrübe

Registriertes Mitglied
Ich habe mich eben mal ein wenig mit der 3D Darstellung dieser Szenerie beschäftigt und hätte da mal einen Wunsch. Könntet ihr eure Ergebnisse jeder Berechnung eventuell in eine CSV Datei gießen (id, x, y, z reichen aus je Objekt), so dass ich eure Ergebnisse (und auch Verbesserung/Verschlechterung selbiger) mal visuell darstellen kann? Danke!

Da ich ein Pathtracing eingebaut habe, würde sogar einer visuellen Gegenüberstellung zur Horizons Datanbank nichts im Wege stehen. Manchmal erkennt das menschliche Auge Fehler anhand von Bewegungen bei übertrieben hohen Geschwindigkeiten. Wenn die Daten also nicht erst live berechnet werden müssen, kann die Darstellung auch viel zügiger ablaufen. Systematisch größer werdende Abweichungen wären nicht mehr nur eine Zahl.

-----

Während der Betrachtung der Pfade ist mir eine Idee gekommen. Mal angenommen, wir würden mit einem Teleskop Langzeitbelichtungen derjenigen Sternsysteme durchführen, in denen wir bereits bestätigte Transits haben. Dann könnten wir mit diesen Aufnahmen einen weit vom Zentralgestirn entfernten Planeten immer dann erkennen, wenn dieser aus unserer Sicht seine horizontale Richtung wechselt. Die Wahrscheinlichkeit sollte ja eigentlich hoch sein, dass die Ebene durch den bestätigten Transit auch für andere Planeten gilt.

Das sähe dann ungefähr so aus:

S... Stern
P... Planetenumkehr


PPP.............................................S.........................................PPP

Bei PPP sollte während der Langzeitbelichtung ein durchaus markanter Punkt zu sehen sein. Hintergründe könnte man ja anhand von Kurzzeitbelichtungen subtrahieren.

Wäre das machbar? Sinnvoll? Wird das eventuell schon gemacht?
 

sirius3100

Registriertes Mitglied
Ich verstehe nicht genau was gemeint ist: Wer soll seine Richtung ändern? Der Transitplanet, oder der Stern, oder der Planet selbst?
Und warum zwangsläufig horizontal und nicht auch vertikal (radial schließe ich mal aus)?
Am ehesten hört sich das ganze ja nach Astrometrie an, aber dafür brauch ich natürlich keinen Transitplaneten in dem System. Das "funktioniert" auch ohne. Astrometrie hat halt den Nachteil das sie explizit von unserem Abstand zum beobachteten Stern abhängt.
Oder soll der Planet an den markierten Punkten der Bahn direkt beobachtet werden?
Direkte Beobachtung ist natürlich umso leichter je heller das Objekt leuchtet und je weiter es vom Zentralgestirn weg ist. Funktioniert aber deshalb eigentlich nur für schwere Gasriesen (andere Objekte sind einfach nicht hell genug um sie direkt zu sehen), und auch bei denen muss man tricksen.
 
Zuletzt bearbeitet:

Bewegt

Registriertes Mitglied
S... Stern
P... Planetenumkehr


PPP.............................................S.........................................PPP

Bei PPP sollte während der Langzeitbelichtung ein durchaus markanter Punkt zu sehen sein. Hintergründe könnte man ja anhand von Kurzzeitbelichtungen subtrahieren.

Wäre das machbar? Sinnvoll? Wird das eventuell schon gemacht?
geniale Idee
die funktioniert um so besser, je mehr sich der Beobachter auf Rotationsebene befindet
 

Runzelrübe

Registriertes Mitglied
Wer soll seine Richtung ändern? Der Transitplanet, oder der Stern, oder der Planet selbst?

Ein (noch nicht entdeckter) Planet innerhalb eines Sternsystems, für das bereits ein Transitplanet bekannt ist.

Dann könnten wir mit diesen Aufnahmen einen weit vom Zentralgestirn entfernten Planeten immer dann erkennen, wenn dieser aus unserer Sicht seine horizontale Richtung wechselt.

Und warum zwangsläufig horizontal und nicht auch vertikal (radial schließe ich mal aus)?

Das habe ich sehr unglücklich ausgedrückt. In Systemen, in denen es aus unserer Sicht heraus Transitplaneten gibt (wir haben nachgewiesen, wie sie vor ihrem Zentralgestirn vorüberziehen), schauen wir aller Wahrscheinlichkeit nach nahezu senkrecht auf den Normalenvektor der Rotationsebene. :)

Dadurch dürften wir eine Planetenbahn als Linie betrachten (und bitte keine Spitzfindigkeiten jetzt!), auf denen der Planet von unserem Standpunkt aus oszilliert. Den Rest der Erklärung kannst Du jetzt gern nochmal durchlesen.
 

sirius3100

Registriertes Mitglied
Hört sich also doch nach Direct Imaging an. Bisher konnten halt nur Exoplaneten mit einem Mindestabstand von ca 8AU (Beta Pictoris b) von ihrem Zentralgestirn direkt abgebildet werden.

Dein Vorschlag läuft also darauf hinaus bei der Auswahl zu beobachtender Sterne solche vorzuziehen bei denen bereits ein Transitplanet entdeckt wurde, da sich damit die Chance erhöht dass ein Planet von uns aus betrachtet weit genug vom Zentralgestirn weg ist um ihn direkt beobachten zu können.

Habe ich das jetzt richtig verstanden? Oder ist ganz was anderes gemeint?
 
Zuletzt bearbeitet:

UMa

Registriertes Mitglied
Hallo Stefan, hallo Bernhard,

ich habe mal mein altes Programm ausgegraben.
Ich habe vom Startpunkt de405iom.pdf Datei genau 360 Tage sowohl relativistisch als auch mit Newton vorrausberechnet. Der Unterschied im Ort der Erde ist 49791 m. Mit Horizons habe ich noch nicht vergleichen können.
Des weiteren könnte es Unterschiede im Zeitablauf geben. Horizons nimmt vermutlich Erdzeit, während ich in der Zeit eines weit entfernten Beobachters rechne (denke ich). Vielleicht erklärt das zusammen die ca. 60km Abweichung.

Nur relativistische Massen oder so einzusetzen verbessert das nicht, ihr müsstet eine bessere Approximation an die ART rechen als Newton z.B. die Einstein-Infeld-Hoffman Gleichungen.
http://en.wikipedia.org/wiki/Einstein–Infeld–Hoffmann_equations
Ergänzung: Die dort verlinkte erste Referenz
Standish, Williams. Orbital Ephemerides of the Sun, Moon, and Planets
düfte eigentlich alles weitere erklären.

Die GM Werte aus der Datei sind vielleicht für die Planeten mit Monden, die aus Horizons ohne Monde?!

Ich habe feste Zeitschritte von 1, 2, 3, 4, 5 und 6 Tagen verwendet. Die Unterschiede sind kleiner als zwischen relativistisch als auch mit Newton, zwischen 1d und 2d weniger als 1 m in der Erdposition.
Die Rechenzeiten für 1000 Jahre (genauer 360000 Tage) liegen bei wenigen Sekunden bis einer Minute je nach Version und Zeitschrittweite. Der Unterschied zwischen 1d und 2d Zeitschrittweite nimmt dabei für der Erde auf über 6km nach ca. 1000Jahren zu.

Ich würde nicht erwarten das RK4 + Extrapolation besser ist als Bulirsch Stoer, könnt ihr aber natürlich selbst ausprobieren.

Grüße

UMa
 
Zuletzt bearbeitet:

Bernhard

Registriertes Mitglied
Ich würde nicht erwarten das RK4 + Extrapolation besser ist als Bulirsch Stoer, könnt ihr aber natürlich selbst ausprobieren.
Hallo UMa, Hallo Stefan,

ich habe jetzt doch mal das Bulirsch-Stoer-Verfahren gemäß UMas Anleitung eingebaut. Man bekommt bei damit bei achtfacher Schrittweite und gleicher Rechenzeit praktisch die gleichen Ergebnisse, wie mit RK4.

BTW: In den numerical recipes wird die Extrapolation etwas anders mit Hilfe von Polynomen, bzw. rationalen Funktionen durchgeführt. UMas Anleitung ist da wesentlich einfacher zu programmieren.

Die EIH-Gleichungen sind noch in Bearbeitung.
MfG
 

Bernhard

Registriertes Mitglied
Nur relativistische Massen oder so einzusetzen verbessert das nicht, ihr müsstet eine bessere Approximation an die ART rechen als Newton z.B. die Einstein-Infeld-Hoffman Gleichungen.
http://en.wikipedia.org/wiki/Einstein–Infeld–Hoffmann_equations
Hallo UMa,

der Tip ist wirklich Spitze. Mit diesen Gleichungen kann man die Abweichung von etwas mehr als 60 km auf etwa 3 km reduzieren und dieser Wert liegt dann größenordungsmäßig in der Größe der Zeitdilatation zwischen Erdzeit und der Zeit im Schwerpunktsystem :cool: .

Ergänzung: Die dort verlinkte erste Referenz
Standish, Williams. Orbital Ephemerides of the Sun, Moon, and Planets
düfte eigentlich alles weitere erklären.
Sehr interessant, dass dort die Bewegung der Sonne aus der numerischen Integrationsroutine herausgenommen wird.

Die GM Werte aus der Datei sind vielleicht für die Planeten mit Monden, die aus Horizons ohne Monde?!
Genau so ist es scheinbar.
Beste Grüße
 

Runzelrübe

Registriertes Mitglied
@SRMeister:

Habe mir eben Deinen Code angeschaut und noch ein paar geschwindigkeitsfördernde Optimierungsvorschläge, die Dir zumindest bei Deiner Umsetzung der relativistischen Berechnung helfen werden.

Aus:

Code:
PSystem::PSystem(double Timestep)
{
	c = 299792458.0;
	Tconst = Timestep;
	sun = 0;
}
....
void PSystem::AddObject(char *Name, double Mass, double x, double y, double z, double vx, double vy, double vz)
{
	...

mach:

Code:
PSystem::PSystem(double Timestep)
{
	c = 299792458.0;
	[B]invC2 = 1.0/(c*c);[/B]
	Tconst = Timestep;
	[B]tConstHalf = 0.5*Tconst;[/B]
	sun = 0;
}
....
void PSystem::AddObject(char *Name, double Mass, double x, double y, double z, double vx, double vy, double vz)
{
	[B]numObjects +=1;[/B]

Dann kannst Du aus:

Code:
/c/c

das hier machen:

Code:
*invC2

und aus:

Code:
Evaluate(Tconst * 0.5,

das hier:

Code:
Evaluate(tConstHalf ,

Einer der wirklichen Geschwindigkeitsbooster, den ich Dir immernoch ans Herz lege, ist die Halbierung Der benötigten Schleifen je Zeitschritt:

Aus:

Code:
for(source = sun; source != 0; source = source->next)
{
	source->der[0].dv = 0;
	for(target = sun->next; target != 0; target = target->next)
	{
		if(source != target)
		{

mach:

Code:
int i,j;
source = sun;

for(i=[B]0[/B]; i<[B]numObjects[/B]; i++)
{
	source->der[0].dv = 0;
	source = source ->next;
}

source = sun;
target = sun->next;

int outerLoopCount = numObjects-1;
for(i=[B]0[/B]; i<[B]outerLoopCount[/B]; i++)
{
	target = source ->next;
	for(j=[B]i+1[/B]; j<[B]numObjects[/B]; j++)
	{
		... Anteil wie bisher auf source anwenden
		... umgekehrter Anteil (andere Masse und Richtung) auf target anwenden
		target = target->next;
	}
	source = source->next;
}

Du hast nur die halbe Anzahl Durchläufe, sparst Dir die n^2 Prüfungen, ob source != target und musst auch nicht den Vektor, den quadratischen Abstand, den Abstand (Wurzelberechnung!) und den inversen Abstand (Division!) doppelt berechnen. Alles, was dazu kommt sind vier Multiplikationen und drei Vorzeichenwechsel (dazu wird ein Bit umgesetzt, ist ein Witz).
 

Bernhard

Registriertes Mitglied
Du hast nur die halbe Anzahl Durchläufe
Hallo Rübe,

ich habe Stefan gerade die neueste Version mit sämtlichen Tips von UMa geschickt. Die Variablen invC2 und tHalfConst habe ich gleich miteingebaut. Die Anzahl der Durchläufe läßt sich IMO nicht reduzieren, weil sämtliche Körper in der Simulation unterschiedliche Massen haben. Man hat damit n^2 verschiedene Beschleunigungen, die immer wieder komplett neu berechnet werden müssen. Eine Umkehrung der Richtung der Beschleunigung alleine reicht dabei nicht aus.

Eventuell hat UMa seine Software noch besser als die neue Version 08 implementiert, weil die scheinbar um einiges schneller läuft.

Nach Einfügen der Körper Pluto, Ceres, Pallas, Vesta und Juno benötigt die Bulirsch-Stoer-Integration der EIH-Gleichungen bei einer Schrittweite von 8640 Sekunden und einer Integration über knapp fünf Jahre, d.h. bis zum 30.03.2016 00:00:00, etwa zweieinhalb Minuten Rechenzeit auf einem aktuellen Standard-PC. Der absolute Fehler in der Position der Erde beträgt im Vergleich zu den Werten aus HORIZONS etwa 53 km.
MfG
 

Runzelrübe

Registriertes Mitglied
Die Anzahl der Durchläufe läßt sich IMO nicht reduzieren, weil sämtliche Körper in der Simulation unterschiedliche Massen haben. Man hat damit n^2 verschiedene Beschleunigungen, die immer wieder komplett neu berechnet werden müssen. Eine Umkehrung der Richtung der Beschleunigung alleine reicht dabei nicht aus.

Hi Bernhard,

geht wohl. Schau Dir mal folgenden funktionierenden C#-Code an (Real ist hierbei der nach Shewchuk implementierte Datentyp - dynamische Präzisionserhaltung garantiert):

Stell doch einfach mal den Code zweier Körper gegenüber. Dann siehst Du, dass Du nur (G)Masse austauschen und den Vektor umdrehen musst.. et voila.

Code:
public void CalculateAccelleration()
{

	Real xx,yy,zz,sqrDistance,distance,invDistance3;
	Real multiplierIII,multiplierJJJ;
	IMovableMass objIII, objJJJ;

	var objects = this.Formation.AsEnumerable().ToArray();
	var objectsCount = objects.Length;
	for (int i = 0; i < objectsCount; i++)
	{
		objects[i].AX = 0;
		objects[i].AY = 0;
		objects[i].AZ = 0;
	}

	var objectsCount1Less = objectsCount - 1;
	for (int i = 0; i < objectsCount1Less; i++)
	{
		obj[B][COLOR="red"]III[/COLOR][/B] = objects[i];
		for (int j = i + 1; j < objectsCount; j++)
		{
			obj[B][COLOR="green"]JJJ[/COLOR][/B] = objects[j];
			xx = obj[B][COLOR="green"]JJJ[/COLOR][/B].X - obj[B][COLOR="red"]III[/COLOR][/B].X;
			yy = obj[B][COLOR="green"]JJJ[/COLOR][/B].Y - obj[B][COLOR="red"]III[/COLOR][/B].Y;
			zz = obj[B][COLOR="green"]JJJ[/COLOR][/B].Z - obj[B][COLOR="red"]III[/COLOR][/B].Z;
			sqrDistance = xx * xx + yy * yy + zz * zz;
			distance = Real.Sqrt(sqrDistance);
			invDistance3 = 1.0 / (sqrDistance * distance);
			multiplier[B][COLOR="red"]III[/COLOR][/B] = obj[B][COLOR="green"]JJJ[/COLOR][/B].Mass * invDistance3;
			multiplier[COLOR="green"][B]JJJ[/B][/COLOR] = obj[COLOR="red"][B]III[/B][/COLOR].Mass * invDistance3;
			obj[B][COLOR="red"]III[/COLOR][/B].AX [B][COLOR="red"]+[/COLOR][/B]= multiplierI * xx;
			obj[B][COLOR="red"]III[/COLOR][/B].AY [B][COLOR="red"]+[/COLOR][/B]= multiplierI * yy;
			obj[B][COLOR="red"]III[/COLOR][/B].AZ [B][COLOR="red"]+[/COLOR][/B]= multiplierI * zz;
			obj[B][COLOR="green"]JJJ[/COLOR][/B].AX [B][COLOR="green"]-[/COLOR][/B]= multiplierJ * xx;
			obj[B][COLOR="green"]JJJ[/COLOR][/B].AY [B][COLOR="green"]-[/COLOR][/B]= multiplierJ * yy;
			obj[B][COLOR="green"]JJJ[/COLOR][/B].AZ [B][COLOR="green"]-[/COLOR][/B]= multiplierJ * zz;
		}
	}
}
 
Zuletzt bearbeitet:

Runzelrübe

Registriertes Mitglied
Ich werde mal versuchen diesen Trick auch auf die Integration der EIH-Gleichung anzuwenden. Die Integration dort läuft aktuell schon etwas langsam.

Hi Bernhard,

viel Erfolg und auch Spaß dabei! :)

Ich hab da noch ein paar Vorschläge. Lass Dich aber bitte nicht gleich davon verwirren. Für den Anfang reicht auch erstmal die Reduktion der Schleifendurchläufe insgesamt, ganz ohne zusätzliche Vorbereitungen.

Ich würde Dir sonst nämlich gern ein mehrstufiges Herangehen empfehlen. Erst die Komponenten errechnen, dann zusammenfügen.

Häufig verwendete Faktoren würde ich mir in lokale Variablen packen.

invC2
invSevenHalfC2

Nach dem Ansatz untere Dreiecksmatrix kannst Du zuerst einmal die benötigten Einzelkomponenten ausrechnen. Ich würde vorschlagen, ein struct dafür einzubauen und dieses je gerichteter Körperkombination in ein lokales Array zu speichern. Wichtig ist, dass diess Array beim Zugriff nach AB und BA unterscheiden kann. Die Indizierung ist Dir selbst überlassen.

Die sofort berechenbaren Einzelkomponenten sind:

nAB, nBA
vA2, vB2
vAB, vBA
rAB, r2AB, invRAB, invR2AB

Ich würde noch zwei weitere Variablen einführen:

pAB = xA-xB (x soll hier jeweils der Ortsvektor sein)
invR3AB = 1.0/(r2AB*rAB)

Somit können die Stellen, an denen "nBA/r2AB" steht, in "pBA*invR3AB" umgewandelt werden und die Stellen, an denen "xA-xB" steht, in "pAB".

Jetzt könnte man hingehen und auch noch die Terme in Variablen packen. Wenn Du das tust, dann lass Dich noch nicht von der Formel irritieren, die nach A, B und C unterscheidet. Du benötigst in den Vorbereitungen die Kombinationen AB und BA. Mehr nicht. Denn auch AC und BC werden in der Schleife irgendwann durchgenommen und aufgesammelt. Erst durch den Zugriff beim Zusammenrechnen greifst Du auf die Logik zurück, jeweils das Körperpaar A, B mit C bis N zu verbinden.

Sieht nach verdammt viel Aufwand aus, aber Deine Zugriffe beim Zusammenrechnen können das um Längen wieder reinholen.

Beim Kombinieren haben wir ja drei ineinander verschachtelte Schleifen.

for i in [0: (N-2)]
for j in [(i+1): (N-1)]
for k in [(j+1): N]

Wenn wir an alle voneinander abhängigen Terme per Zugriff rankommen, schicken wir die Berechnungsliste ja nur noch in die Pipeline. Somit wird dieser Punkt auf einmal der schnellste. Falls Dir das jetzt nicht sofort etwas sagt, hier nochmal ein kleines Beispiel:

Die drei Zeilen:

a = 3*x+1
b = 4*y+1
c = 5*z+1

können gleichzeitig ausgerechnet werden, sogar parallel (wenn man es integriert).

Jedoch:

a = 3*x+1
b = 4*y+a <- Zack, Unterbrechnung, wir benötigen a, um weiterzumachen.


So, da bin ich jetzt tief genug hinabgetaucht. Wie Du es im Endeffekt schneller machst, welche Variablen Du nimmst und wieviel Aufwand Du reinstecken möchtest für das schnellste Ergebnis, das ist Dir überlassen. Und Schritt für Schritt ist da immernoch die beste Lösung. :)
 

Bernhard

Registriertes Mitglied
Hi Rübe,

die nichtrelativistische Rechnung wird durch die Berücksichtigung von \vec{r}_ab = -\vec{r}_ba schon mal um einen Faktor 2 (!) schneller und ist damit ein wirklich wichtiger Hinweis. Die relativistische Rechnung profitiert davon dann auch, weil die Beschleunigungen auf der rechten Seite der EIH-Gleichungen einfach aus der nichtrelativistischen Rechnung genommen werden.

Tricks, wie die Einführung von Konstanten wie SevenHalfC2 bringen ebenfalls Rechenzeit und machen den Code auch lesbarer. Bleibt dann noch die etwas schwierigere Strukturierung der Additionen, damit man auch hier die gleiche Arbeit nicht mehrfach macht.

Man kann mit der aktuellen und noch nicht öffentlich zugänglichen Version mit Runge-Kutta locker mal über 10 Jahre rechnen und das bei einer Schrittweite von 100s. Speziell bei dieser Rechnung waren die Fehler im Vergleich zur Rechnung über ein Jahr wieder deutlich kleiner (?). Auf jeden Fall macht es ziemlich Spaß mit dem Programm und HORIZONS etwas herumzuspielen.
Beste Grüße
 
Zuletzt bearbeitet:

Bernhard

Registriertes Mitglied
Hallo SR und R.Rübe,

Version 09 ist schnell genug, um 3650 (knapp zehn Jahre) Tage mit einer Schrittweite von 8640 Sekunden zu rechnen. Die nichtrelativistische Rechnung nach Bulirsch-Stoer wird bei mir in ca. 12s gerechnet. Die relativistische Rechnung braucht etwa 94s, ist aber leider von den Ergebnissen her schlechter, als die nichtrelativistische Rechnung.

Wer das Programm testen will: In der Datei PSystem.h kann in Zeile 7 ein Define gesetzt werden, das über einige Direktiven (#ifdef, ...) den Code steuert. Mit diesem Define wird relativistisch gerechnet. Wird die Zeile dagegen auskommentiert (und das Programm neu) übersetzt, wird nichtrelativistisch gerechnet. Im Code an SR ist die Zeile auskommentiert.
MfG
 

Bernhard

Registriertes Mitglied
ist aber leider von den Ergebnissen her schlechter, als die nichtrelativistische Rechnung.
Bei einer Rechnung über 20 Jahre ist die relativistische Rechnung wieder deutlich besser.

Es zeigt sich dabei auch, dass die Windows-Kalenderfunktion noch falsch angesteuert wird. Erstens sollten bei den verwendeten Parametern (Vielfache von 86400) keine Uhrzeiten, wie 23:00 entstehen und zweitens kann man das korrekte Datum mittels des julianischen Datums (s. HORIZONS-Ephemeriden) nachprüfen. Bei der Rechnung über zwanzig Jahre beträgt dabei der Fehler vom ersten Hinsehen einen ganzen Tag.

Vorsicht also mit den Kalenderangaben beim Ausprobieren.
MfG
 

SRMeister

Registriertes Mitglied
So ich habs endlich geschafft, eine Windows Version fertigzustellen.
Die zeitkritischen Funktionen (ich sagmal alles was innerhalb von DoSteps() passiert) werden in einer C++ DLL ausgeführt und sind somit so schnell wie die C++ .EXE datei. Die Windowsoberfläche ist C# basiert. Der Wrapperteill (Die Schnittstelle zwischen C++ und .NET) ist in C++.NET geschrieben. Es sind somit 3 Programmiersprachen enthalten.
Das ganze unterteilt sich dann logischerweise in mehrere Projekte:
PSystemDLL <-- enthält alle C++ bestandteile, basierend auf der Version 09 von Bernhard. Es wurde lediglich eine CountPlanets() Fkt. hinzugefügt da sie notwendig wurde.
PSystemWrapper Ist in C++.NET geschrieben und importiert die PSystenDLL in .NET
PSystemWrapperTest Ist eigentlich nicht nötig, aber war zum Testen gedacht. Es ist eine Konsolenanwendung, geschrieben in C# und stellt damit das DotNetGegenstück zur alten PSystemC.exe dar
NewSimulator ist die Windows Anwendung. Sie ist eine .exe und benötigt die PsystemDLL.dll und PSystemWrapper.dll im gleichen Ordner zum starten. Sie besteht wiederum aus einem Hauptprogramm und einem Benutzersteuerelemt "SimViewControl" was die Darstellung übernimmt und somit auch einfach in andere Programme eingebunden werden kann.

VisualStudio 2010 ist für den Quellcode benötigt.

Hier gehts zur 32-bit Binary
Hier zum Quellcode

Have Fun!
 
Zuletzt bearbeitet:
Oben