//---------------------------------------------------------------------------
#include <stdio.h>
#include <math.h>
#define USE_DLL
#include "Swedll.h"
#include "astrolib.h"

#pragma hdrstop

//---------------------------------------------------------------------------

/*

Erzeuge ein komprimiertes Ephemeridenfile, z.B. "j1850.dat"

Datenformat:

HEADER
   jd1
     Erstes Julianisches Datum des Files  (3 Byte Integer)
   jd2
     Letztes Julianisches Datum des Files (3 Byte Integer)
   inc
     Inkrement in Tagen ( 2 Byte Integer )
   a[10][7][2-3]
     Mittlere Bahnelemente für den Gültigkeitsbereich,
     ( Koeffizienten für das Hornerschema 2. Ordnung )
   (min,max)[2][5]
     Kleinste und grösste Werte für Delta r und Delta M
     in 1/10000 und 1/100.000.

DATEN
  n1 Bit für Planet1 am ersten Datum
  n2 Bit für Planet2 am ersten Datum
  ...
  n1 Bit für Planet1 am zweiten Datum
  ...

Platzbedarf für Datenteil, grob geschätzt:
4 Byte * 5 Planeten * 365 = 7,3 KB pro Jahrhundert

 */


#include "coor.c"
#include "kepler.c"

#pragma argsused
/* Reduzierte Massen */
double rm[6] = {0.,0.954786104043041709e-03,
      0.28558373315e-03,
      0.437273164545891804e-04,0.517759138448793616e-04,
      0.277777777777777770e-05};

// Benutzerdialog möglich? (Option -s)
        bool consoleDialog = false;
// Beginn des Anomalienfiles (Option -j)
        int jahr0 = 1850;
// Spanne des Anomalienfiles in Vielfachen des Inkrements (Option -S)
        unsigned int span = 366;
// Inkrement von Zeile zu Zeile in Tagen (Option -i)
        unsigned int increment = 100;
// Koeffizienten für die quadratische Regression
        double a[8][7][3];
// Some globals make life easier
        double jd, jd_start, x[6], y[6];
        double ax,ex,inc,node,omega,n,an;
        int i,rc,se_flag = 0, vorlauf=20, planet = SE_PLUTO;
        char msg[255];
        FILE* daten;


void
// Vergleich mittl. Elemente (ohne m. Anomalien) - "wahre" Werte
// bei Berechnung von Bahnpositionen
     testElements(),
// Vergleich mittl. Elemente - oskulierende Elemente
     testMeanElements(),
// Vergleich mittl. Elemente (mit m. Anomalien) - wahre Werte
     testMeanAnomalies(),
// Berechnen der mittl. Anomaliedifferenzen, im Textformat
// ablegen (Ephem5...Ephem9)
// Voraussetzung: mittl. Bahnelemente wurden bereits berechnet
     computeMeanAnomalies(),
// Reichweite oskulierender Elemente testen
     testOrbitConversion(),
// Oskulierende Bahnelemente in 100d-Abständen berechnen und
// in Elem5..Elem9 ablegen
     computeOsculatingElements(),
// Erzeuge binäres Datenfile aus Ephem5...Ephem9
     compressMeanAnomalies(),
// Eine doppeltgenaue Zahl im Big Endian Format schreiben
     putDouble(FILE* file, double x),
     test();  // für dies und das
void getElementCoefficients();
// Einlesen der durch quadr. Regression ermittelten Bahnelement-Koeffizienten
// Und Auswerten der Regressionspolynome zum Zeitpunkt x
void get_elements(int planet, double x, double *ax, double *ex, double *inc,
                  double *node, double *omega, double *an, double *n);

