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

[aendie]

Using Skyfield 1.54 I came across a somewhat rare warning:

I have reduced the code to an MWE as best I could. I get no information where this fails in my code.
Here is my simplified code to reproduce the error:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from datetime import date, datetime, timedelta
from skyfield import VERSION
from skyfield.api import Loader, wgs84, N, S, E, W
from skyfield import almanac

d = date(2000, 1, 1)        # select DATE       year, month, day
lat = 67.0                  # select latitude
obj = 6                     # select PLANET     0, 1, 2, 3, 4, 5, 6

print("\nSkyfield version = ",VERSION)
load = Loader("./")
ts = load.timescale()
eph = load('de421.bsp')
mercury = eph['mercury']
venus   = eph['venus']
earth   = eph['earth']
mars    = eph['mars']
jupiter = eph['jupiter barycenter']
saturn  = eph['saturn barycenter']
uranus  = eph['uranus barycenter']
neptune = eph['neptune barycenter']

pl_name = ['Mercury', 'Venus', 'Mars', 'Jupiter', 'Saturn', 'Uranus', 'Neptune']
pl_name2 = ['Mercury', 'Venus  ', 'Mars   ', 'Jupiter', 'Saturn ', 'Uranus ', 'Neptune']
pl_eph  = [mercury, venus, mars, jupiter, saturn, uranus, neptune]
planet = pl_eph[obj]

dt = datetime(d.year, d.month, d.day, 0, 0, 0)
dt_start = dt
daysinyear = (date(d.year+1, 1, 1) - date(d.year, 1, 1)).days
latns = 'N' if lat >= 0.0 else 'S'
print("year {}:  {} at latitude {:.1f}°{}".format(d.year,pl_name2[obj],abs(lat),latns))

#---------------------------------------------------------------------------
topos = wgs84.latlon(lat, 0.0 * E, elevation_m=0.0)
observer = earth + topos

print("... raw data from almanac.find_settings and almanac.find_risings:")
for i in range(daysinyear):

    t0 = ts.ut1(dt.year, dt.month, dt.day, dt.hour, dt.minute, dt.second)
    t1 = ts.ut1(dt.year, dt.month, dt.day+1, dt.hour, dt.minute, dt.second)

    planetset,  iS = almanac.find_settings(observer, planet, t0, t1)
    planetrise, iR = almanac.find_risings(observer, planet, t0, t1)
    for ndx, dt0 in enumerate(planetset):
        print("  set ",dt0.utc_iso(' ') ,"   ",iS[ndx])
    for ndx, dt0 in enumerate(planetrise):
        print("  rise",dt0.utc_iso(' '),"   ",iR[ndx])

    dt += timedelta(days = 1)

I am using Python 3.13.9. My code works perfectly, i.e. as expected, in 99.999% of cases and I have since enhanced the code to over double the length above. As far as I can tell, there are TWO ISSUES here:

• the RuntimeWarning from almanc.py in line 339 (Skyfield 1.54)
• almanac.find_settings() or almanac.find_risings() outputs a TIME value that dt0.utc_iso(' ') cannot convert and crashes with the message ValueError: cannot convert float NaN to integer.

Initially I believed that the ValueError is a consequence of the RuntimeWarning. Now I tend to think the reverse is true. I adapted the program below that scans all planets for a given range of latitudes for several years so that the ValueError is trapped: it prints rise ----- NaN error ----- or set ----- NaN error ----- usually without any RuntimeWarning, however runing the program above for the given year/planet/latitude will always print the RuntimeWarning first (on the date on which it fails).

I would REALLY APPRECIATE it if the RuntimeWarning and associated ValueError could be eliminated.
As always.... sincere thanks for all assistance!

UPDATES:
It occurs on 2031-04-18 and you can set for i in range(1): in the last for loop.
It then fails also with a second RuntimeWarning, this time in timelib.py at line 673:

I have now switched to incrementing the latitude in tenths of a degree as an integer, e.g. from 670 to 720, increment 1, and I divide by 10 to get the latitude (degrees) as a float. This way I avoid accumulating the error using a float increment of 0.1 repetitively.

Next failure point I noticed (using a later program version) is with Neptune on 2000-03-18 at 68.7°N:

By comparison, here is correct data from Skyfield 1.49:

Nothing unusual happens on 2000-03-18 ,,,,,,,,
Take a look at my chart for Neptune in 2000 at 68.7°N with Skyfield 1.49:

