Der LowPrecCalculator
Kompakte, niedriggenaue Ephemeriden

Inhalt

• Kurzbeschreibung
• Grundidee der Kompression
• Zur Gewinnung der Binärdateien
• Struktur und Verwendung der Binärdateien
• Die fertigen Dateien zum Herunterladen

Darüberhinaus bietet die Klasse Umrechnungsfunktionen von kartesischen in Polar-, von heliozentrischen in geozentrische und von ekliptischen in äquatoriale Koordinaten (und umgekehrt) an, eine Methode zur Lösung des Keplerproblems (interessant für Asteroidenjäger und alle, denen unser Planetensystem zuwenig Spekulatives enthält und die sich daher Objekten mit ausgedachten Bahnelementen zuwenden), die Delta-T-Funktion nach Stevenson und Morrison (1982) und Methoden für Ekliptikschiefe, Sternzeit und Nutation.

Die Klasse hat schon eine frühere Inkarnation hinter sich. Den Ephemeridenrechner programmierte ich zuerst in C auf dem Macintosh – zu Zeiten, in denen man wirklich noch auf die Kilobytes schaute und Festplatten noch 100 Megabytes hatten (was andererseits auch nicht mehr die Computersteinzeit war). Bei den heutigen Arbeitsspeicher- und Plattenkapazitäten und Verdoppelungszeiten von 18 Monaten für die Speicherdichte ist die Herstellung von möglichst gut komprimierten Ephemeridendateien beinahe ein Luxus, eine Spielerei.

Andererseits haben nicht alle Computerbenutzer Lust, sich jedes Jahr einen neuen Rechner mit noch besseren Eigenschaften zu kaufen. Viele bleiben – wie ich – fünf Jahre und länger ihrem Modell treu. Da ist Platten- und Speicherplatz schon noch ein Thema. Wer eine Webseite unterhält und sich die Monatsauswertungen der User Agents ansieht, wird erstaunt sein, was für Uralt-Browser heutzutage noch im Netz unterwegs sind. Nicht alle diese User sind geneigt, DLL und Datenfiles der Swiss Ephemeris herunterzuladen, bevor sie einen Service mit astrologischen Rechnungen aufrufen können. Andererseits gibt es ein Interesse an Offline-Services, also beispielsweise Applets, die ihren Dienst auch ohne Internetverbindung tun können, nachdem man die Aktion "Webseite offline verfügbar machen" gewählt hat. Intelligente Serverseiten, zum Beispiel Java Server Pages oder Servlets, haben den Vorteil, dass die Ausführung komplizierterer Logik auf dem Server in wohldefinierter Weise erfolgt. Der in ihrer Natur liegende Nachteil ist aber, dass der Offlinebetrieb eben gerade nicht möglich ist.

Wer also ein kleines Applet oder Script mit astrologischen Spezialberechnungen hat, kann den LowPrecCalculator verwenden, um an bogenminutengenaue Planetenstände heranzukommen. Der User kann bei Interesse die Webseite speichern und auch offline verwenden.

So weit, so schön. Was ich bei diesem Ansatz aber unbefriedigend fand, war die direkte Darstellung der gewünschten Planetenkoordinaten als Polynom der Zeit: Interpolation sechsten Grades, ein ganz allgemein verwendbares Verfahren, das hier ohne weitere Fragen auf die Planetenstände angewendet wird. Im Grunde ist es die Holzhammermethode: Man fragt nicht viel nach dem Wesen der zu berechnenden Funktionen, sondern bietet nur eine Tabelle mit den Werten und Interpolationskoeffizienten an.

