Lower Transit in 2019 agrees with Ephem except on 12.01.2019 and 03.04.2019

I am using Skyfield 1.11 and Ephem 3.7.6 with Python 2 mostly. I have produced a Nautical Almanac that agrees or is better than the Pyephem version. However I have a curious error on two dates only within 2019 where the Lower Transit is way off. I cannot explain it and I cannot find a tweak to resolve it.

I have placed a subset of my code in a test file with just four test dates, two of which are incorrect. The result of running it is:

Lower transit on 11.1.2019  from  ephem : 03:53
Lower transit on 11.1.2019 from skyfield: 03:53

Lower transit on 12.1.2019  from  ephem : 04:35
Lower transit on 12.1.2019 from skyfield: 04:16

Lower transit on 13.1.2019  from  ephem : 05:17
Lower transit on 13.1.2019 from skyfield: 05:17

Lower transit on 3.4.2019  from  ephem : 23:12
Lower transit on 3.4.2019 from skyfield: 23:03

A word of explanation: I optimise my search in find_discrete to within one hour (as I already know in which hour the Moon GHA crosses 180°). In order to simplify the supplied test code I have inserted fixed constants for the start and end search hour. The Upper Transit calculation works in my Almanac without error - so you may question my method of implementing Lower Transit.

I would be honoured if you could look at this issue. I am otherwise thinking of implementing my own search for the 180 degree Moon GHA crossover point, which only has to be accurate to the nearest minute (and so will probably be faster).

Last but not least, my sincere thanks for developing Ephem and Skyfield. It's much appreciated!

Andrew Test_01.zip

My first thought would be: is there an official source from which you could learn the real moment of antitransit? Without knowing which library to blame, it's difficult to know which one to spend further time investigating for misbehavior. Could you enter those exact moments above in to HORIZONS and see which times really correspond to anti-transit and which don't? Another check would be: what azimuth Skyfield's altaz() method gives at each moment, to verify that the azimuth is 0° or 180°.

Never having used JPL HORIZONS, I don't see a way to output the Moon GHA for comparison purposes.

I consider the discrepancies I found as follows... every day of the year 2019 agrees perfectly between Ephem and Skyfield except for two dates. I doubt that there's a problem with Skyfield. The problem probably lies with the fact that I call my method "lowertransit" with an argument of 12 or -12. This seems to work in 99.5% of the cases but I suspect that there's a boundary case where it doesn't. Put in other words: I haven't found a perfect solution for the "lowertransit" method. On the other hand, my "uppertransit" method (which takes no argument) works 100% perfectly (so far):

def uppertransit():
    # Build a function of time that returns the moon upper transit time.

    def is_moon_transit_at(t):
        """The function that this returns will expect a single argument that is a 
        :class:`~skyfield.timelib.Time` and will return ``True`` if the moon is up
        or twilight has started, else ``False``."""
        t._nutation_angles = iau2000b(t.tt)
        # Return `True` if the meridian is crossed by time `t`.
        position = earth.at(t).observe(moon)
        ra = position.apparent().radec(epoch='date')[0]
        return t.gast > ra.hours    # works!

    is_moon_transit_at.rough_period = 0.01  # search increment hint
    return is_moon_transit_at

Note also that the choice of the "search increment hint" was pure trial and error on my part.

I would be interested to know if you can code a "lowertransit" method that works correctly. Besides that I think (or hope) that searching for the minute of day that's closest to 180° GHA will provide me with the correct Lower Transit time. (I'm not an astronomer.) Just as searching for the minute of day that's closest to 0° GHA might provide me with the correct Upper Transit time.

Comparing Skyfield and Pyephem as Nautical Almanac generators, I see only two stumbling blocks:

lack of planet magnitudes
lack of a Moon Lower Transit method that's reliable

Regarding performance - a big advantage would be to be able to calculate values faster with lower accuracy, e.g. if Moon Transit times in hours : minutes are sufficient. In Pyephem, I see that mooncolong.c has simplified algorithms compared to those in Astronomical Algorithms, 2nd. Ed. 1998 by Jean Meeus. I am no judge, but I guess this was done to increase performance with minimal loss of accuracy.

Before I dig deeper into what's really going on... here I include a minor (but necessary) correction to Test_01.

I've also added Test_02. This scans all dates from 2010 to 2019 inclusive and displays the Moon Lower Transit inconsistencies (comparing only hours and minutes, rounded) between Ephem and Skyfield using my (imperfect) method for Moon Lower Transit calculation. A discrepancy of 1 minute may, of course, be no error at all. But anything larger is definitely suspicious. The discrepancies found in 10 years are:

