From b91b1d905faa3f377eb9fe4554a04a493e8edc51 Mon Sep 17 00:00:00 2001
From: Don Cross <cosinekitty@gmail.com>
Date: Sat, 1 Oct 2022 21:44:04 -0400
Subject: [PATCH] Python: Reverse chrono search for rise/set, hour angles.

The following Python functions now support searching
in forward or reverse chronological order:

    SearchRiseSet
    SearchAltitude
    SearchHourAngle

Made some minor performance improvements to the
other implementations: return sooner if we
go past time window.
---
 demo/browser/astronomy.browser.js             |   2 +-
 demo/nodejs/astronomy.js                      |   2 +-
 demo/nodejs/calendar/astronomy.ts             |   3 +-
 demo/python/astronomy.py                      | 125 +++++++++++++++---
 generate/dotnet/csharp_test/csharp_test.cs    |  14 +-
 generate/template/astronomy.cs                |   3 +-
 generate/template/astronomy.kt                |   4 +-
 generate/template/astronomy.py                | 125 +++++++++++++++---
 generate/template/astronomy.ts                |   3 +-
 generate/test.py                              |  56 ++++++++
 source/csharp/astronomy.cs                    |   3 +-
 source/js/astronomy.browser.js                |   2 +-
 source/js/astronomy.browser.min.js            |   2 +-
 source/js/astronomy.js                        |   2 +-
 source/js/astronomy.min.js                    |   2 +-
 source/js/astronomy.ts                        |   3 +-
 source/js/esm/astronomy.js                    |   2 +-
 .../github/cosinekitty/astronomy/astronomy.kt |   4 +-
 source/python/README.md                       |   7 +-
 source/python/astronomy/astronomy.py          | 125 +++++++++++++++---
 20 files changed, 412 insertions(+), 77 deletions(-)

diff --git a/demo/browser/astronomy.browser.js b/demo/browser/astronomy.browser.js
index f489a034..39b85b47 100644
--- a/demo/browser/astronomy.browser.js
+++ b/demo/browser/astronomy.browser.js
@@ -4930,9 +4930,9 @@ function ForwardSearchAltitude(body, observer, direction, dateStart, limitDays,
         }
         // If we didn't find the desired event, use time_after to find the next before-event.
         evt_before = SearchHourAngle(body, observer, ha_before, evt_after.time, +1);
-        evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time, +1);
         if (evt_before.time.ut >= time_start.ut + limitDays)
             return null;
+        evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time, +1);
         time_before = evt_before.time;
         error_before = altitude_error(evt_before.time);
         error_after = altitude_error(evt_after.time);
diff --git a/demo/nodejs/astronomy.js b/demo/nodejs/astronomy.js
index 3fdbec6e..0add99fa 100644
--- a/demo/nodejs/astronomy.js
+++ b/demo/nodejs/astronomy.js
@@ -4929,9 +4929,9 @@ function ForwardSearchAltitude(body, observer, direction, dateStart, limitDays,
         }
         // If we didn't find the desired event, use time_after to find the next before-event.
         evt_before = SearchHourAngle(body, observer, ha_before, evt_after.time, +1);
-        evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time, +1);
         if (evt_before.time.ut >= time_start.ut + limitDays)
             return null;
+        evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time, +1);
         time_before = evt_before.time;
         error_before = altitude_error(evt_before.time);
         error_after = altitude_error(evt_after.time);
diff --git a/demo/nodejs/calendar/astronomy.ts b/demo/nodejs/calendar/astronomy.ts
index 1b4d1035..3a9f22ae 100644
--- a/demo/nodejs/calendar/astronomy.ts
+++ b/demo/nodejs/calendar/astronomy.ts
@@ -5480,9 +5480,9 @@ function ForwardSearchAltitude(
 
         // If we didn't find the desired event, use time_after to find the next before-event.
         evt_before = SearchHourAngle(body, observer, ha_before, evt_after.time, +1);
-        evt_after  = SearchHourAngle(body, observer, ha_after, evt_before.time, +1);
         if (evt_before.time.ut >= time_start.ut + limitDays)
             return null;
+        evt_after  = SearchHourAngle(body, observer, ha_after, evt_before.time, +1);
 
         time_before = evt_before.time;
         error_before = altitude_error(evt_before.time);
@@ -5562,7 +5562,6 @@ function BackwardSearchAltitude(
         }
 
         evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time, -1);
-
         if (evt_after.time.ut <= time_start.ut + limitDays)
             return null;
 
diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py
index 22511c94..9a74b0cb 100644
--- a/demo/python/astronomy.py
+++ b/demo/python/astronomy.py
@@ -6332,7 +6332,7 @@ class HourAngleEvent:
     def __repr__(self):
         return 'HourAngleEvent({}, {})'.format(repr(self.time), repr(self.hor))
 
