find_discrete is more reliable than find_risings/settings

[aendie]

I am running initial tests having converted my Nautical Almanac software to use find_risings() and find_settings() when Skyfield VERSION is 1.48 or higher. Lower versions still use find_discrete(). As both are implemented I began comparing them and I came across a discrepancy that I wish to describe in detail here.

To gain confidence in my 'quick and dirty' test program, I quote three test cases: the first two show correct results. I am aware that almanacs print time in UT (=UT1) and my test program simply prints UTC time, which is irrelevant here. Furthermore times in almanacs are rounded to minutes so a 24-hour search begins at 30 seconds before midnight. My test program begins at midnight (again irrelevant here).

TEST CASE 1: 31-May 2023 to 2-June 2023 at latitude 68°

I expect to see these results (my almanac data using find_discrete()) for Wed Thu Fri:

The test Program prints results from find_risings() and find_settings() as 'moonrise' and 'moonset' and results from find_discrete() as 'moonevent' where '0' represents a moonset and '1' a moonrise:

The data looks good!

TEST CASE 2: 31-May 2023 to 2-June 2023 at latitude 72°

Again, the new routines match the data with the legacy find_discrete() logic. Looks good!

TEST CASE 3: 18-Feb 2023 to 20-Feb 2023 at latitude 70°

Here I expected to see this:

but instead saw this:

Note that my old almanac code searches day-by-day forward for the next moonrise or moonset event to determine the current moon state. "All day above horizon" (a white rectangle) is printed because the next event detected is a moonset at 12:05 on Sunday. So what went wrong here? The Test Program shows this:

The "moonrise" at 11:35 on 19-Feb is 'False'. There is no moonrise on 19-Feb. I look into the official HMNAO Nautical Almanac and see a moonrise at 11:28 on Feb 19th. Exactly what find_discrete() prints (when rounded). Why did find_risings() miss it? I use a calculated horizon value but if you omit the last parameter "-horizon" from find_risings(), it makes no difference.

But the horizon value does have an impact on the results. Some reverse-engineering showed that increasing the value of "horizon" by 0.048 degrees not only fixes the problem but gives a result very close to 11:28, i.e. this works:

moonrise, yR = almanac.find_risings(observer, moon, t0, t1, -horizon-0.048)

Is find_risings() using an incorrect horizon value?

Another thing that looks suspicious: on Feb 18th: both "False" events have the same time '2023-02-18 10:34:20Z'.
I understand "False" events to mean that the moon vertical motion changes, even though this is below the horizon, i.e. not an actual moonrise or moonset. But I cannot believe that the moon changes direction twice within a split second.

I cannot dig deeper ....... I hope someone can look into this and either explain or fix the problem.

Here is my Test Program (note I use de421.bsp as my ephemeris):

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

from datetime import date, datetime, timedelta
from math import atan, degrees
from skyfield import VERSION, almanac
from skyfield.api import Loader, Topos, wgs84, N, S, E, W
from skyfield.nutationlib import iau2000b

def getHorizon(t):
    # calculate the angle of the moon below the horizon at moonrise/set

    position = earth.at(t).observe(moon)   # at noontime (for daily average distance)
    distance = position.apparent().radec(epoch='date')[2]
    dist_km = distance.km
    sdm = degrees(atan(1737.4/dist_km))   # volumetric mean radius of moon = 1737.4 km
    horizon = sdm + 0.5666667	# moon's equatorial radius + 34' (atmospheric refraction)
    return horizon

def f_moon(topos, degBelowHorizon):
    # Build a function of time that returns the moon state.
    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,
		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

d = date(2023, 2, 18)       # modify the starting date as required
lat = 70.0                  # modify latitude as required (use float for Topos)

load = Loader("./")
ts = load.timescale()
eph = load('de421.bsp')
earth   = eph['earth']
moon    = eph['moon']
dt = datetime(d.year, d.month, d.day, 0, 0, 0)
print("Skyfield version = {};  Results for latitude {} ...".format(VERSION,lat))

