skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.4k stars 211 forks source link

almanac.find_discrete ERROR: end of civil twilight has no beginning #661

Closed aendie closed 2 years ago

aendie commented 2 years ago

This is a very specific error... 99.9% of the times (at least!) it give me correct values. It appears that this problem has existed for a longer time, but I only hit this odd-ball case just recently. I've created a MWE to reproduce it. I tested this now with Skyfield 1.39 and 1.40.

In my MWE, only the last 40 lines are relevant to the test case, i.e. the daylength() function and below. The rise_set() function is a generic function to return the hours and minutes of the event(s) and has been working fine for months. A few "print" statements show that on 11th July 2021 at 62°N latitude, civil twilight ends at 23:53 and there was no civil twilight start. The correct result should be that there is no Civil Twilight, period. Here's the console output of TEST.py using Python 3.9.4:

image

I noticed this in a Nautical Almanac page I created (the error is highlighted):

image

This error came up while comparing the results that Ephem gave me (Ephem is correct). Thanks in anticipation of you looking into this!! TEST.zip

brandon-rhodes commented 2 years ago

Could you double-check the result against HORIZONS, so we have a value for comparison? You haven't specified how you know the number to be in error. Skyfield’s search routine here is correctly measuring that, at least for Skyfield, the Sun did drop below -6° altitude between 23:53 and 23:54, as we can confirm by adding these lines to your test script:

locn_at = (earth + locn).at
for t in ts.utc(2021, 7, 11, 23, 53), ts.utc(2021, 7, 11, 23, 54):
    alt, az, distance = locn_at(t).observe(sun).apparent().altaz()
    print(alt)

Which prints:

-05deg 59' 56.5"
-06deg 00' 17.6"

Which would be the conclusion of civil twilight, I think, if I understand the context?

aendie commented 2 years ago

Mea Culpa! ... and apologies for wasting your time.

One hour after posting this I realized (in my half-sleep) that I was wrong and Skyfield is correct. I had misread my table: this is the beginning of Civil Twilight, not the end of Civil Twilight. This the marginal condition when Civil Twilight begins for the first time after a pause. Adding the next day shows that on the following day Civil Twilight ends at 00:18 and begins a little earlier at 23:38.

image

And neither is Ephem incorrect: Ephem probably waits another day before beginning Civil Twilight - slightly different data is to be expected from Ephem.

What misled me (apart from some tired gray matter in my head), was Jorrit Visser's website. He also uses Skyfield but posted no beginning or end of Civil Twilight on 11th July 2021 at 62°N latitude. Not only that, but another data item (Sunrise?) was way off. His (very useful and mostly excellent) website is currently down so I can't show you right now. He leaves no contact information, but, like me, he's not infallible.

brandon-rhodes commented 2 years ago

I'm glad that your program turns out to be printing a correct table, at least for the altitudes as returned by Skyfield. I'll go ahead and close the issue, though feel free to follow up with additional comments if you learn more about the reasons behind the diverging data on Visser's site.

aendie commented 2 years ago

Thanks for your message. I am still waiting for the celnav.nl web site to come up again: it closed about a week ago. Maybe it's undergoing maintenance.

But, joking apart, I have another good post to make shortly under Ephem. This time I am triple checking my facts before I make any bold claims. Not having a single Nautical Almanac at home, I find Ephem very useful to compare against. And Ephem is way better than I ever imagined ... except for some unknown quirks that made me open my eyes wide. I'll post once I have a MWE of sorts to demonstrate the discrepancies without my heavy code in between. (I still can't believe no one ever spotted these.)