Lower transit on 20.06.2010  ephem: 06:46  skyfield: 06:47
Lower transit on 13.08.2010  ephem: 02:35  skyfield: --:--
Lower transit on 06.10.2010  ephem: 22:57  skyfield: --:--
Lower transit on 16.04.2011  ephem: 10:23  skyfield: --:--
Lower transit on 24.09.2011  ephem: 21:31  skyfield: 21:32
Lower transit on 05.12.2011  ephem: 07:52  skyfield: 07:53
Lower transit on 05.01.2012  ephem: 08:48  skyfield: 08:49
Lower transit on 22.03.2012  ephem: --:--  skyfield: 13:52
Lower transit on 20.05.2012  ephem: 23:56  skyfield: 23:57
Lower transit on 09.07.2012  ephem: 16:49  skyfield: 16:04
Lower transit on 30.07.2012  ephem: 09:32  skyfield: 09:31
Lower transit on 08.04.2013  ephem: 22:52  skyfield: 22:00
Lower transit on 19.09.2013  ephem: 11:50  skyfield: 11:49
Lower transit on 13.04.2014  ephem: 10:34  skyfield: --:--
Lower transit on 13.08.2014  ephem: 14:33  skyfield: 14:16
Lower transit on 12.01.2015  ephem: 17:26  skyfield: 17:27
Lower transit on 27.05.2015  ephem: 07:13  skyfield: 07:14
Lower transit on 09.06.2015  ephem: 18:09  skyfield: 18:10
Lower transit on 24.06.2015  ephem: 05:51  skyfield: --:--
Lower transit on 07.07.2015  ephem: 16:59  skyfield: 16:54
Lower transit on 16.11.2015  ephem: 03:22  skyfield: 03:21
Lower transit on 19.04.2016  ephem: 10:07  skyfield: --:--
Lower transit on 28.05.2016  ephem: 17:31  skyfield: 17:32
Lower transit on 05.09.2016  ephem: 02:39  skyfield: 02:40
Lower transit on 27.03.2017  ephem: --:--  skyfield: 12:31
Lower transit on 09.04.2017  ephem: 10:49  skyfield: --:--
Lower transit on 24.09.2017  ephem: 02:53  skyfield: 02:52
Lower transit on 22.01.2018  ephem: 03:55  skyfield: 03:27
Lower transit on 08.02.2018  ephem: 18:51  skyfield: 18:50
Lower transit on 17.03.2018  ephem: --:--  skyfield: 16:04
Lower transit on 07.06.2018  ephem: 18:57  skyfield: 18:08
Lower transit on 28.08.2018  ephem: 13:34  skyfield: 13:06
Lower transit on 12.01.2019  ephem: 04:35  skyfield: 04:16
Lower transit on 03.04.2019  ephem: 23:12  skyfield: 23:03

I'll now double-check that Moon Upper Transit times agree. Test_01.zip Test_02.zip

Okay ... now I'm puzzled - because even the Upper Transit method is imperfect. In Test_03 I calculate the Moon Upper Transit for every day from 2010 to 2019 inclusive and report inconsistencies between Ephem and Skyfield. Differences of 1 minute I can ignore because these might be borderline cases where the differences are minute. First the results: EDIT: the text should read "Upper transit" below

Lower transit on 23.09.2010  ephem: --:--  skyfield: 13:07
Lower transit on 02.07.2012  ephem: 23:16  skyfield: 23:15
Lower transit on 23.07.2012  ephem: 15:31  skyfield: 15:30
Lower transit on 02.09.2012  ephem: 00:58  skyfield: 00:57
Lower transit on 09.10.2012  ephem: 06:41  skyfield: 06:40
Lower transit on 22.10.2012  ephem: 18:29  skyfield: 18:28
Lower transit on 14.11.2012  ephem: 12:20  skyfield: 12:19
Lower transit on 19.09.2013  ephem: --:--  skyfield: 19:18
Lower transit on 29.01.2014  ephem: 10:43  skyfield: 10:44
Lower transit on 30.04.2014  ephem: 13:02  skyfield: 13:03
Lower transit on 06.11.2014  ephem: 23:49  skyfield: 23:50
Lower transit on 30.08.2017  ephem: 18:56  skyfield: 18:55
Lower transit on 11.12.2017  ephem: 06:49  skyfield: 06:48
Lower transit on 27.12.2017  ephem: 19:09  skyfield: 19:08
Lower transit on 03.02.2018  ephem: 02:31  skyfield: 02:30
Lower transit on 28.06.2019  ephem: 08:11  skyfield: 08:12
Lower transit on 11.09.2019  ephem: 22:28  skyfield: 22:29
Lower transit on 14.09.2019  ephem: --:--  skyfield: 18:20