Es kam – und kommt – mir geschickter vor, in die Lösung das Wissen hineinzustecken, dass die Planeten sich "im wesentlichen" auf Keplerellipsen bewegen, die Bewegung also nicht durch Polynome, sondern durch geeignete Keplerellipsen zu approximieren. Sicher würde es für diesen Zugang nicht genügen, einfach mittlere Bahnelemente als Funktion der Zeit zu berechnen (zum Beispiel durch lineare Regression der zu jedem Zeitpunkt momentan gültigen Keplerellipse, der sogenannten oskulierenden Bahnelemente) - die Positionen würden für die äusseren Planeten viel zu ungenau. Aber wenn man bei diesen mittleren Bahnelementen ein oder zwei besonders kritische Elemente modifiziert und diese Änderungen tabuliert, bestand die Hoffnung, die tatsächlichen Positionen recht gut zu approximieren.

Nun ist ein Bahnelement, die mittlere Anomalie, für solche Modifikationen besonders gut geeignet. Denn auf einer idealen Keplerellipse verändert sie sich proportional mit der Zeit. Man konnte also hoffen, gut interpolierbare Korrekturen zu erhalten, wenn man die mittlere Anomalie so verschiebt, dass der mit der mittleren Keplerellipse und der korrigierten mittleren Anomalie berechnete Planetenstand den "tatsächlichen" Stand (der für die Erstellung der kompakten Tafeln natürlich benötigt wird und z.B. mit der Swiss Ephemeris berechnet werden kann) möglichst gut approximiert.

Bestärkt wurde ich darin durch die Entdeckung, dass Paul Ahnert in seinen "Astronomisch-chronologischen Tafeln" (1960), offenbar mit einer solchen Methode arbeitete: Man berechnete in diesem Tafelwerk die Mittlere Anomalie von - sagen wir Jupiter - indem man Terme für die Jahreszahl, Tageszahl, Monatszahl und Jahrhundertzahl aus der Tafel entnahm und zusammenzählte. Die Mittlere Anomalie fürs Jahrhundert ist dann logischerweise die einzige Zahl, in die konkrete Informationen zur Störung der Ellipse hineingepackt werden kann (denn die Mittlere Anomalie für Jahr, Monat und Tag muss in jedem Jahrhundert gleich funktionieren und kann daher keine konkrete Information eines bestimmten Jahrhunderts enthalten). Wenn man sich den Verlauf der mittleren Anomalie durch eine kleinen Modelleisenbahn veranschaulicht, die in der Periode des Planeten gleichförmig ihre Kreise zieht, bedeutet diese Methode also, die Eisenbahn bei jedem vollen Jahrhundert von den Gleisen zu nehmen und ein bisschen vorher oder hinterher wieder aufzusetzen. Dabei ergeben sich erstaunlich gute Resultate (der Fehler liegt im Bereich eines Grades). Die Genauigkeit liess sich, wie ich feststellte, auf eine Bogenminute verbessern, wenn man einen solchen Verschub der mittleren Anomalie für jedes einzelne Jahr aufstellte - statt nur jahrhundertweise.

Wie ist nun der Korrekturterm der mittleren Anomalie zu berechnen? Die Abbildung zeigt das Vorgehen. Von der Position des ungestörten Planeten (Pm) auf seiner mittleren Bahn weicht die tatsächliche Position (Pw) um einiges ab. Der Planet kann ausserhalb der Bahnebene liegen, einen anderen Abstand zur Sonne haben und eine andere Anomalie haben als es die Rechnung mit den mittleren Bahnelementen ergibt. Um den Verschub der mittleren Anomalie zu berechnen, gehen wir wie folgt vor: Wir projizieren die tatsächliche Position Pw des Planeten auf seine mittlere Bahnebene. In der Bahnebene berechnen wir die Polarkoordinaten des Projektionspunktes, wobei die Sonne den Mittelpunkt und das Perihel die Nullrichtung kennzeichnen. Die durch den Projektionspunkt gegegene Richtung als wahre Anomalie betrachtend, berechnen wir die zugehörige mittlere Anomalie. Die Differenz zur mittleren Anomalie des mittleren Bahnpunkts ist der zu tabulierende Korrekturterm für "ΔM". Bei ersten Versuchen mit ΔM stellte ich fest, dass mit dem skizzierten Verfahren für Uranus, Neptun und Pluto bereits Bogenminutengenauigkeit erreicht wird. Bei Jupiter und Saturn lagen die Fehler jedoch manchmal auch bei 5' und mehr: Wegen der grösseren Nähe ihrer Bahnen zur Erde führte auch die Abweichung des Radiusvektors zu nicht mehr ignorierbaren Verzerrungen. So ergänzte ich das obige Verfahren um einen Korrekturfaktor für den aus dem Keplerproblem zu berechnenden Radiusvektor: Berechne den Abstand des Projektionspunktes von der Ellipse, gemessen relativ zum Radiusvektor, den der Planet auf der mittleren Ellipse hätte.