for i in range(3):      # sample 3 successive dates

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

    topos = wgs84.latlon(lat, 0.0 * E, elevation_m=0.0)
    observer = earth + topos
    horizon = getHorizon(t0noon)

    moonrise, yR = almanac.find_risings(observer, moon, t0, t1, -horizon)  # -horizon-0.048 fixes problem
    moonset,  yS = almanac.find_settings(observer, moon, t0, t1, -horizon)

    locn = Topos(lat, "0.0 E")
    moonevent, y = almanac.find_discrete(t0, t1, f_moon(locn, horizon))

    print("{}:  {} x moonrise, {} x moonset;  {} x moonevent;  horizon = {:.4f}°".format(dt.strftime("%Y-%m-%d"),len(yR),len(yS),len(y),horizon))
    print("   moonrise ",moonrise.utc_iso(' '),"   ",yR)
    print("   moonset  ",moonset.utc_iso(' ') ,"   ",yS)
    print("   moonevent",moonevent.utc_iso(' '),"   ",y)
    
    dt += timedelta(days=1)

Sincere thanks to whoever looks into this case. I hope I provided a useful test case here.

[brandon-rhodes]

Another thing that looks suspicious: on Feb 18th: both "False" events have the same time 2023-02-18 10:34:20Z. I understand "False" events to mean that the moon vertical motion changes, even though this is below the horizon, i.e. not an actual moonrise or moonset. But I cannot believe that the moon changes direction twice within a split second.

Let's start with this item first, because maybe the documentation could be clearer here. On 02-18 there is no moonrise, so whether you ask for the Moon rise or set, Skyfield instead can only give you the moment of culmination, where the Moon was as close to the horizon as possible — the same moment, whether you are asking about rise or set. Because rise and set only become distinct points when the horizon is intersected.

So we would expect both return values to be exactly the same time and elevation.

How can the documentation be improved so that users are not surprised at both the rise and the set routine returning the same "sorry, it didn't rise or set" value?

Maybe you could clarify what return value you expected in this case. Thanks!

[aendie]

Greetings, Brandon! Thanks for your prompt response.
I see only two scenarios for a 'False' value returned by find_risings() or find_settings() (maybe there are more?). Using the moon as an example:

Scenario 1: almost moonrise

The moon is "all day below the horizon". Its declination is increasing (approaching moonrise); stops before moonrise; then declination decreases again. The moment it stops and changes direction can be labelled "False".

Scenario 2: almost moonset

The moon is "all day above the horizon". Its declination is decreasing (approaching moonset); stops before moonset; then declination increases again. The moment it stops and changes direction can be labelled "False".

The moon cannot be "all day above the horizon" and "all day below the horizon" so you can have only have scenario 1 or 2.

