Fix #998: detect sunrise/moonrise even at 70°N

First, this makes the main search loop more precise by computing the
actual rate at which the target’s HA is changing, instead of assuming
one rotation per day — which worked okay for the Sun and planets, but
less well for the Moon given its high 13°/day velocity.

Second, this bolts on a final step which, a bit expensively, compares
the final two positions computed in the loop to the actual horizon, and
fits a parabola to the altitude and altitude-rate-of-change to estimate
the exact moment at which the object reaches the horizon.

Author: brandon-rhodes

CHANGELOG.rst

@@ -20,6 +20,13 @@ v1.50 — Not yet released
   have to construct the position themselves: ``SSB.at(t)`` returns a
   position whose coordinates and velocity are both zero in the ICRS.
 
+* The routines added last year, :func:`~skyfield.almanac.find_risings()`
+  and :func:`~skyfield.almanac.find_settings()`, could occasionally miss
+  a sunrise or moonrise at high latitudes like 70°N where the Sun or
+  Moon barely crests the horizon.  This has been fixed at a moderate
+  cost to their runtime.
+  `#998 <https://github.com/skyfielders/python-skyfield/issues/998>`_
+
 * Skyfield no longer tries to protect users by raising an exception if,
   contrary to the usual custom in astronomy, they ask for ``ra.degrees``
   or ``dec.hours``.  So users no longer need to add an underscore prefix

skyfield/almanac.py

@@ -4,13 +4,17 @@
 from __future__ import division
 
 import numpy as np
-from numpy import cos, sin, zeros_like
+from numpy import cos, sin, sqrt, zeros_like
 from .constants import pi, tau
 from .framelib import ecliptic_frame
 from .searchlib import find_discrete
 from .nutationlib import iau2000b_radians
 from .units import Angle
 
+_SUN = 10
+_MOON = 301
+_MICROSECOND = 1 / 24.0 / 3600.0 / 1e6
+
 # Not only to support historic code but also for future convenience, let
 # folks import the search routine alongside the almanac routines.
 find_discrete
@@ -330,29 +334,50 @@ def _rising_hour_angle(latitude, declination, altitude_radians):
 def _transit_ha(latitude, declination, altitude_radians):
     return 0.0
 
+def _q(a, b, c, sign):
+    discriminant = np.maximum(b*b - 4*a*c, 0.0)  # avoid tiny negative results
+    return - 2*c / (b + sign * sqrt(discriminant))
+
+def _intersection(y0, y1, v0, v1):
+    # Return x at which a curve reaches y=0, given its position and
+    # velocity y0,v0 at x=0 and y1,v1 at x=1.  For details, see
+    # `design/intersect_function.py` in the Skyfield repository.
+    sign = 1 - 2 * (y0 > y1)
+    return _q(y1 - y0 - v0, v0, y0, sign)
+
 # Per https://aa.usno.navy.mil/faq/RST_defs we estimate 34 arcminutes of
 # atmospheric refraction and 16 arcminutes for the radius of the Sun.
 _sun_horizon_radians = -50.0 / 21600.0 * tau
 _refraction_radians = -34.0 / 21600.0 * tau
 _moon_radius_m = 1.7374e6
+_clip_lower = -1.0
+_clip_upper = +2.0
+
+def build_horizon_function(target):
+    """Build and return a horizon function `h()` for the given `target`.
+
+    The returned function takes a Distance argument giving the distance
+    from the observer to the target, and returns a negative angle in
+    radians giving the altitude which, when reached by the target's
+    center, places it at the moment of rising.
+
+    """
+    target_id = getattr(target, 'target', None)
+    if target_id == _SUN:
+        def h(distance): return _sun_horizon_radians
+    elif target_id == _MOON:
+        def h(distance):
+            return _refraction_radians - _moon_radius_m / distance.m
+    else:
+        def h(distance): return _refraction_radians
+    return h
 
 def _find(observer, target, start_time, end_time, horizon_degrees, f):
-    # Build a function h() that returns the angle above or below the
-    # horizon we are aiming for, in radians.
     if horizon_degrees is None:
-        tt = getattr(target, 'target', None)
-        if tt == 10:
-            horizon_radians = _sun_horizon_radians
-            h = lambda distance: horizon_radians
-        elif tt == 301:
-            horizon_radians = _refraction_radians
-            h = lambda distance: horizon_radians - _moon_radius_m / distance.m
-        else:
-            horizon_radians = _refraction_radians
-            h = lambda distance: horizon_radians
+        h = build_horizon_function(target)
     else:
         horizon_radians = horizon_degrees / 360.0 * tau
-        h = lambda distance: horizon_radians
+        def h(distance): return horizon_radians
 
     geo = observer.vector_functions[-1]  # should we check observer.center?
     latitude = geo.latitude
@@ -384,6 +409,10 @@ def _find(observer, target, start_time, end_time, horizon_degrees, f):
     # throw the first few out and keep only the last one.
     i, = np.nonzero(np.diff(difference) > 0.0)
 
+    # Trim a few arrays down to just the matching elements.
+    old_ha_radians = ha.radians[i]
+    old_t = t[i]
+
     # When might each rising have actually taken place?  Let's
     # interpolate between the two times that bracket each rising.
     a = difference[i]
@@ -392,24 +421,92 @@ def _find(observer, target, start_time, end_time, horizon_degrees, f):
     interpolated_tt = (b * tt[i] + a * tt[i+1]) / (a + b)
     t = ts.tt_jd(interpolated_tt)
 