int main(int argc, char* argv[])
{

        char option, c;
        char *p;

        swe_set_ephe_path("C:\\SWEPH\\EPHE");

// Parameter einlesen
          if ( argc <= 1) {
            printf("Bitte eine Ausführungsart angeben\n");
            exit(EXIT_FAILURE);
            }
          option = ' ';
          while ( --argc > 0 ) {
            if (argv[argc][0] != '-') continue;
            c = argv[argc][1];
            if ( strchr("KaAEoOcCt", c) != NULL ) {
              option = c;
              continue;
              }
            if ( c == 's' ) {
              consoleDialog = true;
              continue;
              }
            if ( c == 'j' ) {
              sscanf(argv[argc]+2,"%d",&jahr0);
              continue;
              }
            if ( c == 'i' ) {
              sscanf(argv[argc]+2,"%d",&increment);
              continue;
              }
            if ( c == 'S' ) {
              sscanf(argv[argc]+2,"%d",&span);
              continue;
              }
            printf("Parameter '%s' kann nicht interpretiert werden\n",argv[argc]);
            exit(EXIT_FAILURE);
            }

       if ( option == ' ') {
          printf("Erlaubte Betriebsarten: KaAEoOcC\n");
          exit(EXIT_FAILURE);
          }

// Startjahr berechnen
       jd_start  = swe_julday(jahr0, 1, 1, 12., SE_GREG_CAL);

       switch(option) {
          case 'K':
            testElements();
            break;
          case 'a':
            computeMeanAnomalies();
            break;
          case 'A':
            testMeanAnomalies();
            break;
          case 'o':
            computeOsculatingElements();
            break;
          case 'E':
            testMeanElements();
            break;
          case 'O':
            testOrbitConversion();
            break;
          case 'C':
            computeMeanAnomalies();
            compressMeanAnomalies();
            break;
          case 'c':
            compressMeanAnomalies();
            break;
          case 't':
            test();
            break;
          }

        c = ' ';
        printf("\nDone.");
        if (consoleDialog)

          while (c!='x') c = getc(stdin);

}


// Dies und das testen
void test() {

}

// Eine doppeltgenaue Zahl im Big Endian Format schreiben
void putDouble(FILE* file, double x) {
  unsigned char *pc, *pc0;
  pc0 = (unsigned char*) &x;
// Pentium speichert Little Endian -> Reihenfolge umkehren für Big Endian
// IEE 754 Datenformat wird auch von Borland verwendet
  for (pc = pc0 + 7; pc >= pc0; pc--) fprintf(file,"%c",(unsigned char) *pc);

}

void testElements() {

    double c = CONV,x0;
    double b,l,r,hb,hl,hr,v,e,man,lmax[5],bmax[5],rmax[5],l1,b1,r1;
    double rd,ax_m,ex_m,inc_m,node_m,omega_m,n_m,an_m;
    double jd0,man0,rd0,jd1,man1,rd1;
    FILE *ephem, *daten;
    char filename[20], test[80];
    long offset;

    daten = fopen("daten.txt","w");

    randomize();
    randomize();
    for (planet = SE_JUPITER; planet <= SE_PLUTO; planet++) {
          lmax[planet-SE_JUPITER] = 0.;
          bmax[planet-SE_JUPITER] = 0.;
          rmax[planet-SE_JUPITER] = 0.;
          for (i = 0; i < 200; i++) {
            jd = random(span)*increment + random(increment) + jd_start;

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

/* Berechnung der exakten heliozentrischen Position */
            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];


/*  Ausdrucken : Differenz der hel. Längen */
            fprintf(daten,"%10.1lf %12.7lf %12.7lf %12.7lf\n",jd,arcdiff(hl,l1),hl,l1);

/* Jetzt Vergleichen !! */
           l1 = fabs(arcdiff(hl,l1));
           if (l1>lmax[planet-SE_JUPITER]) lmax[planet-SE_JUPITER] = l1;
           b = hb - b1;
           if (fabs(b)>bmax[planet-SE_JUPITER]) bmax[planet-SE_JUPITER]=fabs(b);
           r1=fabs(hr-r1);
           if (r1>rmax[planet-SE_JUPITER]) rmax[planet-SE_JUPITER] = r1;
   }
  }

  for (planet = SE_JUPITER; planet <= SE_PLUTO; planet++ ) {
    printf("%10.5lf %10.5lf %10.5lf\n",lmax[planet-SE_JUPITER],
                                       bmax[planet-SE_JUPITER],
                                       rmax[planet-SE_JUPITER]);

    fprintf(daten, "\nStatistik\n%10.5lf %10.5lf %10.5lf\n",lmax[planet-SE_JUPITER],
                                       bmax[planet-SE_JUPITER],
                                       rmax[planet-SE_JUPITER]);
    }
}