Initially I thought that something broke in Skyfield 1.54
Using Windows 11, numpy 2.4.1, and Skyfield 1.54, here are some random test cases using where I detected this issue occurs:

2000-09-17 50.0N Jupiter
2000-12-12 50.0N Neptune
2000-11-17 50.5N Uranus
2000-07-27 50.7N Jupter
2000-04-29 54.9N Neptune
2000-10-26 55.6N Neptune
2000-05-20 56.1N Uranus
2000-10-24 56.8N Neptune
2000-05-29 59.3N Neptune
2000-06-26 60.4N Neptune
2000-05-03 61.1N Neptune
2000-05-03 62.1N Uranus
2000-05-06 63.0N Neptune
2000-04-22 63.8N Uranus
2000-10-10 64.2N Uranus
2000-03-24 64.4N Neptune
2000-05-27 65.6N Uranus
2000-06-08 67.0N Neptune
2000-03-18 68.7N Neptune
2000-10-06 69.3N Neptune
2000-05-13 71.6N Uranus
2001-08-18 50.4N Mars
2001-01-15 51.2N Neptune
2001-10-29 52.3N Uranus
2001-04-29 56.0N Uranus
2001-07-05 56.1N Neptune
2001-11-22 57.3N Jupiter
2001-10-10 57.6N Saturn
2001-11-21 58.1N Jupiter
2001-05-13 58.5N Neptune
2001-01-16 58.8N Jupiter
2001-09-12 58.9M Meptune
2001-06-06 60.6N Neptune
2001-04-13 61.2N Neptune
2001-07-12 63.4N Jupiter
2001-01-02 64.0N Saturn
2001-10-31 64.2N Jupiter
2001-08-31 64.7N Neptune
2001-05-22 65.8N Uranus
2001-05-26 66.0N Neptune
2001-05-14 66.1N Uranus
2001-10-14 67.0N Neptune
2001-06-23 68.1N Uranus
2001-09-20 70.0N Neptune
2001-07-05 70.2N Uranus
2001-10-14 67.0N Neptune
2001-06-23 68.1N Uranus
2001-09-20 70.0N Neptune
2001-07-05 70.2N Uranus
2020-05-23 68.4N Uranus
2021-06-03 67.3N Saturn
2021-06-30 68.7N Uranus
2023-12-31 70.9N Uranus
2025-02-09 69.8N Uranus
2026-02-08 50.7N Uranus
2027-07-16 60.2N Uranus
2028-04-04 56.9N Uranus
2028-07-23 69.3N Saturn
2029-11-04 52.0N Mars
2029-04-23 54.0N Uranus
2029-01-26 62.2N Jupiter
2029-12-23 65.9N Uranus

Initially I suspected this issue is caused by Skyfield 1.54 alone, though I saw no reason for it.

So I decided to test using Skyfield 1.53 (and numpy 2.4.1 in Windows 11):
I tested the years 2027, 2028 and 2029 for latitudes between 50.0°N and 72.0°N in 0.1 degree steps [using integer steps because it is impossible to represent 0.1 (decimal) precisely in binary with a finite number of digits]

Using Windows 11, numpy 2.4.1, and Skyfield 1.53, I found these cases where it fails with the same symptoms:

2027-08-11 Neptune 50.0N
2028-01-20 Jupiter 55.7N
2029-07-14 Saturn 56.8N

So I see that this issue is not limited to Skyfield 1.54 (though 1.54 is more error-prone than 1.53).
(TBD: look at Skyfield 1.49, which I think is 100% free of this problem.)

I switched to another hardware platform: an old ASUS Zenbook with Intel i7-8565U CPU.
Using Windows 11, numpy 2.4.2, and Skyfield 1.54 I got exactly the same four failures for the year 2029 between 50.0°N and 72.0°N as above.

[brandon-rhodes]

Thank you, as usual, for such a thorough report, and for the example script!

Alas, when I run the script as printed above, I get no warning or error:

Skyfield version =  (1, 54)
planet:  Jupiter at latitude 67.49999999999997°N
... raw data from almanac.find_settings and almanac.find_risings:
  set  2031-04-17 04:44:32Z     True
  rise 2031-04-17 03:45:45Z     True
... ditto in chronological order and filtered for True events only:

(Was there supposed to be any output following that second :?)

