/* Zeitrechnung - a class library to determine calendar events Copyright (c) 1984-2023 Ulrich Hilger, http://uhilger.de This program is free software: you can redistribute it and/or modify it under the terms of the GNU Affero General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for more details. You should have received a copy of the GNU Affero General Public License along with this program. If not, see . */ package de.uhilger.zeitrechnung.kalender; import de.uhilger.zeitrechnung.Definition; import de.uhilger.zeitrechnung.Ort; import de.uhilger.zeitrechnung.Zeit; /** * Abstrakte Basisklasse fuer Klassen, die ein Kalendersystem implementieren * * @author Ulrich Hilger * @version 2, 1.10.2022 */ public abstract class BasisKalender implements Zeitrechnung { /* Implementierung der Schnittstelle Zeitrechnung */ @Override public long ganzzahlQuotient(double x, double y) { return (long) Math.floor(x / y); } @Override public long modulo(long x, long y) { return (long) (x - y * Math.floor((double) x / y)); } public double modulo(double x, double y) { return x - y * Math.floor(x / y); } @Override public long tagNach(long datum, int t) { return tagAmOderNach(datum + 7, t); } @Override public long tagAmOderNach(long datum, int t) { return datum - wochentagVonGenerisch(datum - t); } @Override public long wochentagVonGenerisch(long datum) { return modulo(datum, 7); } @Override public long nterTag(int n, int t, long datum) { return n > 0 ? tagVor(datum, t) + 7 * n : tagNach(datum, t) + 7 * n; } @Override public long tagVor(long datum, int t) { return tagAmOderVor(datum - 1, t); } @Override public long tagAmOderVor(long datum, int t) { return datum - wochentagVonGenerisch(datum - t); } @Override public long letzterTag(int t, long datum) { return nterTag(-1, t, datum); } /* ---- Zeit ----- */ public double zuMoment(int stunde, int minute, double sekunde) { return stunde / 24d + minute / (24d * 60) + sekunde / (24d * 60 * 60); } public double zuMoment(Zeit z) { return BasisKalender.this.zuMoment(z.getStunde(), z.getMinute(), z.getSekunde()); } public Zeit vonMoment(double t) { Zeit z = new Zeit(); z.setStunde((int)Math.floor(modulo(t * 24, 24))); z.setMinute((int)Math.floor(modulo(t * 24 * 60, 60))); z.setSekunde(modulo(t * 24 * 60 * 60, 60)); return z; } /* ----------- Mondphase ----------- */ public double mondphase(double t) { return modulo(mondLaenge(t) - solareLaenge(t), 360); } public double mondHoehe(double t, Ort ort) { double phi = ort.getBreite(); double psi = ort.getLaenge(); double varepsilon = obliquity(t); double lambda = mondLaenge(t); double beta = mondBreite(t); double alpha = arcTanGrad( (sinGrad(lambda) * kosGrad(varepsilon) - tanGrad(beta) * sinGrad(varepsilon)) / kosGrad(lambda), (int)ganzzahlQuotient(lambda, (double) (90)) + 1 ); double delta = arcSinGrad(sinGrad(beta) * kosGrad(varepsilon) + kosGrad(beta) * sinGrad(varepsilon) * sinGrad(lambda)); double theta0 = siderischVonMoment(t); double capH = modulo(theta0 + psi - alpha, 360); double altitude = arcSinGrad(sinGrad(phi) * sinGrad(delta) + kosGrad(phi) * kosGrad(delta) * kosGrad(capH)); return modulo(altitude + (double) (180), 360) - (double) (180); } public double mondLaenge(double t) { double c = julJahrhunderte(t); double meanMoon = grad(poly(c, llon.coeffMeanMoon)); double elongation = grad(poly(c, llon.coeffElongation)); double solarAnomaly = grad(poly(c, llon.coeffSolarAnomaly)); double lunarAnomaly = grad(poly(c, llon.coeffLunarAnomaly)); double moonNode = grad(poly(c, llon.coeffMoonNode)); double capE = poly(c, llon.coeffCapE); double sigma = 0; for(int i = 0; i < llon.argsLunarElongation.length; ++i) { double x = llon.argsSolarAnomaly[i]; sigma += llon.sineCoefficients[i] * Math.pow(capE, Math.abs(x)) * sinGrad( llon.argsLunarElongation[i] * elongation + x * solarAnomaly + llon.argsLunarAnomaly[i] * lunarAnomaly + llon.argsMoonFromNode[i] * moonNode); } double correction = ((double) (1) / 1000000) * sigma; double venus = ((double) (3958) / 1000000) * sinGrad(119.75 + c * 131.849); double jupiter = ((double) (318) / 1000000) * sinGrad(53.09 + c * 479264.29); double flatEarth = ((double) (1962) / 1000000) * sinGrad(meanMoon - moonNode); return modulo(meanMoon + correction + venus + jupiter + flatEarth + nutation(t), 360); } private static class llon { private static final double[] coeffMeanMoon = new double[] {218.3164591, 481267.88134236, -0.0013268, 1d/538841, -1d/65194000}; private static final double[] coeffElongation = new double[] {297.8502042, 445267.1115168, -0.00163, 1d/545868, -1d/113065000}; private static final double[] coeffSolarAnomaly = new double[] {357.5291092, 35999.0502909, -0.0001536, 1d/24490000}; private static final double[] coeffLunarAnomaly = new double[] {134.9634114, 477198.8676313, 0.008997, 1d/69699, -1d/14712000}; private static final double[] coeffMoonNode = new double[] {93.2720993, 483202.0175273, -0.0034029, -1d/3526000, 1d/863310000}; private static final double[] coeffCapE = new double[] {1, -0.002516, -0.0000074}; private static final byte[] argsLunarElongation = new byte[] { 0, 2, 2, 0, 0, 0, 2, 2, 2, 2, 0, 1, 0, 2, 0, 0, 4, 0, 4, 2, 2, 1, 1, 2, 2, 4, 2, 0, 2, 2, 1, 2, 0, 0, 2, 2, 2, 4, 0, 3, 2, 4, 0, 2, 2, 2, 4, 0, 4, 1, 2, 0, 1, 3, 4, 2, 0, 1, 2 }; private static final byte[] argsSolarAnomaly = new byte[] { 0, 0, 0, 0, 1, 0, 0, -1, 0, -1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, -1, 0, 0, 0, 1, 0, -1, 0, -2, 1, 2, -2, 0, 0, -1, 0, 0, 1, -1, 2, 2, 1, -1, 0, 0, -1, 0, 1, 0, 1, 0, 0, -1, 2, 1, 0 }; private static final byte[] argsLunarAnomaly = new byte[] { 1, -1, 0, 2, 0, 0, -2, -1, 1, 0, -1, 0, 1, 0, 1, 1, -1, 3, -2, -1, 0, -1, 0, 1, 2, 0, -3, -2, -1, -2, 1, 0, 2, 0, -1, 1, 0, -1, 2, -1, 1, -2, -1, -1, -2, 0, 1, 4, 0, -2, 0, 2, 1, -2, -3, 2, 1, -1, 3 }; private static final byte[] argsMoonFromNode = new byte[] { 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, -2, 2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, -2, 2, 0, 2, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0, -2, -2, 0, 0, 0, 0, 0, 0, 0 }; private static final int[] sineCoefficients = new int[] { 6288774, 1274027, 658314, 213618, -185116, -114332, 58793, 57066, 53322, 45758, -40923, -34720, -30383, 15327, -12528, 10980, 10675, 10034, 8548, -7888, -6766, -5163, 4987, 4036, 3994, 3861, 3665, -2689, -2602, 2390, -2348, 2236, -2120, -2069, 2048, -1773, -1595, 1215, -1110, -892, -810, 759, -713, -700, 691, 596, 549, 537, 520, -487, -399, -381, 351, -340, 330, 327, -323, 299, 294 }; } // aus "Astronomical Algorithms" von Jean Meeus, // Willmann-Bell, Inc., 1998. public double mondBreite(double t) { double c = julJahrhunderte(t); double longitude = grad(poly(c, llat.coeffLongitude)); double elongation = grad(poly(c, llat.coeffElongation)); double solarAnomaly = grad(poly(c, llat.coeffSolarAnomaly)); double lunarAnomaly = grad(poly(c, llat.coeffLunarAnomaly)); double moonNode = grad(poly(c, llat.coeffMoonNode)); double capE = poly(c, llat.coeffCapE); double latitude = 0; for(int i = 0; i < llat.argsLunarElongation.length; ++i) { double x = llat.argsSolarAnomaly[i]; latitude += llat.sineCoefficients[i] * Math.pow(capE, Math.abs(x)) * sinGrad( llat.argsLunarElongation[i] * elongation + x * solarAnomaly + llat.argsLunarAnomaly[i] * lunarAnomaly + llat.argsMoonNode[i] * moonNode); } latitude *= (double) (1) / 1000000; double venus = ((double) (175) / 1000000) * (sinGrad((double) (119.75) + c * 131.849 + moonNode) + sinGrad((double) (119.75) + c * 131.849 - moonNode)); double flatEarth = ((double) (-2235) / 1000000) * sinGrad(longitude) + ((double) (127) / 1000000) * sinGrad(longitude - lunarAnomaly) + ((double) (-115) / 1000000) * sinGrad(longitude + lunarAnomaly); double extra = ((double) (382) / 1000000) * sinGrad((double) (313.45) + c * (double) (481266.484)); return modulo(latitude + venus + flatEarth + extra, 360); } private static class llat { private static final double[] coeffLongitude = new double[] {218.3164591, 481267.88134236, -0.0013268, 1d/538841, -1d/65194000}; private static final double[] coeffElongation = new double[] {297.8502042, 445267.1115168, -0.00163, 1d/545868, -1d/113065000}; private static final double[] coeffSolarAnomaly = new double[] {357.5291092, 35999.0502909, -0.0001536, 1d/24490000}; private static final double[] coeffLunarAnomaly = new double[] {134.9634114, 477198.8676313, 0.008997, 1d/69699, -1d/14712000}; private static final double[] coeffMoonNode = new double[] {93.2720993, 483202.0175273, -0.0034029, -1d/3526000, 1d/863310000}; private static final double[] coeffCapE = new double[] {1, -0.002516, -0.0000074}; private static final byte[] argsLunarElongation = new byte[] { 0, 0, 0, 2, 2, 2, 2, 0, 2, 0, 2, 2, 2, 2, 2, 2, 2, 0, 4, 0, 0, 0, 1, 0, 0, 0, 1, 0, 4, 4, 0, 4, 2, 2, 2, 2, 0, 2, 2, 2, 2, 4, 2, 2, 0, 2, 1, 1, 0, 2, 1, 2, 0, 4, 4, 1, 4, 1, 4, 2}; private static final byte[] argsSolarAnomaly = new byte[] { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, -1, -1, -1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 1, 0, -1, -2, 0, 1, 1, 1, 1, 1, 0, -1, 1, 0, -1, 0, 0, 0, -1, -2}; private static final byte[] argsLunarAnomaly = new byte[] { 0, 1, 1, 0, -1, -1, 0, 2, 1, 2, 0, -2, 1, 0, -1, 0, -1, -1, -1, 0, 0, -1, 0, 1, 1, 0, 0, 3, 0, -1, 1, -2, 0, 2, 1, -2, 3, 2, -3, -1, 0, 0, 1, 0, 1, 1, 0, 0, -2, -1, 1, -2, 2, -2, -1, 1, 1, -2, 0, 0}; private static final byte[] argsMoonNode = new byte[] { 1, 1, -1, -1, 1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, 1, -1, -1, -1, 1, 3, 1, 1, 1, -1, -1, -1, 1, -1, 1, -3, 1, -3, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 3, -1, -1, 1, -1, -1, 1, -1, 1, -1, -1, -1, -1, -1, -1, 1}; private static final int[] sineCoefficients = new int[] { 5128122, 280602, 277693, 173237, 55413, 46271, 32573, 17198, 9266, 8822, 8216, 4324, 4200, -3359, 2463, 2211, 2065, -1870, 1828, -1794, -1749, -1565, -1491, -1475, -1410, -1344, -1335, 1107, 1021, 833, 777, 671, 607, 596, 491, -451, 439, 422, 421, -366, -351, 331, 315, 302, -283, -229, 223, 223, -220, -220, -185, 181, -177, 176, 166, -164, 132, -119, 115, 107}; } public double arcTanGrad(double x, int quad) { double alpha = bogenmassZuGrad(Math.atan(x)); return modulo(quad == 1 || quad == 4 ? alpha : alpha + (double) (180), 360); } public double siderischVonMoment(double t) { double c = (t - j2000()) / 36525; return modulo(poly(c, sfm.siderealCoeff), 360); } private static class sfm { private static final double[] siderealCoeff = new double[] {280.46061837, 36525 * 360.98564736629, 0.000387933, 1d/38710000}; } /* ----------- Sonnenauf- und -untergang ----------- */ public double sonnenaufgang(long date, Ort ort) throws Exception { return morgen(date, ort, alpha(ort)); } public double sonnenuntergang(long date, Ort ort) throws Exception { return abend(date, ort, alpha(ort)); } public double alpha(Ort ort) { double h = Math.max(0, ort.getHoehe()); final double capR = (double) 6.372E+6; double dip = arcKosGrad(capR / (capR + h)); return winkel(0, 50, 0) + dip; } public double morgen(long date, Ort ort, double alpha) throws Exception { double approx; try { approx = zeitVonHorizont(date + .25, ort, alpha); } catch(Exception ex) { approx = date; } double result = zeitVonHorizont(approx, ort, alpha); return standardVonLokal(result, ort); } public double abend(long date, Ort ort, double alpha) throws Exception { double approx; try { approx = zeitVonHorizont(date + .75, ort, alpha); } catch(Exception ex) { approx = date + .99d; } double result = zeitVonHorizont(approx, ort, alpha); return standardVonLokal(result, ort); } 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))); boolean morgen = modulo(approx, 1) < 0.5; double sinusAbstand = tanGrad(phi) * tanGrad(delta) + sinGrad(alpha) / (kosGrad(delta) * kosGrad(phi)); double offset = modulo(0.5 + arcSinGrad(sinusAbstand) / (double) 360, 1) - 0.5; if(Math.abs(sinusAbstand) > 1) { throw new Exception(); } return lokalVonScheinbar(Math.floor(approx) + (morgen ? .25 - offset : .75 + offset)); } public double universalVonLokal(double tl, Ort ort) { return tl - ort.getLaenge() / (double) 360; } public double standardVonLokal(double tl, Ort ort) { return standardVonUniversal(universalVonLokal(tl, ort), ort); } public double lokalVonScheinbar(double t) { return t - zeitgleichung(t); } public double zeitgleichung(double t) { double c = julJahrhunderte(t); double laenge = poly(c, et.koeffLaenge); double anomalie = poly(c, et.koeffAnomalie); double exzentrizitaet = poly(c, et.koeffExzentrizitaet); double varepsilon = obliquity(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)); } private static class et { private static final double[] koeffLaenge = new double[] {280.46645, 36000.76983, 0.0003032}; private static final double[] koeffAnomalie = new double[] {357.52910, 35999.05030, -0.0001559, -0.00000048}; private static final double[] koeffExzentrizitaet = new double[] {0.016708617, -0.000042037, -0.0000001236}; } public double obliquity(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) { if(x < 0) return -1; else if(x > 0) return 1; else return 0; } public double quadrat(double x) { return x * x; } public double kosGrad(double theta) { return Math.cos(gradZuBogenmass(theta)); } public double arcSinGrad(double x) { return bogenmassZuGrad(Math.asin(x)); } public double tanGrad(double theta) { return Math.tan(gradZuBogenmass(theta)); } public double arcKosGrad(double x) { return bogenmassZuGrad(Math.acos(x)); } public double bogenmassZuGrad(double theta) { return grad(theta / Math.PI * 180); } public double winkel(double d, double m, double s) { return d + (m + s / 60) / 60; } /* ---------------- Jahreszeiten ----- */ public static final double TROPISCHES_JAHR = 365.242189; public double standardVonUniversal(double t, Ort ort) { return t + ort.getZeitzone() / 24; } public double solareLaengeNach(double t, double phi) { double varepsilon = 0.00001; double rate = TROPISCHES_JAHR / (double) 360; double tau = t + rate * modulo(phi - solareLaenge(t), 360); double l = Math.max(t, tau - 5); double u = tau + 5; double lo = l, hi = u, x = (hi + lo) / 2; while(hi - lo > varepsilon) { if(modulo(solareLaenge(x) - phi, 360) < (double) 180) hi = x; else lo = x; x = (hi + lo) / 2; } return x; } public double solareLaenge(double t) { double c = julJahrhunderte(t); double sigma = 0; for(int i = 0; i < sLaenge.koeffizienten.length; ++i) { sigma += sLaenge.koeffizienten[i] * sinGrad(sLaenge.multiplikatoren[i] * c + sLaenge.summanden[i]); } double laenge = (double) 282.7771834 + 36000.76953744 * c + 0.000005729577951308232 * sigma; return modulo(laenge + aberration(t) + nutation(t), 360); } public double julJahrhunderte(double t) { return (dynamischVonUniversal(t) - j2000()) / 36525; } public double dynamischVonUniversal(double tee) { return tee + ephemeridenKorrektur(tee); } public double ephemeridenKorrektur(double t) { double[] koeffizient17tes = new double[] {196.58333, -4.0675, 0.0219167}; double[] koeffizient19tes = new double[] {-0.00002, 0.000297, 0.025184, -0.181133, 0.553040, -0.861938, 0.677066, -0.212591}; double[] koeffizient18tes = new double[] {-0.000009, 0.003844, 0.083563, 0.865736, 4.867575, 15.845535, 31.332267, 38.291999, 28.316289, 11.636204, 2.043794}; ISOKalender w = new ISOKalender(); long jahr = w.jahrVonTagen((long)Math.floor(t)); double c = differenz(w.zuTagen(1900, Definition.JANUAR, 1), w.zuTagen(jahr, Definition.JULI, 1)) / 36525d; double result; if(1988 <= jahr && jahr <= 2019) { result = (jahr - 1933) / (24d * 60 * 60); } else if (1900 <= jahr && jahr <= 1987) { result = poly(c, koeffizient19tes); } else if (1800 <= jahr && jahr <= 1899) { result = poly(c, koeffizient18tes); } else if (1620 <= jahr && jahr <= 1799) { result = poly(jahr - 1600, koeffizient17tes) / (24 * 60 * 60); } else { double x = stunde(12) + differenz(w.zuTagen(1810, Definition.JANUAR, 1), w.zuTagen(jahr, Definition.JANUAR, 1)); return (x * x / 41048480 - 15) / (24 * 60 * 60); } return result; } public double nutation(double t) { double[] koeffa = new double[] {124.90, -1934.134, 0.002063}; double[] koeffb = new double[] {201.11, 72001.5377, 0.00057}; double c = julJahrhunderte(t); double capA = poly(c, koeffa); double capB = poly(c, koeffb); return (double) -0.004778 * sinGrad(capA) + (double) -0.0003667 * sinGrad(capB); } public static long differenz(long datum1, long datum2) { return datum2 - datum1; } public double stunde(double x) { return x / 24; } public double j2000() { ISOKalender w = new ISOKalender(); return stunde(12) + w.zuTagen(2000, Definition.JANUAR, 1); } public double sinGrad(double theta) { return Math.sin(gradZuBogenmass(theta)); } public double gradZuBogenmass(double theta) { return grad(theta) * Math.PI / 180; } public double grad(double theta) { return modulo(theta, 360); } public double aberration(double t) { double c = julJahrhunderte(t); return (double) 0.0000974 * gradKosinus((double) 177.63 + (double) 35999.01848 * c) - (double) 0.005575; } public double gradKosinus(double theta) { return Math.cos(gradZuBogenmass(theta)); } public double poly(double x, double[] a) { double ergebnis = a[0]; for(int i = 1; i < a.length; ++i) { ergebnis += a[i] * Math.pow(x, i); } return ergebnis; } private static class sLaenge { private static final int[] koeffizienten = new int[] { 403406, 195207, 119433, 112392, 3891, 2819, 1721, 660, 350, 334, 314, 268, 242, 234, 158, 132, 129, 114, 99, 93, 86, 78, 72, 68, 64, 46, 38, 37, 32, 29, 28, 27, 27, 25, 24, 21, 21, 20, 18, 17, 14, 13, 13, 13, 12, 10, 10, 10, 10 }; private static final double[] multiplikatoren = new double[] { 0.9287892, 35999.1376958, 35999.4089666, 35998.7287385, 71998.20261, 71998.4403, 36000.35726, 71997.4812, 32964.4678, -19.4410, 445267.1117, 45036.8840, 3.1008, 22518.4434, -19.9739, 65928.9345, 9038.0293, 3034.7684, 33718.148, 3034.448, -2280.773, 29929.992, 31556.493, 149.588, 9037.750, 107997.405, -4444.176, 151.771, 67555.316, 31556.080, -4561.540, 107996.706, 1221.655, 62894.167, 31437.369, 14578.298, -31931.757, 34777.243, 1221.999, 62894.511, -4442.039, 107997.909, 119.066, 16859.071, -4.578, 26895.292, -39.127, 12297.536, 90073.778 }; private static final double[] summanden = new double[] { 270.54861, 340.19128, 63.91854, 331.26220, 317.843, 86.631, 240.052, 310.26, 247.23, 260.87, 297.82, 343.14, 166.79, 81.53, 3.50, 132.75, 182.95, 162.03, 29.8, 266.4, 249.2, 157.6, 257.8, 185.1, 69.9, 8.0, 197.1, 250.4, 65.3, 162.7, 341.5, 291.6, 98.5, 146.7, 110.0, 5.2, 342.6, 230.9, 256.1, 45.3, 242.9, 115.2, 151.8, 285.3, 53.3, 126.6, 205.7, 85.9, 146.1 }; } }