-    ha_per_day = tau            # angle the celestrial sphere rotates in 1 day
+    def normalize_zero_to_tau(radians):
+        return radians % tau
 
-    # TODO: How many iterations do we need?  And can we cut down on that
-    # number if we use velocity intelligently?  For now, we experiment
-    # using the ./design/test_sunrise_moonrise.py script in the
-    # repository, that checks both the old Skyfiled routines and this
-    # new one against the USNO.  It suggests that 3 iterations is enough
-    # for the Moon, the fastest-moving Solar System object, to match.
+    def normalize_plus_or_minus_pi(radians):
+        return (radians + pi) % tau - pi
+
+    normalize = normalize_zero_to_tau
+
+    # How did we decide on this many iterations?  We played with the
+    # script ./design/test_sunrise_moonrise.py in the repository.
     for i in 0, 1, 2:
         _fastify(t)
-        ha, dec, distance = observer.at(t).observe(target).apparent().hadec()
+
+        # Expensive: generate true ha/dec at `t`.
+        apparent = observer.at(t).observe(target).apparent()
+        ha, dec, distance = apparent.hadec()
+
+        # Estimate where the horizon-crossing is.
         desired_ha = f(latitude, dec, h(distance))
         ha_adjustment = desired_ha - ha.radians
         ha_adjustment = (ha_adjustment + pi) % tau - pi
+
+        # Figure out how fast the target's HA is changing.
+        ha_diff = normalize(ha.radians - old_ha_radians)
+        t_diff = t - old_t
+        ha_per_day = ha_diff / t_diff
+
+        # Remember this iteration's HA and `t` for the next iteration.
+        old_ha_radians = ha.radians
+        old_t = t
+
+        # The big moment!  Carefully adjust `t` towards intersection.
         timebump = ha_adjustment / ha_per_day
+        timebump[timebump == 0.0] = _MICROSECOND   # avoid divide-by-zero
+        previous_t = t
         t = ts.tt_jd(t.whole, t.tt_fraction + timebump)
 
-    is_above_horizon = (desired_ha % pi != 0.0)
+        # Different `normalize` for all but the first iteration.
+        normalize = normalize_plus_or_minus_pi
+
+    if f is _transit_ha:
+        return t
+
+    # In almost all cases, we are now happy.  But for rising and setting
+    # calculations at high latitudes where the target barely scrapes the
+    # horizon, we might be stuck between two solutions, and need to
+    # interpolate between them.
+
+    # Snag the observer's GeographicPosition and learn the target's
+    # altitude vs the horizon and how fast it's moving vertically at
+    # the second-to-last `t` we computed above.
+    v = observer.vector_functions[-1]
+    altitude0, _, distance0, rate0, _, _ = (
+        apparent.frame_latlon_and_rates(v))
+
+    # Even faster than _fastify(t) is to just assume that nutation
+    # doesn't have much time to move over this short interval.
+    t.M = previous_t.M
+    t._nutation_angles_radians = previous_t._nutation_angles_radians
+
+    # And again, this time with the very final `t` we computed.
+    apparent = observer.at(t).observe(target).apparent()
+    altitude1, _, distance1, rate1, _, _ = (
+        apparent.frame_latlon_and_rates(v))
+
+    # Using the target's altitude and altitude-velocity at the final
+    # two times we computed, compute where it crosses the horizon.
+    tdiff = t - previous_t
+    t_scaled_offset = _intersection(
+        altitude0.radians - h(distance0),
+        altitude1.radians - h(distance1),
+        rate0.radians.per_day * tdiff,
+        rate1.radians.per_day * tdiff,
+    )
+
+    # In case the parabola for some reason goes crazy, don't let our
+    # solution be thrown too far away from our final two times.
+    t_scaled_offset = np.clip(t_scaled_offset, _clip_lower, _clip_upper)
+
+    t = previous_t + t_scaled_offset * tdiff
+
+    is_above_horizon =  (
+        (desired_ha % pi != 0.0)
+        | ((t_scaled_offset > _clip_lower) & (t_scaled_offset < _clip_upper))
+    )
+
     return t, is_above_horizon
 
 def find_risings(observer, target, start_time, end_time, horizon_degrees=None):
@@ -466,5 +563,4 @@ def find_transits(observer, target, start_time, end_time):
     .. versionadded:: 1.47
 
     """
-    t, _ = _find(observer, target, start_time, end_time, 0.0, _transit_ha)
-    return t
+    return _find(observer, target, start_time, end_time, 0.0, _transit_ha)

skyfield/tests/test_almanac.py

@@ -102,6 +102,23 @@ def f(t0, t1, e, topos):
         return t, array([1, 0])
     _sunrise_sunset(f)
 
+def test_new_moonrise_at_70_degrees_north():
+    ts = api.load.timescale()
+    t0 = ts.utc(2023, 2, 19)
+    t1 = ts.utc(2023, 2, 20)
+    e = api.load('de421.bsp')
+    seventy = e['earth'] + api.wgs84.latlon(70, 0)
+
+    t, y = almanac.find_risings(seventy, e['moon'], t0, t1)
+    strings = t.utc_strftime('%Y-%m-%d %H:%M')
+    assert strings == ['2023-02-19 11:28']
+    assert list(y) == [1]
+
+    t, y = almanac.find_settings(seventy, e['moon'], t0, t1)
+    strings = t.utc_strftime('%Y-%m-%d %H:%M')
+    assert strings == ['2023-02-19 12:05']
+    assert list(y) == [1]
+
 def test_dark_twilight_day():
     ts = api.load.timescale()
     t0 = ts.utc(2019, 11, 8, 4)