-def SearchHourAngle(body, observer, hourAngle, startTime):
+def SearchHourAngle(body, observer, hourAngle, startTime, direction = +1):
     """Searches for the time when a celestial body reaches a specified hour angle as seen by an observer on the Earth.
 
     The *hour angle* of a celestial body indicates its position in the sky with respect
@@ -6367,6 +6367,10 @@ def SearchHourAngle(body, observer, hourAngle, startTime):
          body's most recent culmination.
     startTime : Time
          The date and time at which to start the search.
+    direction : int
+        The direction in time to perform the search: a positive value
+        searches forward in time, a negative value searches backward in time.
+        The function throws an exception if `direction` is zero.
 
     Returns
     -------
@@ -6378,6 +6382,9 @@ def SearchHourAngle(body, observer, hourAngle, startTime):
     if hourAngle < 0.0 or hourAngle >= 24.0:
         raise Error('Invalid hour angle.')
 
+    if direction == 0:
+        raise Error('Direction must be positive or negative.')
+
     iter_count = 0
     time = startTime
     while True:
@@ -6390,9 +6397,15 @@ def SearchHourAngle(body, observer, hourAngle, startTime):
         # the hour angle to the desired value.
         delta_sidereal_hours = math.fmod(((hourAngle + ofdate.ra - observer.longitude/15) - gast), 24.0)
         if iter_count == 1:
-            # On the first iteration, always search forward in time.
-            if delta_sidereal_hours < 0.0:
-                delta_sidereal_hours += 24.0
+            # On the first iteration, always search in the requested time direction.
+            if direction > 0:
+                # Search forward in time.
+                if delta_sidereal_hours < 0.0:
+                    delta_sidereal_hours += 24.0
+            else:
+                # Search backward in time.
+                if delta_sidereal_hours > 0.0:
+                    delta_sidereal_hours -= 24.0
         else:
             # On subsequent iterations, we make the smallest possible adjustment,
             # either forward or backward in time.
@@ -6427,8 +6440,10 @@ class Direction(enum.Enum):
     Rise = +1
     Set  = -1
 
+def _ForwardSearchAltitude(body, observer, direction, startTime, limitDays, altitude_error_func, altitude_error_context):
+    if body == Body.Earth:
+        raise EarthNotAllowedError()
 
-def _InternalSearchAltitude(body, observer, direction, startTime, limitDays, altitude_error_func, altitude_error_context):
     if direction == Direction.Rise:
         ha_before = 12.0    # minimum altitude (bottom) happens BEFORE the body rises.
         ha_after  =  0.0    # maximum altitude (culmination) happens AFTER the body rises.
@@ -6438,6 +6453,10 @@ def _InternalSearchAltitude(body, observer, direction, startTime, limitDays, alt
     else:
         raise Error('Invalid value for direction parameter')
 
+    # We cannot possibly satisfy a forward search without a positive time limit.
+    if limitDays <= 0.0:
+        return None
+
     # See if the body is currently above/below the horizon.
     # If we are looking for next rise time and the body is below the horizon,
     # we use the current time as the lower time bound and the next culmination
@@ -6449,7 +6468,7 @@ def _InternalSearchAltitude(body, observer, direction, startTime, limitDays, alt
     alt_before = altitude_error_func(altitude_error_context, time_start)
     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, time_start, +1)
         time_before = evt_before.time
         alt_before = altitude_error_func(altitude_error_context, time_before)
     else:
@@ -6457,7 +6476,7 @@ def _InternalSearchAltitude(body, observer, direction, startTime, limitDays, alt
         # and use the current time as the "before" event.
         time_before = time_start
 
-    evt_after = SearchHourAngle(body, observer, ha_after, time_before)
+    evt_after = SearchHourAngle(body, observer, ha_after, time_before, +1)
     alt_after = altitude_error_func(altitude_error_context, evt_after.time)
 
     while True:
@@ -6465,16 +6484,76 @@ def _InternalSearchAltitude(body, observer, direction, startTime, limitDays, alt
             # Search between the "before time" and the "after time" for the desired event.
             event_time = Search(altitude_error_func, altitude_error_context, time_before, evt_after.time, 1.0)
             if event_time is not None:
+                # If we found the rise/set time, but it falls outside limitDays, fail the search.
+                if event_time.ut > startTime.ut + limitDays:
+                    return None
+                # The search succeeded.
                 return event_time
         # We didn't find the desired event, so use the "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)
+        evt_before = SearchHourAngle(body, observer, ha_before, evt_after.time, +1)
         if evt_before.time.ut >= time_start.ut + limitDays:
             return None
+        evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time, +1)
         time_before = evt_before.time
         alt_before = altitude_error_func(altitude_error_context, evt_before.time)
         alt_after = altitude_error_func(altitude_error_context, evt_after.time)
 
+def _BackwardSearchAltitude(body, observer, direction, startTime, limitDays, altitude_error_func, altitude_error_context):
+    if body == Body.Earth:
+        raise EarthNotAllowedError()
+
+    if direction == Direction.Rise:
+        ha_before = 12.0    # minimum altitude (bottom) happens BEFORE the body rises.
+        ha_after  =  0.0    # maximum altitude (culmination) happens AFTER the body rises.
+    elif direction == Direction.Set:
+        ha_before =  0.0    # culmination happens BEFORE the body sets.
+        ha_after  = 12.0    # bottom happens AFTER the body sets.
+    else:
+        raise Error('Invalid value for direction parameter')
+
+    # We cannot possibly satisfy a backward search without a negative time limit.
+    if limitDays >= 0.0:
+        return None
+
+    # See if the body is currently above/below the horizon.
+    # If we are looking for previous rise time and the body is above the horizon,
+    # we use the current time as the upper time bound and the previous bottom as the lower time bound.
+    # If the body is below the horizon, we search for the previous culmination and use it
+    # as the upper time bound. Then we search for the bottom before that culmination and
+    # use it as the lower time bound.
+    # The same logic applies for finding set times; altitude_error_func and
+    # altitude_error_context ensure that the desired event is represented
+    # by ascending through zero, so the Search function works correctly.
+    time_start = startTime
+    alt_after = altitude_error_func(altitude_error_context, time_start)
+    if alt_after < 0.0:
+        evt_after = SearchHourAngle(body, observer, ha_after, time_start, -1)
+        time_after = evt_after.time
+        alt_after = altitude_error_func(altitude_error_context, time_after)
+    else:
+        time_after = time_start
+
+    evt_before = SearchHourAngle(body, observer, ha_before, time_after, -1)
+    alt_before = altitude_error_func(altitude_error_context, evt_before.time)
+
+    while True:
+        if alt_before <= 0.0 and alt_after > 0.0:
+            # Search between the "before time" and the "after time" for the desired event.
+            event_time = Search(altitude_error_func, altitude_error_context, evt_before.time, time_after, 1.0)
+            if event_time is not None:
+                # If we found the rise/set time, but it falls outside limitDays, fail the search.
+                if event_time.ut < startTime.ut + limitDays:
+                    return None
+                # The search succeeded.
+                return event_time
+        evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time, -1)
+        if evt_after.time.ut <= time_start.ut + limitDays:
+            return None
+        evt_before = SearchHourAngle(body, observer, ha_before, evt_after.time, -1)
+        time_after = evt_after.time
+        alt_before = altitude_error_func(altitude_error_context, evt_before.time)
+        alt_after = altitude_error_func(altitude_error_context, evt_after.time)
+
 
 class _peak_altitude_context:
     def __init__(self, body, direction, observer, body_radius_au):
@@ -6533,11 +6612,14 @@ def SearchRiseSet(body, observer, direction, startTime, limitDays):
     startTime : Time
         The date and time at which to start the search.
     limitDays : float
-        Limits how many days to search for a rise or set time.
+        Limits how many days to search for a rise or set time, and defines
+        the direction in time to search. When `limitDays` is positive, the
+        search is performed into the future, after `startTime`.
+        When negative, the search is performed into the past, before `startTime`.
         To limit a rise or set time to the same day, you can use a value of 1 day.
         In cases where you want to find the next rise or set time no matter how far
-        in the future (for example, for an observer near the south pole), you can pass
-        in a larger value like 365.
+        in the future (for example, for an observer near the south pole), you can
+        pass in a larger value like 365.
 
     Returns
     -------
@@ -6556,7 +6638,9 @@ def SearchRiseSet(body, observer, direction, startTime, limitDays):
         body_radius = 0.0
 
     context = _peak_altitude_context(body, direction, observer, body_radius)
-    return _InternalSearchAltitude(body, observer, direction, startTime, limitDays, _peak_altitude, context)
+    if limitDays < 0.0:
+        return _BackwardSearchAltitude(body, observer, direction, startTime, limitDays, _peak_altitude, context)
+    return _ForwardSearchAltitude(body, observer, direction, startTime, limitDays, _peak_altitude, context)
 
 
 class _altitude_error_context:
@@ -6603,8 +6687,14 @@ def SearchAltitude(body, observer, direction, dateStart, limitDays, altitude):
     startTime : Time
         The date and time at which to start the search.
     limitDays : float
-        The fractional number of days after `dateStart` that limits
-        when the altitude event is to be found. Must be a positive number.
+        Limits how many days to search for the body reaching the altitude angle,
+        and defines the direction in time to search. When `limitDays` is positive, the
+        search is performed into the future, after `startTime`.
+        When negative, the search is performed into the past, before `startTime`.
+        To limit the search to the same day, you can use a value of 1 day.
+        In cases where you want to find the altitude event no matter how far
+        in the future (for example, for an observer near the south pole), you can
+        pass in a larger value like 365.
     altitude : float
         The desired altitude angle of the body's center above (positive)
         or below (negative) the observer's local horizon, expressed in degrees.
@@ -6618,9 +6708,10 @@ def SearchAltitude(body, observer, direction, dateStart, limitDays, altitude):
     """
     if not (-90.0 <= altitude <= +90.0):
         raise Error('Invalid altitude: {}'.format(altitude))
-
     context = _altitude_error_context(body, direction, observer, altitude)