So the error you are seeing might be specific to math on your own processor, that doesn't occur on mine. If you try running Python with the option python -W error, then I think it will turn the RuntimeWarning into an exception, and you'll get a traceback to see where it's happening.

Another solution to your problem of incrementing by 0.1 repeatedly would be to follow up the lat += 0.1 with lat = round(lat, 1), correcting any round-off error and bringing the number even with a single-digit decimal fraction.

I have tried adding this just after your own definition of planet, to try one of the other test cases you mention:

d = date(2000, 3, 18)
lat = 68.7
planet = eph['neptune barycenter']

But the result still showed no error:

Skyfield version =  (1, 54)
planet:  Jupiter at latitude 68.7°N
... raw data from almanac.find_settings and almanac.find_risings:
  set  2000-03-18 11:01:37Z     True
  rise 2000-03-18 06:34:10Z     True
... ditto in chronological order and filtered for True events only:

So I'm not sure yet how to proceed on investigating this.

As usual, your diagrams are very impressive! You have clearly detected and figured out how to display many edge cases over the years. The curve marked "DUSK" in the above diagram is especially elegant.

[aendie]

Greetings !!!! and thanks for your response.

I made a minor modification to my simplified test routine - now in my first comment - to enable selecting a planet using, for example, obj = 6 to select Neptune.

Using Skyfield 1.54, execution begins so ...

and a bit later fails so ...

