Bestimmung der Zeitpunkte von Ereignissen
ulrich
2023-03-23 d20d989f5495492f1258c8313db7c19b429111a3
commit | author | age
66d68b 1 /*
U 2   Zeitrechnung - a class library to determine calendar events
3   Copyright (c) 1984-2023 Ulrich Hilger, http://uhilger.de
4
5   This program is free software: you can redistribute it and/or modify
6   it under the terms of the GNU Affero General Public License as published by
7   the Free Software Foundation, either version 3 of the License, or
8   (at your option) any later version.
9
10   This program is distributed in the hope that it will be useful,
11   but WITHOUT ANY WARRANTY; without even the implied warranty of
12   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13   GNU Affero General Public License for more details.
14
15   You should have received a copy of the GNU Affero General Public License
16   along with this program.  If not, see <http://www.gnu.org/licenses/>.
17  */
18 package de.uhilger.zeitrechnung.kalender;
19
20 import de.uhilger.zeitrechnung.Definition;
21 import de.uhilger.zeitrechnung.Ort;
22 import de.uhilger.zeitrechnung.Zeit;
23
24 /**
d20d98 25  * Abstrakte Basisklasse fuer Klassen, die ein Kalendersystem implementieren.
U 26  * 
27  * Hier sind neben allerlei relevanten Rechenmethoden die grundlegenden 
28  * astronomischen Algorithmen für die Zeit- und Kalenderrechnung implementiert.
66d68b 29  * 
U 30  * @author Ulrich Hilger
31  */
32 public abstract class BasisKalender implements Zeitrechnung {
33   
34   /* Implementierung der Schnittstelle Zeitrechnung */
35   
36   @Override
37   public long ganzzahlQuotient(double x, double y) {
38     return (long) Math.floor(x / y);
39   }
40
41   @Override
42   public long modulo(long x, long y) {
43     return (long) (x - y * Math.floor((double) x / y));
44   }
45   
46     public double modulo(double x, double y) {
47         return x - y * Math.floor(x / y);
48     }
49
50   @Override
51   public long tagNach(long datum, int t) {
52     return tagAmOderNach(datum + 7, t);
53   }
54
55   @Override
56   public long tagAmOderNach(long datum, int t) {
57     return datum - wochentagVonGenerisch(datum - t);
58   }
59
60   @Override
61   public long wochentagVonGenerisch(long datum) {
62     return modulo(datum, 7);
63   }
64     
65   @Override
66     public long nterTag(int n, int t, long datum) {
67         return n > 0 ?
68             tagVor(datum, t) + 7 * n :
69             tagNach(datum, t) + 7 * n;
70     }
71   
72   @Override
73     public long tagVor(long datum, int t) {
74         return tagAmOderVor(datum - 1, t);
75     }
76
77   @Override
78     public long tagAmOderVor(long datum, int t) {
79         return datum - wochentagVonGenerisch(datum - t);
80     }
81
82   @Override
83     public long letzterTag(int t, long datum) {
84         return nterTag(-1, t, datum);
85     }
86   
d20d98 87   /* ----- */
U 88   
89     public double moduloAngepasst(double x, double y) {
90         return y + modulo(x, -y);
91     }
92   
66d68b 93   /* ---- Zeit ----- */
U 94   
95     public double zuMoment(int stunde, int minute, double sekunde) {
96         return stunde / 24d + minute / (24d * 60) + sekunde / (24d * 60 * 60);
97     }
98   
99   public double zuMoment(Zeit z) {
100     return BasisKalender.this.zuMoment(z.getStunde(), z.getMinute(), z.getSekunde());
101   }
102
103     public Zeit vonMoment(double t) {
104     Zeit z = new Zeit();
105         z.setStunde((int)Math.floor(modulo(t * 24, 24)));
106         z.setMinute((int)Math.floor(modulo(t * 24 * 60, 60)));
107         z.setSekunde(modulo(t * 24 * 60 * 60, 60));
108     return z;
109     }  
110     
d20d98 111   /* ----------- Mond ----------- */
U 112   
113   /** durchschnittliche Dauer eines Mondphasenzyklus (synodischer Monat) in Tagen */
114     public static final double MITTLERER_SYNODISCHER_MONAT = 29.530588853;
66d68b 115   
U 116     public double mondphase(double t) {
117         return modulo(mondLaenge(t) - solareLaenge(t), 360);
118     }
119     
120     public double mondHoehe(double t, Ort ort) {
121         double phi = ort.getBreite();
122         double psi = ort.getLaenge();
d20d98 123         double varepsilon = schiefstand(t);
66d68b 124         double lambda = mondLaenge(t);
U 125         double beta = mondBreite(t);
126         double alpha = arcTanGrad(
127             (sinGrad(lambda) * kosGrad(varepsilon) - tanGrad(beta) * sinGrad(varepsilon)) /
128             kosGrad(lambda),
129             (int)ganzzahlQuotient(lambda, (double) (90)) + 1
130         );
131         double delta = arcSinGrad(sinGrad(beta) * kosGrad(varepsilon) +
132             kosGrad(beta) * sinGrad(varepsilon) * sinGrad(lambda));
133         double theta0 = siderischVonMoment(t);
134         double capH = modulo(theta0 + psi - alpha, 360);
135         double altitude = arcSinGrad(sinGrad(phi) * sinGrad(delta) + kosGrad(phi) * kosGrad(delta) * kosGrad(capH));
136         return modulo(altitude + (double) (180), 360) - (double) (180);
137     }
138   
139     public double mondLaenge(double t) {
140         double c = julJahrhunderte(t);
141         double meanMoon = grad(poly(c, llon.coeffMeanMoon));
142         double elongation = grad(poly(c, llon.coeffElongation));
143         double solarAnomaly = grad(poly(c, llon.coeffSolarAnomaly));
144         double lunarAnomaly = grad(poly(c, llon.coeffLunarAnomaly));
145         double moonNode = grad(poly(c, llon.coeffMoonNode));
146         double capE = poly(c, llon.coeffCapE);
147         double sigma = 0;
148         for(int i = 0; i < llon.argsLunarElongation.length; ++i) {
149             double x = llon.argsSolarAnomaly[i];
150             sigma += llon.sineCoefficients[i] *
151                 Math.pow(capE, Math.abs(x)) *
152                 sinGrad(    llon.argsLunarElongation[i] * elongation + 
153                     x * solarAnomaly +
154                     llon.argsLunarAnomaly[i] * lunarAnomaly +
155                     llon.argsMoonFromNode[i] * moonNode);
156         }
157         double correction = ((double) (1) / 1000000) * sigma;
158         double venus = ((double) (3958) / 1000000) * sinGrad(119.75 + c * 131.849);
159         double jupiter = ((double) (318) / 1000000) * sinGrad(53.09 + c * 479264.29);
160         double flatEarth = ((double) (1962) / 1000000) * sinGrad(meanMoon - moonNode);
161         return modulo(meanMoon + correction + venus + jupiter + flatEarth + nutation(t), 360);
162     }
163     private static class llon {
164         private static final double[] coeffMeanMoon = new double[] {218.3164591, 481267.88134236, -0.0013268, 1d/538841, -1d/65194000};
165         private static final double[] coeffElongation = new double[] {297.8502042, 445267.1115168, -0.00163, 1d/545868, -1d/113065000};
166         private static final double[] coeffSolarAnomaly = new double[] {357.5291092, 35999.0502909, -0.0001536, 1d/24490000};
167         private static final double[] coeffLunarAnomaly = new double[] {134.9634114, 477198.8676313, 0.008997, 1d/69699, -1d/14712000};
168         private static final double[] coeffMoonNode = new double[] {93.2720993, 483202.0175273, -0.0034029, -1d/3526000, 1d/863310000};
169         private static final double[] coeffCapE = new double[] {1, -0.002516, -0.0000074};
170         private static final byte[] argsLunarElongation = new byte[] {
171             0, 2, 2, 0, 0, 0, 2, 2, 2, 2, 0, 1, 0, 2, 0, 0, 4, 0, 4, 2, 2, 1,
172             1, 2, 2, 4, 2, 0, 2, 2, 1, 2, 0, 0, 2, 2, 2, 4, 0, 3, 2, 4, 0, 2,
173             2, 2, 4, 0, 4, 1, 2, 0, 1, 3, 4, 2, 0, 1, 2
174         };
175         private static final byte[] argsSolarAnomaly = new byte[] {
176             0, 0, 0, 0, 1, 0, 0, -1, 0, -1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1,
177             0, 1, -1, 0, 0, 0, 1, 0, -1, 0, -2, 1, 2, -2, 0, 0, -1, 0, 0, 1,
178             -1, 2, 2, 1, -1, 0, 0, -1, 0, 1, 0, 1, 0, 0, -1, 2, 1, 0
179         };
180         private static final byte[] argsLunarAnomaly = new byte[] {
181             1, -1, 0, 2, 0, 0, -2, -1, 1, 0, -1, 0, 1, 0, 1, 1, -1, 3, -2,
182             -1, 0, -1, 0, 1, 2, 0, -3, -2, -1, -2, 1, 0, 2, 0, -1, 1, 0,
183             -1, 2, -1, 1, -2, -1, -1, -2, 0, 1, 4, 0, -2, 0, 2, 1, -2, -3,
184             2, 1, -1, 3
185         };
186         private static final byte[] argsMoonFromNode = new byte[] {
187             0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, -2, 2, -2, 0, 0, 0, 0, 0,
188             0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, -2, 2, 0, 2, 0, 0, 0, 0,
189             0, 0, -2, 0, 0, 0, 0, -2, -2, 0, 0, 0, 0, 0, 0, 0
190         };
191         private static final int[] sineCoefficients = new int[] {
192             6288774, 1274027, 658314, 213618, -185116, -114332,
193             58793, 57066, 53322, 45758, -40923, -34720, -30383,
194             15327, -12528, 10980, 10675, 10034, 8548, -7888,
195             -6766, -5163, 4987, 4036, 3994, 3861, 3665, -2689,
196             -2602, 2390, -2348, 2236, -2120, -2069, 2048, -1773,
197             -1595, 1215, -1110, -892, -810, 759, -713, -700, 691,
198             596, 549, 537, 520, -487, -399, -381, 351, -340, 330,
199             327, -323, 299, 294
200         };
201     }
202
203     // aus "Astronomical Algorithms" von Jean Meeus,
204     // Willmann-Bell, Inc., 1998.
205
206     public double mondBreite(double t) {
207         double c = julJahrhunderte(t);
208         double longitude = grad(poly(c, llat.coeffLongitude));
209         double elongation = grad(poly(c, llat.coeffElongation));
210         double solarAnomaly = grad(poly(c, llat.coeffSolarAnomaly));
211         double lunarAnomaly = grad(poly(c, llat.coeffLunarAnomaly));
212         double moonNode = grad(poly(c, llat.coeffMoonNode));
213         double capE = poly(c, llat.coeffCapE);
214         double latitude = 0;
215         for(int i = 0; i < llat.argsLunarElongation.length; ++i) {
216             double x = llat.argsSolarAnomaly[i];
217             latitude += llat.sineCoefficients[i] *
218                 Math.pow(capE, Math.abs(x)) *
219                 sinGrad(    llat.argsLunarElongation[i] * elongation + 
220                     x * solarAnomaly +
221                     llat.argsLunarAnomaly[i] * lunarAnomaly +
222                     llat.argsMoonNode[i] * moonNode);
223         }
224         latitude *= (double) (1) / 1000000;
225         double venus = ((double) (175) / 1000000) * (sinGrad((double) (119.75) + c * 131.849 + moonNode) + sinGrad((double) (119.75) + c * 131.849 - moonNode));
226         double flatEarth = ((double) (-2235) / 1000000) * sinGrad(longitude) +
227             ((double) (127) / 1000000) * sinGrad(longitude - lunarAnomaly) +
228             ((double) (-115) / 1000000) * sinGrad(longitude + lunarAnomaly);
229         double extra = ((double) (382) / 1000000) * sinGrad((double) (313.45) + c * (double) (481266.484));
230         return modulo(latitude + venus + flatEarth + extra, 360);
231     }
232     private static class llat {
233         private static final double[] coeffLongitude = new double[] {218.3164591, 481267.88134236, -0.0013268, 1d/538841, -1d/65194000};
234         private static final double[] coeffElongation = new double[] {297.8502042, 445267.1115168, -0.00163, 1d/545868, -1d/113065000};
235         private static final double[] coeffSolarAnomaly = new double[] {357.5291092, 35999.0502909, -0.0001536, 1d/24490000};
236         private static final double[] coeffLunarAnomaly = new double[] {134.9634114, 477198.8676313, 0.008997, 1d/69699, -1d/14712000};
237         private static final double[] coeffMoonNode = new double[] {93.2720993, 483202.0175273, -0.0034029, -1d/3526000, 1d/863310000};
238         private static final double[] coeffCapE = new double[] {1, -0.002516, -0.0000074};
239         private static final byte[] argsLunarElongation = new byte[] {
240             0, 0, 0, 2, 2, 2, 2, 0, 2, 0, 2, 2, 2, 2, 2, 2, 2, 0, 4, 0, 0, 0,
241             1, 0, 0, 0, 1, 0, 4, 4, 0, 4, 2, 2, 2, 2, 0, 2, 2, 2, 2, 4, 2, 2,
242             0, 2, 1, 1, 0, 2, 1, 2, 0, 4, 4, 1, 4, 1, 4, 2};
243         private static final byte[] argsSolarAnomaly = new byte[] {
244             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, -1, -1, -1, 1, 0, 1,
245             0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 1,
246             0, -1, -2, 0, 1, 1, 1, 1, 1, 0, -1, 1, 0, -1, 0, 0, 0, -1, -2};
247         private static final byte[] argsLunarAnomaly = new byte[] {
248             0, 1, 1, 0, -1, -1, 0, 2, 1, 2, 0, -2, 1, 0, -1, 0, -1, -1, -1,
249             0, 0, -1, 0, 1, 1, 0, 0, 3, 0, -1, 1, -2, 0, 2, 1, -2, 3, 2, -3,
250             -1, 0, 0, 1, 0, 1, 1, 0, 0, -2, -1, 1, -2, 2, -2, -1, 1, 1, -2,
251             0, 0};
252         private static final byte[] argsMoonNode = new byte[] {
253             1, 1, -1, -1, 1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, 1, -1, -1,
254             -1, 1, 3, 1, 1, 1, -1, -1, -1, 1, -1, 1, -3, 1, -3, -1, -1, 1,
255             -1, 1, -1, 1, 1, 1, 1, -1, 3, -1, -1, 1, -1, -1, 1, -1, 1, -1,
256             -1, -1, -1, -1, -1, 1};
257         private static final int[] sineCoefficients = new int[] {
258             5128122, 280602, 277693, 173237, 55413, 46271, 32573,
259             17198, 9266, 8822, 8216, 4324, 4200, -3359, 2463, 2211,
260             2065, -1870, 1828, -1794, -1749, -1565, -1491, -1475,
261             -1410, -1344, -1335, 1107, 1021, 833, 777, 671, 607,
262             596, 491, -451, 439, 422, 421, -366, -351, 331, 315,
263             302, -283, -229, 223, 223, -220, -220, -185, 181,
264             -177, 176, 166, -164, 132, -119, 115, 107};
265     }
266   
267     public double arcTanGrad(double x, int quad) {
268         double alpha = bogenmassZuGrad(Math.atan(x));
269         return modulo(quad == 1 || quad == 4 ? alpha : alpha + (double) (180), 360);
270     }
271   
272     public double siderischVonMoment(double t) {
273         double c = (t - j2000()) / 36525;
274         return modulo(poly(c, sfm.siderealCoeff), 360);
275     }
276     private static class sfm {
277         private static final double[] siderealCoeff = new double[] {280.46061837, 36525 * 360.98564736629, 0.000387933, 1d/38710000};
278     }
279
d20d98 280     public double neumondNach(double tee) {
U 281         return nterNeumond(1 + Math.round(tee / MITTLERER_SYNODISCHER_MONAT - mondphase(tee) / (double) (360)));
282     }
283     
284     public double neumondVor(double tee) {
285         return nterNeumond(Math.round(tee / MITTLERER_SYNODISCHER_MONAT - mondphase(tee) / (double) (360)));
286     }
287
288     public double nterNeumond(long n) {
289         double k = n - 24724;
290         double c = k / 1236.85;
291         double approx = poly(c, nm.coeffApprox);
292         double capE = poly(c, nm.coeffCapE);
293         double solarAnomaly = poly(c, nm.coeffSolarAnomaly);
294         double lunarAnomaly = poly(c, nm.coeffLunarAnomaly);
295         double moonArgument = poly(c, nm.coeffMoonArgument);
296         double capOmega = poly(c, nm.coeffCapOmega);
297         double correction = -0.00017 * sinGrad(capOmega);
298         for(int ix = 0; ix < nm.sineCoeff.length; ++ix) {
299             correction += nm.sineCoeff[ix] * Math.pow(capE, nm.EFactor[ix]) *
300                 sinGrad(nm.solarCoeff[ix] * solarAnomaly +
301                     nm.lunarCoeff[ix] * lunarAnomaly +
302                     nm.moonCoeff[ix] * moonArgument);
303         }
304         double additional = 0;
305         for(int ix = 0; ix < nm.addConst.length; ++ix) {
306             additional += nm.addFactor[ix] *
307                 sinGrad(nm.addConst[ix] + nm.addCoeff[ix] * k);
308         }
309         double extra = 0.000325 * sinGrad(poly(c, nm.extra));
310         return universalVonDynamisch(approx + correction + extra + additional);
311     }
312     private static class nm {
313         private static final double[] coeffApprox = new double[] {730125.59765, MITTLERER_SYNODISCHER_MONAT * 1236.85, 0.0001337, -0.000000150, 0.00000000073};
314         private static final double[] coeffCapE = new double[] {1, -0.002516, -0.0000074};
315         private static final double[] coeffSolarAnomaly = new double[] {2.5534, 29.10535669 * 1236.85, -0.0000218, -0.00000011};
316         private static final double[] coeffLunarAnomaly = new double[] {201.5643, 385.81693528 * 1236.85, 0.0107438, 0.00001239, -0.000000058};
317         private static final double[] coeffMoonArgument = new double[] {160.7108, 390.67050274 * 1236.85, -0.0016341, -0.00000227, 0.000000011};
318         private static final double[] coeffCapOmega = new double[] {124.7746, -1.56375580 * 1236.85, 0.0020691, 0.00000215};
319         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};
320         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};
321         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};
322         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};
323         private static final double[] sineCoeff = new double[] {
324             -0.40720, 0.17241, 0.01608, 0.01039, 0.00739, -0.00514, 0.00208,
325             -0.00111, -0.00057, 0.00056, -0.00042, 0.00042, 0.00038, -0.00024,
326             -0.00007, 0.00004, 0.00004, 0.00003, 0.00003, -0.00003, 0.00003,
327             -0.00002, -0.00002, 0.00002
328         };
329         private static final double[] addConst = new double[] {
330             251.88, 251.83, 349.42, 84.66, 141.74, 207.14, 154.84, 34.52, 207.19,
331             291.34, 161.72, 239.56, 331.55
332         };
333         private static final double[] addCoeff = new double[] {
334             0.016321, 26.641886, 36.412478, 18.206239, 53.303771, 2.453732,
335             7.306860, 27.261239, 0.121824, 1.844379, 24.198154, 25.513099, 3.592518
336         };
337         private static final double[] addFactor = new double[] {
338             0.000165, 0.000164, 0.000126, 0.000110, 0.000062, 0.000060, 0.000056,
339             0.000047, 0.000042, 0.000040, 0.000037, 0.000035, 0.000023
340         };
341         private static final double[] extra = new double[] {
342             299.77, 132.8475848, -0.009173
343         };
344     }
345
346     public double universalVonDynamisch(double tee) {
347         return tee - ephemeridenKorrektur(tee);
348     }
349
350     public double universalVonStandard(double teeS, Ort locale) {
351         return teeS - locale.getZeitzone() / 24;
352     }
353   
66d68b 354   /* ----------- Sonnenauf- und -untergang ----------- */
U 355   
356     public double sonnenaufgang(long date, Ort ort) throws Exception {
357         return morgen(date, ort, alpha(ort));
358     }
359   
360     public double sonnenuntergang(long date, Ort ort)
361         throws Exception
362     {
363         return abend(date, ort, alpha(ort));
364     }
365   
366   public double alpha(Ort ort) {
367         double h = Math.max(0, ort.getHoehe());
368         final double capR = (double) 6.372E+6;
369         double dip = arcKosGrad(capR / (capR + h));
370         return winkel(0, 50, 0) + dip;
371   }
372   
373     public double morgen(long date, Ort ort, double alpha) throws Exception {
374         double approx;
375         try {
376             approx = zeitVonHorizont(date + .25, ort, alpha);
377         } catch(Exception ex) {
378             approx = date;
379         }
380         double result = zeitVonHorizont(approx, ort, alpha);
381         return standardVonLokal(result, ort);
382     }
383     
384     public double abend(long date, Ort ort, double alpha)
385         throws Exception
386     {
387         double approx;
388         try {
389             approx = zeitVonHorizont(date + .75, ort, alpha);
390         } catch(Exception ex) {
391             approx = date + .99d;
392         }
393         double result = zeitVonHorizont(approx, ort, alpha);
394         return standardVonLokal(result, ort);
395     }
396
397     public double zeitVonHorizont(double approx, Ort ort, double alpha) throws Exception {
398         double phi = ort.getBreite();
399         double t = universalVonLokal(approx, ort);
d20d98 400         double delta = arcSinGrad(sinGrad(schiefstand(t)) * sinGrad(solareLaenge(t)));
66d68b 401         boolean morgen = modulo(approx, 1) < 0.5;
U 402         double sinusAbstand = tanGrad(phi) * tanGrad(delta) +
403             sinGrad(alpha) / (kosGrad(delta) * kosGrad(phi));
404         double offset = modulo(0.5 + arcSinGrad(sinusAbstand) / (double) 360, 1) - 0.5;
405         if(Math.abs(sinusAbstand) > 1) {
406             throw new Exception();
407         }
408         return lokalVonScheinbar(Math.floor(approx) + (morgen ? .25 - offset : .75 + offset));
409     }
410     
411     public double universalVonLokal(double tl, Ort ort) {
412         return tl - ort.getLaenge() / (double) 360;
413     }
414
415     public double standardVonLokal(double tl, Ort ort) {
416         return standardVonUniversal(universalVonLokal(tl, ort), ort);
417     }
418     
419     public double lokalVonScheinbar(double t) {
420         return t - zeitgleichung(t);
421     }
422
423     public double zeitgleichung(double t) {
424         double c = julJahrhunderte(t);
425         double laenge = poly(c, et.koeffLaenge);
426         double anomalie = poly(c, et.koeffAnomalie);
427         double exzentrizitaet = poly(c, et.koeffExzentrizitaet);
d20d98 428         double varepsilon = schiefstand(t);
66d68b 429         double y = quadrat(tanGrad(varepsilon / 2));
U 430         double equation = (1d / 2d / Math.PI) * (y * sinGrad(2 * laenge) +
431         -2 * exzentrizitaet * sinGrad(anomalie) +
432         4 * exzentrizitaet * y * sinGrad(anomalie) * kosGrad(2 * laenge) +
433         -0.5 * y * y * sinGrad(4 * laenge) +
434         -1.25 * exzentrizitaet * exzentrizitaet * sinGrad(2 * anomalie));
d20d98 435         return vorzeichen(equation) * Math.min(Math.abs(equation), stunde(12));
66d68b 436     }
U 437     private static class et {
438         private static final double[] koeffLaenge = new double[] {280.46645, 36000.76983, 0.0003032};
439         private static final double[] koeffAnomalie = new double[] {357.52910, 35999.05030, -0.0001559, -0.00000048};
440         private static final double[] koeffExzentrizitaet = new double[] {0.016708617, -0.000042037, -0.0000001236};
441     }
442
d20d98 443     public double schiefstand(double t) {
66d68b 444         double c = julJahrhunderte(t);
U 445         return winkel(23, 26, 21.448) + poly(c, coeffObliquity);
446     }
447     private final double[] coeffObliquity = new double[] {0, winkel(0, 0, -46.8150), winkel(0, 0, -0.00059), winkel(0, 0, 0.001813)};
448   
d20d98 449     public int vorzeichen(double x) {
66d68b 450         if(x < 0)
U 451             return -1;
452         else if(x > 0)
453             return 1;
454         else
455             return 0;
456     }
457     
458     public double quadrat(double x) {
459         return x * x;
460     }
461
462     public double kosGrad(double theta) {
463         return Math.cos(gradZuBogenmass(theta));
464     }
465
466     public double arcSinGrad(double x) {
467         return bogenmassZuGrad(Math.asin(x));
468     }
469     
470     public double tanGrad(double theta) {
471         return Math.tan(gradZuBogenmass(theta));
472     }
473
474     public double arcKosGrad(double x) {
475         return bogenmassZuGrad(Math.acos(x));
476     }
477     
478     public double bogenmassZuGrad(double theta) {
479         return grad(theta / Math.PI * 180);
480     }
481
482     public double winkel(double d, double m, double s) {
483         return d + (m + s / 60) / 60;
484     }
485   
486   /* ---------------- Jahreszeiten ----- */
487
d20d98 488   /** durchschnittliche Dauer eines Umlaufs der Erde um die Sonne in Tagen */
U 489     public static final double MITTLERES_TROPISCHES_JAHR = 365.242189;
66d68b 490     
U 491     public double standardVonUniversal(double t, Ort ort) {
492         return t + ort.getZeitzone() / 24;
493     }
494   
495     public double solareLaengeNach(double t, double phi) {
496         double varepsilon = 0.00001;
d20d98 497     double rate = MITTLERES_TROPISCHES_JAHR / (double) 360;
66d68b 498         double tau = t + rate * modulo(phi - solareLaenge(t), 360);
U 499         double l = Math.max(t, tau - 5);
500         double u = tau + 5;
501         
502         double lo = l, hi = u, x = (hi + lo) / 2;
503         while(hi - lo > varepsilon) {
504             if(modulo(solareLaenge(x) - phi, 360) < (double) 180)
505                 hi = x;
506             else
507                 lo = x;
508
509             x = (hi + lo) / 2;
510         }
511         return x;
512     }
513   
514     public double solareLaenge(double t) {
515         double c = julJahrhunderte(t);
516         double sigma = 0;
517         for(int i = 0; i < sLaenge.koeffizienten.length; ++i) {
518             sigma += sLaenge.koeffizienten[i] * sinGrad(sLaenge.multiplikatoren[i] * c + sLaenge.summanden[i]);
519         }
520     double laenge = (double) 282.7771834 +            
521             36000.76953744 * c +
522             0.000005729577951308232 * sigma;
523         return modulo(laenge + aberration(t) + nutation(t), 360);
524     }
525
d20d98 526     public double geschaetzteSolareLaengeVor(double tee, double phi) {
U 527         double rate = MITTLERES_TROPISCHES_JAHR / (double) (360);
528         double tau = tee - rate * modulo(solareLaenge(tee) - phi, 360);
529         double capDelta = modulo(solareLaenge(tau) - phi + (double) (180), 360) - (double) (180);
530         return Math.min(tee, tau - rate * capDelta);
531     }
532
66d68b 533     public double julJahrhunderte(double t) {
U 534         return (dynamischVonUniversal(t) - j2000()) / 36525;
535     }
536
537     public double dynamischVonUniversal(double tee) {
538         return tee + ephemeridenKorrektur(tee);
539     }
540
541     public double ephemeridenKorrektur(double t) {
542     double[] koeffizient17tes = new double[] {196.58333, -4.0675, 0.0219167};
543         double[] koeffizient19tes = new double[] {-0.00002, 0.000297, 0.025184, -0.181133, 0.553040, -0.861938, 0.677066, -0.212591};
544         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};
545     ISOKalender w = new ISOKalender();
546         long jahr = w.jahrVonTagen((long)Math.floor(t));
547         double c = differenz(w.zuTagen(1900, Definition.JANUAR, 1), 
548             w.zuTagen(jahr, Definition.JULI, 1)) / 36525d;
549         double result;
550         if(1988 <= jahr && jahr <= 2019) {
551             result = (jahr - 1933) / (24d * 60 * 60);
552         } else if (1900 <= jahr && jahr <= 1987) {
553             result = poly(c, koeffizient19tes);
554         } else if (1800 <= jahr && jahr <= 1899) {
555             result = poly(c, koeffizient18tes);
556         } else if (1620 <= jahr && jahr <= 1799) {
557             result = poly(jahr - 1600, koeffizient17tes) / (24 * 60 * 60);
558         } else {
559             double x = stunde(12) + differenz(w.zuTagen(1810, Definition.JANUAR, 1), 
560               w.zuTagen(jahr, Definition.JANUAR, 1));
561             return (x * x / 41048480 - 15) / (24 * 60 * 60);
562         }
563         return result;
564     }
565   
566     public double nutation(double t) {
567         double[] koeffa = new double[] {124.90, -1934.134, 0.002063};
568         double[] koeffb = new double[] {201.11, 72001.5377, 0.00057};
569         double c = julJahrhunderte(t);
570         double capA = poly(c, koeffa);
571         double capB = poly(c, koeffb);
572         return (double) -0.004778 * sinGrad(capA) +
573             (double) -0.0003667 * sinGrad(capB);
574     }
575
576     public static long differenz(long datum1, long datum2) {
577         return datum2 - datum1;
578     }
579   
580     public double stunde(double x) {
581         return x / 24;
582     }
583   
584   public double j2000() {
585     ISOKalender w = new ISOKalender();
586     return stunde(12) + w.zuTagen(2000, Definition.JANUAR, 1);
587   }
588   
589     public double sinGrad(double theta) {
590         return Math.sin(gradZuBogenmass(theta));
591     }
592
593     public double gradZuBogenmass(double theta) {
594         return grad(theta) * Math.PI / 180;
595     }
596
597     public double grad(double theta) {
598         return modulo(theta, 360);
599     }
600
601     public double aberration(double t) {
602         double c = julJahrhunderte(t);
603     return (double) 0.0000974 * gradKosinus((double) 177.63 + (double) 35999.01848 * c) - (double) 0.005575;
604     }
605
606     public double gradKosinus(double theta) {
607         return Math.cos(gradZuBogenmass(theta));
608     }
609
610     public double poly(double x, double[] a) {
611         double ergebnis = a[0];
612         for(int i = 1; i < a.length; ++i) {
613             ergebnis += a[i] * Math.pow(x, i);
614         }
615         return ergebnis;
616     }
617   
618     private static class sLaenge {
619         private static final int[] koeffizienten = new int[] {
620             403406, 195207, 119433, 112392, 3891, 2819, 1721,
621             660, 350, 334, 314, 268, 242, 234, 158, 132, 129, 114,
622             99, 93, 86, 78, 72, 68, 64, 46, 38, 37, 32, 29, 28, 27, 27,
623             25, 24, 21, 21, 20, 18, 17, 14, 13, 13, 13, 12, 10, 10, 10,
624             10
625         };
626         private static final double[] multiplikatoren = new double[] {
627             0.9287892, 35999.1376958, 35999.4089666,
628             35998.7287385, 71998.20261, 71998.4403,
629             36000.35726, 71997.4812, 32964.4678,
630             -19.4410, 445267.1117, 45036.8840, 3.1008,
631             22518.4434, -19.9739, 65928.9345,
632             9038.0293, 3034.7684, 33718.148, 3034.448,
633             -2280.773, 29929.992, 31556.493, 149.588,
634             9037.750, 107997.405, -4444.176, 151.771,
635             67555.316, 31556.080, -4561.540,
636             107996.706, 1221.655, 62894.167,
637             31437.369, 14578.298, -31931.757,
638             34777.243, 1221.999, 62894.511,
639             -4442.039, 107997.909, 119.066, 16859.071,
640             -4.578, 26895.292, -39.127, 12297.536,
641             90073.778
642         };
643         private static final double[] summanden = new double[] {
644             270.54861, 340.19128, 63.91854, 331.26220,
645             317.843, 86.631, 240.052, 310.26, 247.23,
646             260.87, 297.82, 343.14, 166.79, 81.53,
647             3.50, 132.75, 182.95, 162.03, 29.8,
648             266.4, 249.2, 157.6, 257.8, 185.1, 69.9,
649             8.0, 197.1, 250.4, 65.3, 162.7, 341.5,
650             291.6, 98.5, 146.7, 110.0, 5.2, 342.6,
651             230.9, 256.1, 45.3, 242.9, 115.2, 151.8,
652             285.3, 53.3, 126.6, 205.7, 85.9,
653             146.1
654         };
655     } 
656 }