Skip to content

Fix 0/0 = NaN in the rise/set parabola fit - #1140

Merged
brandon-rhodes merged 1 commit into
skyfielders:masterfrom
chaunceygardiner:fix-1114-parabola-nan
Aug 4, 2026
Merged

Fix 0/0 = NaN in the rise/set parabola fit#1140
brandon-rhodes merged 1 commit into
skyfielders:masterfrom
chaunceygardiner:fix-1114-parabola-nan

Conversation

@chaunceygardiner

Copy link
Copy Markdown
Contributor

Fixes #1114. The full analysis is in my comment on that issue; the short version:

By the time _find() reaches its final parabola refinement, the Newton iteration has sometimes already converged so tightly that the last two sample times are numerically adjacent instants, and the recomputed altitude offset from the horizon is bit-identical at both. When that offset happens to land on exactly 0.0, _q() divides 0/0: the numerator -2*c is -0.0 because c = y0 = 0.0, and with y0 == y1 the sign is +1, so the denominator b + sqrt(b*b) is 0.0 whenever b < 0. The resulting NaN survives np.clip() and is returned as the event time, while is_above_horizon stays True — NaN comparisons are False — so callers cannot filter the bad event and crash later in utc_iso()/utc_datetime() with ValueError: cannot convert float NaN to integer.

When c is zero, x = 0 is exactly a root of the quadratic — the curve passes through zero at the first sample point, which is where the iteration converged — so _q() now returns it directly. The errstate guard keeps NumPy quiet while the discarded quotient is computed.

The trigger is last-ULP rounding of the whole trig pipeline, which is why the failing (date, latitude, body) tuples in the issue thread differ across platforms and NumPy builds. The new test therefore skips the ephemeris entirely and feeds _intersection() the exact operands captured from a failing case (a Mars setting at 52.0°N on 2029-11-04, on Linux/aarch64 with NumPy 2.4.1): it computes NaN before this change and 0.0 after.

With the fix, both deterministic cases from the thread return times that agree with an independent find_discrete search over risings_and_settings on the same day:

2029-11-04 mars             52.0°N   set 18:45:31Z   (find_discrete: 18:45:31Z)
2029-10-14 saturn           68.4°N   set 11:15:42Z   (find_discrete: 11:15:42Z)

Test suite: the failing-test set under assay --batch skyfield.tests is identical before and after the change (the failures on my machine are pre-existing NumPy 2.5 deprecation warnings escalated by PYTHONWARNINGS=error, unrelated to this fix), and the new test passes. pyflakes is clean.

When the Newton iteration in _find() happens to land, to the last
bit, exactly on the horizon, the final parabola refinement divided
0/0 and returned NaN as the event time.  Fixes skyfielders#1114.
@brandon-rhodes
brandon-rhodes merged commit 771bb29 into skyfielders:master Aug 4, 2026
2 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

RuntimeWarning from almanac.py line 339 (divide by zero?)

2 participants