C#: Fixed bugs with calendar dates with extreme year values.

Fixed problems converting AstroTime to calendar dates and back.
Also expose struct CalendarDateTime to outside callers,
for convenience dealing with Gregorian calendar dates.

generate/dotnet/csharp_test/csharp_test.cs

@@ -182,7 +182,29 @@ namespace csharp_test
             return line.Split(TokenSeparators, StringSplitOptions.RemoveEmptyEntries);
         }
 
-        static int TestTime()
+        static int CalendarCase(int year, int month, int day, int hour, int minute, double second)
+        {
+            // Convert calendar date to time.
+            var time = new AstroTime(year, month, day, hour, minute, second);
+
+            // Convert time back to calendar date.
+            var cal = time.ToCalendarDateTime();
+
+            // Verify the round-trip was correct.
+            if (cal.year != year || cal.month != month || cal.day != day)
+                return Fail($"{nameof(CalendarCase)}: expected [{year:0000}-{month:00}-{day:00}] but found [{cal.year:0000}-{cal.month:00}-{cal.day:00}]");
+
+            // Be a little more tolerant with time because of roundoff errors.
+            double expectedMillis = 1000.0 * (second + 60*(minute + 60*hour));
+            double calcMillis = 1000.0 * (cal.second + 60*(cal.minute + 60*cal.hour));
+            double diffMillis = abs(calcMillis - expectedMillis);
+            if (diffMillis > 1.0)
+                return Fail($"{nameof(CalendarCase)}: EXCESSIVE time error = {diffMillis} milliseconds.");
+
+            return 0;
+        }
+
+        static int UtDayValueTest()
         {
             const int year = 2018;
             const int month = 12;
@@ -194,13 +216,13 @@ namespace csharp_test
 
             DateTime d = new DateTime(year, month, day, hour, minute, second, milli, DateTimeKind.Utc);
             AstroTime time = new AstroTime(d);
-            Console.WriteLine("C# TestTime: text={0}, ut={1}, tt={2}", time.ToString(), time.ut.ToString("F6"), time.tt.ToString("F6"));
+            Debug("C# UtDayValueTest: text={0}, ut={1}, tt={2}", time, time.ut.ToString("F6"), time.tt.ToString("F6"));
 
             const double expected_ut = 6910.270978506945;
             double diff = time.ut - expected_ut;
             if (abs(diff) > 1.0e-12)
             {
-                Console.WriteLine("C# TestTime: ERROR - excessive UT error {0}", diff);
+                Console.WriteLine("C# UtDayValueTest: ERROR - excessive UT error {0}", diff);
                 return 1;
             }
 
@@ -208,20 +230,33 @@ namespace csharp_test
             diff = time.tt - expected_tt;
             if (abs(diff) > 1.0e-12)
             {
-                Console.WriteLine("C# TestTime: ERROR - excessive TT error {0}", diff);
-                return 1;
-            }
-
-            DateTime utc = time.ToUtcDateTime();
-            if (utc.Year != year || utc.Month != month || utc.Day != day || utc.Hour != hour || utc.Minute != minute || utc.Second != second || utc.Millisecond != milli)
-            {
-                Console.WriteLine("C# TestTime: ERROR - Expected {0:o}, found {1:o}", d, utc);
+                Console.WriteLine("C# UtDayValueTest: ERROR - excessive TT error {0}", diff);
                 return 1;
             }
 
             return 0;
         }
 
