Wir basteln uns eine Simulation

So, da der Blog ja auch die Leute zu tätiger Arbeit animieren soll (auch wenn ich befürchte die meisten sind nur passive Leser) heute ein Thema mit dem man selbst leicht spielen kann:

Wir Simulaieren die Gravitation im Allgemeinen, wenn man es speziell macht, dann kann man Swing-Bys simulieren, das Sonnensystem oder die Entstehung des Sonnensystems oder Sternenhaufen bzw. Galaxien, je nach Lust und Geschwindigkeit des Rechners. Das schöne an der Gravitation ist, das sie so einfach ist:

Die Kraft die auf einen Körper wirkt ist einfach definiert nach:

A = G * Masse / Distanz²

• G ist eine Konstante (6,673×10^-11 m³/kg*s²)
• Masse die Masse des Himmelskörpers in Kilogramm dessen Gravitation wirkt (Erde: 5,976×10²⁴ kg, Sonne 1,989×10³⁰ kg.)
• Distanz die Distanz zwischen dem Körper in Meter (wichtig: Nicht Kilometer!)

Will man die Umlaufbahn der Erde um die Sonne berechnen, so nimmt man die Masse der Sonne, will man die des Mondes um die Erde berechnen, so die der Erde. Kleiner Fallstrick: Die Gravitation ist eine anziehende Kraft, man muss also noch ein Minus als Vorzeichen vor G setzen. Da sich immer zwei Himmelskörper gegenseitig anziehen wird auch der zweite Beschleunigt, man muss nur die Masse austauschen, also bei Erde-Mond-System einmal die Beschleunigung des Mondes durch die Erde und einmal die Beschleunigung der Erde durch den Mond berechnen.

Allgemein gilt für ein System aus N Körpern:

Äußere Schleife I von 1 bis N

Schleifenbeginn

Innere Schleife J von 1 bis N

Schleifenbeginn

Wenn i<>J

dannja

Berechne Gravitationskraft auf Körper[i] mit Masse von Körper[j]

Summiere Gravitationskraft über alle j

dannende

Schleifenende

Berechne neue Position von Körper[i] anhand der Summe der Gravitationskraft

Schleifenende

Es kommt als neues Element nun noch die Positionsberechnung hinzu. Die ist physikalisch so geregelt:

Jeder Körper hat eine Geschwindigkeit v die man in drei Raumvektoren aufteilen kann, genauso wie er eine Position im Raum (Koordinaten der x,y,z-Achse) einnimmt.

Die Distanz ist definiert als die Quadratwurzel aus den Quadraten der Positionsunterschiede zwischen den beiden Himmelskörpern

Distanz = Quadratwurzel aus (Quadrat(Differenz X-koordinaten) + Quadrat(Differenz Y-koordinaten) + Quadrat(Differenz Z-koordinaten))

Nach den Gesetzen der Physik kann man die Beschleunigung A ebenso in Teilbeschleunigungen entlang der Raumachsen teilen indem man folgenden Ansatz macht:

Ax = A * Differenz der X-Koordinaten / Distanz

analog das ganze auch für Ay und Az

Und nach den Gesetzen der Physik erhält man die neue Geschwindigkeit in jeder Raumkoordinate, indem man die Beschleunigung addiert. Wenn man um die Simulation zu beschleunigen nicht 1 Sekunde als Basis nimmt, muss man mit dem Intervall noch multiplizieren:

Vx = Vx + Ax * Intervall

Analog für die anderen zwei Raumachsen

Und die neue Position Px erhält man durch Addition der Geschwindigkeit, auch hier multipliziert mit dem Intervall (für physikalisch fortgebildete: Die Geschwindigkeit ist das Integral der Beschleunigung über die Zeit und der Weg das integral der Geschwindigkeit über die Zeit)

Px = Px + Vx * Intervall

Es gibt nun nur noch eines zu beachten: Wir müssen mit zwei Positionen arbeiten: Einer Position die der Himmelskörper beim Beginn einer jeder Berechnung über alle Körper hat und eine neue die er als Folge der Kräfte einnimmt. Da sich die Distanz durch die Berechnung ändert, muss man überall die gleiche Ausgangsbasis haben, würde man nur eine Position nehmen so würde sie sich als Folge der Berechnungen verändern. Wir müssen also Am Ende jedes Durchgangs nochmals die Positionen umkopieren.

Das war es schon: In Pascal mit einem Rekord pro Himmelskörper (Struct für C/Java Programmierer) sieht dann der Code so aus:

type
  Himmelskoerper=record
    Ax,Ay,Az : Double;
    Vx,Vy,Vz : Double;
    Px,Py,Pz : Double;
    Px2,Py2,Pz2 : Double;
    Radius : Double;
    Masse : Double;
  end;

var
  SSystem : arrayof Himmelskoerper;

const
  Grav=-6.67384E-11;

procedure Step(const T : Double);

var
  I,J : Nativeint;
  Dx,Dy,Dz : Double;
  D,D2 : Double;
  G : Double;

begin
for I:=low(SSystem) to high(SSystem) do
begin
    SSystem[I].Ax:=0;
    SSystem[I].Ay:=0;
    SSystem[I].Az:=0;
    for J:=low(SSystem) to high(SSystem) do
