From 0d24433db35c9643530773652a20c525e9ed75b8 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Fri, 8 Apr 2022 16:51:09 -0400 Subject: [PATCH] Fixed #187 - Seasons() failed for distant years. For years before 1582 or years after 3668, the Seasons functions were unable to find many equinoxes and/or solstices. The problem was that over time, the Earth's axis precesses enough that the calendar dates of these events drifts outside the fixed search ranges I had provided for them. I expanded the search ranges so all season changes can be found for a much wider range of years, as verified by unit tests: C/C++: -2000..9999 C#: 1..9999 JavaScript: -2000..9999 Python: 1..9999 Kotlin: 1..9999 Note: C#, Python, and Kotlin currently do not allow years values below +1. In fact, I discovered we were not noticing when an invalid year was passed into the Kotlin code. I updated that code to throw an exception when the year does not match what was expected. It is disturbing that the GregorianCalendar class silently ignores invalid years! Constricted the search tolerance from 1 second to 0.01 seconds for the seasons search, to ensure more consistent behavior. Fixed a bug in the Kotlin search() function's quadratic interpolation that was causing the convergence to be slower than it should have been. --- README.md | 2 +- demo/browser/astronomy.browser.js | 12 +-- demo/c/test/seasons_correct.txt | 2 +- demo/csharp/test/seasons_correct.txt | 8 +- demo/nodejs/astronomy.js | 12 +-- demo/nodejs/calendar/astronomy.ts | 12 +-- .../nodejs/calendar/test/calendar_correct.txt | 94 +++++++++---------- demo/nodejs/test/seasons_correct.txt | 8 +- demo/python/astronomy.py | 19 ++-- demo/python/test/seasons_correct.txt | 8 +- generate/ctest.c | 39 ++++++++ generate/dotnet/csharp_test/csharp_test.cs | 12 +++ generate/template/astronomy.c | 22 +++-- generate/template/astronomy.cs | 22 +++-- generate/template/astronomy.kt | 18 ++-- generate/template/astronomy.py | 19 ++-- generate/template/astronomy.ts | 12 +-- generate/test.js | 21 +++++ generate/test.py | 15 +++ source/c/astronomy.c | 22 +++-- source/csharp/astronomy.cs | 22 +++-- source/js/astronomy.browser.js | 12 +-- source/js/astronomy.browser.min.js | 12 +-- source/js/astronomy.js | 12 +-- source/js/astronomy.min.js | 7 +- source/js/astronomy.ts | 12 +-- source/js/esm/astronomy.js | 12 +-- .../github/cosinekitty/astronomy/astronomy.kt | 18 ++-- .../io/github/cosinekitty/astronomy/Tests.kt | 9 ++ source/python/astronomy/astronomy.py | 19 ++-- 30 files changed, 334 insertions(+), 180 deletions(-) diff --git a/README.md b/README.md index 86f46a9b..79811759 100644 --- a/README.md +++ b/README.md @@ -168,7 +168,7 @@ of complexity. So I decided to create Astronomy Engine with the following engine - Support JavaScript, C, C#, and Python with the same algorithms, and verify them to produce identical results. - No external dependencies! The code must not require anything outside the standard library for each language. -- Minified JavaScript code less than 120K. (The current size is 108883 bytes.) +- Minified JavaScript code less than 120K. (The current size is 108911 bytes.) - Accuracy always within 1 arcminute of results from NOVAS. - It would be well documented, relatively easy to use, and support a wide variety of common use cases. diff --git a/demo/browser/astronomy.browser.js b/demo/browser/astronomy.browser.js index 82b2da90..55eb9de4 100644 --- a/demo/browser/astronomy.browser.js +++ b/demo/browser/astronomy.browser.js @@ -4022,7 +4022,7 @@ function SearchSunLongitude(targetLon, dateStart, limitDays) { VerifyNumber(limitDays); let t1 = MakeTime(dateStart); let t2 = t1.AddDays(limitDays); - return Search(sun_offset, t1, t2); + return Search(sun_offset, t1, t2, { dt_tolerance_seconds: 0.01 }); } exports.SearchSunLongitude = SearchSunLongitude; /** @@ -4921,7 +4921,7 @@ exports.SeasonInfo = SeasonInfo; function Seasons(year) { function find(targetLon, month, day) { let startDate = new Date(Date.UTC(year, month - 1, day)); - let time = SearchSunLongitude(targetLon, startDate, 4); + let time = SearchSunLongitude(targetLon, startDate, 20); if (!time) throw `Cannot find season change near ${startDate.toISOString()}`; return time; @@ -4930,10 +4930,10 @@ function Seasons(year) { year = year.getUTCFullYear(); if (!Number.isSafeInteger(year)) throw `Cannot calculate seasons because year argument ${year} is neither a Date nor a safe integer.`; - let mar_equinox = find(0, 3, 19); - let jun_solstice = find(90, 6, 19); - let sep_equinox = find(180, 9, 21); - let dec_solstice = find(270, 12, 20); + let mar_equinox = find(0, 3, 10); + let jun_solstice = find(90, 6, 10); + let sep_equinox = find(180, 9, 10); + let dec_solstice = find(270, 12, 10); return new SeasonInfo(mar_equinox, jun_solstice, sep_equinox, dec_solstice); } exports.Seasons = Seasons; diff --git a/demo/c/test/seasons_correct.txt b/demo/c/test/seasons_correct.txt index 2b35489d..c9ea215b 100644 --- a/demo/c/test/seasons_correct.txt +++ b/demo/c/test/seasons_correct.txt @@ -1,4 +1,4 @@ March equinox : 2020-03-20T03:49:57Z June solstice : 2020-06-20T21:43:34Z -September equinox : 2020-09-22T13:30:56Z +September equinox : 2020-09-22T13:30:57Z December solstice : 2020-12-21T10:02:41Z diff --git a/demo/csharp/test/seasons_correct.txt b/demo/csharp/test/seasons_correct.txt index 910bdd28..591050bb 100644 --- a/demo/csharp/test/seasons_correct.txt +++ b/demo/csharp/test/seasons_correct.txt @@ -1,4 +1,4 @@ -March equinox : 2020-03-20T03:49:56.640Z -June solstice : 2020-06-20T21:43:33.612Z -September equinox : 2020-09-22T13:30:56.485Z -December solstice : 2020-12-21T10:02:40.753Z +March equinox : 2020-03-20T03:49:57.177Z +June solstice : 2020-06-20T21:43:33.689Z +September equinox : 2020-09-22T13:30:56.763Z +December solstice : 2020-12-21T10:02:40.981Z diff --git a/demo/nodejs/astronomy.js b/demo/nodejs/astronomy.js index 008a1db9..264f1041 100644 --- a/demo/nodejs/astronomy.js +++ b/demo/nodejs/astronomy.js @@ -4021,7 +4021,7 @@ function SearchSunLongitude(targetLon, dateStart, limitDays) { VerifyNumber(limitDays); let t1 = MakeTime(dateStart); let t2 = t1.AddDays(limitDays); - return Search(sun_offset, t1, t2); + return Search(sun_offset, t1, t2, { dt_tolerance_seconds: 0.01 }); } exports.SearchSunLongitude = SearchSunLongitude; /** @@ -4920,7 +4920,7 @@ exports.SeasonInfo = SeasonInfo; function Seasons(year) { function find(targetLon, month, day) { let startDate = new Date(Date.UTC(year, month - 1, day)); - let time = SearchSunLongitude(targetLon, startDate, 4); + let time = SearchSunLongitude(targetLon, startDate, 20); if (!time) throw `Cannot find season change near ${startDate.toISOString()}`; return time; @@ -4929,10 +4929,10 @@ function Seasons(year) { year = year.getUTCFullYear(); if (!Number.isSafeInteger(year)) throw `Cannot calculate seasons because year argument ${year} is neither a Date nor a safe integer.`; - let mar_equinox = find(0, 3, 19); - let jun_solstice = find(90, 6, 19); - let sep_equinox = find(180, 9, 21); - let dec_solstice = find(270, 12, 20); + let mar_equinox = find(0, 3, 10); + let jun_solstice = find(90, 6, 10); + let sep_equinox = find(180, 9, 10); + let dec_solstice = find(270, 12, 10); return new SeasonInfo(mar_equinox, jun_solstice, sep_equinox, dec_solstice); } exports.Seasons = Seasons; diff --git a/demo/nodejs/calendar/astronomy.ts b/demo/nodejs/calendar/astronomy.ts index 68df29ef..b8e20305 100644 --- a/demo/nodejs/calendar/astronomy.ts +++ b/demo/nodejs/calendar/astronomy.ts @@ -4491,7 +4491,7 @@ export function SearchSunLongitude(targetLon: number, dateStart: FlexibleDateTim VerifyNumber(limitDays); let t1 = MakeTime(dateStart); let t2 = t1.AddDays(limitDays); - return Search(sun_offset, t1, t2); + return Search(sun_offset, t1, t2, {dt_tolerance_seconds: 0.01}); } /** @@ -5439,7 +5439,7 @@ export class SeasonInfo { export function Seasons(year: (number | AstroTime)): SeasonInfo { function find(targetLon: number, month: number, day: number): AstroTime { let startDate = new Date(Date.UTC(year, month-1, day)); - let time = SearchSunLongitude(targetLon, startDate, 4); + let time = SearchSunLongitude(targetLon, startDate, 20); if (!time) throw `Cannot find season change near ${startDate.toISOString()}`; return time; @@ -5451,10 +5451,10 @@ export function Seasons(year: (number | AstroTime)): SeasonInfo { if (!Number.isSafeInteger(year)) throw `Cannot calculate seasons because year argument ${year} is neither a Date nor a safe integer.`; - let mar_equinox = find( 0, 3, 19); - let jun_solstice = find( 90, 6, 19); - let sep_equinox = find(180, 9, 21); - let dec_solstice = find(270, 12, 20); + let mar_equinox = find( 0, 3, 10); + let jun_solstice = find( 90, 6, 10); + let sep_equinox = find(180, 9, 10); + let dec_solstice = find(270, 12, 10); return new SeasonInfo(mar_equinox, jun_solstice, sep_equinox, dec_solstice); } diff --git a/demo/nodejs/calendar/test/calendar_correct.txt b/demo/nodejs/calendar/test/calendar_correct.txt index cbbe91a2..2dc973f6 100644 --- a/demo/nodejs/calendar/test/calendar_correct.txt +++ b/demo/nodejs/calendar/test/calendar_correct.txt @@ -250,7 +250,7 @@ 2021-06-20T20:05:54.827Z moonrise 2021-06-21T00:25:27.970Z sunset 2021-06-21T01:54:48.822Z Moon culminates 49.11 degrees above the horizon -2021-06-21T03:31:54.327Z June solstice +2021-06-21T03:31:54.546Z June solstice 2021-06-21T07:37:22.673Z moonset 2021-06-21T10:27:43.521Z sunrise 2021-06-21T17:26:42.212Z Sun culminates 84.84 degrees above the horizon @@ -846,7 +846,7 @@ 2021-09-22T11:13:22.248Z sunrise 2021-09-22T12:36:15.937Z moonset 2021-09-22T17:17:20.942Z Sun culminates 61.44 degrees above the horizon -2021-09-22T19:20:57.675Z September equinox +2021-09-22T19:20:57.419Z September equinox 2021-09-22T23:20:53.939Z sunset 2021-09-23T00:36:21.105Z moonrise 2021-09-23T07:00:28.922Z Moon culminates 68.71 degrees above the horizon @@ -1419,7 +1419,7 @@ 2021-12-21T07:04:52.754Z Moon culminates 86.91 degrees above the horizon 2021-12-21T12:13:23.197Z sunrise 2021-12-21T14:18:16.236Z moonset -2021-12-21T15:59:16.343Z December solstice +2021-12-21T15:59:16.328Z December solstice 2021-12-21T17:23:03.965Z Sun culminates 37.98 degrees above the horizon 2021-12-21T22:32:44.770Z sunset 2021-12-22T00:44:17.952Z moonrise @@ -1985,7 +1985,7 @@ 2022-03-20T09:32:04.911Z Venus max morning elongation: 46.59 degrees from Sun 2022-03-20T11:28:37.428Z sunrise 2022-03-20T13:01:12.810Z moonset -2022-03-20T15:33:24.790Z March equinox +2022-03-20T15:33:24.472Z March equinox 2022-03-20T17:32:11.528Z Sun culminates 61.44 degrees above the horizon 2022-03-20T23:36:11.734Z sunset 2022-03-21T02:15:18.771Z moonrise @@ -2580,7 +2580,7 @@ 2022-06-21T00:25:25.834Z sunset 2022-06-21T03:11:25.170Z third quarter 2022-06-21T05:45:30.915Z moonrise -2022-06-21T09:13:42.542Z June solstice +2022-06-21T09:13:42.944Z June solstice 2022-06-21T10:27:41.231Z sunrise 2022-06-21T11:48:43.132Z Moon culminates 59.19 degrees above the horizon 2022-06-21T12:27:14.976Z Mars perihelion at 1.3813 AU @@ -3177,7 +3177,7 @@ 2022-09-22T17:17:26.649Z Sun culminates 61.53 degrees above the horizon 2022-09-22T21:53:48.646Z moonset 2022-09-22T23:21:11.783Z sunset -2022-09-23T01:04:04.836Z September equinox +2022-09-23T01:04:04.807Z September equinox 2022-09-23T06:43:36.400Z Mercury inferior conjunction 2022-09-23T08:58:34.119Z moonrise 2022-09-23T11:13:45.681Z sunrise @@ -3746,7 +3746,7 @@ 2022-12-21T15:35:00.203Z Moon culminates 37.55 degrees above the horizon 2022-12-21T17:22:56.536Z Sun culminates 37.98 degrees above the horizon 2022-12-21T20:50:58.914Z moonset -2022-12-21T21:47:58.017Z December solstice +2022-12-21T21:47:58.189Z December solstice 2022-12-21T22:32:37.243Z sunset 2022-12-22T11:25:59.239Z moonrise 2022-12-22T12:13:45.531Z sunrise @@ -4315,7 +4315,7 @@ 2023-03-20T11:28:53.658Z sunrise 2023-03-20T16:46:45.681Z Moon culminates 51.61 degrees above the horizon 2023-03-20T17:32:15.198Z Sun culminates 61.34 degrees above the horizon -2023-03-20T21:24:24.032Z March equinox +2023-03-20T21:24:24.091Z March equinox 2023-03-20T22:41:34.313Z moonset 2023-03-20T23:36:02.848Z sunset 2023-03-21T11:27:44.328Z sunrise @@ -4904,7 +4904,7 @@ 2023-06-21T02:49:28.403Z moonset 2023-06-21T10:27:38.619Z sunrise 2023-06-21T13:18:49.524Z moonrise -2023-06-21T14:57:22.055Z June solstice +2023-06-21T14:57:21.801Z June solstice 2023-06-21T17:26:37.554Z Sun culminates 84.84 degrees above the horizon 2023-06-21T20:26:01.246Z Moon culminates 83.55 degrees above the horizon 2023-06-22T00:25:36.419Z sunset @@ -5503,7 +5503,7 @@ 2023-09-22T23:21:29.277Z sunset 2023-09-22T23:26:33.979Z Moon culminates 32.37 degrees above the horizon 2023-09-23T04:29:55.087Z moonset -2023-09-23T06:50:00.573Z September equinox +2023-09-23T06:50:00.524Z September equinox 2023-09-23T11:13:38.657Z sunrise 2023-09-23T17:17:10.830Z Sun culminates 61.24 degrees above the horizon 2023-09-23T18:01:16.890Z Mercury perihelion at 0.3075 AU @@ -6078,7 +6078,7 @@ 2023-12-21T18:36:31.933Z moonrise 2023-12-21T22:32:31.386Z sunset 2023-12-22T01:10:38.802Z Moon culminates 72.79 degrees above the horizon -2023-12-22T03:27:34.523Z December solstice +2023-12-22T03:27:34.528Z December solstice 2023-12-22T07:52:41.956Z moonset 2023-12-22T12:13:39.811Z sunrise 2023-12-22T17:23:20.528Z Sun culminates 37.98 degrees above the horizon @@ -6646,7 +6646,7 @@ 2024-03-19T22:20:15.419Z Venus aphelion at 0.7282 AU 2024-03-19T23:35:55.569Z sunset 2024-03-20T01:59:43.676Z Moon culminates 85.64 degrees above the horizon -2024-03-20T03:06:23.649Z March equinox +2024-03-20T03:06:23.283Z March equinox 2024-03-20T09:07:26.187Z moonset 2024-03-20T11:28:01.551Z sunrise 2024-03-20T17:32:02.424Z Sun culminates 61.65 degrees above the horizon @@ -7240,7 +7240,7 @@ 2024-06-20T08:51:58.150Z moonset 2024-06-20T10:27:34.413Z sunrise 2024-06-20T17:26:33.246Z Sun culminates 84.84 degrees above the horizon -2024-06-20T20:51:01.134Z June solstice +2024-06-20T20:51:00.974Z June solstice 2024-06-20T23:32:56.409Z moonrise 2024-06-21T00:25:32.168Z sunset 2024-06-21T04:38:39.212Z Moon culminates 32.98 degrees above the horizon @@ -7834,7 +7834,7 @@ 2024-09-22T01:52:58.429Z moonrise 2024-09-22T09:01:49.397Z Moon culminates 85.41 degrees above the horizon 2024-09-22T11:13:30.169Z sunrise -2024-09-22T12:43:52.500Z September equinox +2024-09-22T12:43:52.384Z September equinox 2024-09-22T16:17:28.607Z moonset 2024-09-22T17:17:14.840Z Sun culminates 61.33 degrees above the horizon 2024-09-22T18:34:05.965Z Mercury moves from Leo to Virgo @@ -8403,7 +8403,7 @@ 2024-12-20T17:22:42.423Z Sun culminates 37.98 degrees above the horizon 2024-12-20T22:32:23.166Z sunset 2024-12-21T03:54:09.780Z moonrise -2024-12-21T09:20:17.831Z December solstice +2024-12-21T09:20:17.849Z December solstice 2024-12-21T10:25:14.168Z Moon culminates 68.97 degrees above the horizon 2024-12-21T12:13:31.581Z sunrise 2024-12-21T16:49:16.390Z moonset @@ -8972,7 +8972,7 @@ 2025-03-19T23:31:13.377Z Neptune conjunction 2025-03-19T23:35:45.833Z sunset 2025-03-20T04:39:08.526Z moonrise -2025-03-20T09:01:30.996Z March equinox +2025-03-20T09:01:31.049Z March equinox 2025-03-20T09:47:57.260Z Moon culminates 34.36 degrees above the horizon 2025-03-20T11:28:17.326Z sunrise 2025-03-20T14:53:47.886Z moonset @@ -9567,7 +9567,7 @@ 2025-06-20T17:26:29.253Z Sun culminates 84.84 degrees above the horizon 2025-06-20T19:28:04.410Z moonset 2025-06-21T00:25:28.225Z sunset -2025-06-21T02:42:19.276Z June solstice +2025-06-21T02:42:19.310Z June solstice 2025-06-21T06:50:05.717Z moonrise 2025-06-21T10:27:43.491Z sunrise 2025-06-21T13:39:02.735Z Moon culminates 78.46 degrees above the horizon @@ -10161,7 +10161,7 @@ 2025-09-22T11:52:34.164Z moonrise 2025-09-22T17:17:21.559Z Sun culminates 61.42 degrees above the horizon 2025-09-22T17:52:13.582Z Moon culminates 54.99 degrees above the horizon -2025-09-22T18:19:32.605Z September equinox +2025-09-22T18:19:32.741Z September equinox 2025-09-22T23:20:52.413Z sunset 2025-09-22T23:45:26.827Z moonset 2025-09-23T11:13:54.939Z sunrise @@ -10732,7 +10732,7 @@ 2025-12-20T23:00:13.470Z moonset 2025-12-21T12:13:25.076Z sunrise 2025-12-21T13:41:04.862Z moonrise -2025-12-21T15:03:04.888Z December solstice +2025-12-21T15:03:05.499Z December solstice 2025-12-21T17:23:05.773Z Sun culminates 37.98 degrees above the horizon 2025-12-21T18:48:04.656Z Moon culminates 34.67 degrees above the horizon 2025-12-21T22:32:46.527Z sunset @@ -11302,7 +11302,7 @@ 2026-03-20T00:29:47.663Z moonset 2026-03-20T11:28:34.327Z sunrise 2026-03-20T12:14:56.110Z moonrise -2026-03-20T14:45:39.187Z March equinox +2026-03-20T14:45:39.417Z March equinox 2026-03-20T17:32:10.160Z Sun culminates 61.45 degrees above the horizon 2026-03-20T18:50:57.299Z Moon culminates 73.43 degrees above the horizon 2026-03-20T23:36:12.114Z sunset @@ -11896,7 +11896,7 @@ 2026-06-20T22:41:46.090Z Moon culminates 64.45 degrees above the horizon 2026-06-21T00:25:25.634Z sunset 2026-06-21T04:56:19.042Z moonset -2026-06-21T08:24:58.639Z June solstice +2026-06-21T08:24:58.859Z June solstice 2026-06-21T10:27:41.003Z sunrise 2026-06-21T17:17:48.662Z moonrise 2026-06-21T17:26:39.901Z Sun culminates 84.84 degrees above the horizon @@ -12490,7 +12490,7 @@ 2026-09-22T17:17:26.756Z Sun culminates 61.52 degrees above the horizon 2026-09-22T21:07:14.236Z moonrise 2026-09-22T23:21:09.798Z sunset -2026-09-23T00:05:37.446Z September equinox +2026-09-23T00:05:37.449Z September equinox 2026-09-23T02:36:15.982Z Moon culminates 43.93 degrees above the horizon 2026-09-23T08:10:45.052Z moonset 2026-09-23T11:13:47.742Z sunrise @@ -13058,7 +13058,7 @@ 2026-12-21T12:13:17.718Z sunrise 2026-12-21T17:22:58.529Z Sun culminates 37.98 degrees above the horizon 2026-12-21T20:06:01.240Z moonrise -2026-12-21T20:50:25.605Z December solstice +2026-12-21T20:50:25.604Z December solstice 2026-12-21T22:32:39.271Z sunset 2026-12-22T03:21:43.636Z Moon culminates 87.03 degrees above the horizon 2026-12-22T10:42:37.552Z moonset @@ -13627,7 +13627,7 @@ 2027-03-20T10:11:06.397Z moonset 2027-03-20T11:28:51.943Z sunrise 2027-03-20T17:32:15.602Z Sun culminates 61.36 degrees above the horizon -2027-03-20T20:24:44.239Z March equinox +2027-03-20T20:24:44.049Z March equinox 2027-03-20T22:01:50.267Z moonrise 2027-03-20T23:36:05.352Z sunset 2027-03-21T04:27:32.344Z Moon culminates 65.50 degrees above the horizon @@ -14218,7 +14218,7 @@ 2027-06-21T07:15:26.447Z Moon culminates 38.01 degrees above the horizon 2027-06-21T10:27:38.231Z sunrise 2027-06-21T12:34:12.004Z moonset -2027-06-21T14:10:32.767Z June solstice +2027-06-21T14:10:33.482Z June solstice 2027-06-21T17:26:36.976Z Sun culminates 84.84 degrees above the horizon 2027-06-22T00:25:35.629Z sunset 2027-06-22T02:37:25.069Z moonrise @@ -14814,7 +14814,7 @@ 2027-09-22T17:41:03.231Z moonset 2027-09-22T23:21:26.186Z sunset 2027-09-23T03:58:23.936Z moonrise -2027-09-23T06:01:18.191Z September equinox +2027-09-23T06:01:18.126Z September equinox 2027-09-23T10:21:01.040Z third quarter 2027-09-23T11:13:39.018Z sunrise 2027-09-23T11:19:47.793Z Moon culminates 88.08 degrees above the horizon @@ -15388,7 +15388,7 @@ 2027-12-21T17:22:50.446Z Sun culminates 37.98 degrees above the horizon 2027-12-21T17:53:41.339Z moonset 2027-12-21T22:32:31.332Z sunset -2027-12-22T02:42:19.553Z December solstice +2027-12-22T02:42:19.485Z December solstice 2027-12-22T07:18:06.115Z moonrise 2027-12-22T12:13:39.336Z sunrise 2027-12-22T12:56:38.569Z Moon culminates 45.24 degrees above the horizon @@ -15955,7 +15955,7 @@ 2028-03-19T17:32:19.110Z Sun culminates 61.26 degrees above the horizon 2028-03-19T18:05:12.812Z moonset 2028-03-19T23:35:56.241Z sunset -2028-03-20T02:16:59.224Z March equinox +2028-03-20T02:16:59.678Z March equinox 2028-03-20T08:15:35.560Z moonrise 2028-03-20T11:27:58.715Z sunrise 2028-03-20T13:36:24.296Z Moon culminates 40.33 degrees above the horizon @@ -16547,7 +16547,7 @@ 2028-06-20T10:27:34.632Z sunrise 2028-06-20T15:26:47.813Z Moon culminates 86.02 degrees above the horizon 2028-06-20T17:26:33.215Z Sun culminates 84.84 degrees above the horizon -2028-06-20T20:01:44.806Z June solstice +2028-06-20T20:01:44.926Z June solstice 2028-06-20T22:42:03.418Z moonset 2028-06-21T00:25:31.864Z sunset 2028-06-21T09:07:18.512Z moonrise @@ -17142,7 +17142,7 @@ 2028-09-21T23:21:43.956Z sunset 2028-09-22T01:17:15.345Z moonset 2028-09-22T11:13:32.488Z sunrise -2028-09-22T11:45:20.700Z September equinox +2028-09-22T11:45:21.123Z September equinox 2028-09-22T11:52:55.373Z Venus moves from Cancer to Leo 2028-09-22T15:33:06.268Z moonrise 2028-09-22T17:17:15.125Z Sun culminates 61.32 degrees above the horizon @@ -17712,7 +17712,7 @@ 2028-12-20T21:14:40.661Z Moon culminates 51.50 degrees above the horizon 2028-12-20T22:32:26.175Z sunset 2028-12-21T03:06:21.535Z moonset -2028-12-21T08:20:05.716Z December solstice +2028-12-21T08:20:06.162Z December solstice 2028-12-21T10:30:38.570Z Venus moves from Scorpius to Ophiuchus 2028-12-21T12:13:33.904Z sunrise 2028-12-21T15:58:50.469Z moonrise @@ -18281,7 +18281,7 @@ 2029-03-19T20:42:04.201Z Moon culminates 82.85 degrees above the horizon 2029-03-19T23:35:48.656Z sunset 2029-03-20T03:46:59.967Z moonset -2029-03-20T08:01:42.676Z March equinox +2029-03-20T08:01:43.024Z March equinox 2029-03-20T11:28:15.894Z sunrise 2029-03-20T14:23:12.066Z moonrise 2029-03-20T17:32:06.204Z Sun culminates 61.56 degrees above the horizon @@ -18877,7 +18877,7 @@ 2029-06-20T18:53:05.772Z moonrise 2029-06-21T00:25:28.842Z sunset 2029-06-21T00:38:21.510Z Moon culminates 47.70 degrees above the horizon -2029-06-21T01:48:18.503Z June solstice +2029-06-21T01:48:18.598Z June solstice 2029-06-21T06:17:56.937Z moonset 2029-06-21T10:27:44.966Z sunrise 2029-06-21T17:26:43.367Z Sun culminates 84.84 degrees above the horizon @@ -19470,7 +19470,7 @@ 2029-09-22T11:13:24.531Z sunrise 2029-09-22T16:29:56.389Z full moon 2029-09-22T17:17:19.582Z Sun culminates 61.41 degrees above the horizon -2029-09-22T17:37:48.793Z September equinox +2029-09-22T17:37:48.426Z September equinox 2029-09-22T22:27:07.287Z Mars moves from Libra to Scorpius 2029-09-22T23:08:08.095Z moonrise 2029-09-22T23:20:48.942Z sunset @@ -20049,7 +20049,7 @@ 2029-12-21T05:36:50.259Z Moon culminates 84.04 degrees above the horizon 2029-12-21T12:13:25.508Z sunrise 2029-12-21T12:42:30.714Z moonset -2029-12-21T14:14:11.718Z December solstice +2029-12-21T14:14:11.679Z December solstice 2029-12-21T17:23:06.661Z Sun culminates 37.98 degrees above the horizon 2029-12-21T22:32:47.891Z sunset 2029-12-21T23:24:34.551Z moonrise @@ -20613,7 +20613,7 @@ 2030-03-20T05:49:28.151Z Moon culminates 53.86 degrees above the horizon 2030-03-20T11:28:32.530Z sunrise 2030-03-20T11:42:59.233Z moonset -2030-03-20T13:51:45.836Z March equinox +2030-03-20T13:51:45.343Z March equinox 2030-03-20T17:32:10.274Z Sun culminates 61.47 degrees above the horizon 2030-03-20T23:36:14.122Z sunset 2030-03-21T00:56:06.970Z moonrise @@ -21208,7 +21208,7 @@ 2030-06-20T17:26:28.207Z Sun culminates 84.83 degrees above the horizon 2030-06-21T00:25:26.490Z sunset 2030-06-21T04:16:49.759Z moonrise -2030-06-21T07:31:12.657Z June solstice +2030-06-21T07:31:13.067Z June solstice 2030-06-21T10:22:19.859Z Moon culminates 59.86 degrees above the horizon 2030-06-21T10:27:42.917Z sunrise 2030-06-21T16:33:29.378Z moonset @@ -21807,7 +21807,7 @@ 2030-09-22T17:17:24.626Z Sun culminates 61.51 degrees above the horizon 2030-09-22T20:10:48.664Z moonset 2030-09-22T23:21:06.335Z sunset -2030-09-22T23:27:12.861Z September equinox +2030-09-22T23:27:12.874Z September equinox 2030-09-23T07:28:04.211Z moonrise 2030-09-23T11:13:47.066Z sunrise 2030-09-23T14:10:56.085Z Moon culminates 73.78 degrees above the horizon @@ -22385,7 +22385,7 @@ 2030-12-21T14:14:59.516Z Moon culminates 42.27 degrees above the horizon 2030-12-21T17:22:58.302Z Sun culminates 37.99 degrees above the horizon 2030-12-21T19:43:40.529Z moonset -2030-12-21T20:09:37.646Z December solstice +2030-12-21T20:09:37.885Z December solstice 2030-12-21T22:32:39.648Z sunset 2030-12-22T09:50:07.779Z moonrise 2030-12-22T12:13:46.610Z sunrise @@ -22950,7 +22950,7 @@ 2031-03-20T11:28:49.344Z sunrise 2031-03-20T15:24:22.488Z Moon culminates 53.22 degrees above the horizon 2031-03-20T17:32:14.603Z Sun culminates 61.37 degrees above the horizon -2031-03-20T19:40:52.626Z March equinox +2031-03-20T19:40:52.746Z March equinox 2031-03-20T21:21:53.062Z moonset 2031-03-20T23:36:05.971Z sunset 2031-03-21T10:11:23.400Z moonrise @@ -23543,7 +23543,7 @@ 2031-06-21T10:27:40.559Z sunrise 2031-06-21T11:23:30.922Z lunar apogee at 406441 km 2031-06-21T11:56:48.300Z moonrise -2031-06-21T13:17:00.933Z June solstice +2031-06-21T13:17:00.604Z June solstice 2031-06-21T17:26:38.756Z Sun culminates 84.83 degrees above the horizon 2031-06-21T18:51:13.718Z Moon culminates 79.04 degrees above the horizon 2031-06-22T00:25:36.835Z sunset @@ -24136,7 +24136,7 @@ 2031-09-22T22:14:41.638Z Moon culminates 40.23 degrees above the horizon 2031-09-22T23:21:24.022Z sunset 2031-09-23T03:39:25.429Z moonset -2031-09-23T05:15:16.557Z September equinox +2031-09-23T05:15:16.472Z September equinox 2031-09-23T11:13:40.113Z sunrise 2031-09-23T17:17:08.943Z Sun culminates 61.21 degrees above the horizon 2031-09-23T17:47:20.683Z moonrise @@ -24716,7 +24716,7 @@ 2031-12-21T17:37:59.103Z moonrise 2031-12-21T22:32:33.209Z sunset 2031-12-22T00:04:18.836Z Moon culminates 68.77 degrees above the horizon -2031-12-22T01:55:48.509Z December solstice +2031-12-22T01:55:48.488Z December solstice 2031-12-22T06:36:00.854Z moonset 2031-12-22T12:13:40.006Z sunrise 2031-12-22T17:23:21.580Z Sun culminates 37.99 degrees above the horizon @@ -25283,7 +25283,7 @@ 2032-03-19T17:32:20.307Z Sun culminates 61.28 degrees above the horizon 2032-03-19T23:35:59.418Z sunset 2032-03-20T00:25:33.560Z Moon culminates 79.52 degrees above the horizon -2032-03-20T01:21:49.414Z March equinox +2032-03-20T01:21:49.319Z March equinox 2032-03-20T07:18:31.413Z moonset 2032-03-20T11:27:57.906Z sunrise 2032-03-20T17:32:02.543Z Sun culminates 61.67 degrees above the horizon @@ -25874,7 +25874,7 @@ 2032-06-20T07:55:41.883Z moonset 2032-06-20T10:27:37.386Z sunrise 2032-06-20T17:26:35.344Z Sun culminates 84.83 degrees above the horizon -2032-06-20T19:08:39.855Z June solstice +2032-06-20T19:08:39.587Z June solstice 2032-06-20T21:44:10.492Z moonrise 2032-06-21T00:25:33.343Z sunset 2032-06-21T03:15:26.369Z Moon culminates 42.48 degrees above the horizon @@ -26472,7 +26472,7 @@ 2032-09-21T23:21:40.977Z sunset 2032-09-22T01:06:22.774Z moonrise 2032-09-22T07:49:39.203Z Moon culminates 75.19 degrees above the horizon -2032-09-22T11:10:51.450Z September equinox +2032-09-22T11:10:51.066Z September equinox 2032-09-22T11:13:31.891Z sunrise 2032-09-22T14:37:19.784Z moonset 2032-09-22T17:17:13.269Z Sun culminates 61.31 degrees above the horizon @@ -27046,7 +27046,7 @@ 2032-12-20T22:32:24.336Z sunset 2032-12-20T22:55:55.221Z Uranus opposition 2032-12-21T02:19:57.353Z moonrise -2032-12-21T07:56:05.231Z December solstice +2032-12-21T07:56:05.170Z December solstice 2032-12-21T08:54:51.270Z Moon culminates 70.83 degrees above the horizon 2032-12-21T12:13:30.894Z sunrise 2032-12-21T15:25:18.634Z moonset diff --git a/demo/nodejs/test/seasons_correct.txt b/demo/nodejs/test/seasons_correct.txt index a9773026..7fda0178 100644 --- a/demo/nodejs/test/seasons_correct.txt +++ b/demo/nodejs/test/seasons_correct.txt @@ -1,4 +1,4 @@ -March equinox : 2019-03-20T21:58:19.365Z -June solstice : 2019-06-21T15:54:01.251Z -September equinox : 2019-09-23T07:50:00.100Z -December solstice : 2019-12-22T04:19:32.430Z +March equinox : 2019-03-20T21:58:19.248Z +June solstice : 2019-06-21T15:54:02.000Z +September equinox : 2019-09-23T07:49:59.929Z +December solstice : 2019-12-22T04:19:32.258Z diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py index 616504e7..651f9c99 100644 --- a/demo/python/astronomy.py +++ b/demo/python/astronomy.py @@ -5642,7 +5642,7 @@ def SearchSunLongitude(targetLon, startTime, limitDays): Time or `None` """ t2 = startTime.AddDays(limitDays) - return Search(_sun_offset, targetLon, startTime, t2, 1.0) + return Search(_sun_offset, targetLon, startTime, t2, 0.01) def MoonPhase(time): """Returns the Moon's phase as an angle from 0 to 360 degrees. @@ -6444,7 +6444,7 @@ class SeasonInfo: def _FindSeasonChange(targetLon, year, month, day): startTime = Time.Make(year, month, day, 0, 0, 0) - time = SearchSunLongitude(targetLon, startTime, 4.0) + time = SearchSunLongitude(targetLon, startTime, 20.0) if time is None: # We should always be able to find a season change. raise InternalError() @@ -6486,10 +6486,17 @@ def Seasons(year): ------- SeasonInfo """ - mar_equinox = _FindSeasonChange(0, year, 3, 19) - jun_solstice = _FindSeasonChange(90, year, 6, 19) - sep_equinox = _FindSeasonChange(180, year, 9, 21) - dec_solstice = _FindSeasonChange(270, year, 12, 20) + # https://github.com/cosinekitty/astronomy/issues/187 + # Solstices and equinoxes drift over long spans of time, + # due to precession of the Earth's axis. + # Therefore, we have to search a wider range of time than + # one might expect. It turns out this has very little + # effect on efficiency, thanks to the quick convergence + # of quadratic interpolation inside the `Search` function. + mar_equinox = _FindSeasonChange( 0, year, 3, 10) + jun_solstice = _FindSeasonChange( 90, year, 6, 10) + sep_equinox = _FindSeasonChange(180, year, 9, 10) + dec_solstice = _FindSeasonChange(270, year, 12, 10) return SeasonInfo(mar_equinox, jun_solstice, sep_equinox, dec_solstice) def _MoonDistance(time): diff --git a/demo/python/test/seasons_correct.txt b/demo/python/test/seasons_correct.txt index 52f68f14..15a45ba0 100644 --- a/demo/python/test/seasons_correct.txt +++ b/demo/python/test/seasons_correct.txt @@ -1,4 +1,4 @@ -March equinox : 2019-03-20T21:58:19.366Z -June solstice : 2019-06-21T15:54:01.251Z -September equinox : 2019-09-23T07:50:00.100Z -December solstice : 2019-12-22T04:19:32.431Z +March equinox : 2019-03-20T21:58:19.248Z +June solstice : 2019-06-21T15:54:02.001Z +September equinox : 2019-09-23T07:49:59.929Z +December solstice : 2019-12-22T04:19:32.258Z diff --git a/generate/ctest.c b/generate/ctest.c index 864eaf4f..2d7b5878 100644 --- a/generate/ctest.c +++ b/generate/ctest.c @@ -172,6 +172,7 @@ static int AstroCheck(void); static int Diff(double tolerance, const char *a_filename, const char *b_filename); static int DiffLine(int lnum, const char *aline, const char *bline, maxdiff_column_t column[]); static int SeasonsTest(void); +static int SeasonsIssue187(void); static int MoonPhase(void); static int MoonNodes(void); static int RiseSet(void); @@ -261,6 +262,7 @@ static unit_test_t UnitTests[] = {"riseset", RiseSet}, {"rotation", RotationTest}, {"seasons", SeasonsTest}, + {"seasons187", SeasonsIssue187}, {"sidereal", SiderealTimeTest}, {"time", Test_AstroTime}, {"topostate", TopoStateTest}, @@ -981,6 +983,43 @@ fail: return error; } +static int SeasonsIssue187(void) +{ + int error, year; + astro_seasons_t seasons; + char text[TIME_TEXT_BYTES]; + + // This is a regression test for: + // https://github.com/cosinekitty/astronomy/issues/187 + // For years far from the present, the seasons search was sometimes failing. + + for (year = -2000; year <= +9999; ++year) + { + seasons = Astronomy_Seasons(year); + if (seasons.status != ASTRO_SUCCESS) + FAIL("C SeasonsIssue187: Search error %d for year %d.\n", (int)seasons.status, year); + + if (Verbose && ((year > 0) && (year % 1000 == 999))) + { + printf("C SeasonsIssue187: DEBUG"); + Astronomy_FormatTime(seasons.mar_equinox, TIME_FORMAT_DAY, text, sizeof(text)); + printf(" %s", text); + Astronomy_FormatTime(seasons.jun_solstice, TIME_FORMAT_DAY, text, sizeof(text)); + printf(" %s", text); + Astronomy_FormatTime(seasons.sep_equinox, TIME_FORMAT_DAY, text, sizeof(text)); + printf(" %s", text); + Astronomy_FormatTime(seasons.dec_solstice, TIME_FORMAT_DAY, text, sizeof(text)); + printf(" %s\n", text); + } + } + + printf("C SeasonsIssue187: PASS\n"); + error = 0; + +fail: + return error; +} + /*-----------------------------------------------------------------------------------------------------------*/ static int MoonPhase(void) diff --git a/generate/dotnet/csharp_test/csharp_test.cs b/generate/dotnet/csharp_test/csharp_test.cs index 0f84658b..c0b2c0aa 100644 --- a/generate/dotnet/csharp_test/csharp_test.cs +++ b/generate/dotnet/csharp_test/csharp_test.cs @@ -55,6 +55,7 @@ namespace csharp_test new Test("riseset", RiseSetTest), new Test("rotation", RotationTest), new Test("seasons", SeasonsTest), + new Test("seasons187", SeasonsIssue187), new Test("sidereal", SiderealTimeTest), new Test("transit", TransitTest), new Test("astro_check", AstroCheck), @@ -391,6 +392,17 @@ namespace csharp_test } } + static int SeasonsIssue187() + { + // This is a regression test for: + // https://github.com/cosinekitty/astronomy/issues/187 + // For years far from the present, the seasons search was sometimes failing. + for (int year = 1; year <= 9999; ++year) + Astronomy.Seasons(year); + + return 0; + } + static int MoonPhaseTest() { const string filename = "../../moonphase/moonphases.txt"; diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index c316ba06..66a433ad 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -4442,7 +4442,7 @@ astro_search_result_t Astronomy_SearchSunLongitude( double limitDays) { astro_time_t t2 = Astronomy_AddDays(startTime, limitDays); - return Astronomy_Search(sun_offset, &targetLon, startTime, t2, 1.0); + return Astronomy_Search(sun_offset, &targetLon, startTime, t2, 0.01); } /** @cond DOXYGEN_SKIP */ @@ -4693,7 +4693,7 @@ static int QuadInterp( static astro_status_t FindSeasonChange(double targetLon, int year, int month, int day, astro_time_t *time) { astro_time_t startTime = Astronomy_MakeTime(year, month, day, 0, 0, 0.0); - astro_search_result_t result = Astronomy_SearchSunLongitude(targetLon, startTime, 4.0); + astro_search_result_t result = Astronomy_SearchSunLongitude(targetLon, startTime, 20.0); *time = result.time; return result.status; } @@ -4742,16 +4742,26 @@ astro_seasons_t Astronomy_Seasons(int year) seasons.status = ASTRO_SUCCESS; - status = FindSeasonChange( 0, year, 3, 19, &seasons.mar_equinox); + /* + https://github.com/cosinekitty/astronomy/issues/187 + Solstices and equinoxes drift over long spans of time, + due to precession of the Earth's axis. + Therefore, we have to search a wider range of time than + one might expect. It turns out this has very little + effect on efficiency, thanks to the quick convergence + of quadratic interpolation inside Astronomy_Search(). + */ + + status = FindSeasonChange( 0, year, 3, 10, &seasons.mar_equinox); if (status != ASTRO_SUCCESS) seasons.status = status; - status = FindSeasonChange( 90, year, 6, 19, &seasons.jun_solstice); + status = FindSeasonChange( 90, year, 6, 10, &seasons.jun_solstice); if (status != ASTRO_SUCCESS) seasons.status = status; - status = FindSeasonChange(180, year, 9, 21, &seasons.sep_equinox); + status = FindSeasonChange(180, year, 9, 10, &seasons.sep_equinox); if (status != ASTRO_SUCCESS) seasons.status = status; - status = FindSeasonChange(270, year, 12, 20, &seasons.dec_solstice); + status = FindSeasonChange(270, year, 12, 10, &seasons.dec_solstice); if (status != ASTRO_SUCCESS) seasons.status = status; return seasons; diff --git a/generate/template/astronomy.cs b/generate/template/astronomy.cs index de688a94..3916b105 100644 --- a/generate/template/astronomy.cs +++ b/generate/template/astronomy.cs @@ -4702,19 +4702,27 @@ $ASTRO_IAU_DATA() /// public static SeasonsInfo Seasons(int year) { + // https://github.com/cosinekitty/astronomy/issues/187 + // Solstices and equinoxes drift over long spans of time, + // due to precession of the Earth's axis. + // Therefore, we have to search a wider range of time than + // one might expect. It turns out this has very little + // effect on efficiency, thanks to the quick convergence + // of quadratic interpolation inside the `Search` function. + return new SeasonsInfo( - FindSeasonChange( 0, year, 3, 19), - FindSeasonChange( 90, year, 6, 19), - FindSeasonChange(180, year, 9, 21), - FindSeasonChange(270, year, 12, 20) + FindSeasonChange( 0, year, 3, 10), + FindSeasonChange( 90, year, 6, 10), + FindSeasonChange(180, year, 9, 10), + FindSeasonChange(270, year, 12, 10) ); } private static AstroTime FindSeasonChange(double targetLon, int year, int month, int day) { var startTime = new AstroTime(year, month, day, 0, 0, 0); - return SearchSunLongitude(targetLon, startTime, 4.0) ?? - throw new InternalError($"Cannot find solution for Sun longitude {targetLon}"); + return SearchSunLongitude(targetLon, startTime, 20.0) ?? + throw new InternalError($"Cannot find solution for Sun longitude {targetLon} in year {year}"); } /// @@ -4755,7 +4763,7 @@ $ASTRO_IAU_DATA() { var sun_offset = new SearchContext_SunOffset(targetLon); AstroTime t2 = startTime.AddDays(limitDays); - return Search(sun_offset, startTime, t2, 1.0); + return Search(sun_offset, startTime, t2, 0.01); } /// diff --git a/generate/template/astronomy.kt b/generate/template/astronomy.kt index 7a53a89c..c0c031da 100644 --- a/generate/template/astronomy.kt +++ b/generate/template/astronomy.kt @@ -4233,8 +4233,8 @@ object Astronomy { // (t1,f1), (tmid,fmid), (t2,f2). val tm = tmid.ut val tspan = t2.ut - tmid.ut - val q = (f2 + f1) - fmid - val r = (f2 - f1) / 2.0 + val q = (f2 + f1)/2.0 - fmid + val r = (f2 - f1)/2.0 val s = fmid var foundInterpolation = false var x = Double.NaN @@ -4421,7 +4421,7 @@ object Astronomy { } val context = Context(targetLon) val time2 = startTime.addDays(limitDays) - return search(context, startTime, time2, 1.0) + return search(context, startTime, time2, 0.01) } /** @@ -4460,16 +4460,16 @@ object Astronomy { */ fun seasons(year: Int) = SeasonsInfo( - findSeasonChange( 0.0, year, 3, 19), - findSeasonChange( 90.0, year, 6, 19), - findSeasonChange(180.0, year, 9, 21), - findSeasonChange(270.0, year, 12, 20) + findSeasonChange( 0.0, year, 3, 10), + findSeasonChange( 90.0, year, 6, 10), + findSeasonChange(180.0, year, 9, 10), + findSeasonChange(270.0, year, 12, 10) ) private fun findSeasonChange(targetLon: Double, year: Int, month: Int, day: Int): AstroTime { var startTime = AstroTime(year, month, day, 0, 0, 0.0) - return searchSunLongitude(targetLon, startTime, 4.0) ?: - throw InternalError("Cannot find solution for Sun longitude $targetLon") + return searchSunLongitude(targetLon, startTime, 20.0) ?: + throw InternalError("Cannot find solution for Sun longitude $targetLon for year $year") } /** diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 40713d65..e9bacc8d 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -3600,7 +3600,7 @@ def SearchSunLongitude(targetLon, startTime, limitDays): Time or `None` """ t2 = startTime.AddDays(limitDays) - return Search(_sun_offset, targetLon, startTime, t2, 1.0) + return Search(_sun_offset, targetLon, startTime, t2, 0.01) def MoonPhase(time): """Returns the Moon's phase as an angle from 0 to 360 degrees. @@ -4402,7 +4402,7 @@ class SeasonInfo: def _FindSeasonChange(targetLon, year, month, day): startTime = Time.Make(year, month, day, 0, 0, 0) - time = SearchSunLongitude(targetLon, startTime, 4.0) + time = SearchSunLongitude(targetLon, startTime, 20.0) if time is None: # We should always be able to find a season change. raise InternalError() @@ -4444,10 +4444,17 @@ def Seasons(year): ------- SeasonInfo """ - mar_equinox = _FindSeasonChange(0, year, 3, 19) - jun_solstice = _FindSeasonChange(90, year, 6, 19) - sep_equinox = _FindSeasonChange(180, year, 9, 21) - dec_solstice = _FindSeasonChange(270, year, 12, 20) + # https://github.com/cosinekitty/astronomy/issues/187 + # Solstices and equinoxes drift over long spans of time, + # due to precession of the Earth's axis. + # Therefore, we have to search a wider range of time than + # one might expect. It turns out this has very little + # effect on efficiency, thanks to the quick convergence + # of quadratic interpolation inside the `Search` function. + mar_equinox = _FindSeasonChange( 0, year, 3, 10) + jun_solstice = _FindSeasonChange( 90, year, 6, 10) + sep_equinox = _FindSeasonChange(180, year, 9, 10) + dec_solstice = _FindSeasonChange(270, year, 12, 10) return SeasonInfo(mar_equinox, jun_solstice, sep_equinox, dec_solstice) def _MoonDistance(time): diff --git a/generate/template/astronomy.ts b/generate/template/astronomy.ts index 2cf5d1f8..9e9ed725 100644 --- a/generate/template/astronomy.ts +++ b/generate/template/astronomy.ts @@ -3485,7 +3485,7 @@ export function SearchSunLongitude(targetLon: number, dateStart: FlexibleDateTim VerifyNumber(limitDays); let t1 = MakeTime(dateStart); let t2 = t1.AddDays(limitDays); - return Search(sun_offset, t1, t2); + return Search(sun_offset, t1, t2, {dt_tolerance_seconds: 0.01}); } /** @@ -4433,7 +4433,7 @@ export class SeasonInfo { export function Seasons(year: (number | AstroTime)): SeasonInfo { function find(targetLon: number, month: number, day: number): AstroTime { let startDate = new Date(Date.UTC(year, month-1, day)); - let time = SearchSunLongitude(targetLon, startDate, 4); + let time = SearchSunLongitude(targetLon, startDate, 20); if (!time) throw `Cannot find season change near ${startDate.toISOString()}`; return time; @@ -4445,10 +4445,10 @@ export function Seasons(year: (number | AstroTime)): SeasonInfo { if (!Number.isSafeInteger(year)) throw `Cannot calculate seasons because year argument ${year} is neither a Date nor a safe integer.`; - let mar_equinox = find( 0, 3, 19); - let jun_solstice = find( 90, 6, 19); - let sep_equinox = find(180, 9, 21); - let dec_solstice = find(270, 12, 20); + let mar_equinox = find( 0, 3, 10); + let jun_solstice = find( 90, 6, 10); + let sep_equinox = find(180, 9, 10); + let dec_solstice = find(270, 12, 10); return new SeasonInfo(mar_equinox, jun_solstice, sep_equinox, dec_solstice); } diff --git a/generate/test.js b/generate/test.js index 0860103d..4110a0c3 100644 --- a/generate/test.js +++ b/generate/test.js @@ -1068,6 +1068,26 @@ function Seasons() { } } console.log('JS Seasons: PASS'); + return 0; +} + + +function SeasonsIssue187() { + // This is a regression test for: + // https://github.com/cosinekitty/astronomy/issues/187 + // For years far from the present, the seasons search was sometimes failing. + + for (let year = -2000; year <= +9999; ++year) { + try { + Astronomy.Seasons(year); + } catch (e) { + console.error(`JS SeasonsIssue187: FAIL (year = ${year}): `, e); + return 1; + } + } + + console.log('JS SeasonsIssue187: PASS'); + return 0; } @@ -2743,6 +2763,7 @@ const UnitTests = { rise_set: RiseSet, rotation: Rotation, seasons: Seasons, + seasons187: SeasonsIssue187, sidereal: SiderealTimeTest, topostate: TopoStateTest, transit: Transit, diff --git a/generate/test.py b/generate/test.py index a6777f3b..8a1475ef 100755 --- a/generate/test.py +++ b/generate/test.py @@ -217,6 +217,20 @@ def Seasons(filename = 'seasons/seasons.txt'): print('PY Seasons: Event counts: mar={}, jun={}, sep={}, dec={}'.format(mar_count, jun_count, sep_count, dec_count)) return 0 + +def SeasonsIssue187(): + # This is a regression test for: + # https://github.com/cosinekitty/astronomy/issues/187 + # For years far from the present, the seasons search was sometimes failing. + for year in range(1, 9999, 1): + try: + astronomy.Seasons(year) + except astronomy.InternalError: + print('PY SeasonsIssue187: FAIL - internal error for year {}'.format(year)) + return 1 + print('PY SeasonsIssue187: PASS') + return 0 + #----------------------------------------------------------------------------------------------------------- def MoonPhase(filename = 'moonphase/moonphases.txt'): @@ -2510,6 +2524,7 @@ UnitTests = { 'riseset': RiseSet, 'rotation': Rotation, 'seasons': Seasons, + 'seasons187': SeasonsIssue187, 'sidereal': SiderealTime, 'time': AstroTime, 'topostate': TopoState, diff --git a/source/c/astronomy.c b/source/c/astronomy.c index f04aee29..ff9f89cc 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -5681,7 +5681,7 @@ astro_search_result_t Astronomy_SearchSunLongitude( double limitDays) { astro_time_t t2 = Astronomy_AddDays(startTime, limitDays); - return Astronomy_Search(sun_offset, &targetLon, startTime, t2, 1.0); + return Astronomy_Search(sun_offset, &targetLon, startTime, t2, 0.01); } /** @cond DOXYGEN_SKIP */ @@ -5932,7 +5932,7 @@ static int QuadInterp( static astro_status_t FindSeasonChange(double targetLon, int year, int month, int day, astro_time_t *time) { astro_time_t startTime = Astronomy_MakeTime(year, month, day, 0, 0, 0.0); - astro_search_result_t result = Astronomy_SearchSunLongitude(targetLon, startTime, 4.0); + astro_search_result_t result = Astronomy_SearchSunLongitude(targetLon, startTime, 20.0); *time = result.time; return result.status; } @@ -5981,16 +5981,26 @@ astro_seasons_t Astronomy_Seasons(int year) seasons.status = ASTRO_SUCCESS; - status = FindSeasonChange( 0, year, 3, 19, &seasons.mar_equinox); + /* + https://github.com/cosinekitty/astronomy/issues/187 + Solstices and equinoxes drift over long spans of time, + due to precession of the Earth's axis. + Therefore, we have to search a wider range of time than + one might expect. It turns out this has very little + effect on efficiency, thanks to the quick convergence + of quadratic interpolation inside Astronomy_Search(). + */ + + status = FindSeasonChange( 0, year, 3, 10, &seasons.mar_equinox); if (status != ASTRO_SUCCESS) seasons.status = status; - status = FindSeasonChange( 90, year, 6, 19, &seasons.jun_solstice); + status = FindSeasonChange( 90, year, 6, 10, &seasons.jun_solstice); if (status != ASTRO_SUCCESS) seasons.status = status; - status = FindSeasonChange(180, year, 9, 21, &seasons.sep_equinox); + status = FindSeasonChange(180, year, 9, 10, &seasons.sep_equinox); if (status != ASTRO_SUCCESS) seasons.status = status; - status = FindSeasonChange(270, year, 12, 20, &seasons.dec_solstice); + status = FindSeasonChange(270, year, 12, 10, &seasons.dec_solstice); if (status != ASTRO_SUCCESS) seasons.status = status; return seasons; diff --git a/source/csharp/astronomy.cs b/source/csharp/astronomy.cs index e88dc5e3..fb3edd7e 100644 --- a/source/csharp/astronomy.cs +++ b/source/csharp/astronomy.cs @@ -5914,19 +5914,27 @@ namespace CosineKitty /// public static SeasonsInfo Seasons(int year) { + // https://github.com/cosinekitty/astronomy/issues/187 + // Solstices and equinoxes drift over long spans of time, + // due to precession of the Earth's axis. + // Therefore, we have to search a wider range of time than + // one might expect. It turns out this has very little + // effect on efficiency, thanks to the quick convergence + // of quadratic interpolation inside the `Search` function. + return new SeasonsInfo( - FindSeasonChange( 0, year, 3, 19), - FindSeasonChange( 90, year, 6, 19), - FindSeasonChange(180, year, 9, 21), - FindSeasonChange(270, year, 12, 20) + FindSeasonChange( 0, year, 3, 10), + FindSeasonChange( 90, year, 6, 10), + FindSeasonChange(180, year, 9, 10), + FindSeasonChange(270, year, 12, 10) ); } private static AstroTime FindSeasonChange(double targetLon, int year, int month, int day) { var startTime = new AstroTime(year, month, day, 0, 0, 0); - return SearchSunLongitude(targetLon, startTime, 4.0) ?? - throw new InternalError($"Cannot find solution for Sun longitude {targetLon}"); + return SearchSunLongitude(targetLon, startTime, 20.0) ?? + throw new InternalError($"Cannot find solution for Sun longitude {targetLon} in year {year}"); } /// @@ -5967,7 +5975,7 @@ namespace CosineKitty { var sun_offset = new SearchContext_SunOffset(targetLon); AstroTime t2 = startTime.AddDays(limitDays); - return Search(sun_offset, startTime, t2, 1.0); + return Search(sun_offset, startTime, t2, 0.01); } /// diff --git a/source/js/astronomy.browser.js b/source/js/astronomy.browser.js index 82b2da90..55eb9de4 100644 --- a/source/js/astronomy.browser.js +++ b/source/js/astronomy.browser.js @@ -4022,7 +4022,7 @@ function SearchSunLongitude(targetLon, dateStart, limitDays) { VerifyNumber(limitDays); let t1 = MakeTime(dateStart); let t2 = t1.AddDays(limitDays); - return Search(sun_offset, t1, t2); + return Search(sun_offset, t1, t2, { dt_tolerance_seconds: 0.01 }); } exports.SearchSunLongitude = SearchSunLongitude; /** @@ -4921,7 +4921,7 @@ exports.SeasonInfo = SeasonInfo; function Seasons(year) { function find(targetLon, month, day) { let startDate = new Date(Date.UTC(year, month - 1, day)); - let time = SearchSunLongitude(targetLon, startDate, 4); + let time = SearchSunLongitude(targetLon, startDate, 20); if (!time) throw `Cannot find season change near ${startDate.toISOString()}`; return time; @@ -4930,10 +4930,10 @@ function Seasons(year) { year = year.getUTCFullYear(); if (!Number.isSafeInteger(year)) throw `Cannot calculate seasons because year argument ${year} is neither a Date nor a safe integer.`; - let mar_equinox = find(0, 3, 19); - let jun_solstice = find(90, 6, 19); - let sep_equinox = find(180, 9, 21); - let dec_solstice = find(270, 12, 20); + let mar_equinox = find(0, 3, 10); + let jun_solstice = find(90, 6, 10); + let sep_equinox = find(180, 9, 10); + let dec_solstice = find(270, 12, 10); return new SeasonInfo(mar_equinox, jun_solstice, sep_equinox, dec_solstice); } exports.Seasons = Seasons; diff --git a/source/js/astronomy.browser.min.js b/source/js/astronomy.browser.min.js index 2b53cb4b..b60d92d4 100644 --- a/source/js/astronomy.browser.min.js +++ b/source/js/astronomy.browser.min.js @@ -80,10 +80,10 @@ ha(H[a][2],b.tt/365250,!1):Ma(a,b).Length()}function ca(a,b,c){A(c);b=w(b);if(a= 0,0,0,0,0,b);case q.SSB:return a=new qa(b.tt),new K(-a.Sun.r.x,-a.Sun.r.y,-a.Sun.r.z,-a.Sun.v.x,-a.Sun.v.y,-a.Sun.v.z,b);case q.Mercury:case q.Venus:case q.Earth:case q.Mars:case q.Jupiter:case q.Saturn:case q.Uranus:case q.Neptune:return a=Ka(H[a],b.tt),ra(a,b);case q.Pluto:return wb(b,!0);case q.Moon:case q.EMB:var c=Ka(H.Earth,b.tt);a=a==q.Moon?Ja(b):rb(b);return new K(a.x+c.r.x,a.y+c.r.y,a.z+c.r.z,a.vx+c.v.x,a.vy+c.v.y,a.vz+c.v.z,b);default:throw'HelioState: Unsupported body "'+a+'"';}}function Zc(a, b,c,d,f){var h=(f+c)/2-d;c=(f-c)/2;if(0==h){if(0==c)return null;d=-d/c;if(-1>d||1=d)return null;f=Math.sqrt(d);d=(-c+f)/(2*h);f=(-c-f)/(2*h);if(-1<=d&&1>=d){if(-1<=f&&1>=f)return null}else if(-1<=f&&1>=f)d=f;else return null}return{x:d,t:a+d*b,df_dt:(2*h*d+c)/b}}function I(a,b,c,d){var f=y(d&&d.dt_tolerance_seconds||1);f=Math.abs(f/86400);var h=d&&d.init_f1||a(b),l=d&&d.init_f2||a(c),k=NaN,g=0;d=d&&d.iter_limit||20;for(var m=!0;;){if(++g>d)throw"Excessive iteration in Search()"; var n=new P(b.ut+.5*(c.ut-b.ut)),p=n.ut-b.ut;if(Math.abs(p)(p.ut-b.ut)*(p.ut-c.ut)&&0>(x.ut-b.ut)*(x.ut-c.ut))){t=a(p);var z=a(x);if(0>t&&0<=z){h=t;l=z;b=p;c=x;k=v;m=!1;continue}}}}if(0>h&&0<=k)c=n,l=k;else if(0>k&&0<=l)b=n,h=k;else return null}}function sa(a){for(;-180>=a;)a+=360;for(;180a;)a+=360;for(;360<=a;)a-=360;return a}function gc(a,b,c){y(a);y(c);b=w(b);c=b.AddDays(c);return I(function(d){d=Wb(d);return sa(d.elon-a)},b,c)}function zb(a,b,c){if(a===q.Earth||b===q.Earth)throw"The Earth does not have a longitude as seen from itself.";c=w(c);a=ca(a,c,!1);a=Ha(a);b=ca(b,c,!1);b=Ha(b);return ta(a.elon-b.elon)}function ua(a,b){if(a==q.Earth)throw"The Earth does not have an angle as seen from itself.";var c=w(b);b=ca(q.Sun,c,!0);a=ca(a,c,!0);return O(b, -a)}function ka(a,b){if(a===q.Sun)throw"Cannot calculate heliocentric longitude of the Sun.";a=Ma(a,b);return Ha(a).elon}function hb(a,b){if(a===q.Earth)throw"The illumination of the Earth is not defined.";var c=w(b),d=U(H.Earth,c);if(a===q.Sun){var f=new D(-d.x,-d.y,-d.z,c);b=new D(0,0,0,c);d=0}else a===q.Moon?(f=V(c),b=new D(d.x+f.x,d.y+f.y,d.z+f.z,c)):(b=Ma(a,b),f=new D(b.x-d.x,b.y-d.y,b.z-d.z,c)),d=O(f,b);var h=f.Length(),l=b.Length();if(a===q.Sun)var k=$c+5*Math.log10(h);else if(a===q.Moon)a= -d*e.DEG2RAD,k=a*a,a=-12.717+1.49*Math.abs(a)+.0431*k*k,k=a+=5*Math.log10(h/(385000.6/e.KM_PER_AU)*l);else if(a===q.Saturn){var g=d;a=Ha(f);k=28.06*e.DEG2RAD;var m=e.DEG2RAD*a.elat;a=Math.asin(Math.sin(m)*Math.cos(k)-Math.cos(m)*Math.sin(k)*Math.sin(e.DEG2RAD*a.elon-e.DEG2RAD*(169.51+3.82E-5*c.tt)));k=Math.sin(Math.abs(a));g=-9+.044*g+k*(-2.6+1.2*k)+5*Math.log10(l*h);a*=e.RAD2DEG;k=g;g=a}else{var n=m=k=0;switch(a){case q.Mercury:a=-.6;k=4.98;m=-4.88;n=3.02;break;case q.Venus:163.6>d?(a=-4.47,k=1.03, -m=.57,n=.13):(a=.98,k=-1.02);break;case q.Mars:a=-1.52;k=1.6;break;case q.Jupiter:a=-9.4;k=.5;break;case q.Uranus:a=-7.19;k=.25;break;case q.Neptune:a=-6.87;break;case q.Pluto:a=-1;k=4;break;default:throw"VisualMagnitude: unsupported body "+a;}var p=d/100;k=a+p*(k+p*(m+p*n))+5*Math.log10(l*h)}return new hc(c,k,d,l,h,f,b,g)}function Na(a){if(a===q.Earth)throw"The Earth does not have a synodic period as seen from itself.";if(a===q.Moon)return 29.530588;var b=da[a];if(!b)throw"Not a valid planet name: "+ +360;return a}function ta(a){for(;0>a;)a+=360;for(;360<=a;)a-=360;return a}function gc(a,b,c){y(a);y(c);b=w(b);c=b.AddDays(c);return I(function(d){d=Wb(d);return sa(d.elon-a)},b,c,{dt_tolerance_seconds:.01})}function zb(a,b,c){if(a===q.Earth||b===q.Earth)throw"The Earth does not have a longitude as seen from itself.";c=w(c);a=ca(a,c,!1);a=Ha(a);b=ca(b,c,!1);b=Ha(b);return ta(a.elon-b.elon)}function ua(a,b){if(a==q.Earth)throw"The Earth does not have an angle as seen from itself.";var c=w(b);b=ca(q.Sun, +c,!0);a=ca(a,c,!0);return O(b,a)}function ka(a,b){if(a===q.Sun)throw"Cannot calculate heliocentric longitude of the Sun.";a=Ma(a,b);return Ha(a).elon}function hb(a,b){if(a===q.Earth)throw"The illumination of the Earth is not defined.";var c=w(b),d=U(H.Earth,c);if(a===q.Sun){var f=new D(-d.x,-d.y,-d.z,c);b=new D(0,0,0,c);d=0}else a===q.Moon?(f=V(c),b=new D(d.x+f.x,d.y+f.y,d.z+f.z,c)):(b=Ma(a,b),f=new D(b.x-d.x,b.y-d.y,b.z-d.z,c)),d=O(f,b);var h=f.Length(),l=b.Length();if(a===q.Sun)var k=$c+5*Math.log10(h); +else if(a===q.Moon)a=d*e.DEG2RAD,k=a*a,a=-12.717+1.49*Math.abs(a)+.0431*k*k,k=a+=5*Math.log10(h/(385000.6/e.KM_PER_AU)*l);else if(a===q.Saturn){var g=d;a=Ha(f);k=28.06*e.DEG2RAD;var m=e.DEG2RAD*a.elat;a=Math.asin(Math.sin(m)*Math.cos(k)-Math.cos(m)*Math.sin(k)*Math.sin(e.DEG2RAD*a.elon-e.DEG2RAD*(169.51+3.82E-5*c.tt)));k=Math.sin(Math.abs(a));g=-9+.044*g+k*(-2.6+1.2*k)+5*Math.log10(l*h);a*=e.RAD2DEG;k=g;g=a}else{var n=m=k=0;switch(a){case q.Mercury:a=-.6;k=4.98;m=-4.88;n=3.02;break;case q.Venus:163.6> +d?(a=-4.47,k=1.03,m=.57,n=.13):(a=.98,k=-1.02);break;case q.Mars:a=-1.52;k=1.6;break;case q.Jupiter:a=-9.4;k=.5;break;case q.Uranus:a=-7.19;k=.25;break;case q.Neptune:a=-6.87;break;case q.Pluto:a=-1;k=4;break;default:throw"VisualMagnitude: unsupported body "+a;}var p=d/100;k=a+p*(k+p*(m+p*n))+5*Math.log10(l*h)}return new hc(c,k,d,l,h,f,b,g)}function Na(a){if(a===q.Earth)throw"The Earth does not have a synodic period as seen from itself.";if(a===q.Moon)return 29.530588;var b=da[a];if(!b)throw"Not a valid planet name: "+ a;a=da.Earth.OrbitalPeriod;return Math.abs(a/(a/b.OrbitalPeriod-1))}function va(a,b,c){function d(m){var n=ka(a,m);m=ka(q.Earth,m);return sa(h*(m-n)-b)}y(b);var f=da[a];if(!f)throw"Cannot search relative longitude because body is not a planet: "+a;if(a===q.Earth)throw"Cannot search relative longitude for the Earth (it is always 0)";var h=f.OrbitalPeriod>da.Earth.OrbitalPeriod?1:-1;f=Na(a);c=w(c);var l=d(c);0k;++k){var g=-l/360*f;c=c.AddDays(g);if(1>86400*Math.abs(g))return c; g=l;l=d(c);30>Math.abs(g)&&g!==l&&(g/=g-l,.5g&&(f*=g))}throw"Relative longitude search failed to converge for "+a+" near "+c.toString()+" (error_angle = "+l+").";}function Ab(a){return zb(q.Moon,q.Sun,a)}function Oa(a,b,c){function d(l){l=Ab(l);return sa(l-a)}y(a);y(c);b=w(b);var f=d(b);0c)return null;c=Math.min(c,h+1.5);f=b.AddDays(f);b=b.AddDays(c);return I(d,f,b)}function ic(a){var b=Ab(a);b=(Math.floor(b/90)+1)%4;a=Oa(90*b,a,10);if(!a)throw"Cannot find moon quarter"; return new jc(b,a)}function kc(a,b,c,d,f,h){pa(b);y(f);if(0>=f)throw"Invalid value for limitDays: "+f;if(a===q.Earth)throw"Cannot find altitude event for the Earth.";if(1===c){c=12;var l=0}else if(-1===c)c=0,l=12;else throw"Invalid direction parameter "+c+" -- must be +1 or -1";d=w(d);var k=h(d);var g;if(0=k&&0= @@ -201,8 +201,8 @@ x=(n-d+m*v-p*x)/(1-m*x-p*v);d+=x}while(1E-12<=Math.abs(x));x=Math.cos(d);v=Math. 0,0,0,b);if(a===q.Pluto)return wb(b,!1);var c=new qa(b.tt);switch(a){case q.Sun:return ra(c.Sun,b);case q.Jupiter:return ra(c.Jupiter,b);case q.Saturn:return ra(c.Saturn,b);case q.Uranus:return ra(c.Uranus,b);case q.Neptune:return ra(c.Neptune,b);case q.Moon:case q.EMB:var d=Ka(H[q.Earth],b.tt);a=a===q.Moon?Ja(b):rb(b);return new K(a.x+c.Sun.r.x+d.r.x,a.y+c.Sun.r.y+d.r.y,a.z+c.Sun.r.z+d.r.z,a.vx+c.Sun.v.x+d.v.x,a.vy+c.Sun.v.y+d.v.y,a.vz+c.Sun.v.z+d.v.z,b)}if(a in H)return a=Ka(H[a],b.tt),new K(c.Sun.r.x+ a.r.x,c.Sun.r.y+a.r.y,c.Sun.r.z+a.r.z,c.Sun.v.x+a.v.x,c.Sun.v.y+a.v.y,c.Sun.v.z+a.v.z,b);throw'BaryState: Unsupported body "'+a+'"';};e.HelioState=yb;e.Search=I;e.SearchSunLongitude=gc;e.PairLongitude=zb;e.AngleFromSun=ua;e.EclipticLongitude=ka;var hc=function(a,b,c,d,f,h,l,k){this.time=a;this.mag=b;this.phase_angle=c;this.helio_dist=d;this.geo_dist=f;this.gc=h;this.hc=l;this.ring_tilt=k;this.phase_fraction=(1+Math.cos(e.DEG2RAD*c))/2};e.IlluminationInfo=hc;e.Illumination=hb;e.SearchRelativeLongitude= va;e.MoonPhase=Ab;e.SearchMoonPhase=Oa;var jc=function(a,b){this.quarter=a;this.time=b};e.MoonQuarter=jc;e.SearchMoonQuarter=ic;e.NextMoonQuarter=function(a){a=new Date(a.time.date.getTime()+5184E5);return ic(a)};e.SearchRiseSet=function(a,b,c,d,f){a:switch(a){case q.Sun:var h=gd;break a;case q.Moon:h=hd;break a;default:h=0}return kc(a,b,c,d,f,function(l){var k=Ga(a,l,b,!0,!0);l=Ea(l,b,k.ra,k.dec).altitude+h/k.dist*e.RAD2DEG+id;return c*l})};e.SearchAltitude=function(a,b,c,d,f,h){if(!Number.isFinite(h)|| --90>h||90h||90=++f;){var h=ka(a,b),l=ka(q.Earth,b),k=sa(h-l),g=h=l=void 0;k>=-d.s1&&k<+d.s1?(g=0,l=+d.s1,h=+d.s2):k>=+d.s2||k<-d.s2?(g=0,l=-d.s2,h=-d.s1):0<=k?(g=-Na(a)/4,l=+d.s1,h=+d.s2):(g=-Na(a)/4,l=-d.s2,h=-d.s1);k=b.AddDays(g);l=va(a,l,k);h=va(a,h,l);k=c(l);if(0<=k)throw"SearchMaxElongation: internal error: m1 = "+k;g=c(h);if(0>=g)throw"SearchMaxElongation: internal error: m2 = "+ g;k=I(c,l,h,{init_f1:k,init_f2:g,dt_tolerance_seconds:10});if(!k)throw"SearchMaxElongation: failed search iter "+f+" (t1="+l.toString()+", t2="+h.toString()+")";if(k.tt>=b.tt)return mc(a,k);b=h.AddDays(1)}throw"SearchMaxElongation: failed to find event after 2 tries.";};e.SearchPeakMagnitude=function(a,b){function c(g){var m=g.AddDays(-.005);g=g.AddDays(.005);m=hb(a,m).mag;return(hb(a,g).mag-m)/.01}if(a!==q.Venus)throw"SearchPeakMagnitude currently works for Venus only.";b=w(b);for(var d=0;2>=++d;){var f= ka(a,b),h=ka(q.Earth,b),l=sa(f-h),k=f=h=void 0;-10<=l&&10>l?(k=0,h=10,f=30):30<=l||-30>l?(k=0,h=-30,f=-10):0<=l?(k=-Na(a)/4,h=10,f=30):(k=-Na(a)/4,h=-30,f=-10);l=b.AddDays(k);h=va(a,h,l);f=va(a,f,h);l=c(h);if(0<=l)throw"SearchPeakMagnitude: internal error: m1 = "+l;k=c(f);if(0>=k)throw"SearchPeakMagnitude: internal error: m2 = "+k;l=I(c,h,f,{init_f1:l,init_f2:k,dt_tolerance_seconds:10});if(!l)throw"SearchPeakMagnitude: failed search iter "+d+" (t1="+h.toString()+", t2="+f.toString()+")";if(l.tt>= diff --git a/source/js/astronomy.js b/source/js/astronomy.js index 008a1db9..264f1041 100644 --- a/source/js/astronomy.js +++ b/source/js/astronomy.js @@ -4021,7 +4021,7 @@ function SearchSunLongitude(targetLon, dateStart, limitDays) { VerifyNumber(limitDays); let t1 = MakeTime(dateStart); let t2 = t1.AddDays(limitDays); - return Search(sun_offset, t1, t2); + return Search(sun_offset, t1, t2, { dt_tolerance_seconds: 0.01 }); } exports.SearchSunLongitude = SearchSunLongitude; /** @@ -4920,7 +4920,7 @@ exports.SeasonInfo = SeasonInfo; function Seasons(year) { function find(targetLon, month, day) { let startDate = new Date(Date.UTC(year, month - 1, day)); - let time = SearchSunLongitude(targetLon, startDate, 4); + let time = SearchSunLongitude(targetLon, startDate, 20); if (!time) throw `Cannot find season change near ${startDate.toISOString()}`; return time; @@ -4929,10 +4929,10 @@ function Seasons(year) { year = year.getUTCFullYear(); if (!Number.isSafeInteger(year)) throw `Cannot calculate seasons because year argument ${year} is neither a Date nor a safe integer.`; - let mar_equinox = find(0, 3, 19); - let jun_solstice = find(90, 6, 19); - let sep_equinox = find(180, 9, 21); - let dec_solstice = find(270, 12, 20); + let mar_equinox = find(0, 3, 10); + let jun_solstice = find(90, 6, 10); + let sep_equinox = find(180, 9, 10); + let dec_solstice = find(270, 12, 10); return new SeasonInfo(mar_equinox, jun_solstice, sep_equinox, dec_solstice); } exports.Seasons = Seasons; diff --git a/source/js/astronomy.min.js b/source/js/astronomy.min.js index 0c0a2cd1..ad5a23d8 100644 --- a/source/js/astronomy.min.js +++ b/source/js/astronomy.min.js @@ -204,8 +204,9 @@ b.tt);a=a==Body.Moon?GeoMoonState(b):GeoEmbState(b);return new StateVector(a.x+c function QuadInterp(a,b,c,d,e){var f=(e+c)/2-d;c=(e-c)/2;if(0==f){if(0==c)return null;d=-d/c;if(-1>d||1=d)return null;e=Math.sqrt(d);d=(-c+e)/(2*f);e=(-c-e)/(2*f);if(-1<=d&&1>=d){if(-1<=e&&1>=e)return null}else if(-1<=e&&1>=e)d=e;else return null}return{x:d,t:a+d*b,df_dt:(2*f*d+c)/b}} function Search(a,b,c,d){var e=VerifyNumber(d&&d.dt_tolerance_seconds||1);e=Math.abs(e/SECONDS_PER_DAY);var f=d&&d.init_f1||a(b),g=d&&d.init_f2||a(c),k=NaN,h=0;d=d&&d.iter_limit||20;for(var l=!0;;){if(++h>d)throw"Excessive iteration in Search()";var m=InterpolateTime(b,c,.5),n=m.ut-b.ut;if(Math.abs(n)(n.ut-b.ut)*(n.ut-c.ut)&&0>(r.ut-b.ut)*(r.ut-c.ut))){p=a(n);var t=a(r);if(0>p&&0<=t){f=p;g=t;b=n;c=r;k=q;l=!1;continue}}}}if(0>f&&0<=k)c=m,g=k;else if(0>k&&0<=g)b=m,f=k;else return null}}exports.Search=Search;function LongitudeOffset(a){for(;-180>=a;)a+=360;for(;180a;)a+=360;for(;360<=a;)a-=360;return a} -function SearchSunLongitude(a,b,c){VerifyNumber(a);VerifyNumber(c);b=MakeTime(b);c=b.AddDays(c);return Search(function(d){d=SunPosition(d);return LongitudeOffset(d.elon-a)},b,c)}exports.SearchSunLongitude=SearchSunLongitude;function PairLongitude(a,b,c){if(a===Body.Earth||b===Body.Earth)throw"The Earth does not have a longitude as seen from itself.";c=MakeTime(c);a=GeoVector(a,c,!1);a=Ecliptic(a);b=GeoVector(b,c,!1);b=Ecliptic(b);return NormalizeLongitude(a.elon-b.elon)}exports.PairLongitude=PairLongitude; -function AngleFromSun(a,b){if(a==Body.Earth)throw"The Earth does not have an angle as seen from itself.";var c=MakeTime(b);b=GeoVector(Body.Sun,c,!0);a=GeoVector(a,c,!0);return AngleBetween(b,a)}exports.AngleFromSun=AngleFromSun;function EclipticLongitude(a,b){if(a===Body.Sun)throw"Cannot calculate heliocentric longitude of the Sun.";a=HelioVector(a,b);return Ecliptic(a).elon}exports.EclipticLongitude=EclipticLongitude; +function SearchSunLongitude(a,b,c){VerifyNumber(a);VerifyNumber(c);b=MakeTime(b);c=b.AddDays(c);return Search(function(d){d=SunPosition(d);return LongitudeOffset(d.elon-a)},b,c,{dt_tolerance_seconds:.01})}exports.SearchSunLongitude=SearchSunLongitude; +function PairLongitude(a,b,c){if(a===Body.Earth||b===Body.Earth)throw"The Earth does not have a longitude as seen from itself.";c=MakeTime(c);a=GeoVector(a,c,!1);a=Ecliptic(a);b=GeoVector(b,c,!1);b=Ecliptic(b);return NormalizeLongitude(a.elon-b.elon)}exports.PairLongitude=PairLongitude;function AngleFromSun(a,b){if(a==Body.Earth)throw"The Earth does not have an angle as seen from itself.";var c=MakeTime(b);b=GeoVector(Body.Sun,c,!0);a=GeoVector(a,c,!0);return AngleBetween(b,a)} +exports.AngleFromSun=AngleFromSun;function EclipticLongitude(a,b){if(a===Body.Sun)throw"Cannot calculate heliocentric longitude of the Sun.";a=HelioVector(a,b);return Ecliptic(a).elon}exports.EclipticLongitude=EclipticLongitude; function VisualMagnitude(a,b,c,d){var e=0,f=0,g=0;switch(a){case Body.Mercury:a=-.6;e=4.98;f=-4.88;g=3.02;break;case Body.Venus:163.6>b?(a=-4.47,e=1.03,f=.57,g=.13):(a=.98,e=-1.02);break;case Body.Mars:a=-1.52;e=1.6;break;case Body.Jupiter:a=-9.4;e=.5;break;case Body.Uranus:a=-7.19;e=.25;break;case Body.Neptune:a=-6.87;break;case Body.Pluto:a=-1;e=4;break;default:throw"VisualMagnitude: unsupported body "+a;}b/=100;return a+b*(e+b*(f+b*g))+5*Math.log10(c*d)} function SaturnMagnitude(a,b,c,d,e){d=Ecliptic(d);var f=28.06*exports.DEG2RAD,g=exports.DEG2RAD*d.elat;e=Math.asin(Math.sin(g)*Math.cos(f)-Math.cos(g)*Math.sin(f)*Math.sin(exports.DEG2RAD*d.elon-exports.DEG2RAD*(169.51+3.82E-5*e.tt)));d=Math.sin(Math.abs(e));a=-9+.044*a+d*(-2.6+1.2*d)+5*Math.log10(b*c);return{mag:a,ring_tilt:exports.RAD2DEG*e}}function MoonMagnitude(a,b,c){a*=exports.DEG2RAD;var d=a*a;a=-12.717+1.49*Math.abs(a)+.0431*d*d;return a+=5*Math.log10(c/(385000.6/exports.KM_PER_AU)*b)} var IlluminationInfo=function(a,b,c,d,e,f,g,k){this.time=a;this.mag=b;this.phase_angle=c;this.helio_dist=d;this.geo_dist=e;this.gc=f;this.hc=g;this.ring_tilt=k;this.phase_fraction=(1+Math.cos(exports.DEG2RAD*c))/2};exports.IlluminationInfo=IlluminationInfo; @@ -221,7 +222,7 @@ function InternalSearchAltitude(a,b,c,d,e,f){VerifyObserver(b);VerifyNumber(e);i k=SearchHourAngle(a,b,c,m.time);m=SearchHourAngle(a,b,g,k.time);if(k.time.ut>=d.ut+e)return null;l=k.time;k=f(k.time);h=f(m.time)}}var HourAngleEvent=function(a,b){this.time=a;this.hor=b};exports.HourAngleEvent=HourAngleEvent; function SearchHourAngle(a,b,c,d){VerifyObserver(b);d=MakeTime(d);var e=0;if(a===Body.Earth)throw"Cannot search for hour angle of the Earth.";VerifyNumber(c);if(0>c||24<=c)throw"Invalid hour angle "+c;for(;;){++e;var f=sidereal_time(d),g=Equator(a,d,b,!0,!0);f=(c+g.ra-b.longitude/15-f)%24;1===e?0>f&&(f+=24):-12>f?f+=24:123600*Math.abs(f))return a=Horizon(d,b,g.ra,g.dec,"normal"),new HourAngleEvent(d,a);d=d.AddDays(f/24*SOLAR_DAYS_PER_SIDEREAL_DAY)}}exports.SearchHourAngle=SearchHourAngle; var SeasonInfo=function(a,b,c,d){this.mar_equinox=a;this.jun_solstice=b;this.sep_equinox=c;this.dec_solstice=d};exports.SeasonInfo=SeasonInfo; -function Seasons(a){function b(g,k,h){k=new Date(Date.UTC(a,k-1,h));g=SearchSunLongitude(g,k,4);if(!g)throw"Cannot find season change near "+k.toISOString();return g}a instanceof Date&&Number.isFinite(a.getTime())&&(a=a.getUTCFullYear());if(!Number.isSafeInteger(a))throw"Cannot calculate seasons because year argument "+a+" is neither a Date nor a safe integer.";var c=b(0,3,19),d=b(90,6,19),e=b(180,9,21),f=b(270,12,20);return new SeasonInfo(c,d,e,f)}exports.Seasons=Seasons; +function Seasons(a){function b(g,k,h){k=new Date(Date.UTC(a,k-1,h));g=SearchSunLongitude(g,k,20);if(!g)throw"Cannot find season change near "+k.toISOString();return g}a instanceof Date&&Number.isFinite(a.getTime())&&(a=a.getUTCFullYear());if(!Number.isSafeInteger(a))throw"Cannot calculate seasons because year argument "+a+" is neither a Date nor a safe integer.";var c=b(0,3,10),d=b(90,6,10),e=b(180,9,10),f=b(270,12,10);return new SeasonInfo(c,d,e,f)}exports.Seasons=Seasons; var ElongationEvent=function(a,b,c,d){this.time=a;this.visibility=b;this.elongation=c;this.ecliptic_separation=d};exports.ElongationEvent=ElongationEvent;function Elongation(a,b){b=MakeTime(b);var c=PairLongitude(a,Body.Sun,b);if(180=++e;){var f=EclipticLongitude(a,b),g=EclipticLongitude(Body.Earth,b),k=LongitudeOffset(f-g),h=f=g=void 0;k>=-d.s1&&k<+d.s1?(h=0,g=+d.s1,f=+d.s2):k>=+d.s2||k<-d.s2?(h=0,g=-d.s2,f=-d.s1):0<=k?(h=-SynodicPeriod(a)/ 4,g=+d.s1,f=+d.s2):(h=-SynodicPeriod(a)/4,g=-d.s2,f=-d.s1);k=b.AddDays(h);g=SearchRelativeLongitude(a,g,k);f=SearchRelativeLongitude(a,f,g);k=c(g);if(0<=k)throw"SearchMaxElongation: internal error: m1 = "+k;h=c(f);if(0>=h)throw"SearchMaxElongation: internal error: m2 = "+h;k=Search(c,g,f,{init_f1:k,init_f2:h,dt_tolerance_seconds:10});if(!k)throw"SearchMaxElongation: failed search iter "+e+" (t1="+g.toString()+", t2="+f.toString()+")";if(k.tt>=b.tt)return Elongation(a,k);b=f.AddDays(1)}throw"SearchMaxElongation: failed to find event after 2 tries."; diff --git a/source/js/astronomy.ts b/source/js/astronomy.ts index 68df29ef..b8e20305 100644 --- a/source/js/astronomy.ts +++ b/source/js/astronomy.ts @@ -4491,7 +4491,7 @@ export function SearchSunLongitude(targetLon: number, dateStart: FlexibleDateTim VerifyNumber(limitDays); let t1 = MakeTime(dateStart); let t2 = t1.AddDays(limitDays); - return Search(sun_offset, t1, t2); + return Search(sun_offset, t1, t2, {dt_tolerance_seconds: 0.01}); } /** @@ -5439,7 +5439,7 @@ export class SeasonInfo { export function Seasons(year: (number | AstroTime)): SeasonInfo { function find(targetLon: number, month: number, day: number): AstroTime { let startDate = new Date(Date.UTC(year, month-1, day)); - let time = SearchSunLongitude(targetLon, startDate, 4); + let time = SearchSunLongitude(targetLon, startDate, 20); if (!time) throw `Cannot find season change near ${startDate.toISOString()}`; return time; @@ -5451,10 +5451,10 @@ export function Seasons(year: (number | AstroTime)): SeasonInfo { if (!Number.isSafeInteger(year)) throw `Cannot calculate seasons because year argument ${year} is neither a Date nor a safe integer.`; - let mar_equinox = find( 0, 3, 19); - let jun_solstice = find( 90, 6, 19); - let sep_equinox = find(180, 9, 21); - let dec_solstice = find(270, 12, 20); + let mar_equinox = find( 0, 3, 10); + let jun_solstice = find( 90, 6, 10); + let sep_equinox = find(180, 9, 10); + let dec_solstice = find(270, 12, 10); return new SeasonInfo(mar_equinox, jun_solstice, sep_equinox, dec_solstice); } diff --git a/source/js/esm/astronomy.js b/source/js/esm/astronomy.js index 24cf5219..6377a71e 100644 --- a/source/js/esm/astronomy.js +++ b/source/js/esm/astronomy.js @@ -3978,7 +3978,7 @@ export function SearchSunLongitude(targetLon, dateStart, limitDays) { VerifyNumber(limitDays); let t1 = MakeTime(dateStart); let t2 = t1.AddDays(limitDays); - return Search(sun_offset, t1, t2); + return Search(sun_offset, t1, t2, { dt_tolerance_seconds: 0.01 }); } /** * @brief Returns one body's ecliptic longitude with respect to another, as seen from the Earth. @@ -4860,7 +4860,7 @@ export class SeasonInfo { export function Seasons(year) { function find(targetLon, month, day) { let startDate = new Date(Date.UTC(year, month - 1, day)); - let time = SearchSunLongitude(targetLon, startDate, 4); + let time = SearchSunLongitude(targetLon, startDate, 20); if (!time) throw `Cannot find season change near ${startDate.toISOString()}`; return time; @@ -4869,10 +4869,10 @@ export function Seasons(year) { year = year.getUTCFullYear(); if (!Number.isSafeInteger(year)) throw `Cannot calculate seasons because year argument ${year} is neither a Date nor a safe integer.`; - let mar_equinox = find(0, 3, 19); - let jun_solstice = find(90, 6, 19); - let sep_equinox = find(180, 9, 21); - let dec_solstice = find(270, 12, 20); + let mar_equinox = find(0, 3, 10); + let jun_solstice = find(90, 6, 10); + let sep_equinox = find(180, 9, 10); + let dec_solstice = find(270, 12, 10); return new SeasonInfo(mar_equinox, jun_solstice, sep_equinox, dec_solstice); } /** diff --git a/source/kotlin/src/main/kotlin/io/github/cosinekitty/astronomy/astronomy.kt b/source/kotlin/src/main/kotlin/io/github/cosinekitty/astronomy/astronomy.kt index ab93de3d..43047611 100644 --- a/source/kotlin/src/main/kotlin/io/github/cosinekitty/astronomy/astronomy.kt +++ b/source/kotlin/src/main/kotlin/io/github/cosinekitty/astronomy/astronomy.kt @@ -4233,8 +4233,8 @@ object Astronomy { // (t1,f1), (tmid,fmid), (t2,f2). val tm = tmid.ut val tspan = t2.ut - tmid.ut - val q = (f2 + f1) - fmid - val r = (f2 - f1) / 2.0 + val q = (f2 + f1)/2.0 - fmid + val r = (f2 - f1)/2.0 val s = fmid var foundInterpolation = false var x = Double.NaN @@ -4421,7 +4421,7 @@ object Astronomy { } val context = Context(targetLon) val time2 = startTime.addDays(limitDays) - return search(context, startTime, time2, 1.0) + return search(context, startTime, time2, 0.01) } /** @@ -4460,16 +4460,16 @@ object Astronomy { */ fun seasons(year: Int) = SeasonsInfo( - findSeasonChange( 0.0, year, 3, 19), - findSeasonChange( 90.0, year, 6, 19), - findSeasonChange(180.0, year, 9, 21), - findSeasonChange(270.0, year, 12, 20) + findSeasonChange( 0.0, year, 3, 10), + findSeasonChange( 90.0, year, 6, 10), + findSeasonChange(180.0, year, 9, 10), + findSeasonChange(270.0, year, 12, 10) ) private fun findSeasonChange(targetLon: Double, year: Int, month: Int, day: Int): AstroTime { var startTime = AstroTime(year, month, day, 0, 0, 0.0) - return searchSunLongitude(targetLon, startTime, 4.0) ?: - throw InternalError("Cannot find solution for Sun longitude $targetLon") + return searchSunLongitude(targetLon, startTime, 20.0) ?: + throw InternalError("Cannot find solution for Sun longitude $targetLon for year $year") } /** diff --git a/source/kotlin/src/test/kotlin/io/github/cosinekitty/astronomy/Tests.kt b/source/kotlin/src/test/kotlin/io/github/cosinekitty/astronomy/Tests.kt index 69cad750..45ccdc70 100644 --- a/source/kotlin/src/test/kotlin/io/github/cosinekitty/astronomy/Tests.kt +++ b/source/kotlin/src/test/kotlin/io/github/cosinekitty/astronomy/Tests.kt @@ -832,6 +832,15 @@ class Tests { assertTrue(decemberCount == 301, "decemberCount = $decemberCount") } + @Test + fun `Seasons issue 187`() { + // This is a regression test for: + // https://github.com/cosinekitty/astronomy/issues/187 + // For years far from the present, the seasons search was sometimes failing. + for (year in 1 .. 9999) + Astronomy.seasons(year) + } + //---------------------------------------------------------------------------------------- @Test diff --git a/source/python/astronomy/astronomy.py b/source/python/astronomy/astronomy.py index 616504e7..651f9c99 100644 --- a/source/python/astronomy/astronomy.py +++ b/source/python/astronomy/astronomy.py @@ -5642,7 +5642,7 @@ def SearchSunLongitude(targetLon, startTime, limitDays): Time or `None` """ t2 = startTime.AddDays(limitDays) - return Search(_sun_offset, targetLon, startTime, t2, 1.0) + return Search(_sun_offset, targetLon, startTime, t2, 0.01) def MoonPhase(time): """Returns the Moon's phase as an angle from 0 to 360 degrees. @@ -6444,7 +6444,7 @@ class SeasonInfo: def _FindSeasonChange(targetLon, year, month, day): startTime = Time.Make(year, month, day, 0, 0, 0) - time = SearchSunLongitude(targetLon, startTime, 4.0) + time = SearchSunLongitude(targetLon, startTime, 20.0) if time is None: # We should always be able to find a season change. raise InternalError() @@ -6486,10 +6486,17 @@ def Seasons(year): ------- SeasonInfo """ - mar_equinox = _FindSeasonChange(0, year, 3, 19) - jun_solstice = _FindSeasonChange(90, year, 6, 19) - sep_equinox = _FindSeasonChange(180, year, 9, 21) - dec_solstice = _FindSeasonChange(270, year, 12, 20) + # https://github.com/cosinekitty/astronomy/issues/187 + # Solstices and equinoxes drift over long spans of time, + # due to precession of the Earth's axis. + # Therefore, we have to search a wider range of time than + # one might expect. It turns out this has very little + # effect on efficiency, thanks to the quick convergence + # of quadratic interpolation inside the `Search` function. + mar_equinox = _FindSeasonChange( 0, year, 3, 10) + jun_solstice = _FindSeasonChange( 90, year, 6, 10) + sep_equinox = _FindSeasonChange(180, year, 9, 10) + dec_solstice = _FindSeasonChange(270, year, 12, 10) return SeasonInfo(mar_equinox, jun_solstice, sep_equinox, dec_solstice) def _MoonDistance(time):