+        static int TestTime()
+        {
+            if (0 != UtDayValueTest())
+                return 1;
+
+            for (int year = 2400; year >= -5400; --year)
+                if (0 != CalendarCase(year, 12,  2, 18, 30, 12.543))
+                    return 1;
+
+            for (int year = -999999; year <= -999599; ++year)
+                if (0 != CalendarCase(year, 12,  2, 18, 30, 12.543))
+                    return 1;
+
+            for (int year = +999599; year <= +999999; ++year)
+                if (0 != CalendarCase(year, 12,  2, 18, 30, 12.543))
+                    return 1;
+
+            return Pass(nameof(TestTime));
+        }
+
         static int MoonTest()
         {
             var time = new AstroTime(2019, 6, 24, 15, 45, 37);

generate/template/astronomy.c

@@ -1097,7 +1097,7 @@ astro_utc_t Astronomy_UtcFromTime(astro_time_t time)
     astro_utc_t utc;
     long jd, k, m, n;
     double djd, x;
-    const long c = 2500L;
+    const long c = 2500;
 
     djd = time.ut + 2451545.5;
     jd = (long)floor(djd);
@@ -1120,18 +1120,18 @@ astro_utc_t Astronomy_UtcFromTime(astro_time_t time)
         Any multiple of 400 years has the same number of days,
         because it eliminates all the special cases for leap years.
     */
-    k = jd + (68569L + c*146097L);
-    n = 4L * k / 146097L;
-    k = k - (146097L * n + 3L) / 4L;
-    m = 4000L * (k + 1L) / 1461001L;
-    k = k - 1461L * m / 4L + 31L;
+    k = jd + (68569 + c*146097);
+    n = (4 * k) / 146097;
+    k = k - (146097 * n + 3)/4;
+    m = (4000 * (k+1)) / 1461001;
+    k = k - (1461 * m)/4 + 31;
 
-    utc.month = (int) (80L * k / 2447L);
-    utc.day = (int) (k - 2447L * (long)utc.month / 80L);
-    k = (long) utc.month / 11L;
+    utc.month = (int) ((80 * k) / 2447);
+    utc.day = (int) (k - (2447*utc.month)/80);
+    k = utc.month / 11;
 
-    utc.month = (int) ((long)utc.month + 2L - 12L * k);
-    utc.year = (int) (100L * (n - 49L) + m + k - 400L*c);
+    utc.month = (int) (utc.month + 2 - 12*k);
+    utc.year = (int) (100 * (n - 49) + m + k - 400*c);
 
     return utc;
 }

generate/template/astronomy.cs

@@ -324,23 +324,22 @@ namespace CosineKitty
             return Origin.AddDays(ut).ToUniversalTime();
         }
 
+        /// <summary>
+        /// Converts this object to our custom type #CalendarDateTime.
+        /// </summary>
+        /// <remarks>
+        /// The .NET type `DateTime` can only represent years in the range 0000..9999.
+        /// However, the Astronomy Engine type #CalendarDateTime can represent
+        /// years in the range -999999..+999999. This is a time span of nearly 2 million years.
+        /// This function converts this `AstroTime` object to an equivalent Gregorian calendar representation.
+        /// </remarks>
+        public CalendarDateTime ToCalendarDateTime() => new CalendarDateTime(ut);
+
         /// <summary>
         /// Converts this `AstroTime` to ISO 8601 format, expressed in UTC with millisecond resolution.
         /// </summary>
         /// <returns>Example: "2019-08-30T17:45:22.763Z".</returns>
-        public override string ToString()
-        {
-            var d = new CalendarDateTime(ut);
-            int millis = Math.Max(0, Math.Min(59999, (int)Math.Round(d.second * 1000.0)));
-            string y;
-            if (d.year < 0)
-                y = "-" + (-d.year).ToString("000000");
-            else if (d.year <= 9999)
-                y = d.year.ToString("0000");
-            else
-                y = "+" + d.year.ToString("000000");
-            return $"{y}-{d.month:00}-{d.day:00}T{d.hour:00}:{d.minute:00}:{millis/1000:00}.{millis%1000:000}Z";
-        }
+        public override string ToString() => ToCalendarDateTime().ToString();
 
         private static Regex re = new Regex(
             @"^
@@ -441,19 +440,24 @@ namespace CosineKitty
 
         private static double UniversalTimeFromCalendar(int year, int month, int day, int hour, int minute, double second)
         {
-            // This formula is adapted from NOVAS C 3.1 function julian_date().
-            // It given a Gregorian calendar date/time, it calculates the fractional
-            // number of days since the J2000 epoch.
+            // This formula is adapted from NOVAS C 3.1 function julian_date(),
+            // which in turn comes from Henry F. Fliegel & Thomas C. Van Flendern:
+            // Communications of the ACM, Vol 11, No 10, October 1968, p. 657.
+            // See: https://dl.acm.org/doi/pdf/10.1145/364096.364097
+            //
+            // [Don Cross - 2023-02-25] I modified the formula so that it will
+            // work correctly with years as far back as -999999.
 
             long y = (long)year;
             long m = (long)month;
             long d = (long)day;
+            long f = (14 - m) / 12;
 
             long y2000 = (
-                (d - 2483620L)
-                + 1461L*(y + 4800L - (14L - m) / 12L)/4L
-                + 367L*(m - 2L + (14L - m) / 12L * 12L)/12L
-                - 3L*((y + 4900L - (14L - m) / 12L) / 100L)/4L
+                (d - 365972956)
+                + (1461*(y + 1000000 - f))/4
+                + (367*(m - 2 + 12*f))/12
+                - (3*((y + 1000100 - f) / 100))/4
             );
 
             double ut = (y2000 - 0.5) + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0);
@@ -461,22 +465,45 @@ namespace CosineKitty
         }
     }
 