Alternatively ... and this explanation is independent of the day boundary (let's see what happens within 24 hours beginning at midnight): the declination of an object can only stop under these conditions:

• declination is decreasing - stops - and then increases again. Mathematically this is known as a cusp (ignore the units - I grabbed the graphic from the internet):

• declination is increasing - stops - and then decreases again.

• declination is decreasing - stops - and then decreases again. (Extremely rare but theoretically possible)

• declination is increasing - stops - and then increases again. (Extremely rare but theoretically possible)

Maybe I have overlooked something, but I cannot see how both of the first two conditions can occur simultaneously (at least within a split second) at '2023-02-18 10:34:20Z' to quote my example on above on 2023-02-18 at latitude 70°N.

Regarding the last two conditions, one cannot say that either a moonrise or moonset was almost reached because "the story continues" towards either a moonrise or a moonset. Besides, it is only possible if "stops" has a bandwidth, i.e. if a declination delta less than a given value within a defined time gap counts as "declination stopped increasing or decreasing". In other words ... it will never happen.

And permit me to compliment you on your excellent documentation in Skyfield ... I always have trouble writing documentation because I am never satisfied with it.

[brandon-rhodes]

Could you clarify what return value you are suggesting for the day 2023-02-18, instead of the two False values associated with the time 10:34:20Z? Are you saying that one of the two should be True instead? If so, which one? Or — are you asking that one of them return a different time instead? I am not yet entirely clear on what alternative return value you are proposing. Thanks for clarifying!

[EndlessDex]

There is indeed a bug here. The find_risings and find_settings functions use a shortcut that the older risings_and_settings does not use.

In an attempt to be more efficient, the new functions use an equation to calculate expected hour angle of alt=horizon (and then looks for the time of that hour angle via Newton's method).

See this equation from Fundamental Astronomy:

Where h = hour angle, delta = declination, phi = latitude, and a = altitude:

Unfortunately this equation is undefined for lat - decl > 90 and lat - dec < -90 (plus or minus a bit based on horizon).

Here is a plot of the relevant days:

Because the function used with Newton's Method is not continuous the solution does not converge properly. It gets partway there by pure luck due to one of the sample points being near the rise time.

@brandon-rhodes Proposed solution would be to detect if search function (called _setting_hour_angle in code) is undefined for any of the search times and to fall back to the old method of topos_at(t).observe(target).apparent().altaz()[0].degrees > horizon.

[aendie]

@brandon-rhodes The return values I would expect on 2023-02-18 at latitude 70°N are:

• moonrise ['2023-02-18 10:34:20Z'] [False]
• moonset [] []

to signify that the moon, being below the horizon, rose as high as it could but failed (hence [False]) to reach moonrise.
There is neither a moonset on that day, nor does the moon try to SET by decreasing its declination till it can't any more and then rising again. As it happens the moon is SET anyway, i.e. below the horizon, at the time '2023-02-18 10:34:20Z' when its declination changes from increasing (rising) to decreasing (setting).

Does that make sense?

Please note that I am not interested in the 'False' values in the results from find_risings() and find_settings(). There may be others who are interested in this data and their views are more important than mine. I am not an astronomer and I am willing to be corrected by those who have studied astronomy or are simply experts in this field.

Having mentioned that, I find it awkward to skip the [False] values maintaining chronological order. The data I am interested in is the same as find_discrete() gives me:

• actual moonrise and moonset times
• all sorted by time increasing

It's obvious that a moonrise is always followed by a moonset, and a moonset is always followed by a moonrise - though not necessarily on the same day.

So what I would really appreciate is .....

• either an argument in find_risings() and find_settings() that gives me ONLY the [True] values,
• or some way to remove the [False] time values from a Time object with multiple times, like t.pop(index)

(Neither did I find Skyfield documentation on how to create a Time object with multiple times, by the way. All examples, as far as I am aware, create a Time object with a single time.)

As this request is unrelated to this #998 issue, I could open a new issue if wanted.

[EndlessDex]

some way to remove the [False] time values from a Time object with multiple times, like t.pop(index)

This can be accomplished like so with python semantics:

risings = find_risings(observer, moon, t0, t1)
risings = [time for time, true_event in zip(*risings) if true_event]

>> [<Time tt=2459995.8982223687>, <Time tt=2459996.8716849447>, <Time tt=2459997.8519764347>]

This will return a list of time objects for which the event was true.

As far as accessing the times themselves, its a numpy array in a trench coat. So you can do most anything you can do with one of those.

Here is a way to remove the false events with numpy semantics:

time, true_event = find_risings(observer, moon, t0, t1)
risings = time[true_event == True]

>> <Time tt=[2459995.89822237 2459996.87168494 2459997.85197643]>

This will return a single Time object containing the array of Time's

[EndlessDex]

(Neither did I find Skyfield documentation on how to create a Time object with multiple times, by the way. All examples, as far as I am aware, create a Time object with a single time.)

This can be found here Date Arrays.

Also mentioned here UTC and Your Timezone at the bottom regarding turning a list of datetimes into a date array.

Date Arithmatic might also interest you.

[aendie]

@EndlessDex Many thanks for your explanations above. I will find them very useful!

While scanning for further discrepancies between find_discrete() and find_risings(), find_settings() in 2023 (I only have an official HMNAO Nautical Almanac for 2023), I came across several cases where the latter were more accurate than find_discrete(). I am pleased that the new routines will soon be improved.

I also came across a few test cases which I think belong to the category above that needs a bugfix. I would like to add them here for more comprehensive testing. (I hope more test cases are useful to you.)

NOTE: When comparing results with the HMNAO Nautical Almanac, my Test Program uses a calculated horizon value. For sunrise and sunset, the HMNAO Almanac uses 16' for the semi-diameter and 34' for horizontal refraction but gives no mention what is used for moonrise and moonset ('Explanation', Section 12.) I found that my calculated horizon value gives better agreement with the HMNAO data though I only took a few cases into consideration.

TEST CASE 4: 28-Apr 2023 to 30-Apr 2023 at latitude 72°

HMNAO data: On 28th April the moon is above horizon all day. So the next event must be a moonset (at 7:33) followed by a moonrise (at 7:38) on 29th April. find_settings() does not detect the moonset - only the moonrise is detected (at 7:43). find_discrete() does not detect either. On 30th April all results are in perfect agreement with the HMNAO data.

TEST CASE 5: 3-Jul 2023 to 5-Jul 2023 at latitude 66°

HMNAO data: On 3rd and 4th July the moon is below horizon all day. So the next event must be a moonrise (at 1:46) followed by a moonset (at 1:57) on 5th. July. find_risings() does not detect the moonrise - only the moonset is detected (at 1:58). find_discrete() does not detect either.

TEST CASE 6: 27-Jul 2023 to 29-Jul 2023 at latitude 68°

HMNAO data: On 28th and 29th July the moon is below horizon all day. So the previous event on 27th must be a moonset (at 19:25) preceded by a moonrise (at 18:51) on 27th July. find_settings() does not detect the moonset - only the moonrise is detected (at 18:50). find_discrete() does not detect either.

[aendie]

TAKE TWO: Addressing Issues and Corrective Action

The test cases above (4 to 6) are not in the same category as Test Case 3: here I detected discrepancies not with find_discrete() but with the official HMNAO Nautical Almanac. I can however claim that find_discrete() is more reliable than find_risings()/find_settings() in the sense that I never came across this issue before with find_discrete() no matter for which year (say between 2020 and 2050) I was generating nautical almanac tables.

One point is that using find_discrete(), a moonrise is always followed by a moonset; and a moonset is always followed by a moonrise. Due to the fact that find_risings() and find_settings() are separate functions, one is not aware of the results the other gives.

The test cases above (4 to 6) are marginal cases where it is difficult to detect what the moon state is, and it only fails when given marginal cases at relatively high latitude. Okay, so the onus is on the end user (as ‘step one’) to detect these cases and provide assistance by avoiding a moonrise that follows a moonrise as well as a moonset that follows a moonset. Be aware that a pre-requisite is to filter out all “False” events first (as none of these are candidates for an actual documented moonrise or moonset). This can be done with good accuracy by making a note of the moon state at the time the search begins. This assumes that the moon state can be reliably assessed at the beginning search time, which is not always the case, however it is possibly an improvement rather than discounting it entirely.

We can implement this idea in the Test Program provided to show any ‘True’ moonrise or moonset that was rejected on logical grounds, i.e. “a moonrise cannot follow a moonrise” and “a moonset cannot follow a moonset”. The same test cases as above are still not satisfactory as key events are not detected. (The issue is that BOTH incorrect events ARE detected by find_risings()/find_settings() AND that correct events ARE NOT detected.)

TEST CASE 4: 28-Apr 2023 to 1-May 2023 at latitude 72°

HMNAO data: On 28th April the moon is above horizon all day at latitude 72°N. On 29th the moon sets at 7:33 and rises 5 minutes later. So I am not concerned if we ignore these two events altogether. As you see, the single moonrise has been discarded. A day later there is a moonset at 5:18 and a moonrise at 11:22. Perfect agreement. On 1st May there is a moonset at 3:34 and a moonrise at 13:35. Almost perfect. There is NO ISSUE here - all looks well now.

TEST CASE 5: 3-Jul 2023 to 6-Jul 2023 at latitude 66°

HMNAO data: On 3rd, 4th and 5th July the moon is below horizon all day at latitude 66°N. The incorrect moonset has now been discarded so these days are correct. On 6th June there is a moonrise at 00:13, a moonset at 5:34 and a moonrise at 23:49. Excellent! Perfect agreement now. NO ISSUE here now.

TEST CASE 6: 26-Jul 2023 to 29-Jul 2023 at latitude 68°

HMNAO data: On 26th July at latitude 68°N there is a moonrise at 15:31 and a mooonset at 21:01. Perfect agreement. On 27th there is a moonrise at 18:51 and a moonset at 19:25 (34 minutes later). ONLY THE MOONRISE is detected. No event is detected on the following two days suggesting that moon is above the horizon all day. THIS IS INCORRECT (or to quote Brandon: "Well, Drat!"). HMNAO data shows that on 28th, 29th (and 30th) the moon is BELOW THE HORIZON all day. This discrepancy (lasting several days) is caused by the undetected moonset at 19:25 on 27th. July. I hope a bugfix will cure this mistake too. (It is easy for me to ignore events that must be incorrect; it is impossible for me discover events that are missing using find_risings() and find_settings().

I highly appreciate your time and effort when looking into this issue for me.

The modified Test Program is included here (no points for elegant programming – it merely demonstrates the necessary logic):

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

from datetime import date, datetime, timedelta
from math import atan, degrees
from skyfield import VERSION, almanac
from skyfield.api import Loader, wgs84, N, S, E, W  # Topos is deprecated in Skyfield v1.35!
from skyfield.nutationlib import iau2000b

def initial_moonstate(t, topos, degBelowHorizon):
    # calculate the moonstate at time 't'
    topos_at = (earth + topos).at   # position of this Earth location at time 't'
    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

def getHorizon(t):
    # calculate the angle of the moon below the horizon at moonrise/set

    position = earth.at(t).observe(moon)   # at noontime (for daily average distance)
    distance = position.apparent().radec(epoch='date')[2]
    dist_km = distance.km
    sdm = degrees(atan(1737.4/dist_km))   # volumetric mean radius of moon = 1737.4 km
    horizon = sdm + 0.5666667	# moon's equatorial radius + 34' (atmospheric refraction)

    return horizon

def f_moon(topos, degBelowHorizon):
    # Build a function of time that returns the moon state.
    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,
		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

d = date(2023, 4, 28)       # modify the starting date as required
lat = 72.0                  # modify latitude as required (use float for Topos)

load = Loader("./")
ts = load.timescale()
eph = load('de421.bsp')
earth   = eph['earth']
moon    = eph['moon']
dt = datetime(d.year, d.month, d.day, 0, 0, 0)
print("Skyfield version = {};  Results for latitude {} ...".format(VERSION,lat))
pickRISE = None         # 'True' if we expect a moonrise next (uninitialized)

for i in range(4):      # sample 4 successive dates

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

    topos = wgs84.latlon(lat, 0.0 * E, elevation_m=0.0)
    observer = earth + topos
    horizon = getHorizon(t0noon)
    #horizon = 0.8333333         # 16' + 34' as used in the HMNAO Nautical Almanac

    moonrise, typeR = almanac.find_risings(observer, moon, t0, t1, -horizon)
    moonrise = moonrise[typeR == True]      # discard any 'False' events
    RISEevents = len(moonrise)
    moonset,  typeS = almanac.find_settings(observer, moon, t0, t1, -horizon)
    moonset = moonset[typeS == True]        # discard any 'False' events
    SETevents = len(moonset)
    ALLevents = []
    for t in moonrise:
        ALLevents.append((t,"RISE"))
    for t in moonset:
        ALLevents.append((t,"SET"))
    ALLevents.sort(key=lambda x: x[0])

    moonevent, y = almanac.find_discrete(t0, t1, f_moon(topos, horizon))

    if pickRISE is None:
        pickRISE = not initial_moonstate(t0, topos, horizon)
        if pickRISE: print("begin searching for a moon RISE...")
        else: print("begin searching for a moon SET...")

    print("{}:  {} x moonrise, {} x moonset;  {} x moonevent;  horizon = {:.4f}°".format(dt.strftime("%Y-%m-%d"),len(moonrise),len(moonset),len(y),horizon))

    for t, event_type in ALLevents:
        if event_type == 'RISE':
            msg = "   moonrise  " + t.utc_iso(' ')
            if not pickRISE: msg += "   <= EVENT DISCARDED"
            else: pickRISE = not pickRISE
            print(msg)
        if event_type == 'SET':
            msg = "   moonset   " + t.utc_iso(' ')
            if pickRISE: msg += "   <= EVENT DISCARDED"
            else: pickRISE = not pickRISE
            print(msg)

    print("   moonevent",moonevent.utc_iso(' '),"   ",y)
    
    dt += timedelta(days=1)

[EndlessDex]

@aendie Thanks for all the interesting data!

Skyfield uses this USNO generator as a source of truth for comparison against. This table uses the Refraction + Apparent Radius method. For "normal" locations (ex W112 31, N36 57) Skyfield is accurate to within one minute for all events. Generating a table for (W0, N70) shows that both old and new methods are grossly wrong (as you have noted).

PS. Take a look at the other Data Services of the USNO there's a lot of cool stuff!

A couple requests for you:

1. Can you compare the USNO site against your almanac? Want to make sure we have a common source of truth.
2. Can we agree on a latitude to test at? It does not matter to me which one (68, 70, 72, etc) but changing it in each test case is problematic. By picking a location and then inspecting every event in a year (2023 in this case), we should be able to create a robust solution. My goal would be to get 100% match for a given problematic latitude and then automate a test to cover a variety of other ones.
3. Can you review the error cases you have presented and determine if the moon has an invalid declination? find_risings/find_settings CANNOT return the correct answer if the latitude minus the declination is greater than 90 degrees or if it is even close to 90 for that matter. So when looking for error cases, we need to put all of these cases in a big "Bin 1" to be resolved by fixing the search method. What I am really interested in is "Are there any error cases NOT in Bin 1?".

I am working on resolutions to this and already have a resolution to the same false event being reported for both rise and set. I know you have stated you don't really care about the false events but it is a feature to support regardless.

Side note: I will be using Skyfield's built-in horizon calculation for all testing. (-34' - r_moon/d)

I am looking for both algorithmic (better math) and computational (more CPU time) solutions for the "Bin 1" cases. The latter is easier and I already have a working partial solution but it is significantly more expensive. One of the stated goals of the new method is to be more efficient than find_discrete paired with find_rising_and_settings. Depending on if I can find a better algorithm or not, there may need to be a "fast" segment and a "slow" segment of these functions. Where the "fast" method is used for low latitudes but above a certain point "slow" does some more expensive calculations in the name of accuracy.

I hope we can all agreed that providing WRONG results is worst thing we can do and that it is worth extra time-cost to avoid that.

[brandon-rhodes]

Good afternoon, @aendie and @EndlessDex!

With the holiday weekend upon us, it looks like I have some spare hours for this Skyfield issue. Through the last week, when I’ve had spare time, I have been digging into the behavior of the search routine for this interesting 2023-02-19 70°N case — but haven’t yet sat down to write further comments because of the sheer number of topics that require replies at this point.

So instead of trying to write up a knock-down drag-out comment that addresses them all, I’ll try writing a series of comments through this afternoon that will hopefully get me caught up.

For this first reply, I’ll focus on a few points that I think I can answer quickly:

• Since I got sidetracked a few comments ago by the question of API design and False return values, I should state clearly that I agree that the ‘find’ routines are missing the moonrise on 2023-02-19 for 70°N. And I agree completely that this shortcoming should be fixed!

• No, find_risings() is not using an incorrect horizons value; it’s using the correct one. The problem has to do with the geometry of very slight horizon-grazing moonrises. It was just luck that your twiddling of the horizon value, @aendie, kicked the search routine out of that unhappy edge case by moving the Moon farther from the southern edge of the horizon.

• @aendie, the Skyfield API was designed around the idea that filtering out values in Python is easy, so that callers can easily throw away the False values if they need to. Have you tried out either of the approaches suggested by @EndlessDex for discarding False values and times, to see if they work for your case? It was his second solution that I guessed that folks would most often use, since Skyfield works so closely with NumPy, but the first one — the list comprehension — will also work just fine, at the expense of a bit of speed. (And thank you, @EndlessDex, for jumping in with that quick Python tutorial when I was busy last week! It was a generous contribution to the discussion.)

• Thanks for linking to the ‘Date Arrays’ documentation, @EndlessDex! Let me know if there’s some way I can tweak the docs, @aendie, so that the ‘Date Arrays’ information is easier for folks to discover, since time arrays are the basic return type of pretty much all almanac functions.

• @EndlessDex, I hope you will be encouraged to receive my entire agreement with your point: ‘I hope we can all agreed that providing WRONG results is worst thing we can do and that it is worth extra time-cost to avoid that.’

The other points you both have raised in your comments will require more investigation and/or longer answers, so I’ll pause here and click ‘Comment’, so that these answers can go ahead and become part of the thread while I start work on some others.

[EndlessDex]

@brandon-rhodes Thanks for the updates :) Definitely appreciate that there is a lot of discourse to parse through.

I think the key point for me is that the hour angle equation breaks down for certain times/locations. This leads to the newtons method not converging correctly because it's trying to develop a slope for a non-linear function. Please reference my comment here for plots and equations.

I hope we can get something better than falling back on find_discrete(risings_and_settings) but haven't been able to yet. Doing so seems to fix things for the stated cases. I have not explored the notes @aendie made about find_discrete also being wrong for some cases.

[aendie]

@EndlessDex Thanks for your very prompt response.
In answer to your request #1: I checked all moonrise/moonset data for 2023 at latitude 70°N....

Initial comparison between the official HMNAO Nautical Almanac and the US Navy – Astronomical Applications Department Moonrise/Moonset table data.

The website is interesting and useful. I note that the results provide UT (Universal Time = UT1) just as in the Nautical Almanac. Quickly scanning the dates and times below, they agree very well with differences that occasionally just cross into the next or previous minute. The Nautical Almanac for 2023 was published in 2022 (earlier obviously). The web site might be based on more recent data.

NA = official HMNAO Nautical Almanac 2023
AAD = US Navy - Astronomical Applications Department Moonrise/Moonset table data

All discrepancies in 2023 at latitude 70°N:
2023-02-12 moonset 01:09 (NA), 01:08 (AAD)
2023-02-21 moonset 18:29 (NA), 18:30 (AAD)
2023-02-22 moonset 20:43 (NA), 20:44 (AAD)
2023-03-07 moonset 08:01 (NA), 08:02 (AAD)
2023-03-21 moonset 17:46 (NA), 17:47 (AAD)
2023-04-02 moonrise 12:17 (NA), 12:16 (AAD)
2023-05-15 moonset 14:34 (NA), 14:35 (AAD)
2023-05-28 moonrise 11:09 (NA), 11:08 (AAD)
2023-05-31 moonset 01:29 (NA), 01:30 (AAD)
2023-06-14 moonrise 00:07 (NA), 00:06 (AAD)
2023-06-23 moonset 01:48 (NA), 01:49 (AAD)
2023-07-07 moonset 00:26 (NA), 00:25 (AAD)
2023-07-10 moonset 13:59 (NA), 14:00 (AAD)
2023-07-25 moonset 21:15 (NA), 21:14 (AAD)
2023-09-26 moonrise 19:28 (NA), 19:27 (AAD)
2023-11-12 moonset 13:30 (NA), 13:31 (AAD)
2023-11-24 moonrise 04:31 (NA), 04:32 (AAD)
2023-12-10 moonrise 07:53 (NA), 07:52 (AAD)
2023-12-18 moonrise 13:15 (NA), 13:14 (AAD)

So the data agrees very well and is probably off by only a few seconds (of time).

[brandon-rhodes]

Two quick notes:

The Nautical Almanac for 2023 was published in 2022 (earlier obviously). The web site might be based on more recent data.

Happily, the data source should not be an issue for Moon rise and set, since we have known the Moon's orbit in such exquisite detail for a couple of decades. To double-check, I just ran a full year of 2023 Moon rises and sets from 70° north using ephemerides DE406 (1997), DE421 (2008), and DE440 (2020). All gave exactly the same output. So, for contemporary Moon positions, feel free to use data sources that are 25+ years old, if that's faster or more convenient for you; it shouldn't make a difference.

TEST CASE 4: 28-Apr 2023 to 1-May 2023 at latitude 72° … HMNAO data: …On 29th the moon sets at 7:33 and rises 5 minutes later. So I am not concerned if we ignore these two events altogether.

Note that the USNO does not see the Moon set at 7:33 or rise 5 minutes later; instead, it outputs **** **** for that date, time, and location. So my plan is for Skyfield to continue agreeing with the USNO in that case, and ignore the set and rise reported by the Nautical Almanac.

[brandon-rhodes]

@aendie — Having read back through all of your Test Cases, am I correct that you are happy with the output of find_discrete() and that, at this point, we needn't do any further investigation of its behavior?

I'm actually surprised that you're getting good results from it, because is_moon_up_at.rough_period = 0.5 means that it starts by taking only 24 samples per day, which are each one hour apart. Any rise-and-set sequence (or set-and-rise, if the Moon is mostly above the horizon that day) that is short enough to fit between two of those every-hour sample times would be missed, since find_discrete() only does further investigation if it sees a switch from down-to-up or up-to-down between one sample and the next.

If you want to catch the Moon even when it's only up for one minute, maybe try something like is_moon_up_at.step_days = 1.0 / 24.0 / 60.0 instead of using .rough_period. That will check every minute of the day. And it's far easier to reason about, since step_days directly sets how far apart the samples are, whereas the old deprecated .rough_period attribute is more roundabout.