Bestimmung der Zeitpunkte von Ereignissen
ulrich
2023-03-23 d20d989f5495492f1258c8313db7c19b429111a3
src/de/uhilger/zeitrechnung/kalender/BasisKalender.java
@@ -22,10 +22,12 @@
import de.uhilger.zeitrechnung.Zeit;
/**
 * Abstrakte Basisklasse fuer Klassen, die ein Kalendersystem implementieren
 * Abstrakte Basisklasse fuer Klassen, die ein Kalendersystem implementieren.
 *
 * Hier sind neben allerlei relevanten Rechenmethoden die grundlegenden
 * astronomischen Algorithmen für die Zeit- und Kalenderrechnung implementiert.
 * 
 * @author Ulrich Hilger
 * @version 2, 1.10.2022
 */
public abstract class BasisKalender implements Zeitrechnung {
  
@@ -82,6 +84,12 @@
      return nterTag(-1, t, datum);
   }
  
  /* ----- */
   public double moduloAngepasst(double x, double y) {
      return y + modulo(x, -y);
   }
  /* ---- Zeit ----- */
  
   public double zuMoment(int stunde, int minute, double sekunde) {
@@ -100,7 +108,10 @@
    return z;
   }  
    
  /* ----------- Mondphase ----------- */
  /* ----------- Mond ----------- */
  /** durchschnittliche Dauer eines Mondphasenzyklus (synodischer Monat) in Tagen */
   public static final double MITTLERER_SYNODISCHER_MONAT = 29.530588853;
  
   public double mondphase(double t) {
      return modulo(mondLaenge(t) - solareLaenge(t), 360);
@@ -109,7 +120,7 @@
   public double mondHoehe(double t, Ort ort) {
      double phi = ort.getBreite();
      double psi = ort.getLaenge();
      double varepsilon = obliquity(t);
      double varepsilon = schiefstand(t);
      double lambda = mondLaenge(t);
      double beta = mondBreite(t);
      double alpha = arcTanGrad(
@@ -266,6 +277,80 @@
      private static final double[] siderealCoeff = new double[] {280.46061837, 36525 * 360.98564736629, 0.000387933, 1d/38710000};
   }
   public double neumondNach(double tee) {
      return nterNeumond(1 + Math.round(tee / MITTLERER_SYNODISCHER_MONAT - mondphase(tee) / (double) (360)));
   }
   public double neumondVor(double tee) {
      return nterNeumond(Math.round(tee / MITTLERER_SYNODISCHER_MONAT - mondphase(tee) / (double) (360)));
   }
   public double nterNeumond(long n) {
      double k = n - 24724;
      double c = k / 1236.85;
      double approx = poly(c, nm.coeffApprox);
      double capE = poly(c, nm.coeffCapE);
      double solarAnomaly = poly(c, nm.coeffSolarAnomaly);
      double lunarAnomaly = poly(c, nm.coeffLunarAnomaly);
      double moonArgument = poly(c, nm.coeffMoonArgument);
      double capOmega = poly(c, nm.coeffCapOmega);
      double correction = -0.00017 * sinGrad(capOmega);
      for(int ix = 0; ix < nm.sineCoeff.length; ++ix) {
         correction += nm.sineCoeff[ix] * Math.pow(capE, nm.EFactor[ix]) *
            sinGrad(nm.solarCoeff[ix] * solarAnomaly +
               nm.lunarCoeff[ix] * lunarAnomaly +
               nm.moonCoeff[ix] * moonArgument);
      }
      double additional = 0;
      for(int ix = 0; ix < nm.addConst.length; ++ix) {
         additional += nm.addFactor[ix] *
            sinGrad(nm.addConst[ix] + nm.addCoeff[ix] * k);
      }
      double extra = 0.000325 * sinGrad(poly(c, nm.extra));
      return universalVonDynamisch(approx + correction + extra + additional);
   }
   private static class nm {
      private static final double[] coeffApprox = new double[] {730125.59765, MITTLERER_SYNODISCHER_MONAT * 1236.85, 0.0001337, -0.000000150, 0.00000000073};
      private static final double[] coeffCapE = new double[] {1, -0.002516, -0.0000074};
      private static final double[] coeffSolarAnomaly = new double[] {2.5534, 29.10535669 * 1236.85, -0.0000218, -0.00000011};
      private static final double[] coeffLunarAnomaly = new double[] {201.5643, 385.81693528 * 1236.85, 0.0107438, 0.00001239, -0.000000058};
      private static final double[] coeffMoonArgument = new double[] {160.7108, 390.67050274 * 1236.85, -0.0016341, -0.00000227, 0.000000011};
      private static final double[] coeffCapOmega = new double[] {124.7746, -1.56375580 * 1236.85, 0.0020691, 0.00000215};
      private static final byte[] EFactor = new byte[] {0, 1, 0, 0, 1, 1, 2, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
      private static final byte[] solarCoeff = new byte[] {0, 1, 0, 0, -1, 1, 2, 0, 0, 1, 0, 1, 1, -1, 2, 0, 3, 1, 0, 1, -1, -1, 1, 0};
      private static final byte[] lunarCoeff = new byte[] {1, 0, 2, 0, 1, 1, 0, 1, 1, 2, 3, 0, 0, 2, 1, 2, 0, 1, 2, 1, 1, 1, 3, 4};
      private static final byte[] moonCoeff = new byte[] {0, 0, 0, 2, 0, 0, 0, -2, 2, 0, 0, 2, -2, 0, 0, -2, 0, -2, 2, 2, 2, -2, 0, 0};
      private static final double[] sineCoeff = new double[] {
         -0.40720, 0.17241, 0.01608, 0.01039, 0.00739, -0.00514, 0.00208,
         -0.00111, -0.00057, 0.00056, -0.00042, 0.00042, 0.00038, -0.00024,
         -0.00007, 0.00004, 0.00004, 0.00003, 0.00003, -0.00003, 0.00003,
         -0.00002, -0.00002, 0.00002
      };
      private static final double[] addConst = new double[] {
         251.88, 251.83, 349.42, 84.66, 141.74, 207.14, 154.84, 34.52, 207.19,
         291.34, 161.72, 239.56, 331.55
      };
      private static final double[] addCoeff = new double[] {
         0.016321, 26.641886, 36.412478, 18.206239, 53.303771, 2.453732,
         7.306860, 27.261239, 0.121824, 1.844379, 24.198154, 25.513099, 3.592518
      };
      private static final double[] addFactor = new double[] {
         0.000165, 0.000164, 0.000126, 0.000110, 0.000062, 0.000060, 0.000056,
         0.000047, 0.000042, 0.000040, 0.000037, 0.000035, 0.000023
      };
      private static final double[] extra = new double[] {
         299.77, 132.8475848, -0.009173
      };
   }
   public double universalVonDynamisch(double tee) {
      return tee - ephemeridenKorrektur(tee);
   }
   public double universalVonStandard(double teeS, Ort locale) {
      return teeS - locale.getZeitzone() / 24;
   }
  /* ----------- Sonnenauf- und -untergang ----------- */
  
   public double sonnenaufgang(long date, Ort ort) throws Exception {
@@ -312,7 +397,7 @@
   public double zeitVonHorizont(double approx, Ort ort, double alpha) throws Exception {
      double phi = ort.getBreite();
      double t = universalVonLokal(approx, ort);
      double delta = arcSinGrad(sinGrad(obliquity(t)) * sinGrad(solareLaenge(t)));
      double delta = arcSinGrad(sinGrad(schiefstand(t)) * sinGrad(solareLaenge(t)));
      boolean morgen = modulo(approx, 1) < 0.5;
      double sinusAbstand = tanGrad(phi) * tanGrad(delta) +
         sinGrad(alpha) / (kosGrad(delta) * kosGrad(phi));
@@ -340,14 +425,14 @@
      double laenge = poly(c, et.koeffLaenge);
      double anomalie = poly(c, et.koeffAnomalie);
      double exzentrizitaet = poly(c, et.koeffExzentrizitaet);
      double varepsilon = obliquity(t);
      double varepsilon = schiefstand(t);
      double y = quadrat(tanGrad(varepsilon / 2));
      double equation = (1d / 2d / Math.PI) * (y * sinGrad(2 * laenge) +
      -2 * exzentrizitaet * sinGrad(anomalie) +
      4 * exzentrizitaet * y * sinGrad(anomalie) * kosGrad(2 * laenge) +
      -0.5 * y * y * sinGrad(4 * laenge) +
      -1.25 * exzentrizitaet * exzentrizitaet * sinGrad(2 * anomalie));
      return signum(equation) * Math.min(Math.abs(equation), stunde(12));
      return vorzeichen(equation) * Math.min(Math.abs(equation), stunde(12));
   }
   private static class et {
      private static final double[] koeffLaenge = new double[] {280.46645, 36000.76983, 0.0003032};
@@ -355,13 +440,13 @@
      private static final double[] koeffExzentrizitaet = new double[] {0.016708617, -0.000042037, -0.0000001236};
   }
   public double obliquity(double t) {
   public double schiefstand(double t) {
      double c = julJahrhunderte(t);
      return winkel(23, 26, 21.448) + poly(c, coeffObliquity);
   }
   private final double[] coeffObliquity = new double[] {0, winkel(0, 0, -46.8150), winkel(0, 0, -0.00059), winkel(0, 0, 0.001813)};
  
   public int signum(double x) {
   public int vorzeichen(double x) {
      if(x < 0)
         return -1;
      else if(x > 0)
@@ -400,7 +485,8 @@
  
  /* ---------------- Jahreszeiten ----- */
   public static final double TROPISCHES_JAHR = 365.242189;
  /** durchschnittliche Dauer eines Umlaufs der Erde um die Sonne in Tagen */
   public static final double MITTLERES_TROPISCHES_JAHR = 365.242189;
    
   public double standardVonUniversal(double t, Ort ort) {
      return t + ort.getZeitzone() / 24;
@@ -408,7 +494,7 @@
  
   public double solareLaengeNach(double t, double phi) {
      double varepsilon = 0.00001;
    double rate = TROPISCHES_JAHR / (double) 360;
    double rate = MITTLERES_TROPISCHES_JAHR / (double) 360;
      double tau = t + rate * modulo(phi - solareLaenge(t), 360);
      double l = Math.max(t, tau - 5);
      double u = tau + 5;
@@ -437,6 +523,13 @@
      return modulo(laenge + aberration(t) + nutation(t), 360);
   }
   public double geschaetzteSolareLaengeVor(double tee, double phi) {
      double rate = MITTLERES_TROPISCHES_JAHR / (double) (360);
      double tau = tee - rate * modulo(solareLaenge(tee) - phi, 360);
      double capDelta = modulo(solareLaenge(tau) - phi + (double) (180), 360) - (double) (180);
      return Math.min(tee, tau - rate * capDelta);
   }
   public double julJahrhunderte(double t) {
      return (dynamischVonUniversal(t) - j2000()) / 36525;
   }