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.

[brandon-rhodes]

Skyfield uses this USNO generator as a source of truth for comparison against. This table uses the Refraction + Apparent Radius method. … Generating a table for (W0, N70) shows that both old and new methods are grossly wrong (as you have noted).

@EndlessDex — could you share the table? I would have expected the old method to agree completely with the USNO, and if it doesn't, then that's a separate issue that will need to be fixed. Working code to generate the table would also be helpful. Thanks!

[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?

@brandon-rhodes I agree with your comments. find_discrete() is not perfect, however I was satisfied with the accuracy it gave me. Your idea to define .step_days instead of using .rough_period is also good. I really meant that it never gave me a moonrise following a moonrise, or a moonset following a moonset. So the results were consistent (regarding the rectangles that display 'moon all day above horizon' and 'moon all day below horizon') even if not perfect in marginal cases (typically at higher latitudes). So I would leave it as it is, especially as the new routines should soon be fixed.

[brandon-rhodes]

… even if not perfect in marginal cases …

By "not perfect", I think you mean that it sometimes differs by 1 minute from a USNO moonrise or moonset? Or is there a case where it does worse than that — in which case I would definitely want to get find_discrete() fixed?

[aendie]

By "not perfect", I think you mean that it sometimes differs by 1 minute from a USNO moonrise or moonset? Or is there a case where it does worse than that — in which case I would definitely want to get find_discrete() fixed?

@brandon-rhodes It's well over a year since I was working with the Nautical Almanac tables. I would have to look into that again. I should be able to give you feedback on this but it'll take some time still as I am right in the middle of updating this software that hasn't received an update since April 2023. I am making it compatible with both find_discrete() and the new find_risings() / find_settings() dependent on the Skyfield version the end user has. Also, this needs different logic whether it is running in single-processing mode or multi-processing mode, which the end user can select.

Also .... I was side-tracked by looking into the results that the US Navy - AAD web site gives. I experimented by entering Longitudes that are not 0°E. I can reproduce that results as long as the Time Zone remains at Zero/Zulu. I cannot determine if they use the local elevation (like topos = wgs84.latlon(lat, lon)) or if they employ zero elevation (like topos = wgs84.latlon(lat, lon, elevation_m=0.0)) I expected to see a checkbox where one can select zero elevation. If using a longitude that is not a multiple of 15° I guess it fits the result into the classical 15° time zone boundaries. Anyway ... this is a side-track only.

[aendie]

@EndlessDex Also in answer to your request #1: I checked all moonrise/moonset data for 2023 at latitude 68°N....

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

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

All discrepancies in 2023 at latitude 68°N: 2023-01-10 moonset 11:57 (NA), 11:58 (AAD) 2023-01-12 moonset 11:15 (NA), 11:16 (AAD) 2023-01-16 moonrise 03:19 (NA), 03:18 (AAD) 2023-01-17 moonset 09:13 (NA), 09:14 (AAD) 2023-02-11 moonset 08:39 (NA), 08:40 (AAD) 2023-02-23 moonset 22:43 (NA), 22:44 (AAD) 2023-02-25 moonset 00:49 (NA), 00:50 (AAD) 2023-02-26 moonset 03:16 (NA), 03:17 (AAD) 2023-03-13 moonset 04:32 (NA), 04:33 (AAD) 2023-05-26 moonset 03:33 (NA), 03:34 (AAD) 2023-06-22 moonset 02:10 (NA), 02:11 (AAD) 2023-07-12 moonset 18:23 (NA), 18:24 (AAD) 2023-08-16 moonrise 01:57 (NA), 01:56 (AAD) 2023-08-30 moonrise 20:27 (NA), 20:26 (AAD) 2023-09-02 moonrise 19:18 (NA), 19:17 (AAD) 2023-09-04 moonrise 18:21 (NA), 18:20 (AAD) 2023-09-11 moonrise 23:12 (NA), 23:11 (AAD) 2023-09-25 moonrise 20:30 (NA), 20:29 (AAD) 2023-09-25 moonset 21:29 (NA), 21:30 (AAD) 2023-09-30 moonrise 17:18 (NA), 17:17 (AAD) 2023-10-23 moonset 22:05 (NA), 22:06 (AAD) 2023-11-09 moonrise 02:08 (NA), 02:07 (AAD) 2023-11-22 moonrise 14:22 (NA), 14:21 (AAD) 2023-12-16 moonrise 14:51 (NA), 14:50 (AAD) 2023-12-21 moonrise 11:58 (NA), 11:57 (AAD)

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

[brandon-rhodes]

Happy first day of September!

I’ve generated a few diagrams to help us think through how the search routine works. If you want to see the (work-in-progress) code that produced the diagrams, the script is here:

https://github.com/skyfielders/python-skyfield/blob/investigate-moon-rising/design/illustrate_moonrise_edge_case.py

To keep things simple, I’m going to stick with the Moon as the primary target, and (70°N,0°E) as the location of the observer. (That will let me say simple things like ‘south horizon’, instead of complicated things like ‘the opposite horizon from the hemisphere the observer is in’!)

So let’s dive in.

The search routine’s normal approach is:

1. Generate the target’s position at an arbitrary time to learn its current latitude and HA (hour angle).
2. Use a fast geometry routine to compute the HA at which that line of latitude intersects the horizon.
3. Take the difference between the first HA and the second HA, divide it by the average rate at which the sky rotates, and jump directly to the moment of rising.
4. Just in case, repeat the above steps two more times (in a few paragraphs, we’ll talk about why).

This works stunningly well for the fixed stars. Let’s imagine five stars that range in declination from -20.8° to -21.0°. If we ask for when, say, the -20.85° star rises:

1. The search routine starts Iteration 0, and generates the star’s HA and declination for the starting time of 2023-02-19 00:00:00 UTC. Its starting HA will be something crazy, but we only care about its declination.

2. We predict that, given its declination, the star should rise at 2023-02-19 11:18:57 UTC

3. And it does! The search routine generates a position at that predicted time. It finds the star right where it expects, at Point 1 in the diagram!

4. So the final two iterations produce Point 2 and Point 3 that remain on the same spot.

Here’s a diagram of that process, that leaves out Point 0 (which is way over in some other part of the sky) but shows that the final three iterations all land at exactly the point where the star rises.

If the moment of rising can be found immediately using geometry, why go through with the expense of generating Point 2 and Point 3? The reason is that the rising routine might have been given a moving target, like a planet, instead of a fixed star. In that case the declination will change slightly between the arbitrary 2023-02-19 00:00:00 UTC time of the first data point and the time we computed for rising — thereby throwing off slightly the answer we got from the fast geometry routine.

So Iteration 2 takes the new declination value, and uses it to compute an even better geometric prediction. Repeating that again in Iteration 3 produces excellent results for all the cases that I had tested so far.

But let’s now consider a special case.

But what if the geometry simply doesn’t predict a moment of rising?

What about a star whose declination is so negative that it never rises for our observer at 70°N?

In this case there’s unfortunately no solution to the geometric equation that would provide an intercept between the star’s circle of declination and the horizon. So the search routine does something different in step 2: it falls back to instead computing the moment of transit. Why fall back to the moment of transit?

1. One answer is that the moment of transit should be the moment of the star’s closest approach to the horizon. The caller, if they were interested, could use the time of transit to compute how close the star came closest to rising.

2. But there’s another reason to look for transit: even if a rising is impossible given the target’s initial declination, it’s just barely possible that by the moment of transit, the declination will have shifted enough that the target can peek above the southern horizon.

In the case of our fictional star, whose declination does not, in fact, change over the hours, the search routine immediately finds the moment of transit at Point 1, and Points 2 and 3 simply dwell on the same point:

But for the Moon, there’s trouble!

The Moon moves more rapidly against the background of stars than any of the planets, and thus poses the biggest challenge for the rising and setting routines. And its rising on 2023-02-19 from 70°N is a severe enough event that our search routine, alas, gets stuck!

And it gets stuck in a fairly stable way. Here I’ve expanded the number of search iterations so they go all the way to Iteration 9:

What story is this diagram telling us?

• To start with, let’s note that Point 1 is down to the lower left, far away from the rest of the diagram. Even though the search routine noticed that the Moon is too far south to rise, and aimed at the moment of transit, its first guess fell around 5° short of that goal. Which is to be expected — the Moon moves something like 13° across the sky each day, so the Moon has had the chance to move quite a bit over the eleven hours between 2023-02-19 00:00:00 UTC and this first guess of 2023-02-19 11:09:36 UTC for the moment of transit.

• But the next iteration, which brings us to Point 2, is in one sense a success. The geometry routine again sees that the Moon is too far south, so aims merely for the moment of transit. This time it’s a smaller jump of only 25 minutes, to 2023-02-19 11:34:43 UTC, so the Moon has less time to move, and so Point 2 is much closer to the moment of transit.

• But in another sense Point 2 creates a problem — the Moon is now above the horizon! Suddenly, the geometry routine stops asking for the moment of transit at HA=0.0 and instead can compute a real horizon intercept. Since the Moon is a bit below the -20.9° curve, the geometry routine eyeballs an intercept with the horizon line that’s a bit to the right of the intercept of the -20.9° curve, and chooses 2023-02-19 11:21:43 UTC as the likely moment of moonrise.

• Oh, no! By rewinding the clock back 13 minutes, we have given the Moon enough time to lose a few fractions of a degree of declination, which places it so far back south that it no longer looks like it will reach the horizon at all. The geometry routine, given a declination from which moonrise isn’t even possible, falls back to targeting the moment of transit instead.

• And we’re now stuck. We seem to have found a stable ‘attractor’ that consists of two tight regions, one below the horizon and one above the horizon, between which the search routine will merrily jump forever.

Three quick points:

First, wow, these angles are very small! But they seem to matter. Even though the Moon only rises 0.1° above the horizon, it stays above the horizon for a full 37 minutes! That wasn’t my intuition when I first approached the problem — that such a low rising angle and altitude could produce such a long event. I had assumed that I could ignore tiny-angle cases because they would put a target above the horizon for only a few dozen seconds. But that’s not the case, and we can’t be ignoring events that last 37 minutes, so I agree that this needs to be fixed.

Second, @aendie is particularly unhappy that risings were being missed. It turns out that this is an artifact: it’s the mere result of my happening to choose an odd number of iterations for the search routine loop! Because I stopped after an odd number iterations, I wound up always returning a False point from the lower-left cluster of the attractor. If I had chosen an even number of iterations, like 4, 6, or 8, then I would at least be returning True — though, of course, with a wrong ‘rising time’ that was in fact close to the time of transit.

Third, my vague idea that every possible edge case would somehow just magically resolve itself if I did a few more iterations is clearly false. This case shows that the search routine can become stuck in a very stable behavior from which more iterations will never rescue it.

So some final maneuver is clearly needed here to break out of this cycle — at the very least, one final step that would take the results of the previous two iterations, and jump intelligently between them to reach the horizon.

I’ve made one experiment in that direction already, but I’m not happy with the results, so I’ll play with some more ideas and see what I can come up with.

Finally —

With all of the above as background, I think I can sweep through the rest of the unanswered comments from @EndlessDex and try to provide some kind of response:

Because the function used with Newton's Method is not continuous the solution does not converge properly.

1. I suspect you just typed the wrong word here — in your graph, the function is clearly continuous? I’ll assume that what you mean is that the function doesn’t have a continuous first derivative.

2. In the case of a non-continue first derivative, Newton’s Method would indeed have all kinds of problems. But to guard against that, you’ll see upon re-reading that the search routine doesn’t, in fact, use Newton’s Method — it never tries to compute a slope or rate. Instead, it uses the plain old average speed of the sky’s rotation (360° per day) as its guess about how fast any target will move across the heavens.

This leads to the newtons method not converging correctly because it's trying to develop a slope for a non-linear function.

I’m not sure this quite describes the problem.

First, as noted above, in no case does the search routine try to compute a slope.

Second, even if it did compute a slope, the whole point of Newton’s Method is to explore non-linear functions. If a function is linear, can’t we simply compute its zero-crossing directly, using its slope? My understanding is that Newton’s Method is always used with non-linear functions. In which case, that by itself can’t be the problem here.

I would describe the problem as: the search routine at each step treats the declination as a constant, instead of trying to model it as a moving target with its own rate of motion.

But feel free to follow up with a clarification; maybe you’re getting at something interesting here that I’m not quite catching yet.

It gets partway there by pure luck due to one of the sample points being near the rise time.

I’m not sure luck is involved?

I’m not quite sure what you’re claiming, so please show your work and correct me here if I’m wrong, but as the diagrams above show, when the Moon’s declination is too low for it to rise, the search routine’s behavior is robust and doesn’t depend on the particular initial sample point. In all cases, if the declination is too negative for a moonrise at 70°N, the search routine immediately decides to rotate the sky all the way to the moment of transit to see whether the Moon is above the horizon at transit or not.

In which case the search routine’s behavior will enter the stable two-region attractor, I think? And so the result won’t ever depend on luck?

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 guess is that we shouldn’t pick a single location.

I have, alas, had unfortunate experiences before, where the very tweaks that happen to make a routine work perfectly at one location (like 70°N) promptly wind up breaking the routine for other latitudes — or even for all other latitudes!

This is a specific case of a more general rule with software, where fixing a bug for one specific case often breaks the software for all other cases.

If we want a robust rise-set solution, we should probably generate USNO moonrise and moonset tables for every latitude, and tweak the search routine until it gives good results for all of them.

It would be slow to run so many checks as part of the normal test suite, so maybe the thorough rise-set check could live inside a separate little script that we could run by hand whenever the search routine is being tweaked?

Then we can take each specific edge case that’s uncovered, and promote those to full membership in the normal test suite.

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 … I think the key point for me is that the hour angle equation breaks down for certain times/locations.

Fleshing out your idea a little bit, it sounds like your proposal is to:

1. Go ahead with the first couple of steps of the search routine.

2. Detect, and treat as a special case, the situation where (a) the Moon’s declination is at first so negative that we despair of finding a horizon intercept, (b) but when we look for the moment of transit we notice that the Moon managed to rise anyway.

3. In that special case, launch an expensive-but-reliable find_discrete() binary search for the actual moment of moonrise.

Is that a fair description of the approach you’re thinking of?

Thanks, by the way, for being willing to do such a deep dive on this issue as a first-time contributor! I’ll be interested to hear how your experiments go, and I’ll let you know how mine turn out.

[EndlessDex]

Skyfield uses this USNO generator as a source of truth for comparison against. This table uses the Refraction + Apparent Radius method. … Generating a table for (W0, N70) shows that both old and new methods are grossly wrong (as you have noted).

@EndlessDex — could you share the table? I would have expected the old method to agree completely with the USNO, and if it doesn't, then that's a separate issue that will need to be fixed. Working code to generate the table would also be helpful. Thanks!

I was referencing the output of the design/test_sunrise_moonrise.py test. shows a minute off here and there for Bluffton. And gross errors when I run this same test with the 2023 latitude 70 table. code here

Though for the find_discrete method this failure is likely due to under-sampling.

[aendie]

@brandon-rhodes Using is_moon_up_at.step_days = 1.0 / 24.0 / 60.0 ("New:" below) as opposed to is_moon_up_at.rough_period = 0.5 ("Old:" below) in almanac.find_discrete() did bring some improvements when checking moonrise/moonset in the year 2023.

Most agreed with the HMNAO Nautical Almanac, however I have highlighted those that didn't ("NA:" below) in the list below. So I found some cases where find_discrete() differs in the minute (of time) value, though the actual discrepancy might only be a few seconds that kick it over the minute mark. Thanks for your suggestion to drop .rough_period.

2023-04-14 Lat 66°N Old: moon all day below horizon New: moonrise at 07:02; moonset at 07:34

2023-07-05 Lat 66°N Old: moon all day below horizon New: moonrise at 01:46; moonset at 01:57

2023-07-27 Lat 68°N Old: moon all day below horizon New: moonrise at 18:51; moonset at 19:25

2023-07-30 Lat 62°N Old : moon all day below horizon New: moonrise at 21:56; moonset at 22:43

2023-10-18 Lat 64°N Old: moon all day below horizon New: moonrise at 14:15; moonset at 15:06

2023-10-19 Lat 62°N Old: moon all day below horizon New: moonrise at 15:34; moonset at 15:47 NA: moonrise at 15:35; moonset at 15:46

2023-12-14 Lat 62°N Old: : moon all day below horizon New: moonrise at 13:09; moonset at 13:49 NA: moonrise at 13:09; moonset at 13:48

2023-12-26 Lat 62°N Old: moon all day above horizon New: moonset at 11:08; moonrise at 11:54

P.S. All times in the Nautical Almanac tables I create are rounded to minutes if I print HH:MM or rounded to seconds if I print HH:MM:SS

[EndlessDex]

I suspect you just typed the wrong word here — in your graph, the function is clearly continuous? I’ll assume that what you mean is that the function doesn’t have a continuous first derivative.

Yes, you are correct about my imprecise language above. When referring to "non-continuous" or "non-linear" I was intending to reference the piecewise-like nature of the desired hour angle function where the value is 0 for certain inputs that trigger the clamping while it is a smooth, continuous function for other inputs. I believe the proper mathematical phrasing should been "The desired hour angle function is not a SMOOTH function."

In the case of a non-continue first derivative, Newton’s Method would indeed have all kinds of problems. But to guard against that, you’ll see upon re-reading that the search routine doesn’t, in fact, use Newton’s Method — it never tries to compute a slope or rate. Instead, it uses the plain old average speed of the sky’s rotation (360° per day) as its guess about how fast any target will move across the heavens.

Yes, you are not technically using Newton's method mathematically (except sortof for the first point). It is however similar/equivalent. Your method iteratively approaches a more accurate solution using error and rate of change. Your method IS using a (act-desired)/rate = step format even if the rate here is fixed and not a derivative. So I would say semantics.

I would describe the problem as: the search routine at each step treats the declination as a constant, instead of trying to model it as a moving target with its own rate of motion.

But feel free to follow up with a clarification; maybe you’re getting at something interesting here that I’m not quite catching yet.

Two sides of the same coin I think. The routine treats declination as constant and because the desired hour angle function is non-smooth the model does not converge as the changing declination pushes us from one region to another of the desired hour angle function.

Plot of Desired Hour Angle with Constant Declination


We bounce back and forth over the point/cusp of the desired hour angle function. Which causes the guess not increase in accuracy with each iteration.

But your description does lead more directly to a possible solution. Pre-correcting the step size based on predicted change in declination. My description only leads to the analytical geometry problem of smoothing out the discontinuity from the transformation functions.

I’m not sure luck is involved?

I’m not quite sure what you’re claiming, so please show your work and correct me here if I’m wrong, but as the diagrams above show, when the Moon’s declination is too low for it to rise, the search routine’s behavior is robust and doesn’t depend on the particular initial sample point. In all cases, if the declination is too negative for a moonrise at 70°N, the search routine immediately decides to rotate the sky all the way to the moment of transit to see whether the Moon is above the horizon at transit or not.

This mention of luck came from a misunderstanding on my part about the transit time fallback and how this will always be "best case". I did not know at the time that this was guaranteed.

My guess is that we shouldn’t pick a single location.

I have, alas, had unfortunate experiences before, where the very tweaks that happen to make a routine work perfectly at one location (like 70°N) promptly wind up breaking the routine for other latitudes — or even for all other latitudes!

This is a specific case of a more general rule with software, where fixing a bug for one specific case often breaks the software for all other cases.

I was of course not suggesting to only ever test one location. That would be stupid. However, for the purpose of developing problem analyses and solution proposals, it makes sense to focus on a single error case initially. Once a problem is somewhat understood then you can start exploring how it behaves under different circumstances.

It would be slow to run so many checks as part of the normal test suite, so maybe the thorough rise-set check could live inside a separate little script that we could run by hand whenever the search routine is being tweaked?

Each location takes 1-2 seconds for each table so definitely agreed. I think this already partially exists? design/test_sunrise_moonrise.py

Fleshing out your idea a little bit, it sounds like your proposal is to:

Go ahead with the first couple of steps of the search routine.

Detect, and treat as a special case, the situation where (a) the Moon’s declination is at first so negative that we despair of finding a horizon intercept, (b) but when we look for the moment of transit we notice that the Moon managed to rise anyway.

In that special case, launch an expensive-but-reliable find_discrete() binary search for the actual moment of moonrise.

Basically yes. If the desired_hour_angle function would ever return zero and then non-zero (or vice versa) then we are in the problematic edge case. But I am hoping to find an algorithmic solution instead of having to fall back onfind_discrete.

Thanks, by the way, for being willing to do such a deep dive on this issue as a first-time contributor! I’ll be interested to hear how your experiments go, and I’ll let you know how mine turn out.

You're very welcome! My work deals heavily with programming for satellite propagation and I have used skyfield as an easy server-side TLE propagator for a long time. I was browsing your issues page for my own bug and came across something I could answer. Since then I have been on a bit of a kick haha. Not as familiar with astronomical propagation as I am orbital propagation but it's all very interesting to me. Hope to become a long term contributor actually. Please let me know how I can better help out 😄

(PS take a look at my PR. I really want to be able to use python 3.12 haha)

[brandon-rhodes]

I was referencing the output of the design/test_sunrise_moonrise.py test. shows a minute off here and there for Bluffton. And gross errors when I run this same test with the 2023 latitude 70 table.

Ah, thank you for pointing out the specific script!

Looking back at it, its reporting and logic were a bit fragile, so I've improved them so that missing events, whether on the Skyfield or USNO side, no longer cause it to halt. Instead, it always finishes, having checked every event reported by either library.

Also, it was using time zone of UTC-7, which I think is the source of all the errors you were seeing, since the new table uses UTC instead. Here's an updated copy of the script, where I started with your own version of the script before making any improvements:

https://github.com/EndlessDex/python-skyfield/blob/setting-hour-angle-update/design/test_sunrise_moonrise.py

The script now shows only 1-minute or 2-minute differences of opinion for 2023 moon risings at 70°N using find_discrete():

. . 1 . 1 . . . . . . . . . . . 1 . 1 1 . . -2 . . . . . . 1 . . . . 1 . . . .
. . . . . 1 . . . . . . . . . . . 1 . 1 . . . 1 . . . . . 1 . 2 . . . . . 1 .
. . . . . . 1 . . . . . -1 . . . . . 1 . . . 1 . . . . . . . . 1 1 1 . . 1 . .
. . . . . . . 1 1 . 1 1 1 . -1 . . . . . 1 . . . . . . . . . . . . . 1 1 1 1 1
. . . . . . . . 1 1 . 1 . . . -1 . . . . . 1 .

Which I'm not entirely happy with — I want everything to be within 1 minute of the USNO, or, ideally, a perfect match — but it's at least not the gross errors you feared find_discete() was running into!

And when using find_risings(), it complains about the same missing event that @aendie reported that started this thread:

Skyfield is missing: 2 19 (11, 28)
. . . . . . . . . . . . . . . . . . . . . . skyfield-miss . . . . . . 1 . . .
. . . . . . . . . . . . . . . . . . . . . . . . . 1 . . . . . . . . . 1 . . .
-1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . 1 . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . -1 . . . . . 1 .

So I'm glad to find that find_discrete() is still bulletproof, if given enough starting points that it can find all the events.

[EndlessDex]

@aendie

So I found some cases where find_discrete() differs in the minute (of time) value, though the actual discrepancy might only be a few seconds that kick it over the minute mark.

Try rounding skyfield's output to the nearest minute. It's what the tests do.

How to round skyfield times

```python import skyfield.api as sf def round_to_nearest_minute(sf_time): THIRTY_SECONDS = 1.0 / 24.0 / 60.0 / 2.0 sf_time += THIRTY_SECONDS calendar_tuple = sf_time.utc sf_time = sf_time.ts.utc( calendar_tuple.year, calendar_tuple.month, calendar_tuple.day, calendar_tuple.hour, calendar_tuple.minute, ) return sf_time ts = sf.load.timescale() t = ts.utc(2023, 1, 19, 12, 3, range(0, 60, 15)) print("original", t.utc_iso()) print("rounded ", round_to_nearest_minute(t).utc_iso()) # original ['2023-01-19T12:03:00Z', '2023-01-19T12:03:15Z', '2023-01-19T12:03:30Z', '2023-01-19T12:03:45Z'] # rounded ['2023-01-19T12:03:00Z', '2023-01-19T12:03:00Z', '2023-01-19T12:04:00Z', '2023-01-19T12:04:00Z'] ```

[brandon-rhodes]

(@EndlessDex — Whoops — I just pushed another update to the script, because I had commented out the thirty-second offset to perform rounding, but forgot to uncomment it before committing. I've now updated my previous comment with the result.)

[brandon-rhodes]

@EndlessDex — Oh, I think I have a guess as to why find_discrete() is so much worse than the new routine at 70° N as shown in the output above: it's because it's not using the Moon's dynamic radius, but rather a static estimate for it.

[EndlessDex]

Also, it was using time zone of UTC-7, which I think is the source of all the errors you were seeing, since the new table uses UTC instead. Here's an updated copy of the script, where I started with your own version of the script before making any improvements:

Yes that was a stupid gaff. I found it yesterday and hadn't pushed it yet. Merged your stuff into mine. You missed one timezone value btw on line 236. Either way, it wasn't until I added a f.step_days that I got a well matching solution. As high as 1/24 gets within 2 minutes of truth.

Ooo interesting bit about the apparent radius calc. Would you want to update risings_and_settings to do a similar decision tree to _find for better defaults? Or possible expose the apparent horizon/radius calc functions to give a standard method to users?

[aendie]

@brandon-rhodes

Oh, I think I have a guess as to why find_discrete() is so much worse than the new routine at 70° N as shown in the output above: it's because it's not using the Moon's dynamic radius, but rather a static estimate for it.

I'm not sure what you mean here. My Test Program included in this ISSUE calculates the horizon value to pass on to find_discrete() based on the moon's volumetric mean radius taking it's distance into account plus 34' for atmospheric refraction. I see this as the user's responsibility to pass this on to find_discrete() otherwise it will use the standard value 0.83333 degrees (16' + 34'). Are you referring to another use of find_discrete()?

[EndlessDex]

@brandon-rhodes This works in all cases I've tested but I hate it :( Feel like there should be a mathematical/algorithmic solution but I haven't found one yet. Took a break from that to whip this up.

Regarding the code itself risings_and_settings can't be used due to lack of eph and find_discrete doesn't seem to accept Date Arrays so that has to be serialized. Also not sure on the search time. I landed on +/- 10min and had no issues but that could be a weak point. Ignoring risings/settings since we do not know which we are looking for (also a weak point if the events are very close).

def _find
...
...
    # TODO: How many iterations do we need?  And can we cut down on that
    # number if we use velocity intelligently?  For now, we experiment
    # using the ./design/test_sunrise_moonrise.py script in the
    # repository, that checks both the old Skyfiled routines and this
    # new one against the USNO.  It suggests that 3 iterations is enough
    # for the Moon, the fastest-moving Solar System object, to match.
    last_desired_ha = None
    fallback_mask = np.full_like(t, False)
    for i in 0, 1, 2:
        _fastify(t)
        ha, dec, distance = observer.at(t).observe(target).apparent().hadec()
        desired_ha = f(latitude, dec, h(distance))
        if last_desired_ha is not None:
            # One is zero and one is non zero, ie we are at the edgecase so add to list to redo
            fallback_mask = np.logical_or(fallback_mask, np.logical_xor(last_desired_ha, desired_ha))
        last_desired_ha = desired_ha

        ha_adjustment = desired_ha - ha.radians
        ha_adjustment = (ha_adjustment + pi) % tau - pi
        timebump = ha_adjustment / ha_per_day
        t = ts.tt_jd(t.whole, t.tt_fraction + timebump)

    # find_discrete cannot take time arrays :(
    fallback_times = np.empty_like(fallback_mask)
    for i, f in enumerate(fallback_mask):
        if f:
            fallback_t = t[i]
        else:
            fallback_times[i] = None
            continue
        def is_body_up_at(t_search):
            _fastify(t_search)
            alt, _, dist = observer.at(t_search).observe(target).apparent().altaz()
            return alt.radians > h(dist)
        is_body_up_at.step_days = 1 / 1440.0
        event, _ = find_discrete(fallback_t - (10 / 1440), fallback_t + (10 / 1440), is_body_up_at)
        fallback_times[i] = event[0]

    # Combine t and fallback_times

[aendie]

Initial thoughts from my side (of the 'puddle')...

Stars remain almost constant in the celestial sphere, which can be represented as a dot on a graph of Declination (degrees) versus Sidereal Hour Angle (degrees). This I created long ago:

Well, what does the Moon do, say in February 2023? It moves thus:

I labelled each point with the date (at 00h). The red line is the ecliptic. It more-or-less follows the ecliptic. The Moon's orbital plane is inclined about 5.1° with respect to the ecliptic plane (as seen above). Each calendar month is similar - only the starting position (first day of the month) moves around. The direction of movement is always decreasing SHA (right-to-left). So the Moon's SHA tells us immediately if its Declination is increasing, decreasing or nearly flat at SHA 90° and 270°.

I'm no further than this. It took me 2 days fighting with matplotlib.pyplot.annotate to generate labels with more-or-less equal length arrows :-( But a picture is always great ... methinks.

UPDATE: The matplotlib.pyplot.annotate problem is solved: https://stackoverflow.com/questions/78960741/offset-points-in-matplotlib-pyplot-annotate-gives-unexpected-results