Klasse Elements

package AstroOL;

import java.io.*;
/** Datenklasse für die Berechnung der Bahnelemente.
    Diese Hilfsklasse von {@link LowPrecCalculator LowPrecCalculator}
    berechnet und puffert die angeforderten
    Bahnelemente der Planeten. Zentrale Methode ist
    die Methode {@link #get(double,int,double[]) get()}, die ein Array mit
    den Bahnelementen eines Planeten zu einer gegebenen
    Zeit zurückgibt.

    @see <a href="./doc-files/LowPrecCalculator.html">Datenmodell und die Kompressionsmethode des LowPrecCalculator</a>

    */
public class Elements {

// Verwaltungsdaten
  private static long jd_start = 0;
  private static long jd_end = 0;
  private static short increment;
// Arbeitsbereich mit Koeffizienten für quadratische Interpol.
  private static long jd_act; // Datum für erste Zeile des Arb.ber.
  private static double mAnomCoeff[][],deltaRCoeff[][];
// Koeffizienten für quadratische Regression
  private static double a[][][];
  private static double a0[][][];
// Verwaltungsdaten für Binärdaten mit mittl. Anomalien
  private static short min[][], max[][];
  private static byte size[][];
  private static int recordlength; // in Bit
// Binärdaten mit den mittleren Anomalien
  private static byte data[];

  static void readFile(int jahr0) {
     int bytesAvailable, planet, j;
     short n, m;
     byte i;
     InputStream in;
     DataInputStream din;
     byte int3[] = new byte[3]; // 3-Byte Integer
     String filename = "/rsc/data/J" +
         new Integer(10000 + jahr0).toString().substring(1,5) + ".dat";
// File öffnen, Daten lesen
    try {
      in = Elements.class.getResourceAsStream(filename);
      bytesAvailable = in.available();
//      System.out.println("Grösse:" + bytesAvailable + "\n");
      in.read(int3);
      jd_start = ( (int3[0] & 0xFF) << 16 ) + ( (int3[1] & 0xFF) << 8 ) + (int3[2] & 0xFF);
      in.read(int3);
      jd_end   = ( (int3[0] & 0xFF) << 16 ) + ( (int3[1] & 0xFF) << 8 ) + (int3[2] & 0xFF);
//      System.out.println("jd_start  = " + jd_start);
//      System.out.println("jd_end    = " + jd_end);
      din = new DataInputStream(in);
      increment = din.readShort();
//      System.out.println("increment = " + increment);
// Dann die Koeffizienten
     a0 = new double[3][7][2];
      for (planet = 0; planet < 3; planet++ )
        for (i = 0; i < 7; i++)
          for (j = 0; j < 2; j++) {
            a0[planet][i][j] = din.readDouble();
            }
     a = new double[5][7][3];
      for (planet = 0; planet < 5; planet++ )
        for (i = 0; i < 7; i++)
          for (j = 0; j < 3; j++) {
            a[planet][i][j] = din.readDouble();
            }

// Verwaltungsteil für mittlere Anomalien einlesen
      min = new short[2][5];
      max = new short[2][5];
      size = new byte[2][5];
      recordlength = 0;
      for (planet = 0; planet < 5; planet++ )
        for (j = 0; j < 2; j++ ) {
          min[j][planet] = din.readShort();
          max[j][planet] = din.readShort();
          m = 1;
          i = 0;
          n = (short) (max[j][planet] - min[j][planet]);
          while ( m <= n) { i++; m <<= 1; }
          size[j][planet] = i;
          recordlength += size[j][planet];
          }


// Nun Datenteil für mittl. Anomalien einlesen
      bytesAvailable = in.available();
      data = new byte[bytesAvailable];
      in.read(data);
      }
    catch(Exception e) {
      System.out.println("Fehler beim Oeffnen der Datei" + filename);
      return;
      }
  }


/** Bahnelemente von Planet x für Zeitpunkt jd berechnen.
    Dies ist die zentrale Methode der Klasse Elements. Sie gibt
    Bahnelemente für die Zeit jd und den Planeten planet zurück,
    die eine möglichst gute Positionsberechnung für den gewünschten
    Zeitpunkt gestatten. Die aktuelle Implementierung verwendet
    hierzu eine Binärdatei, die pro Jahrhundert ca. 4.7 KB umfasst.
    Die Datei enthält die linear regredierten Elemente aller Planeten
    sowie für die äusseren Planeten Korrekturterme für Radius
    und mittlere Anomalie. Der Korrekturterm für die mittlere Anomalie
    wird sofort angebracht. Der Korrekturterm für den Radius kann
    erst nach Lösung des Keplerproblems angebracht werden und wird
    daher als "8. Bahnelement" (d.h. an achter Stelle im Rückgabe-Array
    elements[]) zurückgereicht.
    @param jd Julianisches Datum (in Ephemeridenzeit)
    @param planet Nummer des Planeten. Die Numerierung der Planeten
    entspricht ihrer Reihenfolge (mit 0 = SO beginnend) im Attribut
    {@link Calculator#pknam Calculator.pknam}.
    @param elements Rückgabearray mit den Bahnelementen und dem
    Korrekturwert für den Radius. Die Reihenfolge ist
         <table align="center">
         <tr><td>Grosse Halbachse</td><td><code>= elements[0]</code></td></tr>
         <tr><td>Exzentrizität   </td><td><code>= elements[1]</code></td></tr>
         <tr><td>Bahnneigung     </td><td><code>= elements[2]</code></td></tr>
         <tr><td>Aufst. Knoten   </td><td><code>= elements[3]</code></td></tr>
         <tr><td>Perhihellänge   </td><td><code>= elements[4]</code></td></tr>
         <tr><td>Mittl. Anomalie </td><td><code>= elements[5]</code></td></tr>
         <tr><td>Tagesbewegung   </td><td><code>= elements[6]</code></td></tr>
         <tr><td>Korrektur für r </td><td><code>= elements[7]</code></td></tr>
         </table>
         */
public static void get( double jd, int planet,
                        double[] elements ) {
  int jahr0,i = planet - 5, j, k;
  double t;

  if ( (jd_start > jd) || (jd > jd_end) ) {
// Out of Range: Datenfile muss nachgelesen werden
     readFile((Calculator.calendarDate(jd)[0] - 50)/100*100+50);
     }

  t = ( jd - jd_start ) / increment;
  if (planet >= 5) {
    i = planet - 5;
    for (j = 0; j < 7; j++)
      elements[j] = ( a[i][j][2] * t + a[i][j][1] ) * t + a[i][j][0];
    }
  else {
    i = planet - 2;
    for (j = 0; j < 7; j++)
      elements[j] = a0[i][j][1] * ( t * 10 + 1 ) + a0[i][j][0];
    return;  // Keine Anomalie/Radiusdifferenzen verwendet
    }

// Anomalie / deltaR berechnen:

// Ist Arbeitsbereich nicht mehr verwendbar? Dann nachlesen
  if ( ( jd_act > jd ) || ( jd_act + 2*increment < jd) ) {
    get_workarea( jd );
    }

  t = (jd - jd_act ) / increment;
  elements[5] = Calculator.fmod(elements[5]+(( mAnomCoeff[2][i] * t + mAnomCoeff[1][i] )* t + mAnomCoeff[0][i])/1.e3,360.);
  elements[7] =  ( deltaRCoeff[2][i] * t + deltaRCoeff[1][i] )* t + deltaRCoeff[0][i];

}

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;
    }

  }

}