-    internal struct CalendarDateTime
+    /// <summary>
+    /// Represents a Gregorian calendar date and time within plus or minus 1 million years from the year 0.
+    /// </summary>
+    /// <remarks>
+    /// The C# standard type `System.DateTime` only allows years from 0001 to 9999.
+    /// However, the #AstroTime class can represent years in the range -999999 to +999999.
+    /// In order to support formatting dates with extreme year values in an extrapolated
+    /// Gregorian calendar, the `CalendarDateTime` class breaks out the components of
+    /// a date into separate fields.
+    /// </remarks>
+    public struct CalendarDateTime
     {
+        /// <summary>The year value in the range -999999 to +999999.</summary>
         public int year;
+
+        /// <summary>The calendar month in the range 1..12.</summary>
         public int month;
+
+        /// <summary>The day of the month in the reange 1..31.</summary>
         public int day;
+
+        /// <summary>The hour in the range 0..23.</summary>
         public int hour;
+
+        /// <summary>The minute in the range 0..59.</summary>
         public int minute;
+
+        /// <summary>The real-valued second in the half-open range [0, 60).</summary>
         public double second;
 
+        /// <summary>Convert a J2000 day value to a Gregorian calendar date.</summary>
+        /// <param name="ut">The real-valued number of days since the J2000 epoch.</param>
         public CalendarDateTime(double ut)
         {
             // Adapted from the NOVAS C 3.1 function cal_date().
             // Convert fractional days since J2000 into Gregorian calendar date/time.
 
             double djd = ut + 2451545.5;
-            long jd = (long)djd;
+            long jd = (long)Math.Floor(djd);
             double x = 24.0 * (djd % 1.0);
             if (x < 0.0)
                 x += 24.0;
@@ -493,23 +520,43 @@ namespace CosineKitty
             // the calendar date calculations.
             // Any multiple of 400 years has the same number of days,
             // because it eliminates all the special cases for leap years.
-            const long c = 2500L;
+            const long c = 2500;
 
-            long k = jd + (68569L + c*146097L);
-            long n = 4L * k / 146097L;
-            k = k - (146097L * n + 3L) / 4L;
-            long m = 4000L * (k + 1L) / 1461001L;
-            k = k - 1461L * m / 4L + 31L;
+            long k = jd + (68569 + c*146097);
+            long n = (4 * k) / 146097;
+            k = k - (146097*n + 3) / 4;
+            long m = (4000 * (k+1)) / 1461001;
+            k = k - (1461 * m)/4 + 31;
 
-            month = (int) (80L * k / 2447L);
-            day = (int) (k - 2447L * (long)month / 80L);
-            k = (long) month / 11L;
+            month = (int) ((80 * k) / 2447);
+            day = (int) (k - (2447 * month)/80);
+            k = month / 11;
 
-            month = (int) ((long)month + 2L - 12L * k);
-            year = (int) (100L * (n - 49L) + m + k - 400L*c);
+            month = (int) (month + 2 - 12*k);
+            year = (int) (100 * (n - 49) + m + k - 400*c);
 
-            if (month < 1 || day < 1 || year < -999999 || year > +999999)
+            if (year < -999999 || year > +999999)
                 throw new ArgumentOutOfRangeException("The supplied time is too far from the year 2000 to be represented.");
+
+            if (month < 1 || month > 12 || day < 1 || day > 31)
+                throw new InternalError($"Invalid calendar date calculated: month={month}, day={day}.");
+        }
+
+        /// <summary>
+        /// Converts this `CalendarDateTime` to ISO 8601 format, expressed in UTC with millisecond resolution.
+        /// </summary>
+        /// <returns>Example: "2019-08-30T17:45:22.763Z".</returns>
+        public override string ToString()
+        {
+            int millis = Math.Max(0, Math.Min(59999, (int)Math.Round(second * 1000.0)));
+            string y;
+            if (year < 0)
+                y = "-" + (-year).ToString("000000");
+            else if (year <= 9999)
+                y = year.ToString("0000");
+            else
+                y = "+" + year.ToString("000000");
+            return $"{y}-{month:00}-{day:00}T{hour:00}:{minute:00}:{millis/1000:00}.{millis%1000:000}Z";
         }
     }
 

