Change in almanac.find_discrete logic

[aendie]

My Skyfield code was working happily until recently - however I can't pinpoint what changed as previous versions of Skyfield, e.g. 1.15, now also exhibit this (incorrect IMHO) symptom.

The problem is seen only when there is no moonrise or moonset on a particular day at a specific latitude. Take for example 16th. and 17th. February 2020 (Sunday and Monday). I calculated these times:

At latitude 72° N you see a moonrise (05:13) and set (07:02) on Sunday; and none on Monday. Skyfield gives me now this data:

Events on 02/16/2020:
05:13
07:02
True
False

Events on 02/17/2020:
False
False
False
False
False
False
False
False
False
False
False
False
False
False
False
False
False
False
False
False
False
False
False
False

And my complaint is why are there 24 moonsets on Monday (24 x False) ?
The zero moonrise/set event times is correct!

Previously the tuple returned by find_discrete had the same number of elements in each array (event times; event type). And the bug in my code is that I was looking at the number of items in the "event type" array - now 24 moonsets on Monday 02/17/2020 at 72°N.

The code to reproduce this is here:

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

import datetime
from skyfield import almanac
from skyfield.api import Topos, Star, load
from skyfield.nutationlib import iau2000b

def moonrise_set(d, lat, hemisph):
    # returns moonrise and moonset per date and latitude:

    lats = '%3.1f %s' %(abs(lat), hemisph)
    dt = datetime.datetime(d.year, d.month, d.day, 0, 0, 0)

    # Moonrise(s)/Moonset(s) on 1st. day ...
    moonrise, y = search_for_events(dt, lats, hemisph)
    print_events(dt,moonrise,y)

    # Moonrise(s)/Moonset(s) on 2nd. day ...
    dt += datetime.timedelta(days=1)
    moonrise, y = search_for_events(dt, lats, hemisph)
    print_events(dt,moonrise,y)

    return

def search_for_events(dt, lats, hemisph):
    locn = Topos(lats, '0.0 E')
    horizon = 0.8333
    dt -= datetime.timedelta(seconds=30)       # search from 30 seconds before midnight
    t0 = ts.utc(dt.year, dt.month, dt.day, dt.hour, dt.minute, dt.second)
    t1 = ts.utc(dt.year, dt.month, dt.day+1, dt.hour, dt.minute, dt.second)
    moonrise, y = almanac.find_discrete(t0, t1, moonday(locn, horizon))
    return moonrise, y

def print_events(dt, times, y):

    print "\nEvents on %s:" %dt.strftime("%m/%d/%Y")
    for t in times:
        t0 = ts.utc((t.utc_datetime() + datetime.timedelta(seconds=30)).replace(second=0, microsecond=0))
        print t0.utc_iso()[11:16]

    for events in y:
        print events

    return

def moonday(topos, degBelowHorizon):
    # Build a function of time that returns the "moonlight daylength".
    topos_at = (earth + topos).at

    def is_moon_up_at(t):
        """The function that this returns will expect a single argument that is a 
		:class:`~skyfield.timelib.Time` and will return ``True`` if the moon is up
		or twilight has started, else ``False``."""
        t._nutation_angles = iau2000b(t.tt)
        # Return `True` if the moon has risen by time `t`.
        return topos_at(t).observe(moon).apparent().altaz()[0].degrees > -degBelowHorizon

    is_moon_up_at.rough_period = 0.5  # twice a day
    return is_moon_up_at

ts = load.timescale()	# timescale object
eph = load('de421.bsp')	# ephemeris, valid between 1900 and 2050
earth   = eph['earth']
moon    = eph['moon']

d = datetime.date(2020, 2, 16)
moonrise_set(d, 72, "N")

Can we return to the "clean" situation where we have zero event times and zero event types on Monday?

[brandon-rhodes]

Alas! All my own code looks at times before looking at y so the noisy y vector was never getting noticed. Thanks for the full example, which let me easily find the problem. I've fixed it in the repository and added a test so it shouldn't ever happen again. If you'll run:

pip install https://github.com/skyfielders/python-skyfield/archive/master.zip

— then you should find the problem fixed. Please re-open this issue if not!

I will hopefully release a new Skyfield version in a couple of days that includes the fix.

[aendie]

Many thanks for the fast fix - agreed, the problem is fixed in the repository :-)

I was updating my six Repositories on GitHub. I have three versions (each in a Python2 and Python3 variant) that demonstrate my efforts using PyEphem and Skyfield:

• Pyalmanac
  is the fastest with its "somewhat limited" accuracy - but sufficient for nautical navigation though.
• SFalmanac
  is the slowest and most accurate; almost entirely based on Skyfield.
  Pyephem is only used for planet magnitudes (because these are not in Skyfield).
• SkyAlmanac
  is a hybrid version that is significantly faster than SFalmanac.
  Pyephem is used for planet magnitudes and for all planet transits, sunrise and sunset calculations as well as moonrise and moonset.

I would be happy if you take a look yourself at what I've done because it's all based on your hard work. Note that my Nautical Almanacs (the original I forked from Enno Rodegerdts) are published on The Nautical Almanac. These are one version ahead of my GitHub site (they include an image of the moon phase - done with the tikz package).

I'm happy to see that I'm slowly beginning to collect some stars in GitHub :-)