-    return _InternalSearchAltitude(body, observer, direction, dateStart, limitDays, _altitude_error_func, context)
+    if limitDays < 0.0:
+        return _BackwardSearchAltitude(body, observer, direction, dateStart, limitDays, _altitude_error_func, context)
+    return _ForwardSearchAltitude(body, observer, direction, dateStart, limitDays, _altitude_error_func, context)
 
 class SeasonInfo:
     """The dates and times of changes of season for a given calendar year.
diff --git a/generate/dotnet/csharp_test/csharp_test.cs b/generate/dotnet/csharp_test/csharp_test.cs
index 6cb6e238..5db8f114 100644
--- a/generate/dotnet/csharp_test/csharp_test.cs
+++ b/generate/dotnet/csharp_test/csharp_test.cs
@@ -821,6 +821,16 @@ namespace csharp_test
             }
         }
 
+        static Direction Toggle(Direction dir)
+        {
+            switch (dir)
+            {
+            case Direction.Rise: return Direction.Set;
+            case Direction.Set:  return Direction.Rise;
+            default: throw new ArgumentException($"Invalid direction: {dir}");
+            }
+        }
+
         static int RiseSetReverseTest()
         {
             // Verify that the rise/set search works equally well forwards and backwards in time.
@@ -852,7 +862,7 @@ namespace csharp_test
                     if (dt < dtMin) dtMin = dt;
                     if (dt > dtMax) dtMax = dt;
                 }
-                dir = (dir == Direction.Rise) ? Direction.Set : Direction.Rise;
+                dir = Toggle(dir);
                 time = result.AddDays(+nudge);
             }
 
@@ -866,7 +876,7 @@ namespace csharp_test
             // Perform the same search in reverse. Verify we get consistent rise/set times.
             for (int i = nsamples-1; i >= 0; --i)
             {
-                dir = (dir == Direction.Rise) ? Direction.Set : Direction.Rise;
+                dir = Toggle(dir);
                 AstroTime result = Astronomy.SearchRiseSet(Body.Sun, observer, dir, time, -1.0);
                 if (result == null)
                 {
diff --git a/generate/template/astronomy.cs b/generate/template/astronomy.cs
index 810e7611..8e9ede98 100644
--- a/generate/template/astronomy.cs
+++ b/generate/template/astronomy.cs
@@ -5785,11 +5785,12 @@ $ASTRO_IAU_DATA()
 
                 // 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, +1);
-                evt_after  = SearchHourAngle(body, observer, ha_after, evt_before.time, +1);
 
                 if (evt_before.time.ut >= startTime.ut + limitDays)
                     return null;
 
+                evt_after  = SearchHourAngle(body, observer, ha_after, evt_before.time, +1);
+
                 time_before = evt_before.time;
 
                 alt_before = context.Eval(evt_before.time);
diff --git a/generate/template/astronomy.kt b/generate/template/astronomy.kt
index c52ef277..5380ee8b 100644
--- a/generate/template/astronomy.kt
+++ b/generate/template/astronomy.kt
@@ -6121,11 +6121,10 @@ private fun forwardSearchAltitude(
 
         // If we didn't find the desired event, find the next hour angle bracket and try again.
         val evtBefore = searchHourAngle(body, observer, haBefore, timeAfter, +1)
-        val evtAfter  = searchHourAngle(body, observer, haAfter, timeBefore, +1)
-
         if (evtBefore.time.ut >= startTime.ut + limitDays)
             return null
 
+        val evtAfter  = searchHourAngle(body, observer, haAfter, timeBefore, +1)
         timeBefore = evtBefore.time
         timeAfter = evtAfter.time
         altBefore = context.eval(timeBefore)
@@ -6200,7 +6199,6 @@ private fun backwardSearchAltitude(
 
         // If we didn't find the desired event, find the next hour angle bracket and try again.
         val evtAfter  = searchHourAngle(body, observer, haAfter, timeBefore, -1)
-
         if (evtAfter.time.ut <= startTime.ut + limitDays)
             return null
 
diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py
index ac9561be..b304fba9 100644
--- a/generate/template/astronomy.py
+++ b/generate/template/astronomy.py
@@ -4290,7 +4290,7 @@ class HourAngleEvent:
     def __repr__(self):
         return 'HourAngleEvent({}, {})'.format(repr(self.time), repr(self.hor))
 
-def SearchHourAngle(body, observer, hourAngle, startTime):
+def SearchHourAngle(body, observer, hourAngle, startTime, direction = +1):
     """Searches for the time when a celestial body reaches a specified hour angle as seen by an observer on the Earth.
 
     The *hour angle* of a celestial body indicates its position in the sky with respect
@@ -4325,6 +4325,10 @@ def SearchHourAngle(body, observer, hourAngle, startTime):
          body's most recent culmination.
     startTime : Time
          The date and time at which to start the search.
+    direction : int
+        The direction in time to perform the search: a positive value
+        searches forward in time, a negative value searches backward in time.
+        The function throws an exception if `direction` is zero.
 
     Returns
     -------
@@ -4336,6 +4340,9 @@ def SearchHourAngle(body, observer, hourAngle, startTime):
     if hourAngle < 0.0 or hourAngle >= 24.0:
         raise Error('Invalid hour angle.')
 
+    if direction == 0:
+        raise Error('Direction must be positive or negative.')
+
     iter_count = 0
     time = startTime
     while True:
@@ -4348,9 +4355,15 @@ def SearchHourAngle(body, observer, hourAngle, startTime):
         # the hour angle to the desired value.
         delta_sidereal_hours = math.fmod(((hourAngle + ofdate.ra - observer.longitude/15) - gast), 24.0)
         if iter_count == 1:
-            # On the first iteration, always search forward in time.
-            if delta_sidereal_hours < 0.0:
-                delta_sidereal_hours += 24.0
+            # On the first iteration, always search in the requested time direction.
+            if direction > 0:
+                # Search forward in time.
+                if delta_sidereal_hours < 0.0:
+                    delta_sidereal_hours += 24.0
+            else:
+                # Search backward in time.
+                if delta_sidereal_hours > 0.0:
+                    delta_sidereal_hours -= 24.0
         else:
             # On subsequent iterations, we make the smallest possible adjustment,
             # either forward or backward in time.
@@ -4385,8 +4398,10 @@ class Direction(enum.Enum):
     Rise = +1
     Set  = -1
 
+def _ForwardSearchAltitude(body, observer, direction, startTime, limitDays, altitude_error_func, altitude_error_context):
+    if body == Body.Earth:
+        raise EarthNotAllowedError()
 
-def _InternalSearchAltitude(body, observer, direction, startTime, limitDays, altitude_error_func, altitude_error_context):
     if direction == Direction.Rise:
         ha_before = 12.0    # minimum altitude (bottom) happens BEFORE the body rises.
         ha_after  =  0.0    # maximum altitude (culmination) happens AFTER the body rises.