// Wie weit kommt man mit einem mittlere Anomalienfile und
// den linear regredierten Bahnelementen (Routine get_elements)
// Antwort: Ganz schön weit! Bei einem 100 Tage Abstand der Daten
// aus Ephem5...Ephem9 (Differenzen Radiusvektor und mittl. Anomalie)
// bleibt man locker bei einem Fahler unter 0.01° - für alle äußteren
// Planeten
void testMeanAnomalies() {

    double x0;
    double b,r,hb,hl,hr,v,e,man,lmax[5],bmax[5],rmax[5],l1,b1,r1;
    double rd,ax_m,ex_m,inc_m,node_m,omega_m,n_m,an_m;
    double jd0,man0,rd0,jd1,man1,rd1;
    FILE *ephem, *daten;
    char filename[20], test[80];
    long offset;

    daten = fopen("daten.txt","w");

    randomize();
    randomize();
    for (planet = SE_JUPITER; planet <= SE_PLUTO; planet++) {
          sprintf(filename,"Ephem%1d", planet);
          if (planet > SE_JUPITER) fclose(ephem);
          ephem = fopen(filename,"r");
          printf("%s\n",filename);
          lmax[planet-SE_JUPITER] = 0.;
          bmax[planet-SE_JUPITER] = 0.;
          rmax[planet-SE_JUPITER] = 0.;
//          for (i = 0; i < 200; i++) {
//            jd = random(span)*increment + random(increment) + jd_start;
            for (i = 0; i < 1; i++) {
              jd = 2415021;
/* Berechnung der mittleren Bahnelemente */
            get_elements(planet, jd, &ax_m, &ex_m, &inc_m, &node_m, &omega_m, &an_m, &n_m);

/* Einlesen der mittleren Anomalien */
            x0 = (jd - jd_start)/increment;
            offset = floor(x0)*12;
            if (fseek(ephem,offset,SEEK_SET) != 0) {
              printf("Fehler bei fseek auf Position: %ld\n",offset);
              return;
              }

            if (fscanf(ephem,"%lf %lf\n",&man0,&rd0) != 2) {
              printf("Fehler bei fscanf\n");
              return;
              }

            if (fscanf(ephem,"%lf %lf\n",&man1,&rd1) != 2) {
              printf("Fehler bei fscanf\n");
              return;
              }

            jd0 = jd_start + increment*x0;
            jd1 = jd0 + increment;

            if ((jd0 > jd) || (jd > jd1)) {
              printf("Positionsierungsfehler: %lf liegt nicht in [%lf,%lf]\n",
                jd,jd0,jd1);
              return;
              }

            man = man0 + (jd-jd0)/(jd1-jd0)*(man1-man0);
            man /= 1000.;
            rd  = rd0  + (jd-jd0)/(jd1-jd0)*(rd1 -rd0 );

            kepler1(ex_m,ax_m,omega_m,node_m,inc_m,man+an_m,&l1,&b1,&r1);

            printf("ex_m:\t%lf\n",ex_m);
            printf("ax_m:\t%lf\n",ax_m);
            printf("omega_m:\t%lf\n",omega_m);
            printf("node_m:\t%lf\n",node_m);
            printf("inc_m:\t%lf\n",inc_m);
            printf("an_m:\t%lf\n",an_m);
            printf("man\t%lf\n", man);
            printf("rd\t%lf\n", rd);

/* Berechnung der exakten heliozentrischen Position */
            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];


/*  Ausdrucken : M.An. und Radiusdifferenz */
            fprintf(daten,"%10.1lf %12.7lf %12.7lf %12.7lf\n",jd,arcdiff(hl,l1),hl,l1);

/* Jetzt Vergleichen !! */
           l1 = fabs(arcdiff(hl,l1));
           if (l1>lmax[planet-SE_JUPITER]) lmax[planet-SE_JUPITER] = l1;
           b = hb - b1;
           if (fabs(b)>bmax[planet-SE_JUPITER]) bmax[planet-SE_JUPITER]=fabs(b);
           r1=fabs(hr-r1);
           if (r1>rmax[planet-SE_JUPITER]) rmax[planet-SE_JUPITER] = r1;
   }
  }

  for (planet = SE_JUPITER; planet <= SE_PLUTO; planet++ )
    printf("%10.5lf %10.5lf %10.5lf\n",lmax[planet-SE_JUPITER],
                                       bmax[planet-SE_JUPITER],
                                       rmax[planet-SE_JUPITER]);

  }