Take 14.09.2019. Ephem says there's no Upper Transit on that day. Skyfield gives me 18:20 for the Upper Transit time. I tried with different "hint" values when scanning the whole day - that's why "hint" is a parameter. But look at the Skyfield data (below) for the Moon GHA (4th. column) on Saturday (14.09.2019). It starts with 1°09.7 at 0 hours and ends with 336°28.7 at 23 hours. Sunday (15.09.2019) begins with 351°03.6 at 0 hours. SO IT NEVER CROSSES 0° GHA, i.e. Ephem is correct. Skyfield returns 18:20 which is IMPOSSIBLE. Unless I have a bug in Test_03 (which I doubt), this looks likes an issue in Skyfield. I much appreciate you looking into this - THANKS ** 2

SunMoon20190914 Test_03.zip

I have simplified this test case. I ignore Lower Transit for the moment as I am unsure that lowertransit() works correctly. However I think Upper Transit should definitely work. Here I include a test program that compares Skyfield with Skyfield (not Skyfield with Pyephem). I use Skyfield 1.11. In particular it compares:

my uppertransit() method that is used by the find_discrete() method and then rounded to the minute…
with my new find_transit() method that finds the Moon Upper/Lower transit, accurate to the minute, without using find_discrete().

My find_transit() method (described below) produces good and reliable results with data from Skyfield when compared to Pyephem, USNO Data Services and Nautical Almanac data in general. I can now use this as a reference to test against. The incorrect values [from find_discrete()] are not slightly incorrect – they are way off! I have tried various “Search Increment Hint” values [in uppertransit()] and the results are in the attached log files.

This test highlights the discrepancies that I have escalated as an issue. (I will revisit the Lower Transit issue once Upper Transit has been fixed.)

I test all dates in 2019 & 2020 and the Upper Transit discrepancies (using ‘search increment hint’ = 0.01) are listed below. I cannot figure out what is causing these incorrect values.

OLD values use find_discrete(); NEW values use find_transit()
Upper transit on 12.01.2019  OLD: 04:16  NEW: 16:56
   GHA at 16:56 :    359.95432452382823 or -0.04567547617176615
   GHA at 16:57 :    0.19717810522401652
   diff: 0.15150262905225037
Upper transit on 08.02.2019  OLD: 10:34  NEW: 14:54
   GHA at 14:54 :    359.89908665042543 or -0.10091334957456866
   GHA at 14:55 :    0.14198888797329015
   diff: 0.041075538398721495
Upper transit on 01.05.2019  OLD: 06:24  NEW: 09:30
   GHA at 09:29 :    359.81902175692414 or -0.18097824307585597
   GHA at 09:30 :    0.06186491115949183
   diff: -0.11911333191636414
Upper transit on 14.09.2019  OLD: 18:20  NEW: --:--
Upper transit on 12.10.2019  OLD: 00:32  NEW: 23:17
   GHA at 23:16 :    359.78761291850816 or -0.2123870814918405
   GHA at 23:17 :    0.030634201158434293
   diff: -0.1817528803334062
Upper transit on 08.11.2019  OLD: 07:27  NEW: 21:14
   GHA at 21:14 :    359.94291866501754 or -0.057081334982456156
   GHA at 21:15 :    0.1859617542588171
   diff: 0.12888041927636096
Upper transit on 05.12.2019  OLD: 15:14  NEW: 19:10
   GHA at 19:10 :    359.98950135527787 or -0.010498644722133577
   GHA at 19:11 :    0.23257095876829434
   diff: 0.22207231404616076
Upper transit on 29.01.2020  OLD: 07:21  NEW: 15:42
   GHA at 15:41 :    359.80718579418226 or -0.19281420581774
   GHA at 15:42 :    0.05029631850505423
   diff: -0.14251788731268578
Upper transit on 20.04.2020  OLD: 02:34  NEW: 10:18
   GHA at 10:18 :    359.9530184004713 or -0.046981599528692186
   GHA at 10:19 :    0.19616896611178214
   diff: 0.14918736658308995
Upper transit on 11.07.2020  OLD: 00:33  NEW: 04:50
   GHA at 04:49 :    359.844095822866 or -0.15590417713400484
   GHA at 04:50 :    0.08714896602292077
   diff: -0.06875521111108407
