Das Mehrkörperproblem in der Astronomie/ Zusammenfassung wichtiger C-Prozeduren/ Beschleunigung, Ruck und Energie

Eine Mehrkörpersimulation erfordert in jedem Schritt die Berechnung der auf die einzelnen Mitglieder ausgeübten Beschleunigungen auf Grundlage des Newtonschen Gravitationsgesetzes. Werden die nach einem solchen Schritt zu erwartenden neuen Positionen und Geschwindigkeit mittels Hermite-Polynome bestimmt, müssen auch die auf die Körper einwirkenden Rucks ermittelt werden. Für die Berechnung der einzelnen Bahnen wird die Gesamtenergie an für sich nicht benötigt, doch stellt diese aufgrund des Energieerhaltungssatzes ein gutes Kriterium dar, die Stabilität einer Simulation zu überprüfen.

Im Folgenden werden die für die Berechnung der auf individuelle Körper einwirkenden Beschleunigungen und Rucks eingeführten Prozeduren nochmals zusammengestellt, zusammen mit derjenigen für die Bestimmung der Gesamtenergie eines Mehrkörpersystems. In allen Fällen wird davon ausgegangen, dass das Gravitationsgesetz nach der einfachen Methode von Aarseth mittels eines endlichen Plummerradius geglättet wird. Fast alle Variablen sind erneut vom Typ Double, Zähler wiederum Unsigned Integer.

Die hier nochmals skizzierten Prozeduren für die Beschleunigung und den Ruck eines Körpers gelten für Simulationen, bei denen die Beiträge aller übrigen Mitglieder des Ensembles exakt aufsummiert werden (d.h. nicht für hierarchische Verfahren wie den Barnes-Hut-Algorithmus).

Prozedur beschleunigung Bearbeiten

Diese Prozedur ermittelt die auf einen Körper objekt ausgeübte Beschleunigung. Weitere Eingabeparameter sind die Anzahl N der Mitglieder des Systems, der Plummerradius epsilon sowie die Massen m und Positionen r der einzelnen Mitglieder. m ist ein eindimensionales Array mit N Skalaren, r ein zweidimensionales Array mit N aus jeweils drei Komponenten bestehenden Vektoren. Das Resultat wird durch den dreikomponentigen Vektor a[objekt] wiedergegeben (für das ganze System ist a ebenfalls ein zweidimensionales Array mit N Vektoren).

/* Globale Variablen */

unsigned int objekt;
unsigned int N;
double epsilon;
double *m;
double **r;
double **a;

void beschleunigung (unsigned int objekt, unsigned int N, double epsilon, double *m, double **r,
                     double *a)
{
 
/* Lokale Variablen */
 
  unsigned int i,k;
  double dr[3];
  double d3;
 
  for (k = 0;k < 3;k ++)
  a[k] = 0;
 
  for (i = 0;i < N;i ++)
  {
    if (i != objekt)
    { 
      for (k = 0;k < 3;k ++)
      dr[k] = r[i][k] - r[objekt][k];
      d3 = pow (pow (betrag (dr),2) + pow (epsilon,2),1.5);
      for (k = 0;k < 3;k ++)
      a[k] += G * m[i] * dr[k] / d3;
    }
  }
}

Prozedur ruck Bearbeiten

Diese Routine gibt für einen Körper objekt sowohl dessen Beschleunigung als auch Ruck an. Im Vergleich zur Beschleunigungsprozedur tritt als zusätzlicher Eingabeparameter noch die Geschwindigkeit v aller Mitglieder hinzu, welche ein weiteres zweidimensionales Array mit N Vektoren darstellt. Als zweiter Rückgabeparameter taucht der Vektor j[objekt] auf, welcher den Ruck wiedergibt (und wie a für das ganze Ensemble ein zweidiensionales Array mit N Vektoren ist).

/* Globale Variablen */

unsigned int objekt;
unsigned int N;
double epsilon;
double *m;
double **r;
double **v;
double **a;
double **j;

void ruck (unsigned int objekt, unsigned int N, double epsilon, double *m, double **r, double **v,
           double *a, double *j)
{
 
/* Lokale Variablen */
 
  unsigned int i,k;
  double dr[3];
  double d3,d5;
  double dv[3];
  double skalar;

  for (k = 0;k < 3;k ++)
  {
    a[k] = 0;
    j[k] = 0;
  }

  for (i = 0;i < N;i ++)
  {
    if (i != objekt)
    {
      for (k = 0;k < 3;k ++)
      {
        dr[k] = r[i][k] - r[objekt][k];
        dv[k] = v[i][k] - v[objekt][k];
      }
      d3 = pow (pow (betrag (dr),2) + pow (epsilon,2),1.5);
      d5 = pow (pow (betrag (dr),2) + pow (epsilon,2),2.5);
      skalar = skalarprodukt (dr,dv);
      for (k = 0;k < 3;k ++)
      {
        a[k] += G * m[i] * dr[k] / d3;
        j[k] += G * m[i] * (dv[k] / d3 - 3 * skalar * dr[k] / d5);
      }
    }
  }
}

Prozedur energie Bearbeiten

Diese Prozedur dient der Berechnung der kinetischen und potentiellen Energie eines Mehrkörpersystems. Die Eingabeparameter sind weitgehend die gleichen wie für die Ruck-Routine. Da aber nun das ganze Ensemble und nicht ein einzelnes Mitglied betrachtet wird, fällt der Eingabeparameter objekt weg. Zurückgegeben werden die beiden Energieformen durch die Skalare Ekin und Epot.

/* Globale Variablen */

unsigned int N;
double epsilon;
double *m;
double **r;
double **v;
double Ekin,Epot;

void energie (unsigned int N, double epsilon, double *m, double **r, double **v,
              double *Ekin, double *Epot)
{
 
/* Lokale Variablen */
 
  unsigned int i,j,k;
  double dr[3];
  double d;
 
  *Ekin = 0;
  for (i = 0;i < N;i ++)
  *Ekin += m[i] * pow (betrag (v[i]),2) / 2;
 
  *Epot = 0;
  for (i = 0;i < N;i ++)
  {
    for (j = i + 1;j < N;j ++)
    {
      for (k = 0;k < 3;k ++)
      dr[k] = r[j][k] - r[i][k];
      d = betrag (dr);
      *Epot -= G * m[i] * m[j] / epsilon * (pi / 2 - atan (d / epsilon));
    }
  }
}