From 3bed4a9bdcd0b2c34b71e48999851f4e60af6588 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Mon, 26 Sep 2022 21:08:37 -0400 Subject: [PATCH] PY SearchMoonPhase: allow searching backward in time. Enhanced the Python function SearchMoonPhase to allow searching forward in time when the `limitDays` argument is positive, or backward in time when `limitDays` is negative. Added unit test "moon_reverse" to verify this new feature. --- demo/python/astronomy.py | 31 ++++++++++----- generate/template/astronomy.py | 31 ++++++++++----- generate/test.py | 58 ++++++++++++++++++++++++++-- source/python/README.md | 2 +- source/python/astronomy/astronomy.py | 31 ++++++++++----- 5 files changed, 122 insertions(+), 31 deletions(-) diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py index aac0e155..e93851eb 100644 --- a/demo/python/astronomy.py +++ b/demo/python/astronomy.py @@ -5893,7 +5893,9 @@ def SearchMoonPhase(targetLon, startTime, limitDays): startTime : Time The beginning of the time window in which to search for the Moon reaching the specified phase. limitDays : float - The number of days after `startTime` that limits the time window for the search. + The number of days away from `startTime` that limits the time window for the search. + If the value is negative, the search is performed into the past from `startTime`. + Otherwise, the search is performed into the future from `startTime`. Returns ------- @@ -5911,16 +5913,27 @@ def SearchMoonPhase(targetLon, startTime, limitDays): # But we must return None if the final result goes beyond limitDays after startTime. uncertainty = 1.5 ya = _moon_offset(targetLon, startTime) - if ya > 0.0: - ya -= 360.0 # force searching forward in time, not backward - est_dt = -(_MEAN_SYNODIC_MONTH * ya) / 360.0 - dt1 = est_dt - uncertainty - if dt1 > limitDays: - return None # not possible for moon phase to occur within the specified window - dt2 = min(limitDays, est_dt + uncertainty) + if limitDays < 0.0: + # Search backward in time. + if ya < 0.0: + ya += 360.0 # force searching forward in time, not backward + est_dt = -(_MEAN_SYNODIC_MONTH * ya) / 360.0 + dt2 = est_dt + uncertainty + if dt2 < limitDays: + return None # not possible for moon phase to occur within the specified window + dt1 = max(limitDays, est_dt - uncertainty) + else: + # Search forward in time. + if ya > 0.0: + ya -= 360.0 # force searching forward in time, not backward + est_dt = -(_MEAN_SYNODIC_MONTH * ya) / 360.0 + dt1 = est_dt - uncertainty + if dt1 > limitDays: + return None # not possible for moon phase to occur within the specified window + dt2 = min(limitDays, est_dt + uncertainty) t1 = startTime.AddDays(dt1) t2 = startTime.AddDays(dt2) - return Search(_moon_offset, targetLon, t1, t2, 1.0) + return Search(_moon_offset, targetLon, t1, t2, 0.1) class MoonQuarter: """A lunar quarter event along with its date and time. diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index d2f2bb18..afbc646e 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -3851,7 +3851,9 @@ def SearchMoonPhase(targetLon, startTime, limitDays): startTime : Time The beginning of the time window in which to search for the Moon reaching the specified phase. limitDays : float - The number of days after `startTime` that limits the time window for the search. + The number of days away from `startTime` that limits the time window for the search. + If the value is negative, the search is performed into the past from `startTime`. + Otherwise, the search is performed into the future from `startTime`. Returns ------- @@ -3869,16 +3871,27 @@ def SearchMoonPhase(targetLon, startTime, limitDays): # But we must return None if the final result goes beyond limitDays after startTime. uncertainty = 1.5 ya = _moon_offset(targetLon, startTime) - if ya > 0.0: - ya -= 360.0 # force searching forward in time, not backward - est_dt = -(_MEAN_SYNODIC_MONTH * ya) / 360.0 - dt1 = est_dt - uncertainty - if dt1 > limitDays: - return None # not possible for moon phase to occur within the specified window - dt2 = min(limitDays, est_dt + uncertainty) + if limitDays < 0.0: + # Search backward in time. + if ya < 0.0: + ya += 360.0 # force searching forward in time, not backward + est_dt = -(_MEAN_SYNODIC_MONTH * ya) / 360.0 + dt2 = est_dt + uncertainty + if dt2 < limitDays: + return None # not possible for moon phase to occur within the specified window + dt1 = max(limitDays, est_dt - uncertainty) + else: + # Search forward in time. + if ya > 0.0: + ya -= 360.0 # force searching forward in time, not backward + est_dt = -(_MEAN_SYNODIC_MONTH * ya) / 360.0 + dt1 = est_dt - uncertainty + if dt1 > limitDays: + return None # not possible for moon phase to occur within the specified window + dt2 = min(limitDays, est_dt + uncertainty) t1 = startTime.AddDays(dt1) t2 = startTime.AddDays(dt2) - return Search(_moon_offset, targetLon, t1, t2, 1.0) + return Search(_moon_offset, targetLon, t1, t2, 0.1) class MoonQuarter: """A lunar quarter event along with its date and time. diff --git a/generate/test.py b/generate/test.py index 13a1660d..12f3c961 100755 --- a/generate/test.py +++ b/generate/test.py @@ -9,6 +9,7 @@ import astronomy #----------------------------------------------------------------------------------------------------------- Verbose = False +SECONDS_PER_DAY = 86400.0 def Debug(text): if Verbose: @@ -292,7 +293,7 @@ def MoonPhase(filename = 'moonphase/moonphases.txt'): quarter_count += 1 # Make sure the time matches what we expect. - diff_seconds = vabs(mq.time.tt - expected_time.tt) * (24.0 * 3600.0) + diff_seconds = vabs(mq.time.tt - expected_time.tt) * SECONDS_PER_DAY if diff_seconds > threshold_seconds: print('PY MoonPhase({} line {}): excessive time error {:0.3f} seconds.'.format(filename, lnum, diff_seconds)) return 1 @@ -1267,7 +1268,7 @@ def LunarEclipseIssue78(): eclipse = astronomy.SearchLunarEclipse(astronomy.Time.Make(2020, 12, 19, 0, 0, 0)) expected_peak = astronomy.Time.Make(2021, 5, 26, 11, 18, 42) # https://www.timeanddate.com/eclipse/lunar/2021-may-26 - dt = (expected_peak.tt - eclipse.peak.tt) * (24.0 * 3600.0) + dt = (expected_peak.tt - eclipse.peak.tt) * SECONDS_PER_DAY if vabs(dt) > 40.0: print('LunarEclipseIssue78: Excessive prediction error = {} seconds.'.format(dt)) return 1 @@ -2200,7 +2201,7 @@ def Twilight(): for i in range(6): correct = correctTimes[i] calc = calcTimes[i] - diff = 86400.0 * vabs(calc.ut - correct.ut) + diff = SECONDS_PER_DAY * vabs(calc.ut - correct.ut) if diff > tolerance_seconds: print('PY Twilight({} line {}): EXCESSIVE ERROR = {} seconds for case {}'.format(filename, lnum, diff, i)) return 1 @@ -2431,6 +2432,56 @@ def MoonNodes(): #----------------------------------------------------------------------------------------------------------- +def MoonReverse(): + # Verify that SearchMoonPhase works both forward and backward in time. + numNewMoons = 5000 + utList = [] + dtMin = +1000.0 + dtMax = -1000.0 + + # Search forward in time from 1800 to find consecutive new moon events. + time = astronomy.Time.Make(1800, 1, 1, 0, 0, 0.0) + for i in range(numNewMoons): + result = astronomy.SearchMoonPhase(0.0, time, +40.0) + if result is None: + print('PY MoonReverse(i={}): failed to find new moon after {}'.format(i, time)) + return 1 + utList.append(result.ut) + if i > 0: + # Verify that consecutive new moons are reasonably close to the synodic period (29.5 days) apart. + dt = v(utList[i] - utList[i-1]) + if dt < dtMin: + dtMin = dt + if dt > dtMax: + dtMax = dt + time = result.AddDays(+0.1) + + Debug('PY MoonReverse: dtMin={:0.3f} days, dtMax={:0.3f} days.'.format(dtMin, dtMax)) + if (dtMin < 29.273) or (dtMax > 29.832): + print('PY MoonReverse: Time between consecutive new moons is suspicious.') + return 1 + + # Do a reverse chronological search and make sure the results are consistent with the forward search. + time = time.AddDays(20.0) + maxDiff = 0.0 + for i in range(numNewMoons-1, -1, -1): + result = astronomy.SearchMoonPhase(0.0, time, -40.0) + if result is None: + print('PY MoonReverse(i={}): failed to find new moon before {}'.format(i, time)) + return 1 + diff = SECONDS_PER_DAY * abs(result.ut - utList[i]) + if diff > maxDiff: + maxDiff = diff + time = result.AddDays(-0.1) + + print('PY MoonReverse: Maximum discrepancy in reverse search = {:0.3f} seconds.'.format(maxDiff)) + if maxDiff > 0.128: + print('PY MoonReverse: EXCESSIVE DISCREPANCY in reverse search.') + return 1 + return 0 + +#----------------------------------------------------------------------------------------------------------- + class LagrangeFunc: def __init__(self, point, major_body, minor_body): self.point = point @@ -2689,6 +2740,7 @@ UnitTests = { 'magnitude': Magnitude, 'moon': GeoMoon, 'moon_nodes': MoonNodes, + 'moon_reverse': MoonReverse, 'moonphase': MoonPhase, 'planet_apsis': PlanetApsis, 'pluto': PlutoCheck, diff --git a/source/python/README.md b/source/python/README.md index cddfbcff..569cef8d 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -2952,7 +2952,7 @@ This function is useful for finding general phase angles outside those four quar | --- | --- | --- | | `float` | `targetLon` | The difference in geocentric longitude between the Sun and Moon that specifies the lunar phase being sought. This can be any value in the range [0, 360). Certain values have conventional names: 0 = new moon, 90 = first quarter, 180 = full moon, 270 = third quarter. | | [`Time`](#Time) | `startTime` | The beginning of the time window in which to search for the Moon reaching the specified phase. | -| `float` | `limitDays` | The number of days after `startTime` that limits the time window for the search. | +| `float` | `limitDays` | The number of days away from `startTime` that limits the time window for the search. If the value is negative, the search is performed into the past from `startTime`. Otherwise, the search is performed into the future from `startTime`. | **Returns**: [`Time`](#Time) or `None` diff --git a/source/python/astronomy/astronomy.py b/source/python/astronomy/astronomy.py index aac0e155..e93851eb 100644 --- a/source/python/astronomy/astronomy.py +++ b/source/python/astronomy/astronomy.py @@ -5893,7 +5893,9 @@ def SearchMoonPhase(targetLon, startTime, limitDays): startTime : Time The beginning of the time window in which to search for the Moon reaching the specified phase. limitDays : float - The number of days after `startTime` that limits the time window for the search. + The number of days away from `startTime` that limits the time window for the search. + If the value is negative, the search is performed into the past from `startTime`. + Otherwise, the search is performed into the future from `startTime`. Returns ------- @@ -5911,16 +5913,27 @@ def SearchMoonPhase(targetLon, startTime, limitDays): # But we must return None if the final result goes beyond limitDays after startTime. uncertainty = 1.5 ya = _moon_offset(targetLon, startTime) - if ya > 0.0: - ya -= 360.0 # force searching forward in time, not backward - est_dt = -(_MEAN_SYNODIC_MONTH * ya) / 360.0 - dt1 = est_dt - uncertainty - if dt1 > limitDays: - return None # not possible for moon phase to occur within the specified window - dt2 = min(limitDays, est_dt + uncertainty) + if limitDays < 0.0: + # Search backward in time. + if ya < 0.0: + ya += 360.0 # force searching forward in time, not backward + est_dt = -(_MEAN_SYNODIC_MONTH * ya) / 360.0 + dt2 = est_dt + uncertainty + if dt2 < limitDays: + return None # not possible for moon phase to occur within the specified window + dt1 = max(limitDays, est_dt - uncertainty) + else: + # Search forward in time. + if ya > 0.0: + ya -= 360.0 # force searching forward in time, not backward + est_dt = -(_MEAN_SYNODIC_MONTH * ya) / 360.0 + dt1 = est_dt - uncertainty + if dt1 > limitDays: + return None # not possible for moon phase to occur within the specified window + dt2 = min(limitDays, est_dt + uncertainty) t1 = startTime.AddDays(dt1) t2 = startTime.AddDays(dt2) - return Search(_moon_offset, targetLon, t1, t2, 1.0) + return Search(_moon_offset, targetLon, t1, t2, 0.1) class MoonQuarter: """A lunar quarter event along with its date and time.