Ostervollmonde und kirchliches Ostern

#include <stdio.h>
#include <math.h>
#include <stdbool.h>
#include "swephexp.h" // Swiss Ephemeris für die astronomischen Sonnen- und Mondstände

// Berechnet für einen Zeitraum (mit "Jahr von" und "Jahr bis" als Kommandozeilenparameter angegeben)
// 1.) die auf den 1. März folgenden drei Mondzyklen (Voll- und Neomonde)
// 2.) das "astronomische Ostern" (den auf den ersten Frühlingsvollmond folgenden Sonntag)
// 3.) das "kirchliche Ostern" gemäß dem kirchlichen "mittleren Neumond" bzw. "Luna XIV",
//     falls es vom "astronomischen Ostern" abweicht (Basis: gregorianischer Kalender)

// Details: http://ruediger-plantiko.blogspot.com/2019/12/ostertermin-und-osterparadoxa.html

// Datentyp "taylor1" = Zahlenpaar für Funktionswert und 1. Ableitung
typedef struct {
  double y, yp;
  } taylor1;

// Beliebige Funktion, wie sie das Newtonverfahren braucht (mit 1. Ableitung als Ergebnis)
typedef taylor1 function( double x );

// Entscheidung Gregorianischer/Julianischer Kalender aus Julianischer Tageszahl
int is_jd_greg( double jd ) {
  return (jd >= 2299161.30602) ? 1 : 0;
};

// (Grobe) Entscheidung Gregorianischer/Julianischer Kalender aus Jahreszahl
// Wird in der Gaußschen Osterformel verwendet
int is_greg( int year ) {
  return (year >= 1583) ? 1 : 0;
}

// Differenz zweier Winkelwerte, reduziert auf -180°..180°
double arcdiff( double x , double y) {
  double d = x - y;
// Erster Versuch: schnell
  if (d>-180 && d<180) {
    return d;
  }
// Allgemein: reduzieren
  d = d/.36e3 +.5;
  return (d - floor (d ))*.36e3 -.18e3 ;
}

// Newton-Verfahren zur Bestimmung von Nullstellen
double newton( double x0, function *f) {

  const double PREC = 1.e-7; // Genauigkeit
  const int MAXITER = 10000; // Theoretische Obergrenze

  taylor1 t;
  int iter = 0;
  double x = x0;

  for (;;) {

// Funktion evaluieren
    t = f(x);

// Nullstelle gefunden?
    if (fabs(t.y) < PREC) {
      return x; // Ja
    }

// Zur Sicherheit, um Endlosschleifen zu verhindern
// (bei falschen Funktionsdefinitionen, oder bei mangelnder Konvergenz)
    if (++iter > MAXITER) {
      fprintf(stderr,"Maximum number of iterations exceeded, something went wrong\n");
      return 0;
    }

// Division durch Null verhindern
    if (t.yp == 0) {
      fprintf(stderr,"f'(%lf) = 0, Newton algorithm has a problem\n",x);
      return 0;
    }

// Next iteration
    x = x - t.y / t.yp;

  }

}

// Gregorianisches Kalenderdatum aus JD
char* caldate( double jd, char* ds) {
  int year, month, day;
  double hour;
  swe_revjul(jd, is_jd_greg( jd), &year, &month, &day, &hour);
  sprintf(ds,"%2d.%2d.%4d %6.4lf",day,month,year,hour);
  return ds;
}

// Der nachfolgende Sonntag nach einem Ereignis
double sunday_after( double jd ) {
    int jdn = (int) (jd + 0.5); // Julianische Tageszahl, bei 0hET beginnend
    int weekday = (jdn+1)%7; // Sonntag = 0
// Auch wenn der Termin selbst bereits auf einen Sonntag fällt)
    return jdn + (7-weekday)-.5;
}

// Ort (Länge) und Geschwindigkeit eines Planeten berechnen
taylor1 swecalc(int planet, double jd_et) {
    double xx[6];
    int iflag = SEFLG_SWIEPH | SEFLG_SPEED;
    char serr[256];
    serr[0] = '\0';
    int rc = swe_calc( jd_et, planet, iflag, xx, serr);
    if (serr[0] || (rc < 0)) {
      fprintf(stderr,"swe_calc problem: %s\n",serr);
      return (taylor1) {0.,0.};
    }
    return (taylor1) { xx[0], xx[3] };
}

// Frühlingsanfang berechnen
double get_spring( int year ) {
  taylor1 sun_length(double jd_et) {
    taylor1 sun = swecalc(SE_SUN,jd_et);
    return (taylor1) { arcdiff(sun.y,0), sun.yp };
  }
  double jd0 = swe_julday( year, 3, 20, 0., is_greg( year ) );
  return newton(jd0, sun_length);
}