I have a clean Windows 11 25H2 system running on an AMD Ryzen 7 9700X 8-Core Processor with "AGESA PI pre-1.3.0.0 updated" microcode (don't laugh) that was released 9th January 2026. I've seen many microcode releases during 2025 too.

I have two ideas, as I'm sure I can get you to reproduce the error: I'll perform the same test on the same PC using Ubuntu Desktop.

My other idea is that a you run my code that searches for any error in 0.1 degree steps between set latitudes - I choose to scan between 50.0°N and 72.0°N and choose any year within the ephemeris. Here is my all singing-and-dancing scan code, which also includes my code to fix ISSUE 998 finally (needs lots of testing still!):

UPDATE: I have the same as you - it does not fail in Ubuntu. I need to scan the same range of latitudes for all planets for some year also on Linux to see if it fails at all. If so: very strange why it only fails on Windows 11. Let's see....

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from datetime import date, datetime, timedelta
from skyfield import VERSION
from skyfield.api import Angle, Loader, Topos, wgs84, N, S, E, W
from skyfield import almanac
import numpy as np
import sys

year_from = 2029
year_to   = 2029

# increment latitude with an integer step (as 0.1 steps accumulate an error)
lat_from  = 500     # latitude x 10 as integer
lat_to    = 720     # latitude x 10 as integer
lat_step  = 1       # latitude step x 10 as integer

verbose = False
print("Skyfield version = ",VERSION,"\n")
load = Loader("./")
ts = load.timescale()
eph = load('de421.bsp')
mercury = eph['mercury']
venus   = eph['venus']
earth   = eph['earth']
mars    = eph['mars']
jupiter = eph['jupiter barycenter']
saturn  = eph['saturn barycenter']
uranus  = eph['uranus barycenter']
neptune = eph['neptune barycenter']

pl_name = ['Mercury', 'Venus', 'Mars', 'Jupiter', 'Saturn', 'Uranus', 'Neptune']
pl_name2 = ['Mercury', 'Venus  ', 'Mars   ', 'Jupiter', 'Saturn ', 'Uranus ', 'Neptune']
pl_eph  = [mercury, venus, mars, jupiter, saturn, uranus, neptune]
atmosref = 34.0 / 60        # 34' atmospheric refraction in degrees


for year in range(year_from,year_to+1):

    d = date(year, 1, 1)        # select DATE       year, month, day
    dt = datetime(d.year, d.month, d.day, 0, 0, 0)
    dt_start = dt
    daysinyear = (date(d.year+1, 1, 1) - date(d.year, 1, 1)).days

    lat_tenths = lat_from
    while lat_tenths <= lat_to:

        lat = lat_tenths/10
        topos = wgs84.latlon(lat, 0.0 * E, elevation_m=0.0)
        observer = earth + topos

        for obj in range(0,7):      # select PLANET     0, 1, 2, 3, 4, 5, 6
            planet = pl_eph[obj]
            latns = 'N' if lat >= 0.0 else 'S'
            print("year {}:  {} at latitude {}°{}".format(d.year,pl_name2[obj],str(abs(lat)),latns))

            #---------------------------------------------------------------------------
            aboveHorizon_tracked = None     # track 'above horizon' state
            last_true_dt0 = None            # track last true event into next "all day above/below horizon"
            last_true_rs = None
            last_isrise = None
            # a previous false RISE or SET event is tracked until the next (true or false) RISE or SET event
            last_false_rise = None; last_false_set = None

            dt = dt_start
            for i in range(daysinyear):

                t0 = ts.ut1(dt.year, dt.month, dt.day, dt.hour, dt.minute, dt.second)
                t1 = ts.ut1(dt.year, dt.month, dt.day+1, dt.hour, dt.minute, dt.second)

                planetset,  iS = almanac.find_settings(observer, planet, t0, t1)
                planetrise, iR = almanac.find_risings(observer, planet, t0, t1)
                # for ndx, dt0 in enumerate(planetset):
                    # print("  set ",dt0.utc_iso(' '),"   ",iS[ndx])
                # for ndx, dt0 in enumerate(planetrise):
                    # print("  rise",dt0.utc_iso(' '),"   ",iR[ndx])
                # print("... ditto in chronological order:")

                no_events = False
                if len(planetset) == 1 and len(planetrise) == 1:
                    dtS = planetset[0]
                    dtR = planetrise[0]
                    # if RISE == SET and both are False.....
                    if not iS[0] and not iR[0] and dtS == dtR:

                        no_events = True    # no RISE or SET on this day
                        dtX = dtR
                        degBelowHorizon = atmosref

                        # planet_elev = planet elevation in degrees
                        planet_elev = observer.at(dtX).observe(planet).apparent().altaz()[0].degrees
                        # aboveHorizon_tracked = True if planet is above the horizon, else False
                        above_horizon = (planet_elev > -degBelowHorizon)
                        if verbose: print("{} elevation at {} = {:.5f}".format(pl_name[obj],dtX.utc_iso(' '),planet_elev))
                        # print("{} below horizon when elevation < {:f}".format(pl_name[obj],-degBelowHorizon))
                        if verbose: print("{} above horizon: {}".format(pl_name[obj],above_horizon))
                        
                        if aboveHorizon_tracked is None:
                            aboveHorizon_tracked = above_horizon
                        elif aboveHorizon_tracked != above_horizon:
                            if above_horizon:
                                if last_false_rise is not None: # is a recent false RISE available?
                                    print("  rise", last_false_rise.utc_iso(' '),"    False --> True")
                                    last_false_rise = None
                                elif last_true_dt0 is not None and not last_isrise: # no - cancel the last event if SET
                                    print("  set ", last_true_dt0.utc_iso(' '),"    True --> False")
                                    last_true_dt0 = None
                                else:
                                    print("ERROR 1: out of sequence on {}".format(dtX.utc_iso(' '))); sys.exit(0)
                            else:
                                if last_false_set is not None:  # is a recent false SET available?
                                    print("  set ", last_false_set.utc_iso(' '),"    False --> True")
                                    last_false_set = None
                                elif last_true_dt0 is not None and last_isrise:     # no - cancel the last event if RISE
                                    print("  rise", last_true_dt0.utc_iso(' '),"    True --> False")
                                    last_true_dt0 = None
                                else:
                                    print("ERROR 2: out of sequence on {}".format(dtX.utc_iso(' '))); sys.exit(0)
                            aboveHorizon_tracked = above_horizon

                if not no_events:
                    # assemble events in chronological order
                    nR = len(planetrise); nS = len(planetset)
                    ndxR = ndxS = 0
                    riseset_time = []; isrise = []; isTrue = []
                    while ndxR < nR or ndxS < nS:
                        if ndxR < nR and ndxS < nS:
                            if planetrise[ndxR] < planetset[ndxS]:
                                riseset_time.append(planetrise[ndxR]); isrise.append(True); isTrue.append(iR[ndxR])
                                riseset_time.append(planetset[ndxS]); isrise.append(False); isTrue.append(iS[ndxS])
                            else:
                                riseset_time.append(planetset[ndxS]); isrise.append(False); isTrue.append(iS[ndxS])
                                riseset_time.append(planetrise[ndxR]); isrise.append(True); isTrue.append(iR[ndxR])
                            ndxR += 1; ndxS += 1
                        elif ndxR < nR:
                            riseset_time.append(planetrise[ndxR]); isrise.append(True); isTrue.append(iR[ndxR])
                            ndxR += 1
                        elif ndxS < nS:
                            riseset_time.append(planetset[ndxS]); isrise.append(False); isTrue.append(iS[ndxS])
                            ndxS += 1

                    for ndx, dt0 in enumerate(riseset_time):
                        rs = "  rise" if isrise[ndx] else "  set "
                        try:
                            dt0.utc_iso(' ')
                        except ValueError:
                            print(rs," ----- NaN error on {} -----".format(dt.strftime("%Y-%m-%d")))
                        if verbose: print(rs,dt0.utc_iso(' '),"   ",isTrue[ndx])

                        if isTrue[ndx]:     # ---- if it is a TRUE event ----
                            last_true_dt0 = dt0; last_true_rs = rs; last_isrise = isrise[ndx]

                            if aboveHorizon_tracked is None:
                                # initialize aboveHorizon state
                                if isrise[ndx]: # true RISE event
                                    aboveHorizon_tracked = True
                                else:           # true SET event
                                    aboveHorizon_tracked = False

                            elif isrise[ndx]:   # true RISE event
                                last_false_rise = None
                                if not aboveHorizon_tracked:        # if correct...
                                    aboveHorizon_tracked = True     # ...now above horizon
                                elif last_false_set is not None:    # is a recent false SET available?
                                    print("  set ", last_false_set.utc_iso(' '),"    False -> True")
                                    last_false_set = None
                                else:                               # no - then cancel this true RISE event
                                    print(rs, dt0.utc_iso(' '),"    True -> False")

                            else:               # true SET event
                                last_false_set = None
                                if aboveHorizon_tracked:            # if correct...
                                    aboveHorizon_tracked = False    # ...now below horizon
                                elif last_false_rise is not None:   # is a recent false RISE available?
                                    print("  rise", last_false_rise.utc_iso(' '),"    False -> True")
                                    last_false_rise = None
                                else:                               # no - then cancel this true SET event
                                    print(rs, dt0.utc_iso(' '),"    True -> False")

                        else:               # ---- if it is a FALSE event ----
                            if isrise[ndx]:     # false RISE event
                                last_false_rise = dt0
                                last_false_set  = None
                            else:               # false SET event
                                last_false_set  = dt0
                                last_false_rise = None

                dt += timedelta(days = 1)

        # ---------- end of planet FOR loop ----------
        
        lat_tenths += lat_step      # latitude step in tenths of a degree

    # -------------- end of latitude WHILE loop --------------

# ---------- end of year FOR loop ----------

[aendie]

I did manage to reproduce a same failure on Ubuntu (using Python 3.12.3 this time)

I ran the scan for all planets and latitudes between 50.0°N and 72.0°N for the year 2029:

It then failed on Mars at latitude 52.0°N...

So now I run the code in my first comment (with 2029, Mars and 52.0 as data) to list each RISE and SET per day...

And here I see it failed on 2029-11-04 (assuming there is one SET and one RISE per day) - and this agrees with what I detected in Windows 11.

UPDATES:

Please note that this issue does not always cause the "scan all latitudes" program to stop - usually I expect it to detect the runtime warning and continue execution, as in the second case I found for the year 2029 on Ubuntu:

Running the "show each RISE & SET event per day" code then displayed the date it fails on as 2029-10-14:

Please note that the failure points are different to those detected in Windows 11 even running on the same hardware!
On Ubuntu (with numpy 2.1.3) the failures for the year 2029 between latitude 50.0 and 72.0 in 0.1 degree steps are:

2029-11-04 Mars 52.0N
2029-10-14 Saturn 68.4N

I only noticed later that the numpy version was older, and it possibly has an effect on the results.

So, here are the next results from my testing...
On Ubuntu (with numpy 2.4.2) the failures for the years 2027 to 2029 between latitude 50.0 and 72.0 in 0.1 degree steps are:

2027-10-30 Uranus 54.7N
2027-02-10 Neptune 56.9N
2027-07-16 Uranus 60.2N
2028-02-10 Uranus 52.4N
2029-11-04 Mars 52.0N
2029-10-14 Saturn 68.4N

For reasons I cannot explain, these results are different from the consistent results by testing using Windows 11 on two different hardware platforms and CPU architectures [AMD and Intel] ... see my initial comment above.