From f53797453069a98bd69ae1ae6de027adb7c2bea8 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Fri, 24 Feb 2023 20:31:55 -0500 Subject: [PATCH] C: Fixed bugs in Astronomy_UtcFromTime for extreme year values. Addressed limitations of the logic I copied from NOVAS cal_date(). Its formulas did not work for years much before -12000 due to integer division going negative. I figured out how to make the formulas work for plus or minus 1 million years from the present era. --- generate/ctest.c | 8 ++++++++ generate/template/astronomy.c | 15 +++++++++++++-- source/c/astronomy.c | 15 +++++++++++++-- 3 files changed, 34 insertions(+), 4 deletions(-) diff --git a/generate/ctest.c b/generate/ctest.c index ae0df7ce..884e55a9 100644 --- a/generate/ctest.c +++ b/generate/ctest.c @@ -494,6 +494,14 @@ static int Test_AstroTime(void) CHECK(CheckTimeFormat(time, TIME_FORMAT_MINUTE, ASTRO_SUCCESS, "2021-01-01T00:00Z")); CHECK(CheckTimeFormat(time, TIME_FORMAT_DAY, ASTRO_SUCCESS, "2020-12-31")); + // The verification texts are intentionally off by a millisecond or two in the following + // two test cases, because the extreme year values expose floating point precision limitions. + time = Astronomy_MakeTime(-999999, 1, 14, 22, 55, 12.0); + CHECK(CheckTimeFormat(time, TIME_FORMAT_MILLI, ASTRO_SUCCESS, "-999999-01-14T22:55:11.998Z")); + + time = Astronomy_MakeTime(+999999, 11, 30, 8, 15, 45.0); + CHECK(CheckTimeFormat(time, TIME_FORMAT_MILLI, ASTRO_SUCCESS, "+999999-11-30T08:15:44.999Z")); + FPASS(); fail: return error; diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index 062eeecd..9fcf4e1b 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -1088,17 +1088,28 @@ 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; djd = time.ut + 2451545.5; jd = (long)djd; x = 24.0 * fmod(djd, 1.0); + if (x < 0.0) + x += 24.0; utc.hour = (int)x; x = 60.0 * fmod(x, 1.0); utc.minute = (int)x; utc.second = 60.0 * fmod(x, 1.0); - k = jd + 68569L; + // This is my own adjustment to the NOVAS cal_date logic + // so that it can handle dates much farther back in the past. + // I add c*400 years worth of days at the front, + // then subtract c*400 years at the back, + // which avoids negative values in the formulas that mess up + // 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. + k = jd + (68569L + c*146097L); n = 4L * k / 146097L; k = k - (146097L * n + 3L) / 4L; m = 4000L * (k + 1L) / 1461001L; @@ -1109,7 +1120,7 @@ astro_utc_t Astronomy_UtcFromTime(astro_time_t time) k = (long) utc.month / 11L; utc.month = (int) ((long)utc.month + 2L - 12L * k); - utc.year = (int) (100L * (n - 49L) + m + k); + utc.year = (int) (100L * (n - 49L) + m + k - 400L*c); return utc; } diff --git a/source/c/astronomy.c b/source/c/astronomy.c index 43e49a29..8df1cde5 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -1094,17 +1094,28 @@ 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; djd = time.ut + 2451545.5; jd = (long)djd; x = 24.0 * fmod(djd, 1.0); + if (x < 0.0) + x += 24.0; utc.hour = (int)x; x = 60.0 * fmod(x, 1.0); utc.minute = (int)x; utc.second = 60.0 * fmod(x, 1.0); - k = jd + 68569L; + // This is my own adjustment to the NOVAS cal_date logic + // so that it can handle dates much farther back in the past. + // I add c*400 years worth of days at the front, + // then subtract c*400 years at the back, + // which avoids negative values in the formulas that mess up + // 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. + k = jd + (68569L + c*146097L); n = 4L * k / 146097L; k = k - (146097L * n + 3L) / 4L; m = 4000L * (k + 1L) / 1461001L; @@ -1115,7 +1126,7 @@ astro_utc_t Astronomy_UtcFromTime(astro_time_t time) k = (long) utc.month / 11L; utc.month = (int) ((long)utc.month + 2L - 12L * k); - utc.year = (int) (100L * (n - 49L) + m + k); + utc.year = (int) (100L * (n - 49L) + m + k - 400L*c); return utc; }