diff --git a/skyfield/almanac.py b/skyfield/almanac.py index f7fe5959..1c3ed14b 100644 --- a/skyfield/almanac.py +++ b/skyfield/almanac.py @@ -336,7 +336,11 @@ def _transit_ha(latitude, declination, altitude_radians): 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)) + with np.errstate(invalid='ignore', divide='ignore'): + x = - 2*c / (b + sign * sqrt(discriminant)) + # If c is zero, then x=0 is exactly a root, but the quotient above + # can be 0/0 = NaN when the denominator is also zero (see #1114). + return np.where(c == 0.0, 0.0, x) def _intersection(y0, y1, v0, v1): # Return x at which a curve reaches y=0, given its position and diff --git a/skyfield/tests/test_almanac.py b/skyfield/tests/test_almanac.py index aba7725d..5365b40d 100644 --- a/skyfield/tests/test_almanac.py +++ b/skyfield/tests/test_almanac.py @@ -136,6 +136,15 @@ def test_dark_twilight_day(): # Logic. +def test_intersection_when_curve_starts_exactly_at_zero(): + # If the rise/set Newton iteration happens to land, to the last bit, + # exactly on the horizon, then the final parabola fit was dividing + # 0/0 = NaN (#1114). These operands were captured live from a Mars + # setting at latitude 52.0N on 2029-11-04. + v = -1.5990198773309593e-16 + x = almanac._intersection(0.0, 0.0, v, v) + assert x == 0.0 + def test_close_start_and_end(): ts = api.load.timescale() t0 = ts.utc(2018, 9, 23, 1)