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.
This commit is contained in:
Don Cross
2023-02-24 20:31:55 -05:00
parent 88a94b860f
commit f537974530
3 changed files with 34 additions and 4 deletions

View File

@@ -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;

View File

@@ -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;
}

View File

@@ -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;
}