// Korrekturen für Radiusvektor (in Zehntausendsteln) und
// Mittlere Anomalie (in Tausendsteln) berechnen und in Ephem5...Ephem9 ablegen
void computeMeanAnomalies() {
    double c = CONV;
    double b,l,r,hb,hl,hr,v,e,man,lmax[5],bmax[5],rmax[5],l1,b1,r1;
    double ax_m,ex_m,inc_m,node_m,omega_m,n_m,an_m;
    FILE *ephem;
    char filename[20];

    double an_max, r_max, an_min, r_min, dan, dr;
    bool firstRun;

    for (planet = SE_JUPITER; planet <= SE_PLUTO; planet++) {
          sprintf(filename,"Ephem%1d", planet);
          if (planet > SE_JUPITER) fclose(ephem);
          ephem = fopen(filename,"w");
          printf("%s\n",filename);
          lmax[planet-SE_JUPITER] = 0.;
          bmax[planet-SE_JUPITER] = 0.;
          rmax[planet-SE_JUPITER] = 0.;
          firstRun = true;
          for (jd = jd_start; jd <= jd_start + ( span + 2 )* increment; jd += increment) {
/* 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];
/* x/y/z-Koordinaten für Ort und Geschwindigkeit berechnen */
            se_flag = SEFLG_HELCTR + SEFLG_NONUT + SEFLG_SPEED + SEFLG_XYZ;
            rc = swe_calc(jd, planet, se_flag, x, msg);
/* Umrechnung kartesische Koordinaten für Ort und Geschwindigkeit in osk. Bahnelemente */
            car_orbit(x[0],x[1],x[2],x[3],x[4],x[5],
                      rm[planet-4], &ax, &ex, &inc, &node, &omega, &n, &an);

/* 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);

// Für Minimax-Auswertung:
            if (firstRun) {
              an_max = dan;
              an_min = dan;
              r_max  = dr;
              r_min  = dr;
              firstRun = false;
              }
            else {
              if (dan > an_max) an_max = dan;
              if (dan < an_min) an_min = dan;
              if (dr  > r_max)  r_max  = dr;
              if (dr  < r_min)  r_min  = dr;
              }

/* Jetzt Vergleichen !! */
           l1 = fabs(arcdiff(hl,l1));
           if (l1>lmax[planet-SE_JUPITER]) lmax[planet-SE_JUPITER] = l1;
           b = hb - b1;
           if (fabs(b)>bmax[planet-SE_JUPITER]) bmax[planet-SE_JUPITER]=fabs(b);
           r1=fabs(hr-r1);
           if (r1>rmax[planet-SE_JUPITER]) rmax[planet-SE_JUPITER] = r1;
     }

     printf( "%3d: %8.2lf  bis  %8.2lf -> %8.2lf\n", planet, an_min, an_max, an_max - an_min);
     printf( "     %8.2lf  bis  %8.2lf -> %8.2lf\n",         r_min,   r_max,  r_max -  r_min);
  }

  fclose(ephem);

  for (planet = SE_JUPITER; planet <= SE_PLUTO; planet++ )
    printf("%10.5lf %10.5lf %10.5lf\n",lmax[planet-SE_JUPITER],
                                       bmax[planet-SE_JUPITER],
                                       rmax[planet-SE_JUPITER]);
}

