From 409e490728b44d4d0772ef46493fc34736722e57 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Thu, 6 Oct 2022 10:56:17 -0400 Subject: [PATCH] Python: No longer limited to years 0000..9999. I ported the NOVAS C 3.1 functions julian_date and cal_date to Python, and removed the dependence on the standard datetime class for calculating UT. Now we can create Time objects for a much wider range of year values. Simplified the julian_date formula in C and C#. In the Python version, I had to account for a difference in the way integer division works for negative numbers. In Python, integer division always rounds down, not toward zero like it does in C/C#. So I reworked the formulas to avoid dividing a negative integer (month-14), dividing the positive quantity (14-month) instead and toggling addition of the term with subtraction of the term. I use the reworked (14-month) version in C and C# for consistency. Also, the formatting of the formula was wacky and didn't make sense, so now it easier to read and understand. The Python regex for parsing dates has been expanded to allow years before 0 and after 9999. Allow converting Python Time to string for years before 0 and after 9999. --- demo/python/astronomy.py | 54 +++++++++++++++++++++++----- demo/python/correct/lunar_angles.txt | 6 ++-- generate/template/astronomy.c | 20 +++++------ generate/template/astronomy.cs | 18 ++++++---- generate/template/astronomy.py | 54 +++++++++++++++++++++++----- generate/test.py | 34 +++++++++++++++--- source/c/astronomy.c | 20 +++++------ source/csharp/astronomy.cs | 18 ++++++---- source/python/README.md | 2 +- source/python/astronomy/astronomy.py | 54 +++++++++++++++++++++++----- 10 files changed, 210 insertions(+), 70 deletions(-) diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py index 9a74b0cb..7c4650b3 100644 --- a/demo/python/astronomy.py +++ b/demo/python/astronomy.py @@ -577,7 +577,7 @@ def _UniversalTime(tt): return ut dt += err -_TimeRegex = re.compile(r'^([0-9]{1,4})-([0-9]{2})-([0-9]{2})(T([0-9]{2}):([0-9]{2})(:([0-9]{2}(\.[0-9]+)?))?Z)?$') +_TimeRegex = re.compile(r'^([\+\-]?[0-9]{1,5})-([0-9]{2})-([0-9]{2})(T([0-9]{2}):([0-9]{2})(:([0-9]{2}(\.[0-9]+)?))?Z)?$') class Time: """Represents a date and time used for performing astronomy calculations. @@ -705,7 +705,7 @@ class Time: Parameters ---------- year : int - The UTC 4-digit year value, e.g. 2019. + The UTC year value, e.g. 2019. month : int The UTC month in the range 1..12. day : int @@ -721,10 +721,17 @@ class Time: ------- Time """ - micro = round(math.fmod(second, 1.0) * 1000000) - second = math.floor(second - micro/1000000) - d = datetime.datetime(year, month, day, hour, minute, second, micro) - ut = (d - _EPOCH).total_seconds() / 86400 + # This formula is adapted from NOVAS C 3.1 function julian_date(). + y = int(year) + m = int(month) + d = int(day) + y2000 = float( + (d - 2483620) + + 1461 * (y + 4800 - (14 - m) // 12) // 4 + + 367 * (m - 2 + (14 - m) // 12 * 12) // 12 + - 3 * ((y + 4900 - (14 - m) // 12) // 100) // 4 + ) + ut = (y2000 - 0.5) + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0) return Time(ut) @staticmethod @@ -774,9 +781,38 @@ class Time: return 'Time(\'' + str(self) + '\')' def __str__(self): - millis = round(self.ut * 86400000.0) - n = _EPOCH + datetime.timedelta(milliseconds=millis) - return '{:04d}-{:02d}-{:02d}T{:02d}:{:02d}:{:02d}.{:03d}Z'.format(n.year, n.month, n.day, n.hour, n.minute, n.second, math.floor(n.microsecond / 1000)) + # Adapted from the NOVAS C 3.1 function cal_date(). + djd = self.ut + 2451545.5 + jd = int(djd) + x = 24.0 * math.fmod(djd, 1.0) + hour = int(x) + x = 60.0 * math.fmod(x, 1.0) + minute = int(x) + second = round(60.0 * math.fmod(x, 1.0), 3) + if second >= 60.0: + second -= 60.0 + minute += 1 + if minute >= 60: + minute -= 60 + hour += 1 + if hour >= 24: + hour -= 24 + jd += 1 + k = jd + 68569 + n = 4 * k // 146097 + k = k - (146097 * n + 3) // 4 + m = 4000 * (k + 1) // 1461001 + k = k - 1461 * m // 4 + 31 + month = (80 * k) // 2447 + day = k - 2447 * month // 80 + k = month // 11 + month = month + 2 - 12 * k + year = 100 * (n - 49) + m + k + millis = max(0, min(59999, round(1000.0 * second))) + text = '{:04d}-{:02d}-{:02d}T{:02d}:{:02d}:{:02d}.{:03d}Z'.format(year, month, day, hour, minute, millis // 1000, millis % 1000) + if year > 9999: + text = '+' + text + return text def Utc(self): """Returns the UTC date and time as a `datetime` object. diff --git a/demo/python/correct/lunar_angles.txt b/demo/python/correct/lunar_angles.txt index 865eabf5..80b97880 100644 --- a/demo/python/correct/lunar_angles.txt +++ b/demo/python/correct/lunar_angles.txt @@ -1,4 +1,4 @@ -2021-05-15T01:45:15.503Z Jupiter 120 +2021-05-15T01:45:15.502Z Jupiter 120 2021-05-15T17:53:54.557Z Venus 30 2021-05-16T04:23:24.323Z Saturn 150 2021-05-16T05:06:25.856Z Mars 0 @@ -15,7 +15,7 @@ 2021-05-20T23:01:26.198Z Venus 90 2021-05-21T02:57:43.803Z Mars 60 2021-05-21T11:47:54.388Z Mercury 90 -2021-05-22T02:58:57.124Z Jupiter 210 +2021-05-22T02:58:57.125Z Jupiter 210 2021-05-22T03:46:48.175Z Sun 120 2021-05-23T00:10:42.249Z Saturn 240 2021-05-23T06:36:37.117Z Venus 120 @@ -53,7 +53,7 @@ 2021-06-03T21:09:24.870Z Venus 270 2021-06-03T21:22:13.482Z Jupiter 30 2021-06-04T20:23:38.951Z Saturn 60 -2021-06-04T22:38:31.596Z Sun 300 +2021-06-04T22:38:31.597Z Sun 300 2021-06-05T15:57:17.165Z Mercury 300 2021-06-05T22:47:34.963Z Mars 270 2021-06-06T09:32:04.351Z Jupiter 60 diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index e79537f0..35e0fcf7 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -910,19 +910,19 @@ astro_time_t Astronomy_CurrentTime(void) astro_time_t Astronomy_MakeTime(int year, int month, int day, int hour, int minute, double second) { astro_time_t time; - long int jd12h; - long int y2000; + long y = (long)year; + long m = (long)month; + long d = (long)day; /* This formula is adapted from NOVAS C 3.1 function julian_date() */ - jd12h = (long) day - 32075L + 1461L * ((long) year + 4800L - + ((long) month - 14L) / 12L) / 4L - + 367L * ((long) month - 2L - ((long) month - 14L) / 12L * 12L) - / 12L - 3L * (((long) year + 4900L + ((long) month - 14L) / 12L) - / 100L) / 4L; + 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 + ); - y2000 = jd12h - 2451545L; - - time.ut = (double)y2000 - 0.5 + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0); + time.ut = (y2000 - 0.5) + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0); time.tt = TerrestrialTime(time.ut); time.psi = time.eps = time.st = NAN; diff --git a/generate/template/astronomy.cs b/generate/template/astronomy.cs index ef971fae..80881efe 100644 --- a/generate/template/astronomy.cs +++ b/generate/template/astronomy.cs @@ -321,15 +321,19 @@ namespace CosineKitty // It given a Gregorian calendar date/time, it calculates the fractional // number of days since the J2000 epoch. - long jd12h = (long) day - 32075L + 1461L * ((long) year + 4800L - + ((long) month - 14L) / 12L) / 4L - + 367L * ((long) month - 2L - ((long) month - 14L) / 12L * 12L) - / 12L - 3L * (((long) year + 4900L + ((long) month - 14L) / 12L) - / 100L) / 4L; + long y = (long)year; + long m = (long)month; + long d = (long)day; - long y2000 = jd12h - 2451545L; + 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 + ); - return (double)y2000 - 0.5 + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0); + double ut = (y2000 - 0.5) + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0); + return ut; } } diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index b304fba9..99b0a43a 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -577,7 +577,7 @@ def _UniversalTime(tt): return ut dt += err -_TimeRegex = re.compile(r'^([0-9]{1,4})-([0-9]{2})-([0-9]{2})(T([0-9]{2}):([0-9]{2})(:([0-9]{2}(\.[0-9]+)?))?Z)?$') +_TimeRegex = re.compile(r'^([\+\-]?[0-9]{1,5})-([0-9]{2})-([0-9]{2})(T([0-9]{2}):([0-9]{2})(:([0-9]{2}(\.[0-9]+)?))?Z)?$') class Time: """Represents a date and time used for performing astronomy calculations. @@ -705,7 +705,7 @@ class Time: Parameters ---------- year : int - The UTC 4-digit year value, e.g. 2019. + The UTC year value, e.g. 2019. month : int The UTC month in the range 1..12. day : int @@ -721,10 +721,17 @@ class Time: ------- Time """ - micro = round(math.fmod(second, 1.0) * 1000000) - second = math.floor(second - micro/1000000) - d = datetime.datetime(year, month, day, hour, minute, second, micro) - ut = (d - _EPOCH).total_seconds() / 86400 + # This formula is adapted from NOVAS C 3.1 function julian_date(). + y = int(year) + m = int(month) + d = int(day) + y2000 = float( + (d - 2483620) + + 1461 * (y + 4800 - (14 - m) // 12) // 4 + + 367 * (m - 2 + (14 - m) // 12 * 12) // 12 + - 3 * ((y + 4900 - (14 - m) // 12) // 100) // 4 + ) + ut = (y2000 - 0.5) + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0) return Time(ut) @staticmethod @@ -774,9 +781,38 @@ class Time: return 'Time(\'' + str(self) + '\')' def __str__(self): - millis = round(self.ut * 86400000.0) - n = _EPOCH + datetime.timedelta(milliseconds=millis) - return '{:04d}-{:02d}-{:02d}T{:02d}:{:02d}:{:02d}.{:03d}Z'.format(n.year, n.month, n.day, n.hour, n.minute, n.second, math.floor(n.microsecond / 1000)) + # Adapted from the NOVAS C 3.1 function cal_date(). + djd = self.ut + 2451545.5 + jd = int(djd) + x = 24.0 * math.fmod(djd, 1.0) + hour = int(x) + x = 60.0 * math.fmod(x, 1.0) + minute = int(x) + second = round(60.0 * math.fmod(x, 1.0), 3) + if second >= 60.0: + second -= 60.0 + minute += 1 + if minute >= 60: + minute -= 60 + hour += 1 + if hour >= 24: + hour -= 24 + jd += 1 + k = jd + 68569 + n = 4 * k // 146097 + k = k - (146097 * n + 3) // 4 + m = 4000 * (k + 1) // 1461001 + k = k - 1461 * m // 4 + 31 + month = (80 * k) // 2447 + day = k - 2447 * month // 80 + k = month // 11 + month = month + 2 - 12 * k + year = 100 * (n - 49) + m + k + millis = max(0, min(59999, round(1000.0 * second))) + text = '{:04d}-{:02d}-{:02d}T{:02d}:{:02d}:{:02d}.{:03d}Z'.format(year, month, day, hour, minute, millis // 1000, millis % 1000) + if year > 9999: + text = '+' + text + return text def Utc(self): """Returns the UTC date and time as a `datetime` object. diff --git a/generate/test.py b/generate/test.py index 25babef2..89f6cebd 100755 --- a/generate/test.py +++ b/generate/test.py @@ -63,25 +63,25 @@ def AstroTime(): diff = time.ut - expected_ut if vabs(diff) > 1.0e-12: print('PY AstroTime: excessive UT error {}'.format(diff)) - sys.exit(1) + return 1 diff = time.tt - expected_tt if vabs(diff) > 1.0e-12: print('PY AstroTime: excessive TT error {}'.format(diff)) - sys.exit(1) + return 1 s = str(time.Utc()) if s != '2018-12-02 18:30:12.543000': print('PY AstroTime: Utc() returned incorrect string "{}"'.format(s)) - sys.exit(1) + return 1 time = astronomy.Time.Make(2018, 12, 31, 23, 59, 59.9994) s = str(time) if s != '2018-12-31T23:59:59.999Z': print('PY AstroTime: expected 2018-12-31T23:59:59.999Z but found {}'.format(s)) - sys.exit(1) + return 1 time = astronomy.Time.Make(2018, 12, 31, 23, 59, 59.9995) s = str(time) if s != '2019-01-01T00:00:00.000Z': print('PY AstroTime: expected 2019-01-01T00:00:00.000Z but found {}'.format(s)) - sys.exit(1) + return 1 print('PY Current time =', astronomy.Time.Now()) AssertGoodTime('2015-12-31', '2015-12-31T00:00:00.000Z') AssertGoodTime('2015-12-31T23:45Z', '2015-12-31T23:45:00.000Z') @@ -94,6 +94,9 @@ def AstroTime(): AssertBadTime('1971-12-31T23:60:00Z') AssertBadTime('1971-12-31T23:00:60Z') AssertBadTime('1971-03-17T03:30:55.976') + # Extreme year values... + AssertGoodTime('-2300-12-19T16:22:26.325Z', '-2300-12-19T16:22:26.325Z') + AssertGoodTime('+12345-12-11T13:30:10.041Z', '+12345-12-11T13:30:10.041Z') return 0 #----------------------------------------------------------------------------------------------------------- @@ -2851,11 +2854,32 @@ def ArcminVelError(correct, calc): #----------------------------------------------------------------------------------------------------------- +def CheckDecemberSolstice(year, expected): + si = astronomy.Seasons(year) + actual = str(si.dec_solstice) + if actual != expected: + print('PY DatesIssue250: FAIL: year {}, expected [{}], actual [{}]'.format(year, expected, actual)) + return 1 + return 0 + +def DatesIssue250(): + # Make sure we can handle dates outside the range supported by System.DateTime. + # https://github.com/cosinekitty/astronomy/issues/250 + return ( + CheckDecemberSolstice( 2022, "2022-12-21T21:47:58.189Z") or + CheckDecemberSolstice(-2300, "-2300-12-19T16:22:26.325Z") or + CheckDecemberSolstice(12345, "+12345-12-11T13:30:10.041Z") or + Pass('DatesIssue250') + ) + +#----------------------------------------------------------------------------------------------------------- + UnitTests = { 'aberration': Aberration, 'axis': Axis, 'barystate': BaryState, 'constellation': Constellation, + 'dates250': DatesIssue250, 'elongation': Elongation, 'geoid': Geoid, 'global_solar_eclipse': GlobalSolarEclipse, diff --git a/source/c/astronomy.c b/source/c/astronomy.c index a547f5af..4fec4cbc 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -916,19 +916,19 @@ astro_time_t Astronomy_CurrentTime(void) astro_time_t Astronomy_MakeTime(int year, int month, int day, int hour, int minute, double second) { astro_time_t time; - long int jd12h; - long int y2000; + long y = (long)year; + long m = (long)month; + long d = (long)day; /* This formula is adapted from NOVAS C 3.1 function julian_date() */ - jd12h = (long) day - 32075L + 1461L * ((long) year + 4800L - + ((long) month - 14L) / 12L) / 4L - + 367L * ((long) month - 2L - ((long) month - 14L) / 12L * 12L) - / 12L - 3L * (((long) year + 4900L + ((long) month - 14L) / 12L) - / 100L) / 4L; + 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 + ); - y2000 = jd12h - 2451545L; - - time.ut = (double)y2000 - 0.5 + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0); + time.ut = (y2000 - 0.5) + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0); time.tt = TerrestrialTime(time.ut); time.psi = time.eps = time.st = NAN; diff --git a/source/csharp/astronomy.cs b/source/csharp/astronomy.cs index d5222c21..32f42faa 100644 --- a/source/csharp/astronomy.cs +++ b/source/csharp/astronomy.cs @@ -321,15 +321,19 @@ namespace CosineKitty // It given a Gregorian calendar date/time, it calculates the fractional // number of days since the J2000 epoch. - long jd12h = (long) day - 32075L + 1461L * ((long) year + 4800L - + ((long) month - 14L) / 12L) / 4L - + 367L * ((long) month - 2L - ((long) month - 14L) / 12L * 12L) - / 12L - 3L * (((long) year + 4900L + ((long) month - 14L) / 12L) - / 100L) / 4L; + long y = (long)year; + long m = (long)month; + long d = (long)day; - long y2000 = jd12h - 2451545L; + 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 + ); - return (double)y2000 - 0.5 + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0); + double ut = (y2000 - 0.5) + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0); + return ut; } } diff --git a/source/python/README.md b/source/python/README.md index ea6477f2..eecb0ac8 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -917,7 +917,7 @@ The value of the calling object is not modified. This function creates a brand n | Type | Parameter | Description | | --- | --- | --- | -| `int` | `year` | The UTC 4-digit year value, e.g. 2019. | +| `int` | `year` | The UTC year value, e.g. 2019. | | `int` | `month` | The UTC month in the range 1..12. | | `int` | `day` | The UTC day of the month, in the range 1..31. | | `int` | `hour` | The UTC hour, in the range 0..23. | diff --git a/source/python/astronomy/astronomy.py b/source/python/astronomy/astronomy.py index 9a74b0cb..7c4650b3 100644 --- a/source/python/astronomy/astronomy.py +++ b/source/python/astronomy/astronomy.py @@ -577,7 +577,7 @@ def _UniversalTime(tt): return ut dt += err -_TimeRegex = re.compile(r'^([0-9]{1,4})-([0-9]{2})-([0-9]{2})(T([0-9]{2}):([0-9]{2})(:([0-9]{2}(\.[0-9]+)?))?Z)?$') +_TimeRegex = re.compile(r'^([\+\-]?[0-9]{1,5})-([0-9]{2})-([0-9]{2})(T([0-9]{2}):([0-9]{2})(:([0-9]{2}(\.[0-9]+)?))?Z)?$') class Time: """Represents a date and time used for performing astronomy calculations. @@ -705,7 +705,7 @@ class Time: Parameters ---------- year : int - The UTC 4-digit year value, e.g. 2019. + The UTC year value, e.g. 2019. month : int The UTC month in the range 1..12. day : int @@ -721,10 +721,17 @@ class Time: ------- Time """ - micro = round(math.fmod(second, 1.0) * 1000000) - second = math.floor(second - micro/1000000) - d = datetime.datetime(year, month, day, hour, minute, second, micro) - ut = (d - _EPOCH).total_seconds() / 86400 + # This formula is adapted from NOVAS C 3.1 function julian_date(). + y = int(year) + m = int(month) + d = int(day) + y2000 = float( + (d - 2483620) + + 1461 * (y + 4800 - (14 - m) // 12) // 4 + + 367 * (m - 2 + (14 - m) // 12 * 12) // 12 + - 3 * ((y + 4900 - (14 - m) // 12) // 100) // 4 + ) + ut = (y2000 - 0.5) + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0) return Time(ut) @staticmethod @@ -774,9 +781,38 @@ class Time: return 'Time(\'' + str(self) + '\')' def __str__(self): - millis = round(self.ut * 86400000.0) - n = _EPOCH + datetime.timedelta(milliseconds=millis) - return '{:04d}-{:02d}-{:02d}T{:02d}:{:02d}:{:02d}.{:03d}Z'.format(n.year, n.month, n.day, n.hour, n.minute, n.second, math.floor(n.microsecond / 1000)) + # Adapted from the NOVAS C 3.1 function cal_date(). + djd = self.ut + 2451545.5 + jd = int(djd) + x = 24.0 * math.fmod(djd, 1.0) + hour = int(x) + x = 60.0 * math.fmod(x, 1.0) + minute = int(x) + second = round(60.0 * math.fmod(x, 1.0), 3) + if second >= 60.0: + second -= 60.0 + minute += 1 + if minute >= 60: + minute -= 60 + hour += 1 + if hour >= 24: + hour -= 24 + jd += 1 + k = jd + 68569 + n = 4 * k // 146097 + k = k - (146097 * n + 3) // 4 + m = 4000 * (k + 1) // 1461001 + k = k - 1461 * m // 4 + 31 + month = (80 * k) // 2447 + day = k - 2447 * month // 80 + k = month // 11 + month = month + 2 - 12 * k + year = 100 * (n - 49) + m + k + millis = max(0, min(59999, round(1000.0 * second))) + text = '{:04d}-{:02d}-{:02d}T{:02d}:{:02d}:{:02d}.{:03d}Z'.format(year, month, day, hour, minute, millis // 1000, millis % 1000) + if year > 9999: + text = '+' + text + return text def Utc(self): """Returns the UTC date and time as a `datetime` object.