Der LowPrecCalculator
Kompakte, niedriggenaue Ephemeriden

Inhalt


Kurzbeschreibung

Die Klasse LowPrecCalculator ist ein niedriggenauer Ephmeridenrechner auf der Grundlage kleiner, komprimierter Ephemeridendateien. Mond und Merkur bis Pluto werden auf die Bogenminute genau gerechnet, die Sonne mithilfe der Newcombschen Reihe auf die Bogensekunde genau (um auch die Berechnung von Solarhoroskopen zu ermöglichen). Eine Instanz führt stets ein Julianisches Datum, eine Länge und eine Breite sowie die dazugehörigen Planeten- und Hauspositionen nach Placidus (diese werden mit der Klasse Mundan berechnet).

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.


Zum Inhaltsverzeichnis

Grundidee der Kompression

Immer schon am Thema der Planetenberechnung interessiert, besitze ich Publikationen wie die "Planetary Programs and Tables from -4000 to +2800" von Pierre Bretagnon und Jean-Louis Simon (Richmond 1986). Bretagnon und Simon hatten BASIC- und FORTRAN-Programme, mit denen sie die Planetenstände im genannten Zeitraum recht genau berechnen konnten. Während die inneren Planeten Sonne, Merkur, Venus und Mars [der astronomische Begriff des Planeten, den ich nicht verwende, ist mir bekannt] leicht mit je einer einzigen Reihe für Länge, Breite und Radius dargestellt werden konnten, gaben die Verfasser für die äusseren Planeten Polynome sechsten Grades für Länge, Breite und Radius an, deren Koeffizienten sie in Abständen von fünf Jahren auflisteten. Damit hatten sie ein 60mal schmaleres Werk als die berühmten zweibändigen "Tuckerman Tables" geschaffen, das dennoch viel genauere Planetenstände liefert (die "Tuckerman Tables" waren über Jahrzehnte die Referenzephemeride für Historiker. Auch Klöckler empfiehlt sie historisch interessierten Astrologen in seinem "Kursus der Astrologie" (Leipzig 1930, Abschnitt II.A.a.5.)).

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.

Projektion auf die mittlere Bahnebene 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);

Zum Inhaltsverzeichnis

Zur Gewinnung der Binärdateien

Die Herstellung der Binärdateien ging in einem orchestrierten Einsatz eines C-Programms und eines Mathematica-Notebooks vor sich. Da meinem Mathematica-Handbuch keine Informationen darüber zu entlocken waren, wie ein Notebook von der Kommandozeile aufzurufen ist, konnte ich kein separates Skript schreiben, das die Einsätze beider Programme steuert. Ich verwendete daher Mathematica selbst für diese Steuerung und rief die C-Teile in separaten Prozessen auf. Die Kommunikation zwischen Mathematica und C erfolgte über Dateien. Das Vorgehen war im Prinzip folgendes: Hier das Mathematica-Notebook, das den Einsatz steuerte:
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];
]
Um Platz zu sparen, wurden die minimalen und maximalen ΔR- und ΔM-Werte des jeweiligen Jahrhundertfiles ermittelt. Für einen bestimmten Zeitpunkt wurde dann nicht der Wert selbst gespeichert; stattdessen multiplizierte ich die Differenz vom minimalen Wert mit 1000 (für ΔM) bzw. 10.000 (für ΔR) und speicherte die sich ergebende Ganzzahl mit der für ihre Darstellung benötigten Anzahl von Bits. Wenn die ΔM- Werte in einem Jahrhundert beispielsweise zwischen 1,003° und 2,006° variierten, so werden für die einzelnen ΔM-Werte pro Eintrag 10 Bits benötigt (denn es müssen die Zahlen von 0 bis 1003 dargestellt werden können, wofür man 10 Bits benötigt, da die nächstgrössere Zweierpotenz 2^10 = 1024 ist). Zählt man die für jeden einzelnen Planeten benötigten Bits zusammen, so erhält man die Grösse eines Records in Bits. Diese Records werden hintereinander in der Binärdatei gespeichert.

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

Struktur und Verwendung der Binärdateien

Die Binärdateien enthalten am Anfang einen Verwaltungsteil. Diesem folgt der eigentliche Datenteil. Der Verwaltungsteil enthält folgende Zahlen, die in ihren elementaren Datenformaten in der Big Endian Byteordnung geschrieben wurden: Hiermit ist der Verwaltungsteil abgeschlossen. Mit dem folgenden Byte beginnt nun ein binärer Datenstrom der ΔM- und ΔR-Werte. Dabei speichere ich nicht ΔM und ΔR selbst, sondern deren Differenz vom Minimalwert der Periode ab. Die Grösse der einzelnen Records in Bits berechnet sich daher wie folgt:
recordlength = floor(log(dm_ju_max-dm_ju_min)/log(2)+1)
             + ...
             + floor(log(dr_pl_max-dr_pl_min)/log(2)+1)
Das bedeutet, der Offset eines einzelnen Records lässt sich aus den Daten des Verwaltungsteils errechnen, so dass zum Lesen der Elemente für ein bestimmtes Datum ein random file access möglich ist (da die Dateien aber so klein ausfallen, lese ich sie einmalig ein und halte sie dann in aufbereiteter Form im Puffer der Klasse Elements. ). Der folgende Teil des Programms Elements.c ist für die Erzeugung dieses binären Datenstroms zuständig: Eine lange Ganzzahl ohne Vorzeichen (buf) dient als Puffer, in das die Zahlen nacheinander bitweise von rechts nach links hereingeschoben werden. Nach dem Einschieben einer Zahl werden die Bits byteweise von links in die Datei geschrieben und aus dem Puffer gelöscht.
// 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 );
Zur Illustration folgt hier noch die Codestrecke in der Klasse Elements, die für das Einlesen der Binärdaten und die Erzeugung der Elemente zuständig ist. Der Datenteil wird in der Klasse Elements in einem Array von Bytes namens data[] gepuffert. In der Methode get_workarea() wird dieser Byte-Array für ein bestimmtes Julianisches Datum ausgewertet. Es werden drei aufeinanderfolgende Records eingelesen, so dass das angegebene Julianische Datum zwischen die Daten des ersten und des zweiten Records fällt. Aus diesen Records werden dann die Arrays deltaRCoeff[][] und mAnomCoeff[][] quadratischen Interpolationskoeffizienten für ΔR und ΔM gewonnen, die ebenfalls als globale Attribute der Klasse gehalten werden (die "Workarea"). Kommen nun in Folge mehrere nahe beieianderliegende Datumswerte, so kann die Workarea direkt wiederverwendet werden, und die jeweiligen ΔR und ΔM werden durch einfache quadratische Interpolation berechnet.
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;
    }

  }

Zum Inhaltsverzeichnis

Die fertigen Dateien zum Herunterladen

Eine Reihe fertig generierter Dateien zum Herunterladen finden Sie unter http://www.astrotexte.ch/applets/rsc/.
Zum Inhaltsverzeichnis Zur Doku-Startseite Ohne Hintergrund drucken Zurück zur Homepage