Advertisement
rplantiko

Ostervollmonde und kirchliches Ostern

Dec 15th, 2019 (edited)
398
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include <stdbool.h>
  4. #include "swephexp.h" // Swiss Ephemeris für die astronomischen Sonnen- und Mondstände
  5.  
  6. // Berechnet für einen Zeitraum (mit "Jahr von" und "Jahr bis" als Kommandozeilenparameter angegeben)
  7. // 1.) die auf den 1. März folgenden drei Mondzyklen (Voll- und Neomonde)
  8. // 2.) das "astronomische Ostern" (den auf den ersten Frühlingsvollmond folgenden Sonntag)
  9. // 3.) das "kirchliche Ostern" gemäß dem kirchlichen "mittleren Neumond" bzw. "Luna XIV",
  10. //     falls es vom "astronomischen Ostern" abweicht (Basis: gregorianischer Kalender)
  11.  
  12. // Details: http://ruediger-plantiko.blogspot.com/2019/12/ostertermin-und-osterparadoxa.html
  13.  
  14. // Datentyp "taylor1" = Zahlenpaar für Funktionswert und 1. Ableitung
  15. typedef struct {
  16.   double y, yp;
  17.   } taylor1;
  18.  
  19. // Beliebige Funktion, wie sie das Newtonverfahren braucht (mit 1. Ableitung als Ergebnis)
  20. typedef taylor1 function( double x );
  21.  
  22. // Entscheidung Gregorianischer/Julianischer Kalender aus Julianischer Tageszahl
  23. int is_jd_greg( double jd ) {
  24.   return (jd >= 2299161.30602) ? 1 : 0;
  25. };
  26.  
  27. // (Grobe) Entscheidung Gregorianischer/Julianischer Kalender aus Jahreszahl
  28. // Wird in der Gaußschen Osterformel verwendet
  29. int is_greg( int year ) {
  30.   return (year >= 1583) ? 1 : 0;
  31. }
  32.  
  33. // Differenz zweier Winkelwerte, reduziert auf -180°..180°
  34. double arcdiff( double x , double y) {
  35.   double d = x - y;
  36. // Erster Versuch: schnell
  37.   if (d>-180 && d<180) {
  38.     return d;
  39.   }
  40. // Allgemein: reduzieren
  41.   d = d/.36e3 +.5;
  42.   return (d - floor (d ))*.36e3 -.18e3 ;
  43. }
  44.  
  45. // Newton-Verfahren zur Bestimmung von Nullstellen
  46. double newton( double x0, function *f) {
  47.  
  48.   const double PREC = 1.e-7; // Genauigkeit
  49.   const int MAXITER = 10000; // Theoretische Obergrenze
  50.  
  51.   taylor1 t;
  52.   int iter = 0;  
  53.   double x = x0;
  54.  
  55.   for (;;) {
  56.  
  57. // Funktion evaluieren
  58.     t = f(x);    
  59.  
  60. // Nullstelle gefunden?
  61.     if (fabs(t.y) < PREC) {
  62.       return x; // Ja
  63.     }
  64.  
  65. // Zur Sicherheit, um Endlosschleifen zu verhindern
  66. // (bei falschen Funktionsdefinitionen, oder bei mangelnder Konvergenz)
  67.     if (++iter > MAXITER) {
  68.       fprintf(stderr,"Maximum number of iterations exceeded, something went wrong\n");
  69.       return 0;
  70.     }
  71.  
  72. // Division durch Null verhindern
  73.     if (t.yp == 0) {
  74.       fprintf(stderr,"f'(%lf) = 0, Newton algorithm has a problem\n",x);
  75.       return 0;
  76.     }
  77.  
  78. // Next iteration
  79.     x = x - t.y / t.yp;
  80.    
  81.   }
  82.  
  83. }
  84.  
  85. // Gregorianisches Kalenderdatum aus JD
  86. char* caldate( double jd, char* ds) {
  87.   int year, month, day;
  88.   double hour;
  89.   swe_revjul(jd, is_jd_greg( jd), &year, &month, &day, &hour);
  90.   sprintf(ds,"%2d.%2d.%4d %6.4lf",day,month,year,hour);
  91.   return ds;
  92. }
  93.  
  94. // Der nachfolgende Sonntag nach einem Ereignis
  95. double sunday_after( double jd ) {
  96.     int jdn = (int) (jd + 0.5); // Julianische Tageszahl, bei 0hET beginnend
  97.     int weekday = (jdn+1)%7; // Sonntag = 0
  98. // Auch wenn der Termin selbst bereits auf einen Sonntag fällt)
  99.     return jdn + (7-weekday)-.5;
  100. }
  101.  
  102. // Ort (Länge) und Geschwindigkeit eines Planeten berechnen
  103. taylor1 swecalc(int planet, double jd_et) {
  104.     double xx[6];
  105.     int iflag = SEFLG_SWIEPH | SEFLG_SPEED;
  106.     char serr[256];
  107.     serr[0] = '\0';
  108.     int rc = swe_calc( jd_et, planet, iflag, xx, serr);
  109.     if (serr[0] || (rc < 0)) {
  110.       fprintf(stderr,"swe_calc problem: %s\n",serr);
  111.       return (taylor1) {0.,0.};
  112.     }
  113.     return (taylor1) { xx[0], xx[3] };  
  114. }
  115.  
  116. // Frühlingsanfang berechnen
  117. double get_spring( int year ) {
  118.   taylor1 sun_length(double jd_et) {
  119.     taylor1 sun = swecalc(SE_SUN,jd_et);
  120.     return (taylor1) { arcdiff(sun.y,0), sun.yp };
  121.   }
  122.   double jd0 = swe_julday( year, 3, 20, 0., is_greg( year ) );
  123.   return newton(jd0, sun_length);
  124. }
  125.  
  126. // Nächste gewünschte Mondphase nach jd0 berechnen
  127. const bool NEW_MOON = true;
  128. const bool FULL_MOON = false;
  129. const double LUNAR_MONTH = 29.530587981; // aus https://en.wikipedia.org/wiki/Lunar_month
  130. double get_next_syzygy(double jd0,bool new) {
  131.   taylor1 phase(double jd_et) {
  132.     taylor1 sun = swecalc(SE_SUN,jd_et);
  133.     taylor1 moon = swecalc(SE_MOON,jd_et);
  134.     double d = arcdiff(moon.y + (new ? 0. : 180.),sun.y);
  135.     return (taylor1) { d, moon.yp - sun.yp };
  136.   }
  137.  
  138.   taylor1 phase0 = phase(jd0);
  139.   double jd1 = jd0 - phase0.y * LUNAR_MONTH / 360;
  140.  
  141. // Weil das nächste Ereignis nach jd0 gewünscht ist:
  142.   if (phase0.y > 0) jd1 += LUNAR_MONTH;
  143.  
  144.   return newton(jd1,phase);
  145.  
  146. }
  147.  
  148. double easter_date_greg( int year ) {
  149.  
  150. // Gaußsche Osterformel, gregorianische und julianische Version
  151. // siehe https://de.wikipedia.org/wiki/Gau%C3%9Fsche_Osterformel
  152.   int
  153.    a = year % 19,
  154.    b = year % 4,
  155.    c = year % 7,
  156.    M,N;
  157.   if (is_greg(year)) {
  158.    int
  159.      k = (int) (year/100),
  160.      p = (int) ((13+8*k)/25),
  161.      q = (int) (k/4);
  162.      M = (15 - p + k - q) % 30,
  163.      N = (4 + k - q) % 7;
  164.   } else {
  165.      M = 15;
  166.      N = 6;
  167.   }
  168.   int  
  169.    d = (19*a + M) % 30,
  170.    e = (2*b + 4*c + 6*d + N) % 7,
  171.    day,month;  
  172. // Zwei Ausnahmen gegen Ende der erlaubten Ostertermine
  173.    if (is_greg(year) && d == 29 && e == 6) {
  174.      day = 19;
  175.      month = 4;
  176.    } else if (is_greg(year) && d == 28 && e == 6 && ((11*M + 11) % 30) < 19) {
  177.      day = 18;
  178.      month = 4;
  179.    } else if (d+e<=9) {
  180. // März-Termine
  181.      day = 22 + d + e;
  182.      month = 3;
  183.    } else {
  184. // April-Termine
  185.      day = d + e - 9;
  186.      month = 4;
  187.    }
  188.  
  189.    return swe_julday(year,month,day,0.,is_greg(year));
  190. }
  191.  
  192. // Hauptprogramm
  193. int main(int argc, char** argv) {
  194.  
  195. // Wieviele Syzygien ab 1. März berechnet werden sollen
  196.     const int NUMSYZ = 3;
  197.  
  198. // Zwei Aufrufparameter, für Beginn- und Endjahr der Berechnung
  199.   int year0 = 33, year1 = 2500;
  200.   if (argc >= 2) {
  201.     if (sscanf(argv[1],"%d",&year0)!=1) {
  202.       fprintf(stderr,"Argument '%s' ist keine ganze Zahl\n",argv[1]);
  203.       exit(EXIT_FAILURE);
  204.     }
  205.     year1 = year0;
  206.   }
  207.   if (argc >= 3) {
  208.     if (sscanf(argv[2],"%d",&year1)!=1) {
  209.       fprintf(stderr,"Argument '%s' ist keine ganze Zahl\n",argv[2]);
  210.       exit(EXIT_FAILURE);
  211.     }
  212.   }
  213.  
  214.  
  215. // Überschrift
  216.     printf( "%-20s ", "Frühlingsanfang" );
  217.     for (int i=0;i<NUMSYZ;i++) {
  218.       printf( "%-20s", "Neumond" );
  219.       printf( "%-20s", "Vollmond" );
  220.     }
  221.     printf( "%-12s", "A. Ostern" );
  222.     printf( "%-12s", "K. Ostern" );
  223.     printf("\n");
  224.  
  225.  
  226.   for (int year = year0; year<=year1; year++) {
  227.     typedef struct {
  228.        double new, full;
  229.     } syzygy;
  230.  
  231. // Frühlingsanfang berechnen
  232.     double spring = get_spring( year );
  233.  
  234. // Syzygien berechnen
  235.     syzygy s[NUMSYZ];
  236.     double spring_moon = 0.,
  237.            easter_astro = 0.,
  238.            church_moon = 0.,
  239.            easter_church = 0.;
  240.  
  241. // Startdatum 1. März greg. 0:00 ET
  242.     double jd = swe_julday( year, 3, 1, 0., is_greg( year ) );
  243.     for (int j=0;j<2*NUMSYZ;j++) {
  244.       int i = (int) j / 2;
  245.       bool new = (j % 2 == 0);
  246.       jd = get_next_syzygy( jd, new );
  247.       if (new) s[i].new = jd;
  248.       else {
  249.         s[i].full = jd;    
  250. // Astronomischen Frühlingsvollmond merken
  251.         if ((jd >= spring) && (spring_moon == 0)) {
  252.           spring_moon = jd;
  253.         }
  254.       }
  255.     }
  256.     easter_astro = sunday_after( spring_moon );
  257.    
  258.  
  259. // Astronomischen Ostertermin ermitteln
  260.     int jdn = (int) (spring_moon + 0.5); // Julianische Tageszahl, bei 0hET beginnend
  261.     int weekday = (jdn+1)%7; // Sonntag = 0
  262. // Immer den auf den Frühlingsvollmond folgenden Sonntag
  263. // (Auch wenn der Frühlingsvollmond selbst bereits auf einen Sonntag fällt)
  264.     double easter = jdn + (7-weekday)-.5;
  265.    
  266. // Kirchlicher Ostertermin
  267.     easter_church = easter_date_greg( year );
  268.  
  269. // Ausgabe
  270.     char ds[30],ds1[30];
  271.     printf( "%-20s", caldate(spring,ds) );
  272.     for (int i=0;i<NUMSYZ;i++) {
  273.       printf( "%-20s", caldate(s[i].new,ds) );
  274.       printf( "%-20s", caldate(s[i].full,ds) );
  275.     }
  276.     caldate(easter_astro,ds);
  277.     *(ds+11)='\0';
  278.     printf("%-12s",ds);
  279.     caldate(easter_church,ds1);
  280.     *(ds1+11)='\0';
  281.     // Ausgeben, falls abweichend
  282.     if (strcmp(ds,ds1)!=0) {
  283.       printf("%-12s",ds1);
  284.     }
  285.     printf("\n");
  286.   }
  287.  
  288.   return 0;
  289. }
Advertisement
RAW Paste Data Copied
Advertisement