// Nächste gewünschte Mondphase nach jd0 berechnen
const bool NEW_MOON = true;
const bool FULL_MOON = false;
const double LUNAR_MONTH = 29.530587981; // aus https://en.wikipedia.org/wiki/Lunar_month
double get_next_syzygy(double jd0,bool new) {
  taylor1 phase(double jd_et) {
    taylor1 sun = swecalc(SE_SUN,jd_et);
    taylor1 moon = swecalc(SE_MOON,jd_et);
    double d = arcdiff(moon.y + (new ? 0. : 180.),sun.y);
    return (taylor1) { d, moon.yp - sun.yp };
  }

  taylor1 phase0 = phase(jd0);
  double jd1 = jd0 - phase0.y * LUNAR_MONTH / 360;

// Weil das nächste Ereignis nach jd0 gewünscht ist:
  if (phase0.y > 0) jd1 += LUNAR_MONTH;

  return newton(jd1,phase);

}

double easter_date_greg( int year ) {

// Gaußsche Osterformel, gregorianische und julianische Version
// siehe https://de.wikipedia.org/wiki/Gau%C3%9Fsche_Osterformel
  int
   a = year % 19,
   b = year % 4,
   c = year % 7,
   M,N;
  if (is_greg(year)) {
   int
     k = (int) (year/100),
     p = (int) ((13+8*k)/25),
     q = (int) (k/4);
     M = (15 - p + k - q) % 30,
     N = (4 + k - q) % 7;
  } else {
     M = 15;
     N = 6;
  }
  int
   d = (19*a + M) % 30,
   e = (2*b + 4*c + 6*d + N) % 7,
   day,month;
// Zwei Ausnahmen gegen Ende der erlaubten Ostertermine
   if (is_greg(year) && d == 29 && e == 6) {
     day = 19;
     month = 4;
   } else if (is_greg(year) && d == 28 && e == 6 && ((11*M + 11) % 30) < 19) {
     day = 18;
     month = 4;
   } else if (d+e<=9) {
// März-Termine
     day = 22 + d + e;
     month = 3;
   } else {
// April-Termine
     day = d + e - 9;
     month = 4;
   }

   return swe_julday(year,month,day,0.,is_greg(year));
}

// Hauptprogramm
int main(int argc, char** argv) {

// Wieviele Syzygien ab 1. März berechnet werden sollen
    const int NUMSYZ = 3;

// Zwei Aufrufparameter, für Beginn- und Endjahr der Berechnung
  int year0 = 33, year1 = 2500;
  if (argc >= 2) {
    if (sscanf(argv[1],"%d",&year0)!=1) {
      fprintf(stderr,"Argument '%s' ist keine ganze Zahl\n",argv[1]);
      exit(EXIT_FAILURE);
    }
    year1 = year0;
  }
  if (argc >= 3) {
    if (sscanf(argv[2],"%d",&year1)!=1) {
      fprintf(stderr,"Argument '%s' ist keine ganze Zahl\n",argv[2]);
      exit(EXIT_FAILURE);
    }
  }


// Überschrift
    printf( "%-20s ", "Frühlingsanfang" );
    for (int i=0;i<NUMSYZ;i++) {
      printf( "%-20s", "Neumond" );
      printf( "%-20s", "Vollmond" );
    }
    printf( "%-12s", "A. Ostern" );
    printf( "%-12s", "K. Ostern" );
    printf("\n");


  for (int year = year0; year<=year1; year++) {
    typedef struct {
       double new, full;
    } syzygy;

// Frühlingsanfang berechnen
    double spring = get_spring( year );

// Syzygien berechnen
    syzygy s[NUMSYZ];
    double spring_moon = 0.,
           easter_astro = 0.,
           church_moon = 0.,
           easter_church = 0.;

// Startdatum 1. März greg. 0:00 ET
    double jd = swe_julday( year, 3, 1, 0., is_greg( year ) );
    for (int j=0;j<2*NUMSYZ;j++) {
      int i = (int) j / 2;
      bool new = (j % 2 == 0);
      jd = get_next_syzygy( jd, new );
      if (new) s[i].new = jd;
      else {
        s[i].full = jd;
// Astronomischen Frühlingsvollmond merken
        if ((jd >= spring) && (spring_moon == 0)) {
          spring_moon = jd;
        }
      }
    }
    easter_astro = sunday_after( spring_moon );


// Astronomischen Ostertermin ermitteln
    int jdn = (int) (spring_moon + 0.5); // Julianische Tageszahl, bei 0hET beginnend
    int weekday = (jdn+1)%7; // Sonntag = 0
// Immer den auf den Frühlingsvollmond folgenden Sonntag
// (Auch wenn der Frühlingsvollmond selbst bereits auf einen Sonntag fällt)
    double easter = jdn + (7-weekday)-.5;

// Kirchlicher Ostertermin
    easter_church = easter_date_greg( year );

// Ausgabe
    char ds[30],ds1[30];
    printf( "%-20s", caldate(spring,ds) );
    for (int i=0;i<NUMSYZ;i++) {
      printf( "%-20s", caldate(s[i].new,ds) );
      printf( "%-20s", caldate(s[i].full,ds) );
    }
    caldate(easter_astro,ds);
    *(ds+11)='\0';
    printf("%-12s",ds);
    caldate(easter_church,ds1);
    *(ds1+11)='\0';
    // Ausgeben, falls abweichend
    if (strcmp(ds,ds1)!=0) {
      printf("%-12s",ds1);
    }
    printf("\n");
  }

  return 0;
}