Mit diesen beiden Werten, ΔM und ΔR, erreichte ich nun für alle Planeten Bogenminutengenauigkeit. Der Teil meines C-Programms Elements.c, der für diese Umrechnung zuständig ist, ist der folgende:

/* Ort in heliozentrischen Polarkoordinaten berechnen (für den r-Fehler) */
            se_flag = SEFLG_HELCTR + SEFLG_NONUT + SEFLG_SPEED;
            rc = swe_calc(jd, planet, se_flag, x, msg);
            hl = x[0];
            hb = x[1];
            hr = x[2];

/* Berechnung der mittleren Bahnelemente */
            get_elements(planet, jd, &ax_m, &ex_m, &inc_m, &node_m, &omega_m, &an_m, &n_m);

/* Berechnung der Bahnkoordinaten ( l = Differenz Planet - Knoten ) */
            b = asin (-sin(inc_m*c)*cos(hb*c)*sin((hl-node_m)*c)
              + cos(inc_m*c)*sin(hb*c))/c;
            l = atan2( cos(inc_m*c)*cos(hb*c)*sin((hl-node_m)*c)+sin(inc_m*c)
              * sin(hb*c), cos(hb*c)*cos((hl-node_m)*c) )/c;
            if (l < 0) l +=  .36e3;

/* Mittlere Anomalie aus wahrer Anomalie berechnen (Keplerproblem rückwärts!) */
            v = l - omega_m + node_m;
/* Zunächst exzentrische Anomalie: */
            e = 2.*atan(tan(c*v/2.)*sqrt((1.-ex_m)/(1.+ex_m)));
/* Daraus mittlere Anomalie (Keplergleichung) */
            man = (e - ex_m*sin(e))/c;
            if (man<0) man += 360.;

/* Für den r-Fehler: das Keplerproblem zu den mittleren Bahnelementen lösen */
            kepler1(ex_m,ax_m,omega_m,node_m,inc_m,man,&l1,&b1,&r1);

/*  Ausdrucken : M.An. und Radiusdifferenz */
            dan = arcdiff(man,an_m)*1000;
            dr  = (hr-r1)/r1*10000.;
            fprintf(ephem,"%5.0lf %4.0lf\n",dan,dr);

• Zuerst berechnete ich mit Swiss Ephemeris in 100-Tages-Abständen oskulierende Bahnelemente der äusseren Planeten mit Bezug auf die mittlere Ekliptik des Datums und speicherte diese in Dateien.
• Mathematica las diese Dateien ein und führte eine polynomiale Regression dieser Bahnelemente für den gesamten 100-Jahres-Bereich durch. Ergebnis hiervon waren Koeffizienten, durch die eine für das jeweilige Jahrhundert durchschnittlich besonders gut passende Bahn bestimmt war.
• Nun kam wieder das C-Programm an die Reihe und berechnete auf Basis dieser mittleren Bahn jahresweise die Abweichungen ΔM und ΔR, wie oben beschrieben.

