netcore/System.Private.CoreLib/shared/System/Globalization/CalendricalCalculationsHelper.cs

   1 // Licensed to the .NET Foundation under one or more agreements.
   2 // The .NET Foundation licenses this file to you under the MIT license.
   3 // See the LICENSE file in the project root for more information.
   4
   5 using System.Diagnostics;
   6
   7 namespace System.Globalization
   8 {
   9     internal static class CalendricalCalculationsHelper
  10     {
  11         private const double FullCircleOfArc = 360.0; // 360.0;
  12         private const int HalfCircleOfArc = 180;
  13         private const double TwelveHours = 0.5; // half a day
  14         private const double Noon2000Jan01 = 730120.5;
  15         internal const double MeanTropicalYearInDays = 365.242189;
  16         private const double MeanSpeedOfSun = MeanTropicalYearInDays / FullCircleOfArc;
  17         private const double LongitudeSpring = 0.0;
  18         private const double TwoDegreesAfterSpring = 2.0;
  19         private const int SecondsPerDay = 24 * 60 * 60; // 24 hours * 60 minutes * 60 seconds
  20
  21         private const int DaysInUniformLengthCentury = 36525;
  22         private const int SecondsPerMinute = 60;
  23         private const int MinutesPerDegree = 60;
  24
  25         private static readonly long s_startOf1810 = GetNumberOfDays(new DateTime(1810, 1, 1));
  26         private static readonly long s_startOf1900Century = GetNumberOfDays(new DateTime(1900, 1, 1));
  27
  28         private static readonly double[] s_coefficients1900to1987 = new double[] { -0.00002, 0.000297, 0.025184, -0.181133, 0.553040, -0.861938, 0.677066, -0.212591 };
  29         private static readonly double[] s_coefficients1800to1899 = new double[] { -0.000009, 0.003844, 0.083563, 0.865736, 4.867575, 15.845535, 31.332267, 38.291999, 28.316289, 11.636204, 2.043794 };
  30         private static readonly double[] s_coefficients1700to1799 = new double[] { 8.118780842, -0.005092142, 0.003336121, -0.0000266484 };
  31         private static readonly double[] s_coefficients1620to1699 = new double[] { 196.58333, -4.0675, 0.0219167 };
  32         private static readonly double[] s_lambdaCoefficients = new double[] { 280.46645, 36000.76983, 0.0003032 };
  33         private static readonly double[] s_anomalyCoefficients = new double[] { 357.52910, 35999.05030, -0.0001559, -0.00000048 };
  34         private static readonly double[] s_eccentricityCoefficients = new double[] { 0.016708617, -0.000042037, -0.0000001236 };
  35         private static readonly double[] s_coefficients = new double[] { Angle(23, 26, 21.448), Angle(0, 0, -46.8150), Angle(0, 0, -0.00059), Angle(0, 0, 0.001813) };
  36         private static readonly double[] s_coefficientsA = new double[] { 124.90, -1934.134, 0.002063 };
  37         private static readonly double[] s_coefficientsB = new double[] { 201.11, 72001.5377, 0.00057 };
  38
  39         private static double RadiansFromDegrees(double degree)
  40         {
  41             return degree * Math.PI / 180;
  42         }
  43
  44         private static double SinOfDegree(double degree)
  45         {
  46             return Math.Sin(RadiansFromDegrees(degree));
  47         }
  48
  49         private static double CosOfDegree(double degree)
  50         {
  51             return Math.Cos(RadiansFromDegrees(degree));
  52         }
  53         private static double TanOfDegree(double degree)
  54         {
  55             return Math.Tan(RadiansFromDegrees(degree));
  56         }
  57
  58         public static double Angle(int degrees, int minutes, double seconds)
  59         {
  60             return ((seconds / SecondsPerMinute + minutes) / MinutesPerDegree) + degrees;
  61         }
  62
  63         private static double Obliquity(double julianCenturies)
  64         {
  65             return PolynomialSum(s_coefficients, julianCenturies);
  66         }
  67
  68         internal static long GetNumberOfDays(DateTime date)
  69         {
  70             return date.Ticks / GregorianCalendar.TicksPerDay;
  71         }
  72
  73         private static int GetGregorianYear(double numberOfDays)
  74         {
  75             return new DateTime(Math.Min((long)(Math.Floor(numberOfDays) * GregorianCalendar.TicksPerDay), DateTime.MaxValue.Ticks)).Year;
  76         }
  77
  78         private enum CorrectionAlgorithm
  79         {
  80             Default,
  81             Year1988to2019,
  82             Year1900to1987,
  83             Year1800to1899,
  84             Year1700to1799,
  85             Year1620to1699
  86         }
  87
  88         private struct EphemerisCorrectionAlgorithmMap
  89         {
  90             public EphemerisCorrectionAlgorithmMap(int year, CorrectionAlgorithm algorithm)
  91             {
  92                 _lowestYear = year;
  93                 _algorithm = algorithm;
  94             }
  95
  96             internal int _lowestYear;
  97             internal CorrectionAlgorithm _algorithm;
  98         }
  99
 100         private static readonly EphemerisCorrectionAlgorithmMap[] s_ephemerisCorrectionTable = new EphemerisCorrectionAlgorithmMap[]
 101         {
 102             // lowest year that starts algorithm, algorithm to use
 103             new EphemerisCorrectionAlgorithmMap(2020, CorrectionAlgorithm.Default),
 104             new EphemerisCorrectionAlgorithmMap(1988, CorrectionAlgorithm.Year1988to2019),
 105             new EphemerisCorrectionAlgorithmMap(1900, CorrectionAlgorithm.Year1900to1987),
 106             new EphemerisCorrectionAlgorithmMap(1800, CorrectionAlgorithm.Year1800to1899),
 107             new EphemerisCorrectionAlgorithmMap(1700, CorrectionAlgorithm.Year1700to1799),
 108             new EphemerisCorrectionAlgorithmMap(1620, CorrectionAlgorithm.Year1620to1699),
 109             new EphemerisCorrectionAlgorithmMap(int.MinValue, CorrectionAlgorithm.Default) // default must be last
 110         };
 111
 112         private static double Reminder(double divisor, double dividend)
 113         {
 114             double whole = Math.Floor(divisor / dividend);
 115             return divisor - (dividend * whole);
 116         }
 117
 118         private static double NormalizeLongitude(double longitude)
 119         {
 120             longitude = Reminder(longitude, FullCircleOfArc);
 121             if (longitude < 0)
 122             {
 123                 longitude += FullCircleOfArc;
 124             }
 125             return longitude;
 126         }
 127
 128         public static double AsDayFraction(double longitude)
 129         {
 130             return longitude / FullCircleOfArc;
 131         }
 132
 133         private static double PolynomialSum(double[] coefficients, double indeterminate)
 134         {
 135             double sum = coefficients[0];
 136             double indeterminateRaised = 1;
 137             for (int i = 1; i < coefficients.Length; i++)
 138             {
 139                 indeterminateRaised *= indeterminate;
 140                 sum += (coefficients[i] * indeterminateRaised);
 141             }
 142
 143             return sum;
 144         }
 145
 146         private static double CenturiesFrom1900(int gregorianYear)
 147         {
 148             long july1stOfYear = GetNumberOfDays(new DateTime(gregorianYear, 7, 1));
 149             return (double)(july1stOfYear - s_startOf1900Century) / DaysInUniformLengthCentury;
 150         }
 151
 152         // the following formulas defines a polynomial function which gives us the amount that the earth is slowing down for specific year ranges
 153         private static double DefaultEphemerisCorrection(int gregorianYear)
 154         {
 155             Debug.Assert(gregorianYear < 1620 || 2020 <= gregorianYear);
 156             long january1stOfYear = GetNumberOfDays(new DateTime(gregorianYear, 1, 1));
 157             double daysSinceStartOf1810 = january1stOfYear - s_startOf1810;
 158             double x = TwelveHours + daysSinceStartOf1810;
 159             return ((Math.Pow(x, 2) / 41048480) - 15) / SecondsPerDay;
 160         }
 161
 162         private static double EphemerisCorrection1988to2019(int gregorianYear)
 163         {
 164             Debug.Assert(1988 <= gregorianYear && gregorianYear <= 2019);
 165             return (double)(gregorianYear - 1933) / SecondsPerDay;
 166         }
 167
 168         private static double EphemerisCorrection1900to1987(int gregorianYear)
 169         {
 170             Debug.Assert(1900 <= gregorianYear && gregorianYear <= 1987);
 171             double centuriesFrom1900 = CenturiesFrom1900(gregorianYear);
 172             return PolynomialSum(s_coefficients1900to1987, centuriesFrom1900);
 173         }
 174
 175         private static double EphemerisCorrection1800to1899(int gregorianYear)
 176         {
 177             Debug.Assert(1800 <= gregorianYear && gregorianYear <= 1899);
 178             double centuriesFrom1900 = CenturiesFrom1900(gregorianYear);
 179             return PolynomialSum(s_coefficients1800to1899, centuriesFrom1900);
 180         }
 181
 182         private static double EphemerisCorrection1700to1799(int gregorianYear)
 183         {
 184             Debug.Assert(1700 <= gregorianYear && gregorianYear <= 1799);
 185             double yearsSince1700 = gregorianYear - 1700;
 186             return PolynomialSum(s_coefficients1700to1799, yearsSince1700) / SecondsPerDay;
 187         }
 188
 189         private static double EphemerisCorrection1620to1699(int gregorianYear)
 190         {
 191             Debug.Assert(1620 <= gregorianYear && gregorianYear <= 1699);
 192             double yearsSince1600 = gregorianYear - 1600;
 193             return PolynomialSum(s_coefficients1620to1699, yearsSince1600) / SecondsPerDay;
 194         }
 195
 196         // ephemeris-correction: correction to account for the slowing down of the rotation of the earth
 197         private static double EphemerisCorrection(double time)
 198         {
 199             int year = GetGregorianYear(time);
 200             foreach (EphemerisCorrectionAlgorithmMap map in s_ephemerisCorrectionTable)
 201             {
 202                 if (map._lowestYear <= year)
 203                 {
 204                     switch (map._algorithm)
 205                     {
 206                         case CorrectionAlgorithm.Default: return DefaultEphemerisCorrection(year);
 207                         case CorrectionAlgorithm.Year1988to2019: return EphemerisCorrection1988to2019(year);
 208                         case CorrectionAlgorithm.Year1900to1987: return EphemerisCorrection1900to1987(year);
 209                         case CorrectionAlgorithm.Year1800to1899: return EphemerisCorrection1800to1899(year);
 210                         case CorrectionAlgorithm.Year1700to1799: return EphemerisCorrection1700to1799(year);
 211                         case CorrectionAlgorithm.Year1620to1699: return EphemerisCorrection1620to1699(year);
 212                     }
 213
 214                     break; // break the loop and assert eventually
 215                 }
 216             }
 217
 218             Debug.Fail("Not expected to come here");
 219             return DefaultEphemerisCorrection(year);
 220         }
 221
 222         public static double JulianCenturies(double moment)
 223         {
 224             double dynamicalMoment = moment + EphemerisCorrection(moment);
 225             return (dynamicalMoment - Noon2000Jan01) / DaysInUniformLengthCentury;
 226         }
 227
 228         private static bool IsNegative(double value)
 229         {
 230             return Math.Sign(value) == -1;
 231         }
 232
 233         private static double CopySign(double value, double sign)
 234         {
 235             return (IsNegative(value) == IsNegative(sign)) ? value : -value;
 236         }
 237
 238         // equation-of-time; approximate the difference between apparent solar time and mean solar time
 239         // formal definition is EOT = GHA - GMHA
 240         // GHA is the Greenwich Hour Angle of the apparent (actual) Sun
 241         // GMHA is the Greenwich Mean Hour Angle of the mean (fictitious) Sun
 242         // http://www.esrl.noaa.gov/gmd/grad/solcalc/
 243         // http://en.wikipedia.org/wiki/Equation_of_time
 244         private static double EquationOfTime(double time)
 245         {
 246             double julianCenturies = JulianCenturies(time);
 247             double lambda = PolynomialSum(s_lambdaCoefficients, julianCenturies);
 248             double anomaly = PolynomialSum(s_anomalyCoefficients, julianCenturies);
 249             double eccentricity = PolynomialSum(s_eccentricityCoefficients, julianCenturies);
 250
 251             double epsilon = Obliquity(julianCenturies);
 252             double tanHalfEpsilon = TanOfDegree(epsilon / 2);
 253             double y = tanHalfEpsilon * tanHalfEpsilon;
 254
 255             double dividend = ((y * SinOfDegree(2 * lambda))
 256                 - (2 * eccentricity * SinOfDegree(anomaly))
 257                 + (4 * eccentricity * y * SinOfDegree(anomaly) * CosOfDegree(2 * lambda))
 258                 - (0.5 * Math.Pow(y, 2) * SinOfDegree(4 * lambda))
 259                 - (1.25 * Math.Pow(eccentricity, 2) * SinOfDegree(2 * anomaly)));
 260             const double Divisor = 2 * Math.PI;
 261             double equation = dividend / Divisor;
 262
 263             // approximation of equation of time is not valid for dates that are many millennia in the past or future
 264             // thus limited to a half day
 265             return CopySign(Math.Min(Math.Abs(equation), TwelveHours), equation);
 266         }
 267
 268         private static double AsLocalTime(double apparentMidday, double longitude)
 269         {
 270             // slightly inaccurate since equation of time takes mean time not apparent time as its argument, but the difference is negligible
 271             double universalTime = apparentMidday - AsDayFraction(longitude);
 272             return apparentMidday - EquationOfTime(universalTime);
 273         }
 274
 275         // midday
 276         public static double Midday(double date, double longitude)
 277         {
 278             return AsLocalTime(date + TwelveHours, longitude) - AsDayFraction(longitude);
 279         }
 280
 281         private static double InitLongitude(double longitude)
 282         {
 283             return NormalizeLongitude(longitude + HalfCircleOfArc) - HalfCircleOfArc;
 284         }
 285
 286         // midday-in-tehran
 287         public static double MiddayAtPersianObservationSite(double date)
 288         {
 289             return Midday(date, InitLongitude(52.5)); // 52.5 degrees east - longitude of UTC+3:30 which defines Iranian Standard Time
 290         }
 291
 292         private static double PeriodicTerm(double julianCenturies, int x, double y, double z)
 293         {
 294             return x * SinOfDegree(y + z * julianCenturies);
 295         }
 296
 297         private static double SumLongSequenceOfPeriodicTerms(double julianCenturies)
 298         {
 299             double sum = 0.0;
 300             sum += PeriodicTerm(julianCenturies, 403406, 270.54861, 0.9287892);
 301             sum += PeriodicTerm(julianCenturies, 195207, 340.19128, 35999.1376958);
 302             sum += PeriodicTerm(julianCenturies, 119433, 63.91854, 35999.4089666);
 303             sum += PeriodicTerm(julianCenturies, 112392, 331.2622, 35998.7287385);
 304             sum += PeriodicTerm(julianCenturies, 3891, 317.843, 71998.20261);
 305             sum += PeriodicTerm(julianCenturies, 2819, 86.631, 71998.4403);
 306             sum += PeriodicTerm(julianCenturies, 1721, 240.052, 36000.35726);
 307             sum += PeriodicTerm(julianCenturies, 660, 310.26, 71997.4812);
 308             sum += PeriodicTerm(julianCenturies, 350, 247.23, 32964.4678);
 309             sum += PeriodicTerm(julianCenturies, 334, 260.87, -19.441);
 310             sum += PeriodicTerm(julianCenturies, 314, 297.82, 445267.1117);
 311             sum += PeriodicTerm(julianCenturies, 268, 343.14, 45036.884);
 312             sum += PeriodicTerm(julianCenturies, 242, 166.79, 3.1008);
 313             sum += PeriodicTerm(julianCenturies, 234, 81.53, 22518.4434);
 314             sum += PeriodicTerm(julianCenturies, 158, 3.5, -19.9739);
 315             sum += PeriodicTerm(julianCenturies, 132, 132.75, 65928.9345);
 316             sum += PeriodicTerm(julianCenturies, 129, 182.95, 9038.0293);
 317             sum += PeriodicTerm(julianCenturies, 114, 162.03, 3034.7684);
 318             sum += PeriodicTerm(julianCenturies, 99, 29.8, 33718.148);
 319             sum += PeriodicTerm(julianCenturies, 93, 266.4, 3034.448);
 320             sum += PeriodicTerm(julianCenturies, 86, 249.2, -2280.773);
 321             sum += PeriodicTerm(julianCenturies, 78, 157.6, 29929.992);
 322             sum += PeriodicTerm(julianCenturies, 72, 257.8, 31556.493);
 323             sum += PeriodicTerm(julianCenturies, 68, 185.1, 149.588);
 324             sum += PeriodicTerm(julianCenturies, 64, 69.9, 9037.75);
 325             sum += PeriodicTerm(julianCenturies, 46, 8.0, 107997.405);
 326             sum += PeriodicTerm(julianCenturies, 38, 197.1, -4444.176);
 327             sum += PeriodicTerm(julianCenturies, 37, 250.4, 151.771);
 328             sum += PeriodicTerm(julianCenturies, 32, 65.3, 67555.316);
 329             sum += PeriodicTerm(julianCenturies, 29, 162.7, 31556.08);
 330             sum += PeriodicTerm(julianCenturies, 28, 341.5, -4561.54);
 331             sum += PeriodicTerm(julianCenturies, 27, 291.6, 107996.706);
 332             sum += PeriodicTerm(julianCenturies, 27, 98.5, 1221.655);
 333             sum += PeriodicTerm(julianCenturies, 25, 146.7, 62894.167);
 334             sum += PeriodicTerm(julianCenturies, 24, 110.0, 31437.369);
 335             sum += PeriodicTerm(julianCenturies, 21, 5.2, 14578.298);
 336             sum += PeriodicTerm(julianCenturies, 21, 342.6, -31931.757);
 337             sum += PeriodicTerm(julianCenturies, 20, 230.9, 34777.243);
 338             sum += PeriodicTerm(julianCenturies, 18, 256.1, 1221.999);
 339             sum += PeriodicTerm(julianCenturies, 17, 45.3, 62894.511);
 340             sum += PeriodicTerm(julianCenturies, 14, 242.9, -4442.039);
 341             sum += PeriodicTerm(julianCenturies, 13, 115.2, 107997.909);
 342             sum += PeriodicTerm(julianCenturies, 13, 151.8, 119.066);
 343             sum += PeriodicTerm(julianCenturies, 13, 285.3, 16859.071);
 344             sum += PeriodicTerm(julianCenturies, 12, 53.3, -4.578);
 345             sum += PeriodicTerm(julianCenturies, 10, 126.6, 26895.292);
 346             sum += PeriodicTerm(julianCenturies, 10, 205.7, -39.127);
 347             sum += PeriodicTerm(julianCenturies, 10, 85.9, 12297.536);
 348             sum += PeriodicTerm(julianCenturies, 10, 146.1, 90073.778);
 349             return sum;
 350         }
 351
 352         private static double Aberration(double julianCenturies)
 353         {
 354             return (0.0000974 * CosOfDegree(177.63 + (35999.01848 * julianCenturies))) - 0.005575;
 355         }
 356
 357         private static double Nutation(double julianCenturies)
 358         {
 359             double a = PolynomialSum(s_coefficientsA, julianCenturies);
 360             double b = PolynomialSum(s_coefficientsB, julianCenturies);
 361             return (-0.004778 * SinOfDegree(a)) - (0.0003667 * SinOfDegree(b));
 362         }
 363
 364         public static double Compute(double time)
 365         {
 366             double julianCenturies = JulianCenturies(time);
 367             double lambda = 282.7771834
 368                 + (36000.76953744 * julianCenturies)
 369                 + (0.000005729577951308232 * SumLongSequenceOfPeriodicTerms(julianCenturies));
 370
 371             double longitude = lambda + Aberration(julianCenturies) + Nutation(julianCenturies);
 372             return InitLongitude(longitude);
 373         }
 374
 375         public static double AsSeason(double longitude)
 376         {
 377             return (longitude < 0) ? (longitude + FullCircleOfArc) : longitude;
 378         }
 379
 380         private static double EstimatePrior(double longitude, double time)
 381         {
 382             double timeSunLastAtLongitude = time - (MeanSpeedOfSun * AsSeason(InitLongitude(Compute(time) - longitude)));
 383             double longitudeErrorDelta = InitLongitude(Compute(timeSunLastAtLongitude) - longitude);
 384             return Math.Min(time, timeSunLastAtLongitude - (MeanSpeedOfSun * longitudeErrorDelta));
 385         }
 386
 387         // persian-new-year-on-or-before
 388         //  number of days is the absolute date. The absolute date is the number of days from January 1st, 1 A.D.
 389         //  1/1/0001 is absolute date 1.
 390         internal static long PersianNewYearOnOrBefore(long numberOfDays)
 391         {
 392             double date = (double)numberOfDays;
 393
 394             double approx = EstimatePrior(LongitudeSpring, MiddayAtPersianObservationSite(date));
 395             long lowerBoundNewYearDay = (long)Math.Floor(approx) - 1;
 396             long upperBoundNewYearDay = lowerBoundNewYearDay + 3; // estimate is generally within a day of the actual occurrence (at the limits, the error expands, since the calculations rely on the mean tropical year which changes...)
 397             long day = lowerBoundNewYearDay;
 398             for (; day != upperBoundNewYearDay; ++day)
 399             {
 400                 double midday = MiddayAtPersianObservationSite((double)day);
 401                 double l = Compute(midday);
 402                 if ((LongitudeSpring <= l) && (l <= TwoDegreesAfterSpring))
 403                 {
 404                     break;
 405                 }
 406             }
 407             Debug.Assert(day != upperBoundNewYearDay);
 408
 409             return day - 1;
 410         }
 411     }
 412 }