// Die durch Regresseion aus den osk. Bahnelementen gewonnen
// mittleren Bahnelemente werden gegen die tatsächlichen Elemente
// geprüft, wie sie sich aus der Swiss Ephemeris ergeben
void testMeanElements() {

    FILE *daten = fopen("daten.txt", "w");
    double ax1,ex1,inc1,node1,omega1,n1,an1;
    char plname[50];

// Heliocentric without nutation
        se_flag = SEFLG_HELCTR + SEFLG_NONUT + SEFLG_SPEED + SEFLG_XYZ;
        for (planet = SE_JUPITER; planet <= SE_PLUTO; planet++ ) {
          swe_get_planet_name(planet, plname);
          fprintf(daten,"%s\n",plname);
          for (jd = jd_start; jd <= jd_start + span * increment + 1; jd += increment) {
            rc = swe_calc(jd, planet, se_flag, x, msg);
            car_orbit(x[0],x[1],x[2],x[3],x[4],x[5],
                      rm[planet-4], &ax, &ex, &inc, &node, &omega, &n, &an);
            get_elements(planet, jd, &ax1, &ex1, &inc1, &node1, &omega1, &an1, &n1);
            fprintf(daten,"%10.1lf %9.6lf %8.6lf %7.4lf %8.4lf %8.4lf %8.4lf %9.7lf\n",
                  jd,
                  fabs(ax-ax1),
                  fabs(ex-ex1),
                  fabs(arcdiff(inc,inc1)),
                  fabs(arcdiff(node,node1)),
                  fabs(arcdiff(omega,omega1)),
                  fabs(arcdiff(an,an1)),
                  fabs(n-n1));
          }
          }

 }

// Berechne oskulierende Bahnelemente für die äußeren Planeten
// Grundlage sind die präzisen Positionen und Geschwindigkeiten nach
// Swiss Ephemeris
// Ergebnisse werden abgespeichert auf Elem5 bis Elem9.
void computeOsculatingElements() {

       FILE *elem;
       char filename[20];

// Heliozentrisch ohne Nutation
        se_flag = SEFLG_HELCTR + SEFLG_NONUT + SEFLG_SPEED + SEFLG_XYZ;

        for (planet = SE_MERCURY; planet <= SE_PLUTO; planet++ ) {
          sprintf(filename,"Elem%1d", planet);
          if (planet > SE_MERCURY) fclose(elem);
          elem = fopen(filename,"w");
          printf("%s\n",filename);
          for (jd = jd_start; jd < jd_start + ( span + 2 ) * increment + 1;
               jd += (planet > SE_MARS) ? increment : increment / 10) {
            rc = swe_calc(jd, planet, se_flag, x, msg);
            // cartesische Koordinaten (Ort + Geschwindigkeit)
            // in oskulierende Bahnelemente umwandeln
            car_orbit(x[0],x[1],x[2],x[3],x[4],x[5],
                      rm[planet-4], &ax, &ex, &inc, &node, &omega, &n, &an);
            fprintf(elem,"%10.1lf %9.6lf %8.6lf %7.4lf %8.4lf %8.4lf %8.4lf %9.7lf\n",
                  jd,ax,ex,inc,node,omega,an,n);
            }
          }
       fclose(elem);

 }