source/c/astronomy.c

@@ -1103,7 +1103,7 @@ astro_utc_t Astronomy_UtcFromTime(astro_time_t time)
     astro_utc_t utc;
     long jd, k, m, n;
     double djd, x;
-    const long c = 2500L;
+    const long c = 2500;
 
     djd = time.ut + 2451545.5;
     jd = (long)floor(djd);
@@ -1126,18 +1126,18 @@ astro_utc_t Astronomy_UtcFromTime(astro_time_t time)
         Any multiple of 400 years has the same number of days,
         because it eliminates all the special cases for leap years.
     */
-    k = jd + (68569L + c*146097L);
-    n = 4L * k / 146097L;
-    k = k - (146097L * n + 3L) / 4L;
-    m = 4000L * (k + 1L) / 1461001L;
-    k = k - 1461L * m / 4L + 31L;
+    k = jd + (68569 + c*146097);
+    n = (4 * k) / 146097;
+    k = k - (146097 * n + 3)/4;
+    m = (4000 * (k+1)) / 1461001;
+    k = k - (1461 * m)/4 + 31;
 
-    utc.month = (int) (80L * k / 2447L);
-    utc.day = (int) (k - 2447L * (long)utc.month / 80L);
-    k = (long) utc.month / 11L;
+    utc.month = (int) ((80 * k) / 2447);
+    utc.day = (int) (k - (2447*utc.month)/80);
+    k = utc.month / 11;
 
-    utc.month = (int) ((long)utc.month + 2L - 12L * k);
-    utc.year = (int) (100L * (n - 49L) + m + k - 400L*c);
+    utc.month = (int) (utc.month + 2 - 12*k);
+    utc.year = (int) (100 * (n - 49) + m + k - 400*c);
 
     return utc;
 }

source/csharp/README.md

@@ -2589,6 +2589,16 @@ an `AstroTime` value that can be passed to Astronomy Engine functions.
 | --- | --- | --- |
 | `double` | `tt` | The number of days after the J2000 epoch. |
 