Upper transit on 30.09.2020  OLD: 22:27  NEW: 23:21
   GHA at 23:21 :    359.99336763987577 or -0.006632360124228853
   GHA at 23:22 :    0.23641453534512805
   diff: 0.2297821752208992
Upper transit on 28.10.2020  OLD: 04:18  NEW: 22:01
   GHA at 22:01 :    359.9326147295892 or -0.06738527041079578
   GHA at 22:02 :    0.17580003213058748
   diff: 0.1084147617197917
Upper transit on 24.11.2020  OLD: 10:32  NEW: 20:00
   GHA at 20:00 :    359.9307668849618 or -0.06923311503817331
   GHA at 20:01 :    0.173919450166784
   diff: 0.1046863351286107
Upper transit on 21.12.2020  OLD: 18:02  NEW: 17:57
   GHA at 17:56 :    359.79979432794164 or -0.2002056720583596
   GHA at 17:57 :    0.0427791725149973
   diff: -0.1574264995433623

Description of my new find_transit() method:

The find_transit() method detects a transit event (Upper or Lower) in Skyfield without using the find_discrete() method. The result is a time (hh:mm) accurate to the minute (as printed in Nautical Almanacs).

Rounding the time to the nearest minute implies that a transit event (for a particular day) needs to be detected between 23:59:30 the previous day till 23:59:30 the current day. This is because events within 30 seconds before midnight round up to 00:00 on the next day.

The algorithm for detecting an Upper transit is explained as follows. A list of 25 Moon GHA values (in degrees) is calculated for every hour of the day from 0 to 24 inclusive. However the first and last values (Start of Day and End of Day) are slightly different: the first is the GHA at 23:59:30 (i.e. 30 seconds earlier) the previous day; and the last is at 23:59:30 (i.e. 30 seconds earlier) the current day. The 25 GHA values are scanned and when the GHA value drops – the transit event is between the highest and lowest value… so the hour in which the transit occurs is known. (Note: Moon GHA values from 0 to 23 are calculated anyway in a Nautical Almanac.) If the GHA value never drops – there is no transit event on this day.

Now the minutes of the hour detected need to be processed. The GHA is calculated for the minute. Again, when the GHA value drops – the transit event is between the highest and lowest value… so the minute in which the transit occurs is known. And if the lower GHA value (after the event) is less than “360 minus the higher GHA value” (before the event) – it probably happened after the 30 second mark, so the minute can be rounded up. If the values are nearly equal, it’s best to actually calculate the GHA at the 30 second mark to be sure of rounding correctly.

Of course special attention must be paid to the first and last minute of the day. The last minute is easy: events between 23:59:00 and 23:59:30 are never round up. For events between 23:59:30 the previous day and 00:01 the current day – the GHA at 00:00:30 must be calculated to determine if the transit time has to be rounded up to 00:01 (otherwise 00:00 is correct).

This algorithm only needs minor modification to detect a Lower Transit. First store the colongitude of the GHA in list of 25 values (e.g. 270° returns 90°; and 90° returns 270°). Within the find_transit() method, the colongitude of every GHA calculated must be returned. That’s all – and convert the GHA results back taking the colongitude again (if required as output values). In effect we are now detecting when the Moon GHA crosses the 180° mark, but it’s easier to fake it by detecting the crossing from 360° to 0° within the algorithm.

The transit time accuracy is “to the minute”, so the process is fairly fast – but a further optimization is possible: instead of calculating every minute of the hour in question, an educated guess can be made where to start searching from, based purely on the GHA value “at the top of the hour” before the event. (The search only goes forward.) The constants chosen refine the search to beginning only a minute or two before the transit event. The result is a significantly faster algorithm.

Sincere thanks to all who read this! Test_05.zip Test Results 2019 to 2020.zip

I am returned from a trip to a conference and have had a few minutes to look into your problem! Here's a Notebook with some graphs you'll want to examine:

https://github.com/skyfielders/astronomy-notebooks/blob/master/Skyfield-Notes/transit-boundary-problem.ipynb

You will see that your is_moon_transit_at() is not well-behaved at the moments where the gast or the right ascension wrap around from the value 23.99999 back to 0. The test t.gast > ra.hours suddenly becomes the wrong value when one of the values is up near 23.99 and the other is down near 0.

Instead of using >, the standard approach when comparing angles like this is to first compute their difference, then use a modulo operation to map the result into a range like –π,+π for radians or –12,+12 for hours. I give an example expression for doing this in Python down at the bottom.