@@ -4396,6 +4411,10 @@ def _InternalSearchAltitude(body, observer, direction, startTime, limitDays, alt
     else:
         raise Error('Invalid value for direction parameter')
 
+    # We cannot possibly satisfy a forward search without a positive time limit.
+    if limitDays <= 0.0:
+        return None
+
     # See if the body is currently above/below the horizon.
     # If we are looking for next rise time and the body is below the horizon,
     # we use the current time as the lower time bound and the next culmination
@@ -4407,7 +4426,7 @@ def _InternalSearchAltitude(body, observer, direction, startTime, limitDays, alt
     alt_before = altitude_error_func(altitude_error_context, time_start)
     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, time_start, +1)
         time_before = evt_before.time
         alt_before = altitude_error_func(altitude_error_context, time_before)
     else:
@@ -4415,7 +4434,7 @@ def _InternalSearchAltitude(body, observer, direction, startTime, limitDays, alt
         # and use the current time as the "before" event.
         time_before = time_start
 
-    evt_after = SearchHourAngle(body, observer, ha_after, time_before)
+    evt_after = SearchHourAngle(body, observer, ha_after, time_before, +1)
     alt_after = altitude_error_func(altitude_error_context, evt_after.time)
 
     while True:
@@ -4423,16 +4442,76 @@ def _InternalSearchAltitude(body, observer, direction, startTime, limitDays, alt
             # Search between the "before time" and the "after time" for the desired event.
             event_time = Search(altitude_error_func, altitude_error_context, time_before, evt_after.time, 1.0)
             if event_time is not None:
+                # If we found the rise/set time, but it falls outside limitDays, fail the search.
+                if event_time.ut > startTime.ut + limitDays:
+                    return None
+                # The search succeeded.
                 return event_time
         # We didn't find the desired event, so use the "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)
+        evt_before = SearchHourAngle(body, observer, ha_before, evt_after.time, +1)
         if evt_before.time.ut >= time_start.ut + limitDays:
             return None
+        evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time, +1)
         time_before = evt_before.time
         alt_before = altitude_error_func(altitude_error_context, evt_before.time)
         alt_after = altitude_error_func(altitude_error_context, evt_after.time)
 
+def _BackwardSearchAltitude(body, observer, direction, startTime, limitDays, altitude_error_func, altitude_error_context):
+    if body == Body.Earth:
+        raise EarthNotAllowedError()
+
+    if direction == Direction.Rise:
+        ha_before = 12.0    # minimum altitude (bottom) happens BEFORE the body rises.
+        ha_after  =  0.0    # maximum altitude (culmination) happens AFTER the body rises.
+    elif direction == Direction.Set:
+        ha_before =  0.0    # culmination happens BEFORE the body sets.
+        ha_after  = 12.0    # bottom happens AFTER the body sets.
+    else:
+        raise Error('Invalid value for direction parameter')
+
+    # We cannot possibly satisfy a backward search without a negative time limit.
+    if limitDays >= 0.0:
+        return None
+
+    # See if the body is currently above/below the horizon.
+    # If we are looking for previous rise time and the body is above the horizon,
+    # we use the current time as the upper time bound and the previous bottom as the lower time bound.
+    # If the body is below the horizon, we search for the previous culmination and use it
+    # as the upper time bound. Then we search for the bottom before that culmination and
+    # use it as the lower time bound.
+    # The same logic applies for finding set times; altitude_error_func and
+    # altitude_error_context ensure that the desired event is represented
+    # by ascending through zero, so the Search function works correctly.
+    time_start = startTime
+    alt_after = altitude_error_func(altitude_error_context, time_start)
+    if alt_after < 0.0:
+        evt_after = SearchHourAngle(body, observer, ha_after, time_start, -1)
+        time_after = evt_after.time
+        alt_after = altitude_error_func(altitude_error_context, time_after)
+    else:
+        time_after = time_start
+
+    evt_before = SearchHourAngle(body, observer, ha_before, time_after, -1)
+    alt_before = altitude_error_func(altitude_error_context, evt_before.time)
+
+    while True:
+        if alt_before <= 0.0 and alt_after > 0.0:
+            # Search between the "before time" and the "after time" for the desired event.
+            event_time = Search(altitude_error_func, altitude_error_context, evt_before.time, time_after, 1.0)
+            if event_time is not None:
+                # If we found the rise/set time, but it falls outside limitDays, fail the search.
+                if event_time.ut < startTime.ut + limitDays:
+                    return None
+                # The search succeeded.
+                return event_time
+        evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time, -1)
+        if evt_after.time.ut <= time_start.ut + limitDays:
+            return None
+        evt_before = SearchHourAngle(body, observer, ha_before, evt_after.time, -1)
+        time_after = evt_after.time
+        alt_before = altitude_error_func(altitude_error_context, evt_before.time)
+        alt_after = altitude_error_func(altitude_error_context, evt_after.time)
+
 
 class _peak_altitude_context:
     def __init__(self, body, direction, observer, body_radius_au):
@@ -4491,11 +4570,14 @@ def SearchRiseSet(body, observer, direction, startTime, limitDays):
     startTime : Time
         The date and time at which to start the search.
     limitDays : float
-        Limits how many days to search for a rise or set time.
+        Limits how many days to search for a rise or set time, and defines
+        the direction in time to search. When `limitDays` is positive, the
+        search is performed into the future, after `startTime`.
+        When negative, the search is performed into the past, before `startTime`.
         To limit a rise or set time to the same day, you can use a value of 1 day.
         In cases where you want to find the next rise or set time no matter how far
-        in the future (for example, for an observer near the south pole), you can pass
-        in a larger value like 365.
+        in the future (for example, for an observer near the south pole), you can
+        pass in a larger value like 365.
 
     Returns
     -------
@@ -4514,7 +4596,9 @@ def SearchRiseSet(body, observer, direction, startTime, limitDays):
         body_radius = 0.0
 
     context = _peak_altitude_context(body, direction, observer, body_radius)
-    return _InternalSearchAltitude(body, observer, direction, startTime, limitDays, _peak_altitude, context)
+    if limitDays < 0.0:
+        return _BackwardSearchAltitude(body, observer, direction, startTime, limitDays, _peak_altitude, context)
+    return _ForwardSearchAltitude(body, observer, direction, startTime, limitDays, _peak_altitude, context)
 
 
 class _altitude_error_context:
@@ -4561,8 +4645,14 @@ def SearchAltitude(body, observer, direction, dateStart, limitDays, altitude):
     startTime : Time
         The date and time at which to start the search.
     limitDays : float
-        The fractional number of days after `dateStart` that limits
-        when the altitude event is to be found. Must be a positive number.
+        Limits how many days to search for the body reaching the altitude angle,
+        and defines the direction in time to search. When `limitDays` is positive, the
+        search is performed into the future, after `startTime`.
+        When negative, the search is performed into the past, before `startTime`.
+        To limit the search to the same day, you can use a value of 1 day.
+        In cases where you want to find the altitude event no matter how far
+        in the future (for example, for an observer near the south pole), you can
+        pass in a larger value like 365.
     altitude : float
         The desired altitude angle of the body's center above (positive)
         or below (negative) the observer's local horizon, expressed in degrees.
@@ -4576,9 +4666,10 @@ def SearchAltitude(body, observer, direction, dateStart, limitDays, altitude):
     """
     if not (-90.0 <= altitude <= +90.0):
         raise Error('Invalid altitude: {}'.format(altitude))
-
     context = _altitude_error_context(body, direction, observer, altitude)
-    return _InternalSearchAltitude(body, observer, direction, dateStart, limitDays, _altitude_error_func, context)
+    if limitDays < 0.0:
+        return _BackwardSearchAltitude(body, observer, direction, dateStart, limitDays, _altitude_error_func, context)
+    return _ForwardSearchAltitude(body, observer, direction, dateStart, limitDays, _altitude_error_func, context)
 
 class SeasonInfo:
     """The dates and times of changes of season for a given calendar year.
diff --git a/generate/template/astronomy.ts b/generate/template/astronomy.ts
index 927dd15a..ca230329 100644
--- a/generate/template/astronomy.ts
+++ b/generate/template/astronomy.ts
@@ -4474,9 +4474,9 @@ function ForwardSearchAltitude(
 
         // If we didn't find the desired event, use time_after to find the next before-event.
         evt_before = SearchHourAngle(body, observer, ha_before, evt_after.time, +1);
-        evt_after  = SearchHourAngle(body, observer, ha_after, evt_before.time, +1);
         if (evt_before.time.ut >= time_start.ut + limitDays)
             return null;
+        evt_after  = SearchHourAngle(body, observer, ha_after, evt_before.time, +1);
 
         time_before = evt_before.time;
         error_before = altitude_error(evt_before.time);
@@ -4556,7 +4556,6 @@ function BackwardSearchAltitude(
         }
 
         evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time, -1);
-
         if (evt_after.time.ut <= time_start.ut + limitDays)
             return null;
 
diff --git a/generate/test.py b/generate/test.py
index 12f3c961..66401261 100755
--- a/generate/test.py
+++ b/generate/test.py
@@ -648,6 +648,61 @@ def Magnitude():
 
 #-----------------------------------------------------------------------------------------------------------
 
+def ToggleDir(dir):
+    return astronomy.Direction(-dir.value)
+
+def RiseSetReverse():
+    nsamples = 5000
+    nudge = 0.1
+    utList = []
+    observer = astronomy.Observer(30.5, -90.7, 0.0)
+    dtMin = +1000.0
+    dtMax = -1000.0
+    maxDiff = 0.0
+
+    # Find alternating sunrise/sunset events in forward chronological order.
+    dir = astronomy.Direction.Rise
+    time = astronomy.Time.Make(2022, 1, 1, 0, 0, 0)
+    for i in range(nsamples):
+        result = astronomy.SearchRiseSet(astronomy.Body.Sun, observer, dir, time, +1.0)
+        if not result:
+            print('PY RiseSetReverse: cannot find {} event after {}'.format(dir, time))
+            return 1
+        utList.append(result.ut)
+        if i > 0:
+            # Check the time between consecutive sunrise/sunset events.
+            # These will vary considerably with the seasons, so just make sure we don't miss any entirely.
+            dt = v(utList[i] - utList[i-1])
+            dtMin = min(dtMin, dt)
+            dtMax = max(dtMax, dt)
+        dir = ToggleDir(dir)
+        time = result.AddDays(+nudge)
+
+    print('PY RiseSetReverse: dtMin={:0.6f} days, dtMax={:0.6f} days.'.format(dtMin, dtMax))
+    if (dtMin < 0.411) or (dtMax > 0.589):
+        print('PY RiseSetReverse: Invalid intervals between sunrise/sunset.')
+        return 1
+
+    # Perform the same search in reverse. Verify we get consistent rise/set times.
+    for i in range(nsamples-1, -1, -1):
+        dir = ToggleDir(dir)
+        result = astronomy.SearchRiseSet(astronomy.Body.Sun, observer, dir, time, -1.0)
+        if not result:
+            print('PY RiseSetReverse: cannot find {] event before {}.'.format(dir, time))
+            return 1
+        diff = SECONDS_PER_DAY * vabs(utList[i] - result.ut)
+        maxDiff = max(maxDiff, diff)
+        time = result.AddDays(-nudge)
+
+    if maxDiff > 0.982:
+        print('PY RiseSetReverse: EXCESSIVE discrepancy = {:0.6f} seconds.'.format(maxDiff))
+        return 1
+
+    print('PY RiseSetReverse: PASS - max diff = {:0.6f} seconds.'.format(maxDiff))
+    return 0
+
+#-----------------------------------------------------------------------------------------------------------
+
 def RiseSet(filename = 'riseset/riseset.txt'):
     sum_minutes = 0.0
     max_minutes = 0.0
@@ -2747,6 +2802,7 @@ UnitTests = {
     'refraction':               Refraction,
     'repr':                     Repr,
     'riseset':                  RiseSet,
+    'riseset_reverse':          RiseSetReverse,
     'rotation':                 Rotation,
     'seasons':                  Seasons,
     'seasons187':               SeasonsIssue187,
diff --git a/source/csharp/astronomy.cs b/source/csharp/astronomy.cs
index 8752a873..88568186 100644
--- a/source/csharp/astronomy.cs
+++ b/source/csharp/astronomy.cs
@@ -6997,11 +6997,12 @@ namespace CosineKitty
 
                 // 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, +1);
-                evt_after  = SearchHourAngle(body, observer, ha_after, evt_before.time, +1);
 
                 if (evt_before.time.ut >= startTime.ut + limitDays)
                     return null;
 
+                evt_after  = SearchHourAngle(body, observer, ha_after, evt_before.time, +1);
+
                 time_before = evt_before.time;
 
                 alt_before = context.Eval(evt_before.time);
diff --git a/source/js/astronomy.browser.js b/source/js/astronomy.browser.js
index f489a034..39b85b47 100644
--- a/source/js/astronomy.browser.js
+++ b/source/js/astronomy.browser.js
@@ -4930,9 +4930,9 @@ function ForwardSearchAltitude(body, observer, direction, dateStart, limitDays,
         }
         // If we didn't find the desired event, use time_after to find the next before-event.
         evt_before = SearchHourAngle(body, observer, ha_before, evt_after.time, +1);
-        evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time, +1);
         if (evt_before.time.ut >= time_start.ut + limitDays)
             return null;
+        evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time, +1);
         time_before = evt_before.time;
         error_before = altitude_error(evt_before.time);
         error_after = altitude_error(evt_after.time);