+<a name="AstroTime.ToCalendarDateTime"></a>
+### AstroTime.ToCalendarDateTime() &#8658; [`CalendarDateTime`](#CalendarDateTime)
+
+**Converts this object to our custom type [`CalendarDateTime`](#CalendarDateTime).**
+
+The .NET type `DateTime` can only represent years in the range 0000..9999.
+However, the Astronomy Engine type [`CalendarDateTime`](#CalendarDateTime) can represent
+years in the range -999999..+999999. This is a time span of nearly 2 million years.
+This function converts this `AstroTime` object to an equivalent Gregorian calendar representation.
+
 <a name="AstroTime.ToString"></a>
 ### AstroTime.ToString() &#8658; `string`
 
@@ -2719,6 +2729,50 @@ It is expressed in the J2000 mean equator system (EQJ).
 
 ---
 
+<a name="CalendarDateTime"></a>
+## `struct CalendarDateTime`
+
+**Represents a Gregorian calendar date and time within plus or minus 1 million years from the year 0.**
+
+The C# standard type `System.DateTime` only allows years from 0001 to 9999.
+However, the [`AstroTime`](#AstroTime) class can represent years in the range -999999 to +999999.
+In order to support formatting dates with extreme year values in an extrapolated
+Gregorian calendar, the `CalendarDateTime` class breaks out the components of
+a date into separate fields.
+
+### constructors
+
+### new CalendarDateTime(ut)
+
+**Convert a J2000 day value to a Gregorian calendar date.**
+
+| Type | Parameter | Description |
+| --- | --- | --- |
+| `double` | `ut` | The real-valued number of days since the J2000 epoch. |
+
+
+### member variables
+
+| Type | Name | Description |
+| --- | --- | --- |
+| `int` | `year` | The year value in the range -999999 to +999999. |
+| `int` | `month` | The calendar month in the range 1..12. |
+| `int` | `day` | The day of the month in the reange 1..31. |
+| `int` | `hour` | The hour in the range 0..23. |
+| `int` | `minute` | The minute in the range 0..59. |
+| `double` | `second` | The real-valued second in the half-open range [0, 60). |
+
+### member functions
+
+<a name="CalendarDateTime.ToString"></a>
+### CalendarDateTime.ToString() &#8658; `string`
+
+**Converts this `CalendarDateTime` to ISO 8601 format, expressed in UTC with millisecond resolution.**
+
+**Returns:** Example: "2019-08-30T17:45:22.763Z".
+
+---
+
 <a name="ConstellationInfo"></a>
 ## `struct ConstellationInfo`
 

source/csharp/astronomy.cs

@@ -324,23 +324,22 @@ namespace CosineKitty
             return Origin.AddDays(ut).ToUniversalTime();
         }
 
+        /// <summary>
+        /// Converts this object to our custom type #CalendarDateTime.
+        /// </summary>
+        /// <remarks>
+        /// The .NET type `DateTime` can only represent years in the range 0000..9999.
+        /// However, the Astronomy Engine type #CalendarDateTime can represent
+        /// years in the range -999999..+999999. This is a time span of nearly 2 million years.
+        /// This function converts this `AstroTime` object to an equivalent Gregorian calendar representation.
+        /// </remarks>
+        public CalendarDateTime ToCalendarDateTime() => new CalendarDateTime(ut);
+
         /// <summary>
         /// Converts this `AstroTime` to ISO 8601 format, expressed in UTC with millisecond resolution.
         /// </summary>
         /// <returns>Example: "2019-08-30T17:45:22.763Z".</returns>
-        public override string ToString()
-        {
-            var d = new CalendarDateTime(ut);
-            int millis = Math.Max(0, Math.Min(59999, (int)Math.Round(d.second * 1000.0)));
-            string y;
-            if (d.year < 0)
-                y = "-" + (-d.year).ToString("000000");
-            else if (d.year <= 9999)
-                y = d.year.ToString("0000");
-            else
-                y = "+" + d.year.ToString("000000");
-            return $"{y}-{d.month:00}-{d.day:00}T{d.hour:00}:{d.minute:00}:{millis/1000:00}.{millis%1000:000}Z";
-        }
+        public override string ToString() => ToCalendarDateTime().ToString();
 
         private static Regex re = new Regex(
             @"^
@@ -441,19 +440,24 @@ namespace CosineKitty
 
         private static double UniversalTimeFromCalendar(int year, int month, int day, int hour, int minute, double second)
         {
-            // This formula is adapted from NOVAS C 3.1 function julian_date().
-            // It given a Gregorian calendar date/time, it calculates the fractional
-            // number of days since the J2000 epoch.
+            // This formula is adapted from NOVAS C 3.1 function julian_date(),
+            // which in turn comes from Henry F. Fliegel & Thomas C. Van Flendern:
+            // Communications of the ACM, Vol 11, No 10, October 1968, p. 657.
+            // See: https://dl.acm.org/doi/pdf/10.1145/364096.364097
+            //
+            // [Don Cross - 2023-02-25] I modified the formula so that it will
+            // work correctly with years as far back as -999999.
 
             long y = (long)year;
             long m = (long)month;
             long d = (long)day;
+            long f = (14 - m) / 12;
 
             long y2000 = (
-                (d - 2483620L)
-                + 1461L*(y + 4800L - (14L - m) / 12L)/4L
-                + 367L*(m - 2L + (14L - m) / 12L * 12L)/12L
-                - 3L*((y + 4900L - (14L - m) / 12L) / 100L)/4L
+                (d - 365972956)
+                + (1461*(y + 1000000 - f))/4
+                + (367*(m - 2 + 12*f))/12
+                - (3*((y + 1000100 - f) / 100))/4
             );
 
             double ut = (y2000 - 0.5) + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0);
@@ -461,22 +465,45 @@ namespace CosineKitty
         }
     }
 
-    internal struct CalendarDateTime
+    /// <summary>
+    /// Represents a Gregorian calendar date and time within plus or minus 1 million years from the year 0.
+    /// </summary>
+    /// <remarks>
+    /// The C# standard type `System.DateTime` only allows years from 0001 to 9999.
+    /// However, the #AstroTime class can represent years in the range -999999 to +999999.
+    /// In order to support formatting dates with extreme year values in an extrapolated
+    /// Gregorian calendar, the `CalendarDateTime` class breaks out the components of
+    /// a date into separate fields.
+    /// </remarks>
+    public struct CalendarDateTime
     {
+        /// <summary>The year value in the range -999999 to +999999.</summary>
         public int year;
+
+        /// <summary>The calendar month in the range 1..12.</summary>
         public int month;
+
+        /// <summary>The day of the month in the reange 1..31.</summary>
         public int day;
+
+        /// <summary>The hour in the range 0..23.</summary>
         public int hour;
+
+        /// <summary>The minute in the range 0..59.</summary>
         public int minute;
+
+        /// <summary>The real-valued second in the half-open range [0, 60).</summary>
         public double second;
 
+        /// <summary>Convert a J2000 day value to a Gregorian calendar date.</summary>
+        /// <param name="ut">The real-valued number of days since the J2000 epoch.</param>
         public CalendarDateTime(double ut)
         {
             // Adapted from the NOVAS C 3.1 function cal_date().
             // Convert fractional days since J2000 into Gregorian calendar date/time.
 
             double djd = ut + 2451545.5;
-            long jd = (long)djd;
+            long jd = (long)Math.Floor(djd);
             double x = 24.0 * (djd % 1.0);
             if (x < 0.0)
                 x += 24.0;
@@ -493,23 +520,43 @@ namespace CosineKitty
             // the calendar date calculations.
             // Any multiple of 400 years has the same number of days,
             // because it eliminates all the special cases for leap years.
-            const long c = 2500L;
+            const long c = 2500;
 
-            long k = jd + (68569L + c*146097L);
-            long n = 4L * k / 146097L;
-            k = k - (146097L * n + 3L) / 4L;
-            long m = 4000L * (k + 1L) / 1461001L;
-            k = k - 1461L * m / 4L + 31L;
+            long k = jd + (68569 + c*146097);
+            long n = (4 * k) / 146097;
+            k = k - (146097*n + 3) / 4;
+            long m = (4000 * (k+1)) / 1461001;
+            k = k - (1461 * m)/4 + 31;
 
-            month = (int) (80L * k / 2447L);
-            day = (int) (k - 2447L * (long)month / 80L);
-            k = (long) month / 11L;
+            month = (int) ((80 * k) / 2447);
+            day = (int) (k - (2447 * month)/80);
+            k = month / 11;
 
-            month = (int) ((long)month + 2L - 12L * k);
-            year = (int) (100L * (n - 49L) + m + k - 400L*c);
+            month = (int) (month + 2 - 12*k);
+            year = (int) (100 * (n - 49) + m + k - 400*c);
 
-            if (month < 1 || day < 1 || year < -999999 || year > +999999)
+            if (year < -999999 || year > +999999)
                 throw new ArgumentOutOfRangeException("The supplied time is too far from the year 2000 to be represented.");
+
+            if (month < 1 || month > 12 || day < 1 || day > 31)
+                throw new InternalError($"Invalid calendar date calculated: month={month}, day={day}.");
+        }
+
+        /// <summary>
+        /// Converts this `CalendarDateTime` to ISO 8601 format, expressed in UTC with millisecond resolution.
+        /// </summary>
+        /// <returns>Example: "2019-08-30T17:45:22.763Z".</returns>
+        public override string ToString()
+        {
+            int millis = Math.Max(0, Math.Min(59999, (int)Math.Round(second * 1000.0)));
+            string y;
+            if (year < 0)
+                y = "-" + (-year).ToString("000000");
+            else if (year <= 9999)
+                y = year.ToString("0000");
+            else
+                y = "+" + year.ToString("000000");
+            return $"{y}-{month:00}-{day:00}T{hour:00}:{minute:00}:{millis/1000:00}.{millis%1000:000}Z";
         }
     }