Almanac.find_discrete returning inconsistent results

[pczhang-tony]

Hi Brandon,

After re-installing Anaconda on my PC, I found that calculating the moon phases using the almanac.find_discrete routine gives inconsistent results.

I'm using an excerpt of the DE431 ephemeris, Skyfield 1.16, Spyder 4.0.1, Python 3.7.6, Numpy 1.18.1.

For a test case, I calculated the moon phases between -1150-01-01T00:00:00 and -1149-01-01T00:00:00, in both the proleptic Gregorian calendar and the proleptic Julian calendar.

Case one, Gregorian calendar

T0 = Timescale.ut1_jd(1301031.5)
T1 = Timescale.ut1_jd(1301396.5)
t,y = almanac.find_discrete(T0,T1,almanac.moon_phases(e))
ans = t[y==1] # Times of the first quarter moon

ans.ut1 has the following values:
0 1301046.12572536
1 1301075.42562022
2 1301104.65598818
3 1301133.87319827
4 1301163.13824748
5 1301192.50759745
6 1301222.02205536
7 1301251.69407743
8 1301281.49975272
9 1301311.38248919
10 1301341.26599025
11 1301371.06986399
12 1301378.01611522

If converted to the proleptic Julian calendar, the dates are:
0 -1150-01-26T15:01:03
1 -1150-02-24T22:12:54
2 -1150-03-26T03:44:37
3 -1150-04-24T08:57:24
4 -1150-05-23T15:19:05
5 -1150-06-22T00:10:56
6 -1150-07-21T12:31:46
7 -1150-08-20T04:39:28
8 -1150-09-18T23:59:39
9 -1150-10-18T21:10:47
10 -1150-11-17T18:23:02
11 -1150-12-17T13:40:36
12 -1150-12-24T12:23:12

From the results above it is already clear that ans[11] and ans[12] cannot both be correct, as they are only seven days apart. Setting this problem aside, let us examine case two.

Case two, Julian calendar

T0 = Timescale.ut1_jd(1301020.5)
T1 = Timescale.ut1_jd(1301385.5)
t,y = almanac.find_discrete(T0,T1,almanac.moon_phases(e))
ans2 = t[y==1] # Times of the first quarter moon

ans2.ut1 has the following values:
0 1301046.12572536
1 1301075.42562022
2 1301082.25535569
3 1301104.65598818
4 1301133.87319827
5 1301163.13824748
6 1301192.50759745
7 1301222.02205536
8 1301251.69407743
9 1301281.49975272
10 1301311.38248919
11 1301341.26599025
12 1301371.06986399
13 1301378.01611522

If converted to the proleptic Julian calendar, the dates are:
0 -1150-01-26T15:01:03
1 -1150-02-24T22:12:54
2 -1150-03-03T18:07:43
3 -1150-03-26T03:44:37
4 -1150-04-24T08:57:24
5 -1150-05-23T15:19:05
6 -1150-06-22T00:10:56
7 -1150-07-21T12:31:46
8 -1150-08-20T04:39:28
9 -1150-09-18T23:59:39
10 -1150-10-18T21:10:47
11 -1150-11-17T18:23:02
12 -1150-12-17T13:40:36
13 -1150-12-24T12:23:12

Compared to ans, ans2 has one extra value (-1150-03-03T18:07:43), which cannot be a true solution as it is only about seven days after -1150-02-24T22:12:54.

In the two examples above, the start and end time points are offset by only 11 days, basically covering the same time period, but produces different results.

I never had this issue before. I don't know whether it's an issue with the new version of Skyfield, or is it something about Numpy, Spyder, or Python 3.7.6. Can anyone else reproduce the issue? I hope it isn't something peculiar about my computer...

[brandon-rhodes]

I had not seen this issue either! But I can confirm that it happens to me as well. It looks like the search is so exact that it is finding moments when the various motions involved — the Earth's, the planet — are "ticking" at different rates, are ratcheting forward at different moments, so that the following sequence occurs:

1. It's First Quarter.
2. One of the bodies moves forward, making it Full Moon.
3. But the other body, whose motion opposes slightly the achievement, ratchets and now it's back to First Quarter because the angle is slightly shy of Full Moon.
4. The first body ratchets forward again and it's now Full Moon again.

I have, indeed, here made a fundamental error: just because I stop searching once the answer is within epsilon of correct doesn't mean that I don't accidentally, sometimes, find that moment of backwards-ratcheting if one of my data points happens accidentally to straddle it.

Interestingly, I just fixed this in a maxima-finding routine in almanac.py. I'll go in and apply the same logic to discrete searching.

Thanks for letting me know about this — I'll let you know when I have a fixed version to try out! Here's the code from the maxima-finding routine that I'll try moving over to the discrete-finding routine:

            mask = concatenate(((True,), diff(ends) > epsilon))
            ends = ends[mask]
            y = y[mask]

[pczhang-tony]

What's really weird about this is I've been using the same code for over a year, and I only had this issue starting yesterday. I tried downgrading to Skyfield 1.10, numpy 1.15.1, Python 3.6, to try to see if using earlier versions fixes the issue, but so far no luck.

[brandon-rhodes]

Earlier versions, alas, will not help. The phenomenon should be very sensitive to exactly the dates you search for, and whether those dates happen to land one of their test times at exactly the moment that a jagged curve happens to pass a boundary twice.

[pczhang-tony]

Earlier versions, alas, will not help. The phenomenon should be very sensitive to exactly the dates you search for, and whether those dates happen to land one of their test times at exactly the moment that a jagged curve happens to pass a boundary twice.

Hi Brandon, I was able to resolve the problem by downgrading:

jplephem 2.9
skyfield 1.10
spyder 3.3.6
numpy 1.16.5
python 3.7.4

Part of the issue might be in jplephem. It seems that my code broke after updating jplephem to 2.12. Did somethings change in the latest update?

[brandon-rhodes]

Yes, jplephem has changes with every release. Here are the changes between those two releases:

brandon-rhodes/python-jplephem@v2.9...v2.12

I made several tweaks to make the library produce numbers closer (at very high precision) to those produced by NASA's SPICE library. It could be that slight difference in the last decimal place of each number that randomly triggers an extra data point in one case but not another.

[pczhang-tony]

I tried the two test cases again under the downgraded conditions, and the results are consistent:

The universal times are:
1301046.12572536
1301075.42562022
1301104.65598818
1301133.87319827
1301163.13824748
1301192.50759745
1301222.02205536
1301251.69407743
1301281.49975272
1301311.38248919
1301341.26599025
1301371.06986399

Corresponding to the dates in the proleptic Julian calendar:
-1150-01-26T15:01:03
-1150-02-24T22:12:54
-1150-03-26T03:44:37
-1150-04-24T08:57:24
-1150-05-23T15:19:05
-1150-06-22T00:10:56
-1150-07-21T12:31:46
-1150-08-20T04:39:28
-1150-09-18T23:59:39
-1150-10-18T21:10:47
-1150-11-17T18:23:02
-1150-12-17T13:40:36

I suspect something in jplephem 2.12 must have triggered a change in behavior of skyfield, hopefully it will be helpful for your troubleshooting.

Thanks for looking into this! For the time being, I'll stick with my current environment.

[brandon-rhodes]

There — I was able to reproduce the behavior in a test case and then, hopefully, resolve it. I plan to release a new Skyfield version later today, so trying again tomorrow should give you a good shot at seeing this behavior go away.

Thanks again for reporting it!