// ----------------------------------------------------------------------
// Ein Planet wird mit hochgenauer Geschwindigkeit mittels Swiss Ephemeris
// berechnet (cartesisch, ergibt x[0]...x[5])
// Ort und Geschwindigkeit werden unter Verwendung der reduzierten Massen
// (rm[]) in Bahnelemente umgerechnet
// Schließlich wird im Intervall [jd-1000,jd+1000] die Position des
// Planeten einerseits aus Swiss Ephemeris, andererseits mittels der
// Bahnelemente zum Zeitpunkt jd berechnet. Die Abweichungen in Länge,
// Breite und Radiusvektor werden auf das File daten.txt geschrieben.
// Ergebnisse: Jupiter und Saturn reproduzieren die Pos. noch +/- 500d
// auf 0.01 Grad. Für Uranus, Neptun und Pluto sind die Zahlen noch weit
// besser (+/- 1000 Tage reichen noch für 0.002°!).
void testOrbitConversion() {

daten = fopen("daten.txt","w");

        for (planet = SE_JUPITER; planet <= SE_PLUTO; planet++) {
          fprintf(daten,"Planet: %3d\n",planet);
// Heliocentric without nutation
        se_flag = SEFLG_HELCTR + SEFLG_NONUT + SEFLG_SPEED + SEFLG_XYZ;
        rc = swe_calc(jd_start, planet, se_flag, x, msg);
//        if ((rc < 0)||(rc != se_flag)) printf("rc = %d\nse = %d\n%s\n",rc,se_flag, msg);
        car_orbit(x[0],x[1],x[2],x[3],x[4],x[5],
                  rm[planet-4], &ax, &ex, &inc, &node, &omega, &n, &an);

        for (vorlauf = -1000; vorlauf < 1001; vorlauf++) {
          for (i=0;i<6;i++) x[i] = 0.;
          kepler1 (ex,ax,omega,node,inc,an+vorlauf*n,
                   &x[0],&x[1],&x[2]);
          se_flag = SEFLG_HELCTR + SEFLG_NONUT + SEFLG_SPEED;
          rc = swe_calc(jd_start+vorlauf, planet, se_flag, y, msg);
 //         if ((rc < 0)||(rc != se_flag)) printf("rc = %d\nse = %d\n%s\n",rc,se_flag, msg);
          fprintf(daten,"%5d",vorlauf);
          for (i=0; i<3; i++) fprintf(daten,"%15.10lf",arcdiff(x[i],y[i]));
          fprintf(daten,"\n");
          }

        }

}
//---------------------------------------------------------------------------
// Unterprogramm zur Berechnung mittlerer Bahnelemente, die zuvor durch
// lineare Regression (mittels eines MATHEMATICA Notebooks) aus Listen der
// oskulierenden Bahnelemente gewonnen waren.
void get_elements(int iplanet, double jd, double *ax, double *ex, double *inc,
             double *node, double *omega, double *an, double *n)
{

static bool dataRead = false;

double t = ( jd - jd_start ) / 100;
int ipl = iplanet - SE_MERCURY;

if (!dataRead) {
  getElementCoefficients();
  dataRead = true;
}


  *ax    = a[ipl][0][0] + ( a[ipl][0][1] + a[ipl][0][2] * t ) * t;
  *ex    = a[ipl][1][0] + ( a[ipl][1][1] + a[ipl][1][2] * t ) * t;
  *inc   = a[ipl][2][0] + ( a[ipl][2][1] + a[ipl][2][2] * t ) * t;
  *node  = a[ipl][3][0] + ( a[ipl][3][1] + a[ipl][3][2] * t ) * t;
  *omega = a[ipl][4][0] + ( a[ipl][4][1] + a[ipl][4][2] * t ) * t;
  *an    = a[ipl][5][0] + ( a[ipl][5][1] + a[ipl][5][2] * t ) * t;
  *n     = a[ipl][6][0] + ( a[ipl][6][1] + a[ipl][6][2] * t ) * t;

  *an = frac(*an);

  return; // Hier reicht's

}

void getElementCoefficients() {
  char filename[20];
  FILE *coefficients;

  sprintf(filename,"c%04d\.DAT", jahr0);
  coefficients = fopen(filename,"r");
  for (int iplanet = SE_MERCURY; iplanet <= SE_PLUTO; iplanet++)
    for (i = 0; i < 7; i++)
      if (iplanet <= SE_MARS) {
        fscanf(coefficients,"%lf %lf", &(a[iplanet-SE_MERCURY][i][0]),
                                       &(a[iplanet-SE_MERCURY][i][1]));
        a[iplanet-SE_MERCURY][i][2] = 0.;
        }
      else
        fscanf(coefficients,"%lf %lf %lf", &(a[iplanet-SE_MERCURY][i][0]),
                                           &(a[iplanet-SE_MERCURY][i][1]),
                                           &(a[iplanet-SE_MERCURY][i][2]));
  fclose(coefficients);
}