diff --git a/source/js/astronomy.browser.min.js b/source/js/astronomy.browser.min.js
index e3a89706..e186cc4b 100644
--- a/source/js/astronomy.browser.min.js
+++ b/source/js/astronomy.browser.min.js
@@ -94,7 +94,7 @@ Math.cos(l)-Math.cos(m)*Math.sin(l)*Math.sin(e.DEG2RAD*a.elon-e.DEG2RAD*(169.51+
 a;}var r=d/100;l=a+r*(l+r*(m+r*n))+5*Math.log10(k*g)}return new pc(c,l,d,k,g,f,b,h)}function Ua(a){if(a===p.Earth)throw"The Earth does not have a synodic period as seen from itself.";if(a===p.Moon)return 29.530588;var b=ia[a];if(!b)throw"Not a valid planet name: "+a;a=ia.Earth.OrbitalPeriod;return Math.abs(a/(a/b.OrbitalPeriod-1))}function Ca(a,b,c){function d(m){var n=qa(a,m);m=qa(p.Earth,m);return za(g*(m-n)-b)}x(b);var f=ia[a];if(!f)throw"Cannot search relative longitude because body is not a planet: "+
 a;if(a===p.Earth)throw"Cannot search relative longitude for the Earth (it is always 0)";var g=f.OrbitalPeriod>ia.Earth.OrbitalPeriod?1:-1;f=Ua(a);c=v(c);var k=d(c);0<k&&(k-=360);for(var l=0;100>l;++l){var h=-k/360*f;c=c.AddDays(h);if(1>86400*Math.abs(h))return c;h=k;k=d(c);30>Math.abs(h)&&h!==k&&(h/=h-k,.5<h&&2>h&&(f*=h))}throw"Relative longitude search failed to converge for "+a+" near "+c.toString()+" (error_angle = "+k+").";}function Ib(a){return Hb(p.Moon,p.Sun,a)}function Va(a,b,c){function d(l){l=
 Ib(l);return za(l-a)}x(a);x(c);b=v(b);var f=d(b);if(0>c){0>f&&(f+=360);var g=-(29.530588*f)/360;f=g+1.5;if(f<c)return null;var k=Math.max(c,g-1.5)}else{0<f&&(f-=360);g=-(29.530588*f)/360;k=g-1.5;if(k>c)return null;f=Math.min(c,g+1.5)}c=b.AddDays(k);b=b.AddDays(f);return I(d,c,b,{dt_tolerance_seconds:.1})}function qc(a){var b=Ib(a);b=(Math.floor(b/90)+1)%4;a=Va(90*b,a,10);if(!a)throw"Cannot find moon quarter";return new rc(b,a)}function sc(a,b,c,d,f,g){ka(b);x(f);if(a===p.Earth)throw"Cannot find altitude event for the Earth.";
-if(1===c){c=12;var k=0}else if(-1===c)c=0,k=12;else throw"Invalid direction parameter "+c+" -- must be +1 or -1";if(0>=f)return null;d=v(d);var l=g(d);var h;if(0<l){l=ba(a,b,c,d,1);var m=l.time;l=g(m)}else m=d;var n=ba(a,b,k,m,1);for(h=g(n.time);;){if(0>=l&&0<h&&(m=I(g,m,n.time,{init_f1:l,init_f2:h})))return m.ut>d.ut+f?null:m;l=ba(a,b,c,n.time,1);n=ba(a,b,k,l.time,1);if(l.time.ut>=d.ut+f)return null;m=l.time;l=g(l.time);h=g(n.time)}}function tc(a,b,c,d,f,g){ka(b);x(f);if(a===p.Earth)throw"Cannot find altitude event for the Earth.";
+if(1===c){c=12;var k=0}else if(-1===c)c=0,k=12;else throw"Invalid direction parameter "+c+" -- must be +1 or -1";if(0>=f)return null;d=v(d);var l=g(d);var h;if(0<l){l=ba(a,b,c,d,1);var m=l.time;l=g(m)}else m=d;var n=ba(a,b,k,m,1);for(h=g(n.time);;){if(0>=l&&0<h&&(m=I(g,m,n.time,{init_f1:l,init_f2:h})))return m.ut>d.ut+f?null:m;l=ba(a,b,c,n.time,1);if(l.time.ut>=d.ut+f)return null;n=ba(a,b,k,l.time,1);m=l.time;l=g(l.time);h=g(n.time)}}function tc(a,b,c,d,f,g){ka(b);x(f);if(a===p.Earth)throw"Cannot find altitude event for the Earth.";
 if(1===c){c=12;var k=0}else if(-1===c)c=0,k=12;else throw"Invalid direction parameter "+c+" -- must be +1 or -1";if(0<=f)return null;d=v(d);var l=g(d);var h;if(0>l){l=ba(a,b,k,d,-1);var m=l.time;l=g(m)}else m=d;var n=ba(a,b,c,m,-1);for(h=g(n.time);;){if(0>=h&&0<l&&(m=I(g,n.time,m,{init_f1:h,init_f2:l})))return m.ut<d.ut+f?null:m;l=ba(a,b,k,n.time,-1);if(l.time.ut<=d.ut+f)return null;n=ba(a,b,c,l.time,-1);m=n.time;h=g(n.time);l=g(l.time)}}function ba(a,b,c,d,f){f=void 0===f?1:f;ka(b);d=v(d);var g=
 0;if(a===p.Earth)throw"Cannot search for hour angle of the Earth.";x(c);if(0>c||24<=c)throw"Invalid hour angle "+c;x(f);if(0===f)throw"Direction must be positive or negative.";for(;;){++g;var k=Z(d),l=Na(a,d,b,!0,!0);k=(c+l.ra-b.longitude/15-k)%24;1===g?0<f?0>k&&(k+=24):0<k&&(k-=24):-12>k?k+=24:12<k&&(k-=24);if(.1>3600*Math.abs(k))return a=La(d,b,l.ra,l.dec,"normal"),new uc(d,a);d=d.AddDays(k/24*.9972695717592592)}}function vc(a,b){b=v(b);var c=Hb(a,p.Sun,b);if(180<c){var d="morning";c=360-c}else d=
 "evening";a=Ba(a,b);return new wc(b,d,a,c)}function xc(a){function b(l){var h=l.AddDays(-5E-4);l=l.AddDays(5E-4);h=ca(h).distance_au;return(ca(l).distance_au-h)/.001}function c(l){return-b(l)}a=v(a);for(var d=b(a),f=0;59.061176>5*f;++f){var g=a.AddDays(5),k=b(g);if(0>=d*k){if(0>d||0<k){a=I(b,a,g,{init_f1:d,init_f2:k});if(!a)throw"SearchLunarApsis INTERNAL ERROR: perigee search failed!";d=ca(a).distance_au;return new Wa(a,0,d)}if(0<d||0>k){a=I(c,a,g,{init_f1:-d,init_f2:-k});if(!a)throw"SearchLunarApsis INTERNAL ERROR: apogee search failed!";
diff --git a/source/js/astronomy.js b/source/js/astronomy.js
index 3fdbec6e..0add99fa 100644
--- a/source/js/astronomy.js
+++ b/source/js/astronomy.js
@@ -4929,9 +4929,9 @@ function ForwardSearchAltitude(body, observer, direction, dateStart, limitDays,
         }
         // If we didn't find the desired event, use time_after to find the next before-event.
         evt_before = SearchHourAngle(body, observer, ha_before, evt_after.time, +1);
-        evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time, +1);
         if (evt_before.time.ut >= time_start.ut + limitDays)
             return null;
+        evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time, +1);
         time_before = evt_before.time;
         error_before = altitude_error(evt_before.time);
         error_after = altitude_error(evt_after.time);
diff --git a/source/js/astronomy.min.js b/source/js/astronomy.min.js
index 3ff8de9e..f5c7731b 100644
--- a/source/js/astronomy.min.js
+++ b/source/js/astronomy.min.js
@@ -227,7 +227,7 @@ var MoonQuarter=function(a,b){this.quarter=a;this.time=b};exports.MoonQuarter=Mo
 function BodyRadiusAu(a){switch(a){case Body.Sun:return SUN_RADIUS_AU;case Body.Moon:return MOON_EQUATORIAL_RADIUS_AU;default:return 0}}function SearchRiseSet(a,b,c,d,e){function f(k){var h=Equator(a,k,b,!0,!0);k=Horizon(k,b,h.ra,h.dec).altitude+g/h.dist*exports.RAD2DEG+REFRACTION_NEAR_HORIZON;return c*k}var g=BodyRadiusAu(a);return 0>e?BackwardSearchAltitude(a,b,c,d,e,f):ForwardSearchAltitude(a,b,c,d,e,f)}exports.SearchRiseSet=SearchRiseSet;
 function SearchAltitude(a,b,c,d,e,f){function g(k){var h=Equator(a,k,b,!0,!0);k=Horizon(k,b,h.ra,h.dec);return c*(k.altitude-f)}if(!Number.isFinite(f)||-90>f||90<f)throw"Invalid altitude angle: "+f;return 0>e?BackwardSearchAltitude(a,b,c,d,e,g):ForwardSearchAltitude(a,b,c,d,e,g)}exports.SearchAltitude=SearchAltitude;
 function ForwardSearchAltitude(a,b,c,d,e,f){VerifyObserver(b);VerifyNumber(e);if(a===Body.Earth)throw"Cannot find altitude event for the Earth.";if(1===c){c=12;var g=0}else if(-1===c)c=0,g=12;else throw"Invalid direction parameter "+c+" -- must be +1 or -1";if(0>=e)return null;d=MakeTime(d);var k=f(d);var h;if(0<k){k=SearchHourAngle(a,b,c,d,1);var l=k.time;k=f(l)}else l=d;var m=SearchHourAngle(a,b,g,l,1);for(h=f(m.time);;){if(0>=k&&0<h&&(l=Search(f,l,m.time,{init_f1:k,init_f2:h})))return l.ut>d.ut+
-e?null:l;k=SearchHourAngle(a,b,c,m.time,1);m=SearchHourAngle(a,b,g,k.time,1);if(k.time.ut>=d.ut+e)return null;l=k.time;k=f(k.time);h=f(m.time)}}
+e?null:l;k=SearchHourAngle(a,b,c,m.time,1);if(k.time.ut>=d.ut+e)return null;m=SearchHourAngle(a,b,g,k.time,1);l=k.time;k=f(k.time);h=f(m.time)}}
 function BackwardSearchAltitude(a,b,c,d,e,f){VerifyObserver(b);VerifyNumber(e);if(a===Body.Earth)throw"Cannot find altitude event for the Earth.";if(1===c){c=12;var g=0}else if(-1===c)c=0,g=12;else throw"Invalid direction parameter "+c+" -- must be +1 or -1";if(0<=e)return null;d=MakeTime(d);var k=f(d);var h;if(0>k){k=SearchHourAngle(a,b,g,d,-1);var l=k.time;k=f(l)}else l=d;var m=SearchHourAngle(a,b,c,l,-1);for(h=f(m.time);;){if(0>=h&&0<k&&(l=Search(f,m.time,l,{init_f1:h,init_f2:k})))return l.ut<
 d.ut+e?null:l;k=SearchHourAngle(a,b,g,m.time,-1);if(k.time.ut<=d.ut+e)return null;m=SearchHourAngle(a,b,c,k.time,-1);l=m.time;h=f(m.time);k=f(k.time)}}var HourAngleEvent=function(a,b){this.time=a;this.hor=b};exports.HourAngleEvent=HourAngleEvent;
 function SearchHourAngle(a,b,c,d,e){e=void 0===e?1:e;VerifyObserver(b);d=MakeTime(d);var f=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;VerifyNumber(e);if(0===e)throw"Direction must be positive or negative.";for(;;){++f;var g=sidereal_time(d),k=Equator(a,d,b,!0,!0);g=(c+k.ra-b.longitude/15-g)%24;1===f?0<e?0>g&&(g+=24):0<g&&(g-=24):-12>g?g+=24:12<g&&(g-=24);if(.1>3600*Math.abs(g))return a=Horizon(d,b,k.ra,k.dec,"normal"),
diff --git a/source/js/astronomy.ts b/source/js/astronomy.ts
index 1b4d1035..3a9f22ae 100644
--- a/source/js/astronomy.ts
+++ b/source/js/astronomy.ts
@@ -5480,9 +5480,9 @@ function ForwardSearchAltitude(
 
         // If we didn't find the desired event, use time_after to find the next before-event.
         evt_before = SearchHourAngle(body, observer, ha_before, evt_after.time, +1);
-        evt_after  = SearchHourAngle(body, observer, ha_after, evt_before.time, +1);
         if (evt_before.time.ut >= time_start.ut + limitDays)
             return null;
+        evt_after  = SearchHourAngle(body, observer, ha_after, evt_before.time, +1);
 
         time_before = evt_before.time;
         error_before = altitude_error(evt_before.time);
@@ -5562,7 +5562,6 @@ function BackwardSearchAltitude(
         }
 
         evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time, -1);
-
         if (evt_after.time.ut <= time_start.ut + limitDays)
             return null;
 
diff --git a/source/js/esm/astronomy.js b/source/js/esm/astronomy.js
index 96639739..7004c4c7 100644
--- a/source/js/esm/astronomy.js
+++ b/source/js/esm/astronomy.js
@@ -4869,9 +4869,9 @@ function ForwardSearchAltitude(body, observer, direction, dateStart, limitDays,
         }
         // If we didn't find the desired event, use time_after to find the next before-event.
         evt_before = SearchHourAngle(body, observer, ha_before, evt_after.time, +1);
-        evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time, +1);
         if (evt_before.time.ut >= time_start.ut + limitDays)
             return null;
+        evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time, +1);
         time_before = evt_before.time;
         error_before = altitude_error(evt_before.time);
         error_after = altitude_error(evt_after.time);
