From 0c4e3dfec3de5ce5b2a9dc0ffb267001c5bae8a7 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Fri, 8 Apr 2022 18:18:45 -0400 Subject: [PATCH] Fixed #187 - Seasons() fixes from kotlin branch. Backported fixes to the Seasons functions in C, C#, Python, and JavaScript. They were failing to find equinoxes and/or solstices for distant year values. Also brought over some other minor code cleanup. --- 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 | 49 +++++++++- generate/dotnet/csharp_test/csharp_test.cs | 51 +++++++--- generate/template/astronomy.c | 22 +++-- generate/template/astronomy.cs | 69 +++++++------- generate/template/astronomy.py | 19 ++-- generate/template/astronomy.ts | 12 +-- generate/test.js | 31 +++++- generate/test.py | 25 ++++- source/c/astronomy.c | 22 +++-- source/csharp/README.md | 4 +- source/csharp/astronomy.cs | 69 +++++++------- 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 +-- source/python/astronomy/astronomy.py | 19 ++-- 28 files changed, 391 insertions(+), 245 deletions(-) diff --git a/README.md b/README.md index 5d66d06a..fdb21d84 100644 --- a/README.md +++ b/README.md @@ -151,7 +151,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 92b3fd14..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) @@ -997,7 +1036,7 @@ static int MoonPhase(void) astro_angle_result_t result; double degree_error, arcmin, max_arcmin = 0.0; double diff_seconds, maxdiff = 0.0; - const double threshold_seconds = 120.0; /* max tolerable prediction error in seconds */ + const double threshold_seconds = 90.0; /* max tolerable prediction error in seconds */ astro_moon_quarter_t mq; char line[200]; @@ -6187,13 +6226,13 @@ static int SiderealTimeTest(void) int error; astro_time_t time; double gast, diff; - const double correct = 140.975528 / 15; /* https://eco.mtk.nao.ac.jp/cgi-bin/koyomi/cande/gst_en.cgi */ + const double correct = 9.398368460418821; time = Astronomy_MakeTime(2022, 3, 15, 21, 50, 0.0); gast = Astronomy_SiderealTime(&time); - diff = 3.6e+6 * ABS(gast - correct); /* calculate error in milliseconds */ - printf("C SiderealTimeTest: gast=%0.10lf, correct=%0.10lf, diff=%0.3lf milliseconds.\n", gast, correct, diff); - if (diff > 0.263) + diff = ABS(gast - correct); + printf("C SiderealTimeTest: gast=%0.15f, correct=%0.15lf, diff=%0.3le hours.\n", gast, correct, diff); + if (diff > 1.0e-15) FAIL("C SiderealTimeTest: EXCESSIVE ERROR\n"); printf("C SiderealTimeTest: PASS\n"); diff --git a/generate/dotnet/csharp_test/csharp_test.cs b/generate/dotnet/csharp_test/csharp_test.cs index f30c2fb3..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,12 +392,23 @@ 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"; using (StreamReader infile = File.OpenText(filename)) { - const double threshold_seconds = 120.0; + const double threshold_seconds = 90.0; int lnum = 0; string line; double max_arcmin = 0.0; @@ -626,8 +638,8 @@ namespace csharp_test Observer observer = new Observer(); bool foundObserver = false; AstroTime r_search_date = null, s_search_date = null; - AstroTime r_evt = null, s_evt = null; /* rise event, set event: search results */ - AstroTime a_evt = null, b_evt = null; /* chronologically first and second events */ + AstroTime r_evt = null, s_evt = null; // rise event, set event: search results + AstroTime a_evt = null, b_evt = null; // chronologically first and second events Direction a_dir = Direction.Rise, b_dir = Direction.Rise; const double nudge_days = 0.01; double sum_minutes = 0.0; @@ -656,8 +668,8 @@ namespace csharp_test Direction direction = (m.Groups[9].Value == "r") ? Direction.Rise : Direction.Set; var correct_date = new AstroTime(year, month, day, hour, minute, 0); - /* Every time we see a new geographic location or body, start a new iteration */ - /* of finding all rise/set times for that UTC calendar year. */ + // Every time we see a new geographic location or body, start a new iteration + // of finding all rise/set times for that UTC calendar year. if (!foundObserver || observer.latitude != latitude || observer.longitude != longitude || current_body != body) { current_body = body; @@ -670,6 +682,9 @@ namespace csharp_test if (b_evt != null) { + // The previous iteration found two events. + // We already processed the earlier event (a_evt). + // Now it is time to process the later event (b_evt). a_evt = b_evt; a_dir = b_dir; b_evt = null; @@ -686,11 +701,13 @@ namespace csharp_test s_evt = Astronomy.SearchRiseSet(body, observer, Direction.Set, s_search_date, 366.0); if (s_evt == null) { - Console.WriteLine("C# RiseSetTest({0} line {1}): Did not find {2} rise event.", filename, lnum, body); + Console.WriteLine("C# RiseSetTest({0} line {1}): Did not find {2} set event.", filename, lnum, body); return 1; } - /* Expect the current event to match the earlier of the found dates. */ + // Sort the two events chronologically. + // We will check the earlier event in this iteration, + // and check the later event in the next iteration. if (r_evt.tt < s_evt.tt) { a_evt = r_evt; @@ -706,16 +723,20 @@ namespace csharp_test b_dir = Direction.Rise; } - /* Nudge the event times forward a tiny amount. */ + // Nudge the event times forward a tiny amount. + // This prevents us from getting stuck in a loop, finding the same event repeatedly. r_search_date = r_evt.AddDays(nudge_days); s_search_date = s_evt.AddDays(nudge_days); } + // Expect the current search result to match the earlier of the found dates. + if (a_dir != direction) { Console.WriteLine("C# RiseSetTest({0} line {1}): expected dir={2} but found {3}", filename, lnum, a_dir, direction); return 1; } + double error_minutes = (24.0 * 60.0) * abs(a_evt.tt - correct_date.tt); sum_minutes += error_minutes * error_minutes; if (error_minutes > max_minutes) @@ -3164,9 +3185,9 @@ namespace csharp_test Astronomy.SearchAltitude(Body.Sun, observer, Direction.Rise, searchDate, 1.0, -18.0), // astronomical dawn Astronomy.SearchAltitude(Body.Sun, observer, Direction.Rise, searchDate, 1.0, -12.0), // nautical dawn Astronomy.SearchAltitude(Body.Sun, observer, Direction.Rise, searchDate, 1.0, -6.0), // civil dawn - Astronomy.SearchAltitude(Body.Sun, observer, Direction.Set, searchDate, 1.0, -6.0), // civil dawn - Astronomy.SearchAltitude(Body.Sun, observer, Direction.Set, searchDate, 1.0, -12.0), // nautical dawn - Astronomy.SearchAltitude(Body.Sun, observer, Direction.Set, searchDate, 1.0, -18.0), // astronomical dawn + Astronomy.SearchAltitude(Body.Sun, observer, Direction.Set, searchDate, 1.0, -6.0), // civil dawn + Astronomy.SearchAltitude(Body.Sun, observer, Direction.Set, searchDate, 1.0, -12.0), // nautical dawn + Astronomy.SearchAltitude(Body.Sun, observer, Direction.Set, searchDate, 1.0, -18.0), // astronomical dawn }; for (int i = 0; i < 6; ++i) @@ -3520,12 +3541,12 @@ namespace csharp_test static int SiderealTimeTest() { - const double correct = 140.975528 / 15; // https://eco.mtk.nao.ac.jp/cgi-bin/koyomi/cande/gst_en.cgi + const double correct = 9.398368460418821; var time = new AstroTime(2022, 3, 15, 21, 50, 0); double gast = Astronomy.SiderealTime(time); - double diff = 3.6e+6 * abs(gast - correct); // calculate error in milliseconds - Console.WriteLine($"C# SiderealTimeTest: gast={gast:F10}, correct={correct:F10}, diff={diff:F3} milliseconds."); - if (diff > 0.263) + double diff = abs(gast - correct); + Console.WriteLine($"C# SiderealTimeTest: gast={gast:F10}, correct={correct:F10}, diff={diff:E3}."); + if (diff > 1.0e-15) { Console.WriteLine("C# SiderealTimeTest: EXCESSIVE ERROR"); return 1; 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 25d0ab5b..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}"); } /// @@ -4729,7 +4737,7 @@ $ASTRO_IAU_DATA() /// to calculate all equinoxes and solstices for a given calendar year. /// /// The function searches the window of time specified by `startTime` and `startTime+limitDays`. - /// The search will return an error if the Sun never reaches the longitude `targetLon` or + /// The search will return `null` if the Sun never reaches the longitude `targetLon` or /// if the window is so large that the longitude ranges more than 180 degrees within it. /// It is recommended to keep the window smaller than 10 days when possible. /// @@ -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); } /// @@ -4779,8 +4787,6 @@ $ASTRO_IAU_DATA() /// thus the difference would equal zero at the moment in time the planet reaches the /// desired longitude. /// - /// Every call to `func.Eval` must either return a valid #AstroTime or throw an exception. - /// /// The search calls `func.Eval` repeatedly to rapidly narrow in on any ascending /// root within the time window specified by `t1` and `t2`. The search never /// reports a solution outside this time window. @@ -4856,8 +4862,8 @@ $ASTRO_IAU_DATA() /* Try to find a parabola that passes through the 3 points we have sampled: */ /* (t1,f1), (tmid,fmid), (t2,f2) */ - double q_x, q_ut, q_df_dt; - if (QuadInterp(tmid.ut, t2.ut - tmid.ut, f1, fmid, f2, out q_x, out q_ut, out q_df_dt)) + double q_ut, q_df_dt; + if (QuadInterp(tmid.ut, t2.ut - tmid.ut, f1, fmid, f2, out q_ut, out q_df_dt)) { var tq = new AstroTime(q_ut); double fq = func.Eval(tq); @@ -4880,9 +4886,8 @@ $ASTRO_IAU_DATA() { if ((tright.ut - t1.ut)*(tright.ut - t2.ut) < 0.0) { - double fleft, fright; - fleft = func.Eval(tleft); - fright = func.Eval(tright); + double fleft = func.Eval(tleft); + double fright = func.Eval(tright); if (fleft<0.0 && fright>=0.0) { f1 = fleft; @@ -4923,12 +4928,13 @@ $ASTRO_IAU_DATA() private static bool QuadInterp( double tm, double dt, double fa, double fm, double fb, - out double out_x, out double out_t, out double out_df_dt) + out double out_t, out double out_df_dt) { double Q, R, S; - double u, ru, x1, x2; + double u, ru; + double x, x1, x2; - out_x = out_t = out_df_dt = 0.0; + out_t = out_df_dt = 0.0; Q = (fb + fa)/2.0 - fm; R = (fb - fa)/2.0; @@ -4939,8 +4945,8 @@ $ASTRO_IAU_DATA() /* This is a line, not a parabola. */ if (R == 0.0) return false; /* This is a HORIZONTAL line... can't make progress! */ - out_x = -S / R; - if (out_x < -1.0 || out_x > +1.0) + x = -S / R; + if (x < -1.0 || x > +1.0) return false; /* out of bounds */ } else @@ -4957,16 +4963,16 @@ $ASTRO_IAU_DATA() { if (-1.0 <= x2 && x2 <= +1.0) return false; /* two roots are within bounds; we require a unique zero-crossing. */ - out_x = x1; + x = x1; } else if (-1.0 <= x2 && x2 <= +1.0) - out_x = x2; + x = x2; else return false; /* neither root is within bounds */ } - out_t = tm + out_x*dt; - out_df_dt = (2*Q*out_x + R) / dt; + out_t = tm + x*dt; + out_df_dt = (2*Q*x + R) / dt; return true; /* success */ } @@ -5184,13 +5190,12 @@ $ASTRO_IAU_DATA() */ HourAngleInfo evt_before, evt_after; - AstroTime time_start = startTime; - double alt_before = context.Eval(time_start); + double alt_before = context.Eval(startTime); AstroTime time_before; if (alt_before > 0.0) { /* We are past the sought event, so we have to wait for the next "before" event (culm/bottom). */ - evt_before = SearchHourAngle(body, observer, ha_before, time_start); + evt_before = SearchHourAngle(body, observer, ha_before, startTime); time_before = evt_before.time; alt_before = context.Eval(time_before); } @@ -5198,7 +5203,7 @@ $ASTRO_IAU_DATA() { /* We are before or at the sought event, so we find the next "after" event (bottom/culm), */ /* and use the current time as the "before" event. */ - time_before = time_start; + time_before = startTime; } evt_after = SearchHourAngle(body, observer, ha_after, time_before); @@ -5209,16 +5214,16 @@ $ASTRO_IAU_DATA() if (alt_before <= 0.0 && alt_after > 0.0) { /* Search between evt_before and evt_after for the desired event. */ - AstroTime result = Search(context, time_before, evt_after.time, 1.0); - if (result != null) - return result; + AstroTime time = Search(context, time_before, evt_after.time, 1.0); + if (time != null) + return time; } /* If we didn't find the desired event, use evt_after.time to find the next before-event. */ evt_before = SearchHourAngle(body, observer, ha_before, evt_after.time); evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time); - if (evt_before.time.ut >= time_start.ut + limitDays) + if (evt_before.time.ut >= startTime.ut + limitDays) return null; time_before = evt_before.time; 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 98e08069..4110a0c3 100644 --- a/generate/test.js +++ b/generate/test.js @@ -141,7 +141,7 @@ function MoonPhase() { function SearchYear(year, data, index) { const millis_per_minute = 60*1000; - const threshold_minutes = 2; // max tolerable prediction error in minutes + const threshold_minutes = 1.5; // max tolerable prediction error in minutes let date = new Date(Date.UTC(year, 0, 1)); let maxdiff = 0; let count = 0; @@ -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; } @@ -2705,10 +2725,10 @@ function LagrangeTest() { function SiderealTimeTest() { const date = new Date('2022-03-15T21:50:00Z'); const gast = Astronomy.SiderealTime(date); - const correct = 140.975528 / 15; // https://eco.mtk.nao.ac.jp/cgi-bin/koyomi/cande/gst_en.cgi - const diff_ms = 3.6e+6 * abs(gast - correct); // calculate time error in milliseconds - console.log(`JS SiderealTimeTest: gast=${gast.toFixed(10)}, correct=${correct.toFixed(10)}, diff=${diff_ms.toFixed(3)} milliseconds.`); - if (diff_ms > 0.263) { + const correct = 9.398368460418821; + const diff = abs(gast - correct); + console.log(`JS SiderealTimeTest: gast=${gast.toFixed(15)}, correct=${correct.toFixed(15)}, diff=${diff.toExponential(3)}`); + if (diff > 1.0e-15) { console.error('JS SiderealTimeTest: FAIL - excessive error.'); return 1; } @@ -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 8bc7819b..8a1475ef 100755 --- a/generate/test.py +++ b/generate/test.py @@ -217,10 +217,24 @@ 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'): - threshold_seconds = 120.0 # max tolerable prediction error in seconds + threshold_seconds = 90.0 # max tolerable prediction error in seconds max_arcmin = 0.0 maxdiff = 0.0 quarter_count = 0 @@ -2424,12 +2438,12 @@ def Lagrange(): #----------------------------------------------------------------------------------------------------------- def SiderealTime(): - correct = 140.975528 / 15 # https://eco.mtk.nao.ac.jp/cgi-bin/koyomi/cande/gst_en.cgi + correct = 9.398368460418821 time = astronomy.Time.Make(2022, 3, 15, 21, 50, 0) gast = astronomy.SiderealTime(time) - diff = 3.6e+6 * abs(gast - correct) # calculate error in millseconds - print('PY SiderealTime: gast={:0.10f}, correct={:0.10f}, diff={:0.3f} milliseconds.'.format(gast, correct, diff)) - if diff > 0.263: + diff = abs(gast - correct) + print('PY SiderealTime: gast={:0.15f}, correct={:0.15f}, diff={:0.3e}'.format(gast, correct, diff)) + if diff > 1.0e-15: print('PY SiderealTime: EXCESSIVE ERROR') return 1 print('PY SiderealTime: PASS') @@ -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/README.md b/source/csharp/README.md index edda7ed9..509a3ae4 100644 --- a/source/csharp/README.md +++ b/source/csharp/README.md @@ -1501,8 +1501,6 @@ It could subtract the target longitude from the actual longitude at a given time thus the difference would equal zero at the moment in time the planet reaches the desired longitude. -Every call to `func.Eval` must either return a valid [`AstroTime`](#AstroTime) or throw an exception. - The search calls `func.Eval` repeatedly to rapidly narrow in on any ascending root within the time window specified by `t1` and `t2`. The search never reports a solution outside this time window. @@ -1896,7 +1894,7 @@ However, it is usually more convenient and efficient to call [`Astronomy.Seasons to calculate all equinoxes and solstices for a given calendar year. The function searches the window of time specified by `startTime` and `startTime+limitDays`. -The search will return an error if the Sun never reaches the longitude `targetLon` or +The search will return `null` if the Sun never reaches the longitude `targetLon` or if the window is so large that the longitude ranges more than 180 degrees within it. It is recommended to keep the window smaller than 10 days when possible. diff --git a/source/csharp/astronomy.cs b/source/csharp/astronomy.cs index 54f90446..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}"); } /// @@ -5941,7 +5949,7 @@ namespace CosineKitty /// to calculate all equinoxes and solstices for a given calendar year. /// /// The function searches the window of time specified by `startTime` and `startTime+limitDays`. - /// The search will return an error if the Sun never reaches the longitude `targetLon` or + /// The search will return `null` if the Sun never reaches the longitude `targetLon` or /// if the window is so large that the longitude ranges more than 180 degrees within it. /// It is recommended to keep the window smaller than 10 days when possible. /// @@ -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); } /// @@ -5991,8 +5999,6 @@ namespace CosineKitty /// thus the difference would equal zero at the moment in time the planet reaches the /// desired longitude. /// - /// Every call to `func.Eval` must either return a valid #AstroTime or throw an exception. - /// /// The search calls `func.Eval` repeatedly to rapidly narrow in on any ascending /// root within the time window specified by `t1` and `t2`. The search never /// reports a solution outside this time window. @@ -6068,8 +6074,8 @@ namespace CosineKitty /* Try to find a parabola that passes through the 3 points we have sampled: */ /* (t1,f1), (tmid,fmid), (t2,f2) */ - double q_x, q_ut, q_df_dt; - if (QuadInterp(tmid.ut, t2.ut - tmid.ut, f1, fmid, f2, out q_x, out q_ut, out q_df_dt)) + double q_ut, q_df_dt; + if (QuadInterp(tmid.ut, t2.ut - tmid.ut, f1, fmid, f2, out q_ut, out q_df_dt)) { var tq = new AstroTime(q_ut); double fq = func.Eval(tq); @@ -6092,9 +6098,8 @@ namespace CosineKitty { if ((tright.ut - t1.ut)*(tright.ut - t2.ut) < 0.0) { - double fleft, fright; - fleft = func.Eval(tleft); - fright = func.Eval(tright); + double fleft = func.Eval(tleft); + double fright = func.Eval(tright); if (fleft<0.0 && fright>=0.0) { f1 = fleft; @@ -6135,12 +6140,13 @@ namespace CosineKitty private static bool QuadInterp( double tm, double dt, double fa, double fm, double fb, - out double out_x, out double out_t, out double out_df_dt) + out double out_t, out double out_df_dt) { double Q, R, S; - double u, ru, x1, x2; + double u, ru; + double x, x1, x2; - out_x = out_t = out_df_dt = 0.0; + out_t = out_df_dt = 0.0; Q = (fb + fa)/2.0 - fm; R = (fb - fa)/2.0; @@ -6151,8 +6157,8 @@ namespace CosineKitty /* This is a line, not a parabola. */ if (R == 0.0) return false; /* This is a HORIZONTAL line... can't make progress! */ - out_x = -S / R; - if (out_x < -1.0 || out_x > +1.0) + x = -S / R; + if (x < -1.0 || x > +1.0) return false; /* out of bounds */ } else @@ -6169,16 +6175,16 @@ namespace CosineKitty { if (-1.0 <= x2 && x2 <= +1.0) return false; /* two roots are within bounds; we require a unique zero-crossing. */ - out_x = x1; + x = x1; } else if (-1.0 <= x2 && x2 <= +1.0) - out_x = x2; + x = x2; else return false; /* neither root is within bounds */ } - out_t = tm + out_x*dt; - out_df_dt = (2*Q*out_x + R) / dt; + out_t = tm + x*dt; + out_df_dt = (2*Q*x + R) / dt; return true; /* success */ } @@ -6396,13 +6402,12 @@ namespace CosineKitty */ HourAngleInfo evt_before, evt_after; - AstroTime time_start = startTime; - double alt_before = context.Eval(time_start); + double alt_before = context.Eval(startTime); AstroTime time_before; if (alt_before > 0.0) { /* We are past the sought event, so we have to wait for the next "before" event (culm/bottom). */ - evt_before = SearchHourAngle(body, observer, ha_before, time_start); + evt_before = SearchHourAngle(body, observer, ha_before, startTime); time_before = evt_before.time; alt_before = context.Eval(time_before); } @@ -6410,7 +6415,7 @@ namespace CosineKitty { /* We are before or at the sought event, so we find the next "after" event (bottom/culm), */ /* and use the current time as the "before" event. */ - time_before = time_start; + time_before = startTime; } evt_after = SearchHourAngle(body, observer, ha_after, time_before); @@ -6421,16 +6426,16 @@ namespace CosineKitty if (alt_before <= 0.0 && alt_after > 0.0) { /* Search between evt_before and evt_after for the desired event. */ - AstroTime result = Search(context, time_before, evt_after.time, 1.0); - if (result != null) - return result; + AstroTime time = Search(context, time_before, evt_after.time, 1.0); + if (time != null) + return time; } /* If we didn't find the desired event, use evt_after.time to find the next before-event. */ evt_before = SearchHourAngle(body, observer, ha_before, evt_after.time); evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time); - if (evt_before.time.ut >= time_start.ut + limitDays) + if (evt_before.time.ut >= startTime.ut + limitDays) return null; time_before = evt_before.time; 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/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):