void compressMeanAnomalies() {

    FILE *ephem, *compress;
    char filename[20], lsb, msb;
    int max[2][5], min[2][5], interval[2][5], size[2][5], mask[2][5] ;
    int x[2][5][1600],j, k, bitzeiger;
    unsigned long buf;

// Zum Wegschreiben im Dateiheader
    unsigned long jdstart = (int) jd_start,
                  jdend = (int) ( jd_start + span * increment ) ;

// Einlesen
    for (planet = SE_JUPITER; planet <= SE_PLUTO; planet++ ) {
      sprintf(filename,"Ephem%d",planet);
      ephem = fopen(filename,"r");
      for (i = 0; i <= span + 2; i++) {
        fscanf(ephem,"%d %d\n", &x[0][planet-SE_JUPITER][i],
                                &x[1][planet-SE_JUPITER][i]);

        if (i == 0) {
          min[0][planet-SE_JUPITER] = x[0][planet-SE_JUPITER][i];
          min[1][planet-SE_JUPITER] = x[1][planet-SE_JUPITER][i];
          max[0][planet-SE_JUPITER] = x[0][planet-SE_JUPITER][i];
          max[1][planet-SE_JUPITER] = x[1][planet-SE_JUPITER][i];
          }
        else {
          if (min[0][planet-SE_JUPITER] > x[0][planet-SE_JUPITER][i])
              min[0][planet-SE_JUPITER] = x[0][planet-SE_JUPITER][i];
          if (min[1][planet-SE_JUPITER] > x[1][planet-SE_JUPITER][i])
              min[1][planet-SE_JUPITER] = x[1][planet-SE_JUPITER][i];
          if (max[0][planet-SE_JUPITER] < x[0][planet-SE_JUPITER][i])
              max[0][planet-SE_JUPITER] = x[0][planet-SE_JUPITER][i];
          if (max[1][planet-SE_JUPITER] < x[1][planet-SE_JUPITER][i])
              max[1][planet-SE_JUPITER] = x[1][planet-SE_JUPITER][i];
          }
        }
      fclose(ephem);
    }

// Auch die mittleren Bahnelemente
   getElementCoefficients();

// Daten sind jetzt da --> Comprimieren und wechschreiben
    sprintf(filename,"J%04d.dat",jahr0);
    compress = fopen(filename,"wb");
// Zunächst die Giltigkeit (Anfang, Ende, Inkrement)
    fprintf(compress,"%c%c%c", ( jdstart >> 16 ) & 0xFF,
                                 ( jdstart >> 8 ) & 0xFF, jdstart & 0xFF );
    fprintf(compress,"%c%c%c", ( jdend >> 16 ) & 0xFF,
                                 ( jdend >> 8 ) & 0xFF, jdend & 0xFF );
    fprintf(compress,"%c%c", increment >> 8, increment & 0xFF);

// Dann die Koeffizienten
   for (planet = SE_MERCURY; planet <= SE_PLUTO ; planet++ )
     for (i = 0; i < 7; i++)
       for (j = 0; j < ((planet <= SE_MARS) ? 2:3); j++)
         putDouble( compress, a[planet-SE_MERCURY][i][j] );

// Sodann die Minima und Maxima
    for (i = 0; i <= (SE_PLUTO - SE_JUPITER); i++)
      for (j = 0; j<2 ; j++) {
       fprintf(compress,"%c%c",(min[j][i] & 0xFF00) / 256, min[j][i] & 255 );
       fprintf(compress,"%c%c",(max[j][i] & 0xFF00) / 256, max[j][i] & 255 );
       interval[j][i] = max[j][i] - min[j][i];
       size[j][i] =  log(interval[j][i])/log(2) + 1;
       mask[j][i] =  ( 1 << size[j][i] ) - 1;
       printf("%d %d %d %d\n", min[j][i], max[j][i], interval[j][i], size[j][i]);
       }

// 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 );
    fclose(compress);
}