Skip to content

Fix NaN in rising/setting search near a crossing (#1128) - #1145

Open
karimosama7 wants to merge 1 commit into
skyfielders:masterfrom
karimosama7:fix-1128-almanac-nan
Open

karimosama7 wants to merge 1 commit into
skyfielders:masterfrom
karimosama7:fix-1128-almanac-nan

Conversation

@karimosama7

Copy link
Copy Markdown

This fixes the almanac half of #1128: find_settings() / find_risings() raising a misleading EphemerisRangeError when start_time falls just before a crossing.

Cause. In @ejeschke's reproduction, t0 is itself almost exactly a sunset, so the Sun sets again only ~16 µs later and difference[0] in _find() is just 3.58e-10 radians. The interpolated first guess was built as a single float Julian date, whose resolution near JD 2.46 million is ~40 µs, so it rounded back to exactly tt[i]. That made t_diff zero, ha_per_day = 0/0 = NaN, and the NaN time then reached jplephem. It reproduces with NumPy 2.4.4 as well as 2.5.3 (on Windows), so it isn't NumPy-specific.

Fix.

  • Build the interpolated time as ts.tt_jd(tt[i], offset), so a tiny offset survives in the fraction.
  • If t_diff is still exactly zero (a sample landing exactly on a crossing), assume about one turn of the sky per day for that step instead of dividing 0/0.

Behavior change. For the reproduction, find_settings() now returns two settings, 2026-08-08T05:01:51.54Z (16 µs after t0) and 2026-08-09T05:01:15Z, because the first one is inside [t0, t1].

Testing.

  • New test_setting_a_few_microseconds_after_start_time uses the inputs from the issue. It fails on master with RuntimeWarning: invalid value encountered in divide and passes with this change.
  • With only the zero-t_diff guard applied, and without the interpolation change, the reproduction also passes and returns the same two settings to within nanoseconds.
  • I compared old and new code over all of 2026: 8 latitudes from −60° to 70°, Sun and Moon, risings and settings. Event counts and y flags are identical. Times differ by under 0.03 µs, except up to a few µs at 66°–70°, mostly on y=False events.
  • Caveat: I ran the tests with pytest on Windows, because assay needs fcntl. All almanac tests pass on NumPy 2.4.4 and 2.5.3. The rest of the suite has the same results before and after this change.

The other half of #1128 (.apparent() on a star with NaN coordinates) is not touched here; I described it in the issue.

🤖 Generated with Claude Code

If `start_time` fell just microseconds before a rising or setting, the
first interpolated guess in `_find()` was built as a single float Julian
date, which rounded back to exactly the sample time and made the hour
angle rate 0/0 = NaN.  The NaN time then reached the ephemeris, which
raised a misleading EphemerisRangeError.  Now the interpolated offset is
kept in the fraction, and a step that doesn't move `t` assumes about one
turn of the sky per day instead of dividing by zero.

Addresses the almanac half of skyfielders#1128.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@brandon-rhodes

Copy link
Copy Markdown
Member

I will have time for more open source work in November, and hopefully will look at this in the first week or two of the month. It will be nice to eliminate another edge case with the rising / setting finder.

@cstolting-collab

Copy link
Copy Markdown

Thanks for the careful split-time fix, @karimosama7. I tested e2cbfe4 on Linux; your original reproducer passes. I have published two small follow-up commits directly on top of your PR for you to incorporate if useful:

  1. Retain a usable hour-angle rate. On a later refinement, t_diff can be nonzero while ha_diff rounds to zero. The resulting zero rate makes the next time adjustment divide by zero. Initializing the rate to tau, then updating it only when both differences are nonzero, keeps your initial fallback and retains the last measured rate. This is an existing edge case, not a regression introduced by your PR.
  2. Preserve split-time precision throughout sampling and interpolation. Building the initial grid from start_time.tt can move a boundary by microseconds. Using end_time - start_time and offsets from the original Time also prevents a 20-microsecond window from losing an event. The commit includes a rise/set/transit regression with containing, preceding, following, and zero-width windows.

A small reproduction for the first case:

from skyfield import almanac, api

ts = api.load.timescale()
eph = api.load('de421.bsp')
observer = eph['earth'] + api.wgs84.latlon(0, -155.4761, 4139)
start = ts.tt_jd(2461260.709134074)
t, y = almanac.find_settings(
    observer, eph['sun'], start, start + 1.05,
    horizon_degrees=-7.779891398269016)
print(t.utc_iso(), y)

The rate addition recovers 23 failing exact-crossing cases in the original 532-case set. The precision addition resolves another 29 failing scenarios in the expanded set. The final candidate passes all 1,420 scenarios on NumPy 2.3.5 and 2.5.3. Across 960 ordinary/holdout cases, event counts and flags are unchanged; a separate bisection check of 96 events agrees within 1 microsecond using the same astronomical model.

The published evidence package includes patches, inputs, observations, scripts, and limitations. The combined native suite has one unchanged upstream datetime-list casting failure, so I am not presenting it as a completely green upstream release. The separate NaN-star change is in #1156.

The optional review harness uses my real CML runtime to enforce explicit checks; the numerical diagnosis and corrections were prepared with Codex assistance. Skyfield gains no CML dependency. I am happy to share the CML engine as well if it would help you or the maintainers. Thank you for the original investigation and fix.

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.

3 participants