arcdiff[x_,y_] := x-y-Floor[(x-y)/360 + 0.5]*360;
For[jd0 = 150, jd0 <= 950, jd0+=100,
Print[jd0];
(* Oskulierende Elemente in 100d-Schritten erzeugen  *)
Command =
    ToString[StringForm[
        "C:\Programme\Borland\CBuilder5\Projects\Elements\Elements.exe -o \
         -S366 -i100 -j``",jd0]];
Run[Command];
(* Datei für Koeff. der mittleren Elemente öffnen  *)
stMeanElements =
    OpenWrite[
      ToString["C:\Programme\Borland\CBuilder5\Projects\Elements\c" <>
          StringTake[ToString[jd0+10000],{2,5}] <> ".dat"],
      FormatType->OutputForm];
For[planet = 2, planet <=9, planet++,
(* Oskulierende Elemente einlesen (Piping des kleinen Mannes) *)
elem = ReadList[
        ToString[
          StringForm[
            "C:\Programme\Borland\CBuilder5\Projects\Elements\Elem``",
            planet]],{Number,Number,Number,Number,Number, Number, Number,
          Number}];
(* Die winkelförmigen Elemente hochrechnen (keine 360°-Sprünge!) *)
size=Dimensions[elem][[1]];
For[j=5,j<=7,j++,
		For[ i=2, i <= size, i++,
        elem[[i,j]] = elem[[i-1,j]] + arcdiff[elem[[i,j]],elem[[i-1,j]]]]];
(* Regression 2. Ordnung durchführen *)
(* Zeilenstruktur: jd,ax,ex,inc,node,omega,an,n *)
For[i=2,i<=8,i++,
		elemfit =
        N[Fit[Transpose[elem][[i]],If[ planet<=4,{1,t-1},{1,t-1,(t-1)^2}],t],
          10];
(* Koeffizienten fortschreiben *)
	 Write[stMeanElements,
        TableForm[
          Table[CForm[Coefficient[elemfit,t,k]],{k,0,If[planet<=4,1,2]}],
          TableDirections -> {Row}]]
		]
];
Close[stMeanElements];
(* Binärdatei jnnnn.dat mit ΔM und ΔR errechnen: *)
Command =
    ToString[StringForm[
        "C:\Programme\Borland\CBuilder5\Projects\Elements\Elements.exe -C \
         -S366 -i100 -j``",jd0]];
Run[Command];
]

Wenn nun in der Java-Klasse LowPrecCalculator die Position eines (äusseren) Planeten berechnet werden soll, werden die folgenden Schritte durchgeführt:

• Von der Klasse Elements (die für das Lesen und Auswerten der binären Datenfiles zuständig ist) werden die mittleren Bahnelemente und die beiden Korrekturterme ΔM und ΔR angefordert.
• ΔM wird sofort in die mittlere Anomalie eingerechnet
• Nun wird das Keplerproblem mit den mittleren Bahnelementen zur so modifizierten mittleren Anomalie gelöst.
• Nach Berechnung der heliozentrischen Polarkoordinaten wird der Abstand von der Sonne mit dem Wert 1 + ΔR multipliziert (Korrektur des Radiusvektors).
• Das Ergebnis ist die (niedriggenaue) heliozentrische Position des Planeten. Diese wird nun unter Verwendung der Sonnenposition in geozentrische Koordinaten umgerechnet.

