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):
#for i in range(5):

    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)

dt = dt_start
print("... ditto in chronological order and filtered for True events only:")
for i in range(daysinyear):
#for i in range(5):

    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)

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

    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)
• and most likely as a consequence of the RuntimeWarning: almanac.find_settings() or almanac.find_risings() can output NaN instead of an event TIME in the first output parameter, which should never be the case - this causes my code to crash

I would REALLY APPRECIATE it if this abrupt warning 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
Here are some random test cases using Skyfield 1.54 where I detected this issue occurs:

2000-06-08 Neptune 67.0N
2000-03-18 Neptune 68.7N
2000-10-06 Neptune 69.3N
2000-05-13 Uranus 71.6N
2001-10-14 Neptune 67.0N
2001-06-23 Uranus 68.1N
2001-09-20 Neptune 70.0N
2001-07-05 Uranus 70.2N
2020-05-23 Uranus 68.4N
2021-06-03 Saturn 67.3N
2021-06-30 Uranus 68.7N
2023-12-31 Uranus 70.9N
2025-02-09 Uranus 69.8N
2026-02-08 Uranus 50.7N
2027-07-16 Uranus 60.2N
2028-04-04 Uranus 56.9N
2028-07-23 Saturn 69.3N
2029-11-04 Mars 52.0N
2029-04-23 Uranus 54.0N
2029-01-26 Jupiter 62.2N
2029-12-23 Uranus 65.9N

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.
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]

I found these cases where Skyfield 1.53 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).
I now want to look at Skyfield 1.49, which I thought was 100% free of this problem.

[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  = 670     # latitude x 10 as integer
lat_from  = 500     # latitude x 10 as integer
lat_to    = 720     # latitude x 10 as integer
#lat_to    = 670     # 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 "
                        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.

So I think you should be able to reproduce this error.