From 503538da12cdc345fc6e5724733e141da76ff010 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Sat, 25 Feb 2023 22:37:58 -0500 Subject: [PATCH] PY: Fixed calendar/time conversion functions for extreme year values. Applying the same recent fixes to C and C# to the Python code. I'm also changing my philosophy of representing times. From now on, they will be truncated to the floor millisecond, not rounded to the nearest millisecond. This means we don't reach another calendar date until we have had 60 full seconds after the last minute. Otherwise there is too much nasty logic for rounding up calendar dates. I will follow suit across all languages. --- demo/python/astronomy.py | 81 ++++++++++++++------------- demo/python/correct/constellation.txt | 20 +++---- demo/python/correct/culminate.txt | 8 +-- demo/python/correct/lunar_angles.txt | 78 +++++++++++++------------- demo/python/correct/lunar_eclipse.txt | 40 ++++++------- demo/python/correct/moonphase.txt | 12 ++-- demo/python/correct/riseset.txt | 4 +- demo/python/correct/seasons.txt | 4 +- generate/template/astronomy.py | 81 ++++++++++++++------------- generate/test.py | 40 +++++++------ source/python/README.md | 10 ++++ source/python/astronomy/astronomy.py | 81 ++++++++++++++------------- 12 files changed, 242 insertions(+), 217 deletions(-) diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py index e4d030a6..9ca546f7 100644 --- a/demo/python/astronomy.py +++ b/demo/python/astronomy.py @@ -48,6 +48,11 @@ def _cbrt(x: float) -> float: return y raise InternalError() +def _cdiv(numer: int, denom: int) -> int: + '''Divide negative numbers the same way C does: rounding toward zero, not down.''' + sign = -1 if (numer * denom) < 0 else +1 + return sign * (abs(numer) // abs(denom)) + KM_PER_AU = 1.4959787069098932e+8 # The number of kilometers per astronomical unit. C_AUDAY = 173.1446326846693 # The speed of light expressed in astronomical units per day. AU_PER_LY = 63241.07708807546 # The number of astronomical units in one light-year. @@ -275,11 +280,6 @@ def _UniversalTime(tt: float) -> float: _TimeRegex = re.compile(r'^([\+\-]?[0-9]+)-([0-9]{2})-([0-9]{2})(T([0-9]{2}):([0-9]{2})(:([0-9]{2}(\.[0-9]+)?))?Z)?$') -def _cdiv(numer: int, denom: int) -> int: - '''Divide negative numbers the same way C does: rounding toward zero, not down.''' - sign = +1 if (numer * denom) >= 0 else -1 - return sign * (abs(numer) // abs(denom)) - class Time: """Represents a date and time used for performing astronomy calculations. @@ -428,10 +428,10 @@ class Time: d = int(day) f = (14 - m) // 12 y2000 = ( - (d - 2483620) - + _cdiv(1461 * (y + 4800 - f), 4) + (d - 365972956) + + _cdiv(1461 * (y + 1000000 - f), 4) + _cdiv(367 * (m - 2 + f*12), 12) - - _cdiv(3 * _cdiv(y + 4900 - f, 100), 4) + - _cdiv(3 * _cdiv(y + 1000100 - f, 100), 4) ) ut = (y2000 - 0.5) + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0) return Time(ut) @@ -483,37 +483,8 @@ class Time: return 'Time(\'' + str(self) + '\')' def __str__(self) -> str: - # 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) - if x < 0.0: - x += 24.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 - c = 2500 - 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 - month = (80 * k) // 2447 - day = k - 2447 * month // 80 - k = month // 11 - month = month + 2 - 12 * k - year = 100 * (n - 49) + m + k - c*400 - millis = max(0, min(59999, round(1000.0 * second))) + (year, month, day, hour, minute, second) = self.Calendar() + millis = max(0, min(59999, int(math.floor(1000.0 * second)))) if year < 0: text = '-{:06d}'.format(-year) elif year <= 9999: @@ -523,6 +494,38 @@ class Time: text += '-{:02d}-{:02d}T{:02d}:{:02d}:{:02d}.{:03d}Z'.format(month, day, hour, minute, millis // 1000, millis % 1000) return text + def Calendar(self) -> Tuple[int, int, int, int, int, float]: + """Returns a tuple of the form (year, month, day, hour, minute, second). + + This is a convenience method for converting a `Time` value into + its Gregorian calendar date/time representation. + Unlike the built-in `datetime` class, this method can represent + dates over a nearly 2 million year range: the years -999999 to +999999. + """ + # Adapted from the NOVAS C 3.1 function cal_date(). + # [Don Cross - 2023-02-25] Fixed to handle a much wider range of years. + djd = self.ut + 2451545.5 + jd = int(math.floor(djd)) + x = 24.0 * math.fmod(djd, 1.0) + if x < 0.0: + x += 24.0 + hour = int(x) + x = 60.0 * math.fmod(x, 1.0) + minute = int(x) + second = 60.0 * math.fmod(x, 1.0) + c = 2500 + k = jd + (68569 + c*146097) + n = _cdiv(4*k, 146097) + k = k - _cdiv(146097*n + 3, 4) + m = _cdiv(4000*(k+1), 1461001) + k = k - _cdiv(1461 * m, 4) + 31 + month = _cdiv(80 * k, 2447) + day = k - _cdiv(2447 * month, 80) + k = _cdiv(month, 11) + month = month + 2 - 12 * k + year = 100 * (n - 49) + m + k - c*400 + return (year, month, day, hour, minute, second) + def Utc(self) -> datetime.datetime: """Returns the UTC date and time as a `datetime` object. diff --git a/demo/python/correct/constellation.txt b/demo/python/correct/constellation.txt index eb9765ce..fb7295a5 100644 --- a/demo/python/correct/constellation.txt +++ b/demo/python/correct/constellation.txt @@ -1,18 +1,18 @@ -2021-06-01T01:43:20.607Z : Moon leaves Capricornus and enters Aquarius. +2021-06-01T01:43:20.606Z : Moon leaves Capricornus and enters Aquarius. 2021-06-03T12:20:46.904Z : Moon leaves Aquarius and enters Pisces. 2021-06-04T04:34:27.528Z : Moon leaves Pisces and enters Cetus. -2021-06-05T03:50:55.644Z : Moon leaves Cetus and enters Pisces. -2021-06-06T11:54:31.671Z : Moon leaves Pisces and enters Cetus. -2021-06-06T17:34:49.729Z : Moon leaves Cetus and enters Aries. +2021-06-05T03:50:55.643Z : Moon leaves Cetus and enters Pisces. +2021-06-06T11:54:31.670Z : Moon leaves Pisces and enters Cetus. +2021-06-06T17:34:49.728Z : Moon leaves Cetus and enters Aries. 2021-06-08T05:07:03.380Z : Moon leaves Aries and enters Taurus. 2021-06-11T08:17:23.569Z : Moon leaves Taurus and enters Gemini. -2021-06-13T13:37:29.742Z : Moon leaves Gemini and enters Cancer. +2021-06-13T13:37:29.741Z : Moon leaves Gemini and enters Cancer. 2021-06-15T02:27:52.921Z : Moon leaves Cancer and enters Leo. 2021-06-17T18:40:24.949Z : Moon leaves Leo and enters Virgo. -2021-06-20T23:40:23.974Z : Moon leaves Virgo and enters Libra. -2021-06-22T17:24:01.187Z : Moon leaves Libra and enters Scorpius. -2021-06-23T01:36:00.307Z : Moon leaves Scorpius and enters Ophiuchus. -2021-06-24T07:36:08.596Z : Moon leaves Ophiuchus and enters Sagittarius. -2021-06-26T12:43:42.908Z : Moon leaves Sagittarius and enters Capricornus. +2021-06-20T23:40:23.973Z : Moon leaves Virgo and enters Libra. +2021-06-22T17:24:01.186Z : Moon leaves Libra and enters Scorpius. +2021-06-23T01:36:00.306Z : Moon leaves Scorpius and enters Ophiuchus. +2021-06-24T07:36:08.595Z : Moon leaves Ophiuchus and enters Sagittarius. +2021-06-26T12:43:42.907Z : Moon leaves Sagittarius and enters Capricornus. 2021-06-28T10:40:55.380Z : Moon leaves Capricornus and enters Aquarius. 2021-06-30T19:53:26.746Z : Moon leaves Aquarius and enters Pisces. diff --git a/demo/python/correct/culminate.txt b/demo/python/correct/culminate.txt index e08f992b..b8e68c72 100644 --- a/demo/python/correct/culminate.txt +++ b/demo/python/correct/culminate.txt @@ -1,11 +1,11 @@ search : 2015-02-28T00:00:00.000Z Sun : 2015-02-28T18:12:31.044Z altitude= 52.13 azimuth= 180.00 -Moon : 2015-02-28T01:58:15.630Z altitude= 77.88 azimuth= 180.00 +Moon : 2015-02-28T01:58:15.629Z altitude= 77.88 azimuth= 180.00 Mercury : 2015-02-28T16:31:18.939Z altitude= 42.57 azimuth= 180.00 Venus : 2015-02-28T20:03:56.711Z altitude= 63.24 azimuth= 180.00 -Mars : 2015-02-28T19:52:21.644Z altitude= 62.24 azimuth= 180.00 +Mars : 2015-02-28T19:52:21.643Z altitude= 62.24 azimuth= 180.00 Jupiter : 2015-02-28T04:40:04.003Z altitude= 77.29 azimuth= 180.00 -Saturn : 2015-02-28T11:40:53.878Z altitude= 40.96 azimuth= 180.00 +Saturn : 2015-02-28T11:40:53.877Z altitude= 40.96 azimuth= 180.00 Uranus : 2015-02-28T20:20:51.246Z altitude= 65.12 azimuth= 180.00 Neptune : 2015-02-28T18:04:23.941Z altitude= 50.53 azimuth= 180.00 -Pluto : 2015-02-28T14:31:39.864Z altitude= 39.51 azimuth= 180.00 +Pluto : 2015-02-28T14:31:39.863Z altitude= 39.51 azimuth= 180.00 diff --git a/demo/python/correct/lunar_angles.txt b/demo/python/correct/lunar_angles.txt index 667856c6..83779a62 100644 --- a/demo/python/correct/lunar_angles.txt +++ b/demo/python/correct/lunar_angles.txt @@ -1,77 +1,77 @@ -2021-05-15T01:45:15.706Z Jupiter 120 +2021-05-15T01:45:15.705Z Jupiter 120 2021-05-15T17:53:54.314Z Venus 30 2021-05-16T04:23:24.170Z Saturn 150 2021-05-16T05:06:25.682Z Mars 0 -2021-05-16T12:51:07.810Z Mercury 30 +2021-05-16T12:51:07.809Z Mercury 30 2021-05-17T06:05:45.825Z Sun 60 -2021-05-17T13:28:33.770Z Jupiter 150 +2021-05-17T13:28:33.769Z Jupiter 150 2021-05-18T10:43:22.784Z Venus 60 2021-05-18T14:28:30.696Z Saturn 180 -2021-05-18T18:01:39.264Z Mars 30 +2021-05-18T18:01:39.263Z Mars 30 2021-05-19T02:55:46.814Z Mercury 60 2021-05-19T19:13:15.714Z Sun 90 2021-05-19T22:06:37.540Z Jupiter 180 -2021-05-20T21:07:42.837Z Saturn 210 -2021-05-20T23:01:24.412Z Venus 90 -2021-05-21T02:57:42.297Z Mars 60 +2021-05-20T21:07:42.836Z Saturn 210 +2021-05-20T23:01:24.411Z Venus 90 +2021-05-21T02:57:42.296Z Mars 60 2021-05-21T11:47:52.705Z Mercury 90 2021-05-22T02:58:55.870Z Jupiter 210 2021-05-22T03:46:46.605Z Sun 120 -2021-05-23T00:10:41.162Z Saturn 240 -2021-05-23T06:36:35.859Z Venus 120 +2021-05-23T00:10:41.161Z Saturn 240 +2021-05-23T06:36:35.858Z Venus 120 2021-05-23T07:51:24.443Z Mars 90 -2021-05-23T15:59:54.578Z Mercury 120 +2021-05-23T15:59:54.577Z Mercury 120 2021-05-24T04:37:02.655Z Jupiter 240 -2021-05-24T08:28:11.867Z Sun 150 +2021-05-24T08:28:11.866Z Sun 150 2021-05-25T00:34:55.447Z Saturn 270 2021-05-25T10:01:10.768Z Mars 120 2021-05-25T11:01:14.390Z Venus 150 2021-05-25T17:13:39.027Z Mercury 150 -2021-05-26T04:29:33.608Z Jupiter 270 -2021-05-26T11:14:24.797Z Sun 180 -2021-05-27T00:00:58.651Z Saturn 300 +2021-05-26T04:29:33.607Z Jupiter 270 +2021-05-26T11:14:24.796Z Sun 180 +2021-05-27T00:00:58.650Z Saturn 300 2021-05-27T11:31:14.024Z Mars 150 2021-05-27T14:43:34.090Z Venus 180 2021-05-27T17:36:05.802Z Mercury 180 2021-05-28T04:30:48.680Z Jupiter 300 2021-05-28T14:32:15.449Z Sun 210 -2021-05-29T00:23:12.823Z Saturn 330 -2021-05-29T14:36:18.888Z Mars 180 -2021-05-29T19:07:05.905Z Mercury 210 -2021-05-29T20:23:22.440Z Venus 210 -2021-05-30T06:34:01.580Z Jupiter 330 -2021-05-30T20:43:27.790Z Sun 240 -2021-05-31T03:24:48.177Z Saturn 0 -2021-05-31T21:10:58.832Z Mars 210 -2021-05-31T23:13:56.300Z Mercury 240 -2021-06-01T06:14:13.103Z Venus 240 +2021-05-29T00:23:12.822Z Saturn 330 +2021-05-29T14:36:18.887Z Mars 180 +2021-05-29T19:07:05.904Z Mercury 210 +2021-05-29T20:23:22.439Z Venus 210 +2021-05-30T06:34:01.579Z Jupiter 330 +2021-05-30T20:43:27.789Z Sun 240 +2021-05-31T03:24:48.176Z Saturn 0 +2021-05-31T21:10:58.831Z Mars 210 +2021-05-31T23:13:56.299Z Mercury 240 +2021-06-01T06:14:13.102Z Venus 240 2021-06-01T12:03:45.645Z Jupiter 0 -2021-06-02T07:25:02.187Z Sun 270 -2021-06-02T10:07:59.646Z Saturn 30 -2021-06-03T06:24:39.897Z Mercury 270 +2021-06-02T07:25:02.186Z Sun 270 +2021-06-02T10:07:59.645Z Saturn 30 +2021-06-03T06:24:39.896Z Mercury 270 2021-06-03T08:08:56.940Z Mars 240 2021-06-03T21:09:23.029Z Venus 270 2021-06-03T21:22:12.060Z Jupiter 30 -2021-06-04T20:23:37.640Z Saturn 60 -2021-06-04T22:38:30.011Z Sun 300 -2021-06-05T15:57:16.003Z Mercury 300 +2021-06-04T20:23:37.639Z Saturn 60 +2021-06-04T22:38:30.010Z Sun 300 +2021-06-05T15:57:16.002Z Mercury 300 2021-06-05T22:47:34.042Z Mars 270 -2021-06-06T09:32:03.690Z Jupiter 60 +2021-06-06T09:32:03.689Z Jupiter 60 2021-06-06T15:57:31.165Z Venus 300 -2021-06-07T08:51:04.462Z Saturn 90 +2021-06-07T08:51:04.461Z Saturn 90 2021-06-07T16:37:11.096Z Sun 330 -2021-06-08T02:27:50.634Z Mercury 330 -2021-06-08T15:07:18.273Z Mars 300 +2021-06-08T02:27:50.633Z Mercury 330 +2021-06-08T15:07:18.272Z Mars 300 2021-06-08T22:47:04.179Z Jupiter 90 -2021-06-09T12:02:10.110Z Venus 330 -2021-06-09T21:43:15.531Z Saturn 120 +2021-06-09T12:02:10.109Z Venus 330 +2021-06-09T21:43:15.530Z Saturn 120 2021-06-10T10:53:20.717Z Sun 0 2021-06-10T12:38:05.109Z Mercury 0 2021-06-11T07:02:57.745Z Mars 330 -2021-06-11T11:27:40.329Z Jupiter 120 +2021-06-11T11:27:40.328Z Jupiter 120 2021-06-12T06:59:51.450Z Venus 0 2021-06-12T09:34:14.672Z Saturn 150 2021-06-12T21:38:18.839Z Mercury 30 -2021-06-13T03:33:35.931Z Sun 30 +2021-06-13T03:33:35.930Z Sun 30 2021-06-13T21:08:06.168Z Mars 0 -2021-06-13T22:26:25.894Z Jupiter 150 +2021-06-13T22:26:25.893Z Jupiter 150 diff --git a/demo/python/correct/lunar_eclipse.txt b/demo/python/correct/lunar_eclipse.txt index 96df4170..e0404be9 100644 --- a/demo/python/correct/lunar_eclipse.txt +++ b/demo/python/correct/lunar_eclipse.txt @@ -2,49 +2,49 @@ 1988-03-03T16:12:44.158Z Peak of partial eclipse. 1988-03-03T16:20:26.766Z Partial eclipse ends. -1988-08-27T10:07:28.658Z Partial eclipse begins. -1988-08-27T11:04:30.915Z Peak of partial eclipse. -1988-08-27T12:01:33.172Z Partial eclipse ends. +1988-08-27T10:07:28.657Z Partial eclipse begins. +1988-08-27T11:04:30.914Z Peak of partial eclipse. +1988-08-27T12:01:33.171Z Partial eclipse ends. -1989-02-20T13:43:24.718Z Partial eclipse begins. +1989-02-20T13:43:24.717Z Partial eclipse begins. 1989-02-20T14:55:34.839Z Total eclipse begins. -1989-02-20T15:35:19.956Z Peak of total eclipse. -1989-02-20T16:15:05.073Z Total eclipse ends. -1989-02-20T17:27:15.194Z Partial eclipse ends. +1989-02-20T15:35:19.955Z Peak of total eclipse. +1989-02-20T16:15:05.072Z Total eclipse ends. +1989-02-20T17:27:15.193Z Partial eclipse ends. 1989-08-17T01:20:43.868Z Partial eclipse begins. -1989-08-17T02:19:56.976Z Total eclipse begins. +1989-08-17T02:19:56.975Z Total eclipse begins. 1989-08-17T03:08:10.847Z Peak of total eclipse. -1989-08-17T03:56:24.719Z Total eclipse ends. +1989-08-17T03:56:24.718Z Total eclipse ends. 1989-08-17T04:55:37.826Z Partial eclipse ends. 1990-02-09T17:28:36.915Z Partial eclipse begins. 1990-02-09T18:49:12.803Z Total eclipse begins. 1990-02-09T19:11:06.009Z Peak of total eclipse. -1990-02-09T19:32:59.215Z Total eclipse ends. -1990-02-09T20:53:35.103Z Partial eclipse ends. +1990-02-09T19:32:59.214Z Total eclipse ends. +1990-02-09T20:53:35.102Z Partial eclipse ends. -1990-08-06T12:44:10.959Z Partial eclipse begins. -1990-08-06T14:12:20.199Z Peak of partial eclipse. +1990-08-06T12:44:10.958Z Partial eclipse begins. +1990-08-06T14:12:20.198Z Peak of partial eclipse. 1990-08-06T15:40:29.439Z Partial eclipse ends. 1991-12-21T10:00:02.660Z Partial eclipse begins. -1991-12-21T10:33:04.084Z Peak of partial eclipse. +1991-12-21T10:33:04.083Z Peak of partial eclipse. 1991-12-21T11:06:05.507Z Partial eclipse ends. 1992-06-15T03:26:36.541Z Partial eclipse begins. -1992-06-15T04:56:56.392Z Peak of partial eclipse. +1992-06-15T04:56:56.391Z Peak of partial eclipse. 1992-06-15T06:27:16.242Z Partial eclipse ends. -1992-12-09T21:59:22.080Z Partial eclipse begins. +1992-12-09T21:59:22.079Z Partial eclipse begins. 1992-12-09T23:06:41.497Z Total eclipse begins. 1992-12-09T23:44:04.198Z Peak of total eclipse. -1992-12-10T00:21:26.899Z Total eclipse ends. -1992-12-10T01:28:46.317Z Partial eclipse ends. +1992-12-10T00:21:26.898Z Total eclipse ends. +1992-12-10T01:28:46.316Z Partial eclipse ends. -1993-06-04T11:11:10.310Z Partial eclipse begins. +1993-06-04T11:11:10.309Z Partial eclipse begins. 1993-06-04T12:12:10.636Z Total eclipse begins. -1993-06-04T13:00:24.284Z Peak of total eclipse. +1993-06-04T13:00:24.283Z Peak of total eclipse. 1993-06-04T13:48:37.931Z Total eclipse ends. 1993-06-04T14:49:38.257Z Partial eclipse ends. diff --git a/demo/python/correct/moonphase.txt b/demo/python/correct/moonphase.txt index 38091930..ce44af00 100644 --- a/demo/python/correct/moonphase.txt +++ b/demo/python/correct/moonphase.txt @@ -2,13 +2,13 @@ 2019-06-15T09:15:32.987Z : Moon's illuminated fraction = 95.63%. The next 10 lunar quarters are: -2019-06-17T08:31:17.943Z : Full Moon -2019-06-25T09:47:05.867Z : Third Quarter -2019-07-02T19:16:46.600Z : New Moon +2019-06-17T08:31:17.942Z : Full Moon +2019-06-25T09:47:05.866Z : Third Quarter +2019-07-02T19:16:46.599Z : New Moon 2019-07-09T10:55:27.398Z : First Quarter -2019-07-16T21:38:53.593Z : Full Moon +2019-07-16T21:38:53.592Z : Full Moon 2019-07-25T01:18:41.310Z : Third Quarter 2019-08-01T03:12:25.793Z : New Moon 2019-08-07T17:31:34.799Z : First Quarter -2019-08-15T12:29:56.277Z : Full Moon -2019-08-23T14:56:45.685Z : Third Quarter +2019-08-15T12:29:56.276Z : Full Moon +2019-08-23T14:56:45.684Z : Third Quarter diff --git a/demo/python/correct/riseset.txt b/demo/python/correct/riseset.txt index da63a58b..a31a01b4 100644 --- a/demo/python/correct/riseset.txt +++ b/demo/python/correct/riseset.txt @@ -1,5 +1,5 @@ search : 2018-11-30T17:55:07.234Z sunrise : 2018-12-01T13:22:51.098Z -sunset : 2018-11-30T22:21:03.107Z +sunset : 2018-11-30T22:21:03.106Z moonrise : 2018-12-01T06:49:32.983Z -moonset : 2018-11-30T19:22:02.736Z +moonset : 2018-11-30T19:22:02.735Z diff --git a/demo/python/correct/seasons.txt b/demo/python/correct/seasons.txt index 6abbe056..0f818869 100644 --- a/demo/python/correct/seasons.txt +++ b/demo/python/correct/seasons.txt @@ -1,4 +1,4 @@ -March equinox : 2019-03-20T21:58:18.532Z -June solstice : 2019-06-21T15:54:02.462Z +March equinox : 2019-03-20T21:58:18.531Z +June solstice : 2019-06-21T15:54:02.461Z September equinox : 2019-09-23T07:49:57.766Z December solstice : 2019-12-22T04:19:32.601Z diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 264073bc..35440d1a 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -48,6 +48,11 @@ def _cbrt(x: float) -> float: return y raise InternalError() +def _cdiv(numer: int, denom: int) -> int: + '''Divide negative numbers the same way C does: rounding toward zero, not down.''' + sign = -1 if (numer * denom) < 0 else +1 + return sign * (abs(numer) // abs(denom)) + KM_PER_AU = 1.4959787069098932e+8 # The number of kilometers per astronomical unit. C_AUDAY = 173.1446326846693 # The speed of light expressed in astronomical units per day. AU_PER_LY = 63241.07708807546 # The number of astronomical units in one light-year. @@ -275,11 +280,6 @@ def _UniversalTime(tt: float) -> float: _TimeRegex = re.compile(r'^([\+\-]?[0-9]+)-([0-9]{2})-([0-9]{2})(T([0-9]{2}):([0-9]{2})(:([0-9]{2}(\.[0-9]+)?))?Z)?$') -def _cdiv(numer: int, denom: int) -> int: - '''Divide negative numbers the same way C does: rounding toward zero, not down.''' - sign = +1 if (numer * denom) >= 0 else -1 - return sign * (abs(numer) // abs(denom)) - class Time: """Represents a date and time used for performing astronomy calculations. @@ -428,10 +428,10 @@ class Time: d = int(day) f = (14 - m) // 12 y2000 = ( - (d - 2483620) - + _cdiv(1461 * (y + 4800 - f), 4) + (d - 365972956) + + _cdiv(1461 * (y + 1000000 - f), 4) + _cdiv(367 * (m - 2 + f*12), 12) - - _cdiv(3 * _cdiv(y + 4900 - f, 100), 4) + - _cdiv(3 * _cdiv(y + 1000100 - f, 100), 4) ) ut = (y2000 - 0.5) + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0) return Time(ut) @@ -483,37 +483,8 @@ class Time: return 'Time(\'' + str(self) + '\')' def __str__(self) -> str: - # 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) - if x < 0.0: - x += 24.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 - c = 2500 - 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 - month = (80 * k) // 2447 - day = k - 2447 * month // 80 - k = month // 11 - month = month + 2 - 12 * k - year = 100 * (n - 49) + m + k - c*400 - millis = max(0, min(59999, round(1000.0 * second))) + (year, month, day, hour, minute, second) = self.Calendar() + millis = max(0, min(59999, int(math.floor(1000.0 * second)))) if year < 0: text = '-{:06d}'.format(-year) elif year <= 9999: @@ -523,6 +494,38 @@ class Time: text += '-{:02d}-{:02d}T{:02d}:{:02d}:{:02d}.{:03d}Z'.format(month, day, hour, minute, millis // 1000, millis % 1000) return text + def Calendar(self) -> Tuple[int, int, int, int, int, float]: + """Returns a tuple of the form (year, month, day, hour, minute, second). + + This is a convenience method for converting a `Time` value into + its Gregorian calendar date/time representation. + Unlike the built-in `datetime` class, this method can represent + dates over a nearly 2 million year range: the years -999999 to +999999. + """ + # Adapted from the NOVAS C 3.1 function cal_date(). + # [Don Cross - 2023-02-25] Fixed to handle a much wider range of years. + djd = self.ut + 2451545.5 + jd = int(math.floor(djd)) + x = 24.0 * math.fmod(djd, 1.0) + if x < 0.0: + x += 24.0 + hour = int(x) + x = 60.0 * math.fmod(x, 1.0) + minute = int(x) + second = 60.0 * math.fmod(x, 1.0) + c = 2500 + k = jd + (68569 + c*146097) + n = _cdiv(4*k, 146097) + k = k - _cdiv(146097*n + 3, 4) + m = _cdiv(4000*(k+1), 1461001) + k = k - _cdiv(1461 * m, 4) + 31 + month = _cdiv(80 * k, 2447) + day = k - _cdiv(2447 * month, 80) + k = _cdiv(month, 11) + month = month + 2 - 12 * k + year = 100 * (n - 49) + m + k - c*400 + return (year, month, day, hour, minute, second) + def Utc(self) -> datetime.datetime: """Returns the UTC date and time as a `datetime` object. diff --git a/generate/test.py b/generate/test.py index a167af63..4fcf31e9 100755 --- a/generate/test.py +++ b/generate/test.py @@ -77,18 +77,13 @@ def AstroTime(): if s != '2018-12-02 18:30:12.543000': print('PY AstroTime: Utc() returned incorrect string "{}"'.format(s)) return 1 - time = astronomy.Time.Make(2018, 12, 31, 23, 59, 59.9994) + time = astronomy.Time.Make(2018, 12, 31, 23, 59, 59.9999) + expected = '2018-12-31T23:59:59.999Z' 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)) - 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)) + if s != expected: + print('PY AstroTime: expected {} but found {}'.format(expected, s)) 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') AssertGoodTime('2015-01-02T23:45:17Z', '2015-01-02T23:45:17.000Z') AssertGoodTime('1971-03-17T03:30:55.976Z', '1971-03-17T03:30:55.976Z') @@ -100,14 +95,25 @@ def AstroTime(): AssertBadTime('1971-12-31T23:00:60Z') AssertBadTime('1971-03-17T03:30:55.976') # Extreme year values... + AssertGoodTime('-4172-12-02T14:30:45.123Z', '-004172-12-02T14:30:45.123Z') + AssertGoodTime('-4173-12-02T14:30:45.123Z', '-004173-12-02T14:30:45.123Z') + AssertGoodTime('-4174-12-02T14:30:45.123Z', '-004174-12-02T14:30:45.123Z') + AssertGoodTime('-4175-12-02T14:30:45.123Z', '-004175-12-02T14:30:45.123Z') + AssertGoodTime('-4176-12-02T14:30:45.123Z', '-004176-12-02T14:30:45.123Z') AssertGoodTime('-2300-12-19T16:22:26.325Z', '-002300-12-19T16:22:26.325Z') - AssertGoodTime('+12345-12-11T13:30:10.041Z', '+012345-12-11T13:30:10.041Z') - AssertGoodTime('+12345-12-11T13:30:10.041Z', '+012345-12-11T13:30:10.041Z') - AssertGoodTime('+12345-12-11T13:30:10.041Z', '+012345-12-11T13:30:10.041Z') - AssertGoodTime('-123456-01-14T22:55:12.000Z', '-123456-01-14T22:55:12.000Z') - AssertGoodTime('+123456-01-14T22:55:12.000Z', '+123456-01-14T22:55:12.000Z') - AssertGoodTime('-999999-01-14T22:55:12.000Z', '-999999-01-14T22:55:11.999Z') - AssertGoodTime('+999999-01-14T22:55:12.000Z', '+999999-01-14T22:55:11.999Z') + AssertGoodTime('-2300-12-19T16:22:26.325Z', '-002300-12-19T16:22:26.325Z') + AssertGoodTime('+12345-12-11T13:30:10.041Z', '+012345-12-11T13:30:10.040Z') + AssertGoodTime('+12346-12-11T13:30:10.041Z', '+012346-12-11T13:30:10.040Z') + AssertGoodTime('+12347-12-11T13:30:10.041Z', '+012347-12-11T13:30:10.040Z') + AssertGoodTime('+12348-12-11T13:30:10.041Z', '+012348-12-11T13:30:10.040Z') + AssertGoodTime('-123456-01-14T22:55:12.000Z', '-123456-01-14T22:55:11.999Z') + AssertGoodTime('+123456-01-14T22:55:12.000Z', '+123456-01-14T22:55:11.999Z') + AssertGoodTime('-999995-01-14T22:55:12.297Z', '-999995-01-14T22:55:12.297Z') + AssertGoodTime('-999996-01-14T22:55:12.297Z', '-999996-01-14T22:55:12.297Z') + AssertGoodTime('-999997-01-14T22:55:12.297Z', '-999997-01-14T22:55:12.297Z') + AssertGoodTime('-999998-01-14T22:55:12.297Z', '-999998-01-14T22:55:12.297Z') + AssertGoodTime('-999999-01-14T22:55:12.000Z', '-999999-01-14T22:55:11.998Z') + AssertGoodTime('+999999-01-14T22:55:12.000Z', '+999999-01-14T22:55:11.998Z') return 0 #----------------------------------------------------------------------------------------------------------- @@ -3033,7 +3039,7 @@ 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:54.456Z") or + CheckDecemberSolstice( 2022, "2022-12-21T21:47:54.455Z") or CheckDecemberSolstice(-2300, "-002300-12-19T16:22:27.929Z") or CheckDecemberSolstice(12345, "+012345-12-11T13:30:10.276Z") or Pass('DatesIssue250') diff --git a/source/python/README.md b/source/python/README.md index 1db05228..73011fc2 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -947,6 +947,16 @@ The value of the calling object is not modified. This function creates a brand n **Returns**: [`Time`](#Time) + +### Time.Calendar(self) → Tuple\[int, int, int, int, int, float\] + +**Returns a tuple of the form (year, month, day, hour, minute, second).** + +This is a convenience method for converting a `Time` value into +its Gregorian calendar date/time representation. +Unlike the built-in `datetime` class, this method can represent +dates over a nearly 2 million year range: the years -999999 to +999999. + ### Time.FromTerrestrialTime(tt: float) → [`Time`](#Time) diff --git a/source/python/astronomy/astronomy.py b/source/python/astronomy/astronomy.py index e4d030a6..9ca546f7 100644 --- a/source/python/astronomy/astronomy.py +++ b/source/python/astronomy/astronomy.py @@ -48,6 +48,11 @@ def _cbrt(x: float) -> float: return y raise InternalError() +def _cdiv(numer: int, denom: int) -> int: + '''Divide negative numbers the same way C does: rounding toward zero, not down.''' + sign = -1 if (numer * denom) < 0 else +1 + return sign * (abs(numer) // abs(denom)) + KM_PER_AU = 1.4959787069098932e+8 # The number of kilometers per astronomical unit. C_AUDAY = 173.1446326846693 # The speed of light expressed in astronomical units per day. AU_PER_LY = 63241.07708807546 # The number of astronomical units in one light-year. @@ -275,11 +280,6 @@ def _UniversalTime(tt: float) -> float: _TimeRegex = re.compile(r'^([\+\-]?[0-9]+)-([0-9]{2})-([0-9]{2})(T([0-9]{2}):([0-9]{2})(:([0-9]{2}(\.[0-9]+)?))?Z)?$') -def _cdiv(numer: int, denom: int) -> int: - '''Divide negative numbers the same way C does: rounding toward zero, not down.''' - sign = +1 if (numer * denom) >= 0 else -1 - return sign * (abs(numer) // abs(denom)) - class Time: """Represents a date and time used for performing astronomy calculations. @@ -428,10 +428,10 @@ class Time: d = int(day) f = (14 - m) // 12 y2000 = ( - (d - 2483620) - + _cdiv(1461 * (y + 4800 - f), 4) + (d - 365972956) + + _cdiv(1461 * (y + 1000000 - f), 4) + _cdiv(367 * (m - 2 + f*12), 12) - - _cdiv(3 * _cdiv(y + 4900 - f, 100), 4) + - _cdiv(3 * _cdiv(y + 1000100 - f, 100), 4) ) ut = (y2000 - 0.5) + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0) return Time(ut) @@ -483,37 +483,8 @@ class Time: return 'Time(\'' + str(self) + '\')' def __str__(self) -> str: - # 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) - if x < 0.0: - x += 24.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 - c = 2500 - 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 - month = (80 * k) // 2447 - day = k - 2447 * month // 80 - k = month // 11 - month = month + 2 - 12 * k - year = 100 * (n - 49) + m + k - c*400 - millis = max(0, min(59999, round(1000.0 * second))) + (year, month, day, hour, minute, second) = self.Calendar() + millis = max(0, min(59999, int(math.floor(1000.0 * second)))) if year < 0: text = '-{:06d}'.format(-year) elif year <= 9999: @@ -523,6 +494,38 @@ class Time: text += '-{:02d}-{:02d}T{:02d}:{:02d}:{:02d}.{:03d}Z'.format(month, day, hour, minute, millis // 1000, millis % 1000) return text + def Calendar(self) -> Tuple[int, int, int, int, int, float]: + """Returns a tuple of the form (year, month, day, hour, minute, second). + + This is a convenience method for converting a `Time` value into + its Gregorian calendar date/time representation. + Unlike the built-in `datetime` class, this method can represent + dates over a nearly 2 million year range: the years -999999 to +999999. + """ + # Adapted from the NOVAS C 3.1 function cal_date(). + # [Don Cross - 2023-02-25] Fixed to handle a much wider range of years. + djd = self.ut + 2451545.5 + jd = int(math.floor(djd)) + x = 24.0 * math.fmod(djd, 1.0) + if x < 0.0: + x += 24.0 + hour = int(x) + x = 60.0 * math.fmod(x, 1.0) + minute = int(x) + second = 60.0 * math.fmod(x, 1.0) + c = 2500 + k = jd + (68569 + c*146097) + n = _cdiv(4*k, 146097) + k = k - _cdiv(146097*n + 3, 4) + m = _cdiv(4000*(k+1), 1461001) + k = k - _cdiv(1461 * m, 4) + 31 + month = _cdiv(80 * k, 2447) + day = k - _cdiv(2447 * month, 80) + k = _cdiv(month, 11) + month = month + 2 - 12 * k + year = 100 * (n - 49) + m + k - c*400 + return (year, month, day, hour, minute, second) + def Utc(self) -> datetime.datetime: """Returns the UTC date and time as a `datetime` object.