• Das julianische Startdatum als Drei-Byte-Integer. Für das File jnnnn.dat habe ich dabei den 1.1. des Jahres nnnn, 12 Uhr Ephemeridenzeit als Startdatum gewählt (was immer eine ganze Zahl für das Julianische Datum ergibt).
• Das julianische Endedatum (Datum des letzten Records) als Drei-Byte-Integer. Das Endedatum ergibt sich aus dem Startdatum durch die Formel jdend = jdstart + span * increment . Dabei ist span die Anzahl der Records und increment das Inkrement in Ephemeridentagen von einem Record zum nächsten.
• Das Inkrement in Tagen von einem Record zum nächsten, gespeichert als Zwei-Byte-Integer.
• Es folgen die Koeffizienten der sieben Bahnelemente a, e, i, Ω, ω, M, n der acht Planeten Merkur, Venus, Mars, Jupiter, Saturn, Uranus, Neptun und Pluto. Für die inneren Planeten Merkur, Venus und Mars wurde nur lineare Regression im Zeitraum ausgeführt, für die äusseren Regression mit Parabeln. Für die inneren Planeten braucht man daher pro Element zwei Koeffizienten, für die äusseren drei. Auch diese Zahlen werden in der Big Endian Byteordnung gespeichert, als 8 Byte Gleitkommazahlen gemäss der IEE 754 Spezifikation. Sie werden in folgender Reihenfolge in die Datei geschrieben:

  a_me_0 a_me_1 e_me_0 e_me_1 i_me_0 i_me_1 Ω_me_0 Ω_me_1 ω_me_0 ω_me_1 M_me_0 M_me_1 n_me_0 n_me_1
  a_ve_0 a_ve_1 e_ve_0 e_ve_1 i_ve_0 i_ve_1 Ω_ve_0 Ω_ve_1 ω_ve_0 ω_ve_1 M_ve_0 M_ve_1 n_ve_0 n_ve_1
  a_ma_0 a_ma_1 e_ma_0 e_ma_1 i_ma_0 i_ma_1 Ω_ma_0 Ω_ma_1 ω_ma_0 ω_ma_1 M_ma_0 M_ma_1 n_ma_0 n_ma_1
  a_ju_0 a_ju_1 a_ju_2 e_ju_0 e_ju_1 e_ju_2 i_ju_0 i_ju_1 i_ju_2 Ω_ju_0 Ω_ju_1 Ω_ju_2 ω_ju_0 ω_ju_1
  ω_ju_2 M_ju_0 M_ju_1 M_ju_2 n_ju_0 n_ju_1 n_ju_2 a_sa_0 a_sa_1 a_sa_2 e_sa_0 e_sa_1 e_sa_2 i_sa_0
  i_sa_1 i_sa_2 Ω_sa_0 Ω_sa_1 Ω_sa_2 ω_sa_0 ω_sa_1 ω_sa_2 M_sa_0 M_sa_1 M_sa_2 n_sa_0 n_sa_1 n_sa_2
  a_ur_0 a_ur_1 a_ur_2 e_ur_0 e_ur_1 e_ur_2 i_ur_0 i_ur_1 i_ur_2 Ω_ur_0 Ω_ur_1 Ω_ur_2 ω_ur_0 ω_ur_1
  ω_ur_2 M_ur_0 M_ur_1 M_ur_2 n_ur_0 n_ur_1 n_ur_2 a_ne_0 a_ne_1 a_ne_2 e_ne_0 e_ne_1 e_ne_2 i_ne_0
  i_ne_1 i_ne_2 Ω_ne_0 Ω_ne_1 Ω_ne_2 ω_ne_0 ω_ne_1 ω_ne_2 M_ne_0 M_ne_1 M_ne_2 n_ne_0 n_ne_1 n_ne_2
  a_pl_0 a_pl_1 a_pl_2 e_pl_0 e_pl_1 e_pl_2 i_pl_0 i_pl_1 i_pl_2 Ω_pl_0 Ω_pl_1 Ω_pl_2 ω_pl_0 ω_pl_1
  ω_pl_2 M_pl_0 M_pl_1 M_pl_2 n_pl_0 n_pl_1 n_pl_2
  

• Die Minima und Maxima für ΔM und ΔR für die Planeten Jupiter, Saturn, Uranus, Neptun und Pluto. Jede dieser Zahlen wird als 2-Byte-Integer fortgeschrieben. Das ist möglich, da es keine ΔM-Werte gibt, die 32.767° überschreiten, und erst recht keine ΔR Korrekturfaktoren, die den Wert 3.2767 übersteigen! Die Reihenfolge ist (dm steht für ΔM, dr für ΔR):

  dm_ju_min dr_ju_min dm_ju_max dr_ju_max ... dr_pl_max