Let me know if this helps fix your script!

Brilliant answer! I tested every date from 2010 to 2029 inclusive are no discrepancy was found. I can now state that the following method correctly finds the Upper Transit (assuming of course that my find_transit() method also works correctly, which appears to be the case):

def uppertransit():
    # Build a function of time that returns the moon upper transit time.

    def is_moon_transit_at(t):
        """The function that this returns will expect a single argument that is a 
        :class:`~skyfield.timelib.Time` and will return ``True`` if the moon is up
        or twilight has started, else ``False``."""
        t._nutation_angles = iau2000b(t.tt)
        # Return `True` if the meridian is crossed by time `t`.
        position = earth.at(t).observe(moon)
        ra = position.apparent().radec(epoch='date')[0]
        return (t.gast-ra.hours+12) % 24 - 12 > 0

    is_moon_transit_at.rough_period = 0.01  # search increment hint
    return is_moon_transit_at

I'm baffled why and how this solution works - I would have never discovered this. Sincere thanks for your cooperation (from an Englishman in Munich).

Here is the fixed test program (in case anyone is interested): Test_05a.zip

In a couple of days I'll take a look at the problems I had with Lower Transit ... I guess it was the same problem here. Then I'll close this officially.

Sincere thanks for your cooperation

I'm happy to have been able to help!

I'm baffled why and how this solution works - I would have never discovered this.

It has to do with the dance the numbers do right at that dip in the graph where your original routine goes from True to False. You could try this: given the time on 2019-1-12 at which the notch happens in the graph, choose 2 times from earlier that day, 2 times during the notch, and 2 times afterward. For each, print out the actual values of gast; of the right ascension; and of the true/false for your inequality. That should help you get a better intuitive sense of where the numbers are going at the moment that your function gives the opposite answer to what you wanted!

Complete success! I like the simile regarding the "dance of numbers". Normally one's invited to dance. I'll call this issue the "Woodstock Effect" - get up and dance, just "do your thing" - and don't care what the number next to you is doing.

For anyone who's interested, here's the initial Moon Lower Transit test that fails with lots of discrepancies: Test_06.zip

I thought my earlier Lower Transit code may be faulty, so I was utterly surprised when it worked. Here's the corrected test program: Test_06a.zip

I used a "delta" parameter that's either 12 or -12:

def lowertransit(delta):
    # Build a function of time that returns the moon lower transit time.

    def is_moon_transit_at(t):
        t._nutation_angles = iau2000b(t.tt)
        # Return `True` if the antimeridian is crossed by time `t`.
        position = earth.at(t).observe(moon)
        ra = position.apparent().radec(epoch='date')[0]
        #return t.gast > ra.hours + delta
        return (t.gast-ra.hours+12-delta) % 24 - 12 > 0

    is_moon_transit_at.rough_period = 0.01  # search increment hint
    return is_moon_transit_at

... and "delta" is set by this logic:

    if tfr.gast < ra.hours:
        transit_low, y = almanac.find_discrete(tfr, tto, lowertransit(-12))
    else:
        transit_low, y = almanac.find_discrete(tfr, tto, lowertransit(12))

So then I thought if there's a way to do away with the parameter. Here's the next version of the test program: Test_06b.zip

The lowertransit method now looks like this:

def lowertransit():
    # Build a function of time that returns the moon lower transit time.

    def is_moon_transit_at(t):
        t._nutation_angles = iau2000b(t.tt)
        # Return `True` if the antimeridian is crossed by time `t`.
        position = earth.at(t).observe(moon)
        ra = position.apparent().radec(epoch='date')[0]
        #return t.gast > ra.hours + delta
        return (t.gast-ra.hours) % 24 - 12 > 0

    is_moon_transit_at.rough_period = 0.01  # search increment hint
    return is_moon_transit_at

... and it's called simply so:

    transit_low, y = almanac.find_discrete(tfr, tto, lowertransit())

All dates from 2010 to 2029 give correct results for the Moon Lower Transit! And it's a tad faster. My find_discrete method wins on speed but it only finds the Moon Lower Transit accurate to the minute.

So this completes the discussion and resolution of the issue I raised. Many thanks for your assistance!

I'm glad we were able to fix the edge case that was troubling your code! This will be a helpful issue to point folks at in the future who might suffer from a similar edge case.

skyfielders / python-skyfield

Lower Transit in 2019 agrees with Ephem except on 12.01.2019 and 03.04.2019 #280