diff --git a/source/kotlin/src/main/kotlin/io/github/cosinekitty/astronomy/astronomy.kt b/source/kotlin/src/main/kotlin/io/github/cosinekitty/astronomy/astronomy.kt
index e50c3408..eb93a469 100644
--- a/source/kotlin/src/main/kotlin/io/github/cosinekitty/astronomy/astronomy.kt
+++ b/source/kotlin/src/main/kotlin/io/github/cosinekitty/astronomy/astronomy.kt
@@ -6121,11 +6121,10 @@ private fun forwardSearchAltitude(
 
         // If we didn't find the desired event, find the next hour angle bracket and try again.
         val evtBefore = searchHourAngle(body, observer, haBefore, timeAfter, +1)
-        val evtAfter  = searchHourAngle(body, observer, haAfter, timeBefore, +1)
-
         if (evtBefore.time.ut >= startTime.ut + limitDays)
             return null
 
+        val evtAfter  = searchHourAngle(body, observer, haAfter, timeBefore, +1)
         timeBefore = evtBefore.time
         timeAfter = evtAfter.time
         altBefore = context.eval(timeBefore)
@@ -6200,7 +6199,6 @@ private fun backwardSearchAltitude(
 
         // If we didn't find the desired event, find the next hour angle bracket and try again.
         val evtAfter  = searchHourAngle(body, observer, haAfter, timeBefore, -1)
-
         if (evtAfter.time.ut <= startTime.ut + limitDays)
             return null
 
diff --git a/source/python/README.md b/source/python/README.md
index 569cef8d..ea6477f2 100644
--- a/source/python/README.md
+++ b/source/python/README.md
@@ -2757,7 +2757,7 @@ Astronomical twilight uses -18 degrees as the `altitude` value.
 | [`Observer`](#Observer) | `observer` | The location where observation takes place. |
 | [`Direction`](#Direction) | `direction` | Either `Direction.Rise` to find an ascending altitude event or `Direction.Set` to find a descending altitude event. |
 | [`Time`](#Time) | `startTime` | The date and time at which to start the search. |
-| `float` | `limitDays` | The fractional number of days after `dateStart` that limits when the altitude event is to be found. Must be a positive number. |
+| `float` | `limitDays` | Limits how many days to search for the body reaching the altitude angle, and defines the direction in time to search. When `limitDays` is positive, the search is performed into the future, after `startTime`. When negative, the search is performed into the past, before `startTime`. To limit the search to the same day, you can use a value of 1 day. In cases where you want to find the altitude event no matter how far in the future (for example, for an observer near the south pole), you can pass in a larger value like 365. |
 | `float` | `altitude` | The desired altitude angle of the body's center above (positive) or below (negative) the observer's local horizon, expressed in degrees. Must be in the range [-90, +90]. |
 
 **Returns**: [`Time`](#Time) or `None`
@@ -2787,7 +2787,7 @@ passing in the `peak` value returned from the previous call.
 ---
 
 <a name="SearchHourAngle"></a>
-### SearchHourAngle(body, observer, hourAngle, startTime)
+### SearchHourAngle(body, observer, hourAngle, startTime, direction=1)
 
 **Searches for the time when a celestial body reaches a specified hour angle as seen by an observer on the Earth.**
 
@@ -2815,6 +2815,7 @@ of the body at that time, as seen by the given observer.
 | [`Observer`](#Observer) | `observer` | Indicates a location on or near the surface of the Earth where the observer is located. |
 | `float` | `hourAngle` | An hour angle value in the range [0.0, 24.0) indicating the number of sidereal hours after the body's most recent culmination. |
 | [`Time`](#Time) | `startTime` | The date and time at which to start the search. |
+| `int` | `direction` | The direction in time to perform the search: a positive value searches forward in time, a negative value searches backward in time. The function throws an exception if `direction` is zero. |
 
 **Returns**: [`HourAngleEvent`](#HourAngleEvent)
 
@@ -3097,7 +3098,7 @@ Therefore callers must not assume that the function will always succeed.
 | [`Observer`](#Observer) | `observer` | The location where observation takes place. |
 | [`Direction`](#Direction) | `direction` | Either `Direction.Rise` to find a rise time or `Direction.Set` to find a set time. |
 | [`Time`](#Time) | `startTime` | The date and time at which to start the search. |
-| `float` | `limitDays` | Limits how many days to search for a rise or set time. To limit a rise or set time to the same day, you can use a value of 1 day. In cases where you want to find the next rise or set time no matter how far in the future (for example, for an observer near the south pole), you can pass in a larger value like 365. |
+| `float` | `limitDays` | Limits how many days to search for a rise or set time, and defines the direction in time to search. When `limitDays` is positive, the search is performed into the future, after `startTime`. When negative, the search is performed into the past, before `startTime`. To limit a rise or set time to the same day, you can use a value of 1 day. In cases where you want to find the next rise or set time no matter how far in the future (for example, for an observer near the south pole), you can pass in a larger value like 365. |
 
 **Returns**: [`Time`](#Time) or `None`
 If the rise or set time is found within the specified time window,
diff --git a/source/python/astronomy/astronomy.py b/source/python/astronomy/astronomy.py
index 22511c94..9a74b0cb 100644
--- a/source/python/astronomy/astronomy.py
+++ b/source/python/astronomy/astronomy.py
@@ -6332,7 +6332,7 @@ class HourAngleEvent:
     def __repr__(self):
         return 'HourAngleEvent({}, {})'.format(repr(self.time), repr(self.hor))
 
-def SearchHourAngle(body, observer, hourAngle, startTime):
+def SearchHourAngle(body, observer, hourAngle, startTime, direction = +1):
     """Searches for the time when a celestial body reaches a specified hour angle as seen by an observer on the Earth.
 
     The *hour angle* of a celestial body indicates its position in the sky with respect
@@ -6367,6 +6367,10 @@ def SearchHourAngle(body, observer, hourAngle, startTime):
          body's most recent culmination.
     startTime : Time
          The date and time at which to start the search.
+    direction : int
+        The direction in time to perform the search: a positive value
+        searches forward in time, a negative value searches backward in time.
+        The function throws an exception if `direction` is zero.
 
     Returns
     -------
@@ -6378,6 +6382,9 @@ def SearchHourAngle(body, observer, hourAngle, startTime):
     if hourAngle < 0.0 or hourAngle >= 24.0:
         raise Error('Invalid hour angle.')
 
+    if direction == 0:
+        raise Error('Direction must be positive or negative.')
+
     iter_count = 0
     time = startTime
     while True:
@@ -6390,9 +6397,15 @@ def SearchHourAngle(body, observer, hourAngle, startTime):
         # the hour angle to the desired value.
         delta_sidereal_hours = math.fmod(((hourAngle + ofdate.ra - observer.longitude/15) - gast), 24.0)
         if iter_count == 1:
-            # On the first iteration, always search forward in time.
-            if delta_sidereal_hours < 0.0:
-                delta_sidereal_hours += 24.0
+            # On the first iteration, always search in the requested time direction.
+            if direction > 0:
+                # Search forward in time.
+                if delta_sidereal_hours < 0.0:
+                    delta_sidereal_hours += 24.0
+            else:
+                # Search backward in time.
+                if delta_sidereal_hours > 0.0:
+                    delta_sidereal_hours -= 24.0
         else:
             # On subsequent iterations, we make the smallest possible adjustment,
             # either forward or backward in time.
@@ -6427,8 +6440,10 @@ class Direction(enum.Enum):
     Rise = +1
     Set  = -1
 
+def _ForwardSearchAltitude(body, observer, direction, startTime, limitDays, altitude_error_func, altitude_error_context):
+    if body == Body.Earth:
+        raise EarthNotAllowedError()
 
-def _InternalSearchAltitude(body, observer, direction, startTime, limitDays, altitude_error_func, altitude_error_context):
     if direction == Direction.Rise:
         ha_before = 12.0    # minimum altitude (bottom) happens BEFORE the body rises.
         ha_after  =  0.0    # maximum altitude (culmination) happens AFTER the body rises.
@@ -6438,6 +6453,10 @@ def _InternalSearchAltitude(body, observer, direction, startTime, limitDays, alt
     else:
         raise Error('Invalid value for direction parameter')
 
+    # We cannot possibly satisfy a forward search without a positive time limit.
+    if limitDays <= 0.0:
+        return None
+
     # See if the body is currently above/below the horizon.
     # If we are looking for next rise time and the body is below the horizon,
     # we use the current time as the lower time bound and the next culmination
@@ -6449,7 +6468,7 @@ def _InternalSearchAltitude(body, observer, direction, startTime, limitDays, alt
     alt_before = altitude_error_func(altitude_error_context, time_start)
     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, time_start, +1)
         time_before = evt_before.time
         alt_before = altitude_error_func(altitude_error_context, time_before)
     else:
@@ -6457,7 +6476,7 @@ def _InternalSearchAltitude(body, observer, direction, startTime, limitDays, alt
         # and use the current time as the "before" event.
         time_before = time_start
 
-    evt_after = SearchHourAngle(body, observer, ha_after, time_before)
+    evt_after = SearchHourAngle(body, observer, ha_after, time_before, +1)
     alt_after = altitude_error_func(altitude_error_context, evt_after.time)
 
     while True:
@@ -6465,16 +6484,76 @@ def _InternalSearchAltitude(body, observer, direction, startTime, limitDays, alt
             # Search between the "before time" and the "after time" for the desired event.
             event_time = Search(altitude_error_func, altitude_error_context, time_before, evt_after.time, 1.0)
             if event_time is not None:
+                # If we found the rise/set time, but it falls outside limitDays, fail the search.
+                if event_time.ut > startTime.ut + limitDays:
+                    return None
+                # The search succeeded.
                 return event_time
         # We didn't find the desired event, so use the "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)
+        evt_before = SearchHourAngle(body, observer, ha_before, evt_after.time, +1)
         if evt_before.time.ut >= time_start.ut + limitDays:
             return None
+        evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time, +1)
         time_before = evt_before.time
         alt_before = altitude_error_func(altitude_error_context, evt_before.time)
         alt_after = altitude_error_func(altitude_error_context, evt_after.time)
 
+def _BackwardSearchAltitude(body, observer, direction, startTime, limitDays, altitude_error_func, altitude_error_context):
+    if body == Body.Earth:
+        raise EarthNotAllowedError()
+
+    if direction == Direction.Rise:
+        ha_before = 12.0    # minimum altitude (bottom) happens BEFORE the body rises.
+        ha_after  =  0.0    # maximum altitude (culmination) happens AFTER the body rises.
+    elif direction == Direction.Set:
+        ha_before =  0.0    # culmination happens BEFORE the body sets.
+        ha_after  = 12.0    # bottom happens AFTER the body sets.
+    else:
+        raise Error('Invalid value for direction parameter')
+
+    # We cannot possibly satisfy a backward search without a negative time limit.
+    if limitDays >= 0.0:
+        return None
+
+    # See if the body is currently above/below the horizon.
+    # If we are looking for previous rise time and the body is above the horizon,
+    # we use the current time as the upper time bound and the previous bottom as the lower time bound.
+    # If the body is below the horizon, we search for the previous culmination and use it
+    # as the upper time bound. Then we search for the bottom before that culmination and
+    # use it as the lower time bound.
+    # The same logic applies for finding set times; altitude_error_func and
+    # altitude_error_context ensure that the desired event is represented
+    # by ascending through zero, so the Search function works correctly.
+    time_start = startTime
+    alt_after = altitude_error_func(altitude_error_context, time_start)
+    if alt_after < 0.0:
+        evt_after = SearchHourAngle(body, observer, ha_after, time_start, -1)
+        time_after = evt_after.time
+        alt_after = altitude_error_func(altitude_error_context, time_after)
+    else:
+        time_after = time_start
+
+    evt_before = SearchHourAngle(body, observer, ha_before, time_after, -1)
+    alt_before = altitude_error_func(altitude_error_context, evt_before.time)
+
+    while True:
+        if alt_before <= 0.0 and alt_after > 0.0:
+            # Search between the "before time" and the "after time" for the desired event.
+            event_time = Search(altitude_error_func, altitude_error_context, evt_before.time, time_after, 1.0)
+            if event_time is not None:
+                # If we found the rise/set time, but it falls outside limitDays, fail the search.
+                if event_time.ut < startTime.ut + limitDays:
+                    return None
+                # The search succeeded.
+                return event_time
+        evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time, -1)
+        if evt_after.time.ut <= time_start.ut + limitDays:
+            return None
+        evt_before = SearchHourAngle(body, observer, ha_before, evt_after.time, -1)
+        time_after = evt_after.time
+        alt_before = altitude_error_func(altitude_error_context, evt_before.time)
+        alt_after = altitude_error_func(altitude_error_context, evt_after.time)
+
 
 class _peak_altitude_context:
     def __init__(self, body, direction, observer, body_radius_au):
@@ -6533,11 +6612,14 @@ def SearchRiseSet(body, observer, direction, startTime, limitDays):
     startTime : Time
         The date and time at which to start the search.
     limitDays : float
-        Limits how many days to search for a rise or set time.
+        Limits how many days to search for a rise or set time, and defines
+        the direction in time to search. When `limitDays` is positive, the
+        search is performed into the future, after `startTime`.
+        When negative, the search is performed into the past, before `startTime`.
         To limit a rise or set time to the same day, you can use a value of 1 day.
         In cases where you want to find the next rise or set time no matter how far
-        in the future (for example, for an observer near the south pole), you can pass
-        in a larger value like 365.
+        in the future (for example, for an observer near the south pole), you can
+        pass in a larger value like 365.
 
     Returns
     -------
@@ -6556,7 +6638,9 @@ def SearchRiseSet(body, observer, direction, startTime, limitDays):
         body_radius = 0.0
 
     context = _peak_altitude_context(body, direction, observer, body_radius)
-    return _InternalSearchAltitude(body, observer, direction, startTime, limitDays, _peak_altitude, context)
+    if limitDays < 0.0:
+        return _BackwardSearchAltitude(body, observer, direction, startTime, limitDays, _peak_altitude, context)
+    return _ForwardSearchAltitude(body, observer, direction, startTime, limitDays, _peak_altitude, context)
 
 
 class _altitude_error_context:
@@ -6603,8 +6687,14 @@ def SearchAltitude(body, observer, direction, dateStart, limitDays, altitude):
     startTime : Time
         The date and time at which to start the search.
     limitDays : float
-        The fractional number of days after `dateStart` that limits
-        when the altitude event is to be found. Must be a positive number.
+        Limits how many days to search for the body reaching the altitude angle,
+        and defines the direction in time to search. When `limitDays` is positive, the
+        search is performed into the future, after `startTime`.
+        When negative, the search is performed into the past, before `startTime`.
+        To limit the search to the same day, you can use a value of 1 day.
+        In cases where you want to find the altitude event no matter how far
+        in the future (for example, for an observer near the south pole), you can
+        pass in a larger value like 365.
     altitude : float
         The desired altitude angle of the body's center above (positive)
         or below (negative) the observer's local horizon, expressed in degrees.
@@ -6618,9 +6708,10 @@ def SearchAltitude(body, observer, direction, dateStart, limitDays, altitude):
     """
     if not (-90.0 <= altitude <= +90.0):
         raise Error('Invalid altitude: {}'.format(altitude))
-
     context = _altitude_error_context(body, direction, observer, altitude)
-    return _InternalSearchAltitude(body, observer, direction, dateStart, limitDays, _altitude_error_func, context)
+    if limitDays < 0.0:
+        return _BackwardSearchAltitude(body, observer, direction, dateStart, limitDays, _altitude_error_func, context)
+    return _ForwardSearchAltitude(body, observer, direction, dateStart, limitDays, _altitude_error_func, context)
 
 class SeasonInfo:
     """The dates and times of changes of season for a given calendar year.