recordlength = floor(log(dm_ju_max-dm_ju_min)/log(2)+1)
             + ...
             + floor(log(dr_pl_max-dr_pl_min)/log(2)+1)

// Und schliesslich die Daten selbst
    bitzeiger = 0;  // Höchste relevante Bitposition in buf
    buf = 0;        // Bitpuffer
    for (k = 0; k <= span + 2; k++) {
      for (i = 0; i <= (SE_PLUTO - SE_JUPITER); i++)
        for (j = 0; j<2 ; j++) {
          buf <<= size[j][i]; // Platz machen
          buf |= (( x[j][i][k] - min[j][i] ) & mask[j][i]);  // Wert einfügen
          bitzeiger += size[j][i];  // Zeiger inkrementieren
          while (bitzeiger > 8) {
            fprintf(compress,"%c", ( buf >> (bitzeiger - 8) ) & 0xFF);
            buf = ( ( 1 << ( bitzeiger - 8 ) ) - 1 ) & buf;
            bitzeiger -= 8;
            }
         }
       }
    if (bitzeiger > 0) fprintf(compress, "%c", buf & 0xff );

private static void get_workarea( double jd ) {

  long offset, buffer = 0;
  byte size0, x, bits, i, k, pl;
  double x1,x2;

  if (mAnomCoeff  == null) mAnomCoeff = new double[3][5];
  if (deltaRCoeff == null) deltaRCoeff = new double[3][5];
// Arbeitsbereiche aus binärem Datenstrom einlesen
  offset = (long) ((jd - jd_start) / increment);
  jd_act = jd_start + increment * offset;
  offset *= recordlength;  // Bit-Offset
  for (i = 0; i < 3; i++ ) {       // 3 Zeiten
    for ( pl = 0; pl < 5; pl++ )   // 5 Planeten
      for ( k = 0; k < 2; k++ )  {     // 2 Zahlen ( Anomalie, deltaR )
        size0  = size[k][pl]; // so viele Bits müssen eingelesen werden
        buffer = 0;
        while ( size0 > 0 ) {
          x = (byte) ( data[(int) (offset / 8)] & ( 0xff >> ( offset % 8 ) ) ) ;
          bits = (byte) ( 8 - (offset % 8) ) ; // so viele Bits sind "eingelesen"
          if ( size0 < bits ) {
            x = (byte) ((x & 0xff) >> ( bits - size0 ));
            bits = size0;
            }
          buffer = ( buffer << bits ) + ( x & 0xff);
          size0 -= bits;
          offset += bits;
          }
// Nun sind size0 bits eingelesen, Ergebnis steht in buffer
        switch(k) {
          case 0:
            mAnomCoeff[i][pl] =  (double) min[k][pl] + buffer;
            break;
          case 1:
            deltaRCoeff[i][pl] = (double) min[k][pl] + buffer;
            break;
          }
        }
      }

// Koeffizienten für quadr. Polynom aufbereiten
// Funktion f an Stützstellen 0, 1, 2 durch Parabel interpolieren:
//
// q(x) =        f(0)
//      +  x   * ( 2*f(1) - (f(2)+3*f(0))/2 )
//      +  x^2 * ( -f(1)  + (f(2)+f(0))/2   )
//

  for (pl = 0; pl < 5; pl++) {
    x1 = mAnomCoeff[1][pl];
    x2 = mAnomCoeff[2][pl];
    mAnomCoeff[1][pl] = 2*x1 - ( x2 + 3*mAnomCoeff[0][pl] )/2;
    mAnomCoeff[2][pl] = ( x2 + mAnomCoeff[0][pl] ) / 2 - x1;
    x1 = deltaRCoeff[1][pl];
    x2 = deltaRCoeff[2][pl];
    deltaRCoeff[1][pl] = 2*x1 - ( x2 + 3*deltaRCoeff[0][pl] )/2;
    deltaRCoeff[2][pl] = ( x2 + deltaRCoeff[0][pl] ) / 2 - x1;
    }

  }