begin
if I<>J then
begin
        Dx:=SSystem[I].Px-SSystem[J].Px;
        Dy:=SSystem[I].Py-SSystem[J].Py;
        Dz:=SSystem[I].Pz-SSystem[J].Pz;
        D2:=Sqr(Dx)+Sqr(Dy)+Sqr(Dz);
        D:=Sqrt(D2);
        G:=(SSystem[J].Masse*Grav)/D2;
        SSystem[I].Ax:=SSystem[I].Ax+(G*Dx/D);
        SSystem[I].Ay:=SSystem[I].Ay+(G*Dy/D);
        SSystem[I].Az:=SSystem[I].Az+(G*Dz/D);
      end;
    end;
    SSystem[I].Vx:=SSystem[I].Vx+(SSystem[I].Ax*T);
    SSystem[I].Vy:=SSystem[I].Vy+(SSystem[I].Ay*T);
    SSystem[I].Vz:=SSystem[I].Vz+(SSystem[I].Az*T);
    SSystem[I].Px2:=SSystem[I].Px+(SSystem[I].Vx*T);
    SSystem[I].Py2:=SSystem[I].Py+(SSystem[I].Vy*T);
    SSystem[I].Pz2:=SSystem[I].Pz+(SSystem[I].Vz*T);
  end;
  for I:=low(SSystem) to high(SSystem) do
begin
    SSystem[I].Px:=SSystem[I].Px2;
    SSystem[I].Py:=SSystem[I].Py2;
    SSystem[I].Pz:=SSystem[I].Pz2;
  end;
end;

Radius wird hier nicht benutzt kann aber genutzt werden, um eine Kollision zu erkennen. Ein System beliebiger Art kann nun einfach simuliert werden indem man dem dynamischen Array SSystem Himmelskörper mit den entsprechenden physikalischen Parametern zuweist und in einer Schleife immer wieder die neue Position vermisst und sie plottet. Ich habe hier zwei Grafiken mit einer Berechnung pro Tag und einer Laufzeit von 1 Million Tagen, also rund 3000 Jahren wiedergegeben.

• Der Blaue Kreis ist die Erde
• Der Rote Kreis ist Jupiter. Der schwarz eingezeichnete Körper ist ein Asteroid bei Simulationsbeginn in 380 Millionen km Entfernung
• Fall 1: Er hat eine Geschwindigkeit die 1000 m/s kleiner als die Kreisbahngeschwindigkeit ist (Das Perihel liegt näher als 380 Millionen km)
• Fall 2:Er hat eine Geschwindigkeit die 3000 m/s höher als die Kreisbahngeschwindigkeit ist (Das Perihel liegt höher als 380 Millionen km)

Man sieht in beiden Fällen wie die Gravitation, vor allem von Jupiter, die Bahn verzerrt. Der Effekt wird immer größer je näher der Körper Jupiter bekommt. Da er im zweiten Fall sehr bald die jupiterbahn kreuzt ist es nru eine Frage der Zeit bis er auf Jupiter stürzt oder aus dem Sonnensystem herauskatapultiert wird (3000 Jahre sind in geologischen Maßstäben keine lange Zeit). Auch die Position der Erde schwankt wie man an den etwas dickeren Linien sieht, hier sind es aber nur wenige Millionen Kilometer (1 Pixel = 1 Million km). Wenn jemand so etwas plottet sollte er darauf achten vor dem Plotten die Position der Sonne abzuziehen, da diese auch angezogen wird und sich bewegen wird und damit das ganze System.

Hier das ganze als Delphi Projekt, es sollte mit kleinen Anpassungen auch bei Lazarus funktionieren. Der EnhancedEdit fehlt im Standard Delphi, das ist ein Edit Button für Zahlen. der eingetippte Wert steht als Integer/Double in Value/Valueint. Man kann das ganze aber auch leicht mit einem normalen Edit machen und dann eben den String in eine Zahl konvertieren. In diesem Beispiel sind die Himmelskörper hart im Quelltext codiert. Wer will kann es noch um eine Maske für die Eingabe derer Werte ergänzen.

Dann noch eine kleine Nachricht: Die Neuauflage der europäischen Trägerraketen ist nun auch als ebook erschienen. Anders als die alte Auflage nun zu besseren Konditionen, sprich günstigerem Preis. Für 4 Wochen wird sie nochmals günstiger, denn die ersten 4 Wochen gibt es sie für 13,99 Euro anstatt 17,99 Euro. (Die Printausgabe kostet 34,99, man spart also eine Menge und das ganze auch noch in Farbe …). Meine gedruckte Ausgabe habe ich inzwischen erhalten. Ich hatte zuerst die Befürchtung Michel Vans Covergrafiken sind verpixelt, aber das sind sie nicht. Sie sind nur etwas klein, aber das liegt eben an dem Format das dreimal breiter als hoch ist. Der Druck ist aber okay, nicht zu dunkel wie ich auch schon Fälle hatte.

Inzwischen ist es mir auch gelungen für Createspace ein sauberes PDF zu erzeugen, plötzlich ging es wieder, fragt nicht wieso. Wenn ich noch das Cover von Michel habe, dann könnt ihr die Gesamtausgabe mit 540 Seiten im A-4 Format für einen Hammerpreis erwerben. Leider nur über Amazon.

Inzwischen bin ich schon an dem nächsten Buch und arbeite mich durch die Neuauflage der Internationalen Trägerraketen. Das wird aber nur eine moderate Ergänzung sein, also neue Träger, aktualisierte Startzahlen und Texte und soweit möglich auch noch die Triebwerksdaten / Startprofile. Es wird aber sicher nicht das Mammutwerk wie die US-Trägerraketen werden.