Fix twilights API to return correct results

[andris-zalitis]

Here is a code that I used to print astronomical dawn and dusk times for each day in 2017 for a location in Paris:

    func calculateTwilights() {
        let paris = GeographicCoordinates(CLLocation(latitude: 48.866667, longitude: 2.350000)) 

        let gregorianCalendar = Calendar(identifier: .gregorian)
        var date = JulianDay(year: 2017, month: 1, day: 1).date

        for _ in 1...365 {
            let earth = Earth(julianDay: JulianDay(date))

            let twilights = earth.twilights(forSunAltitude: TwilightSunAltitude.astronomical.rawValue, coordinates: paris)

            if let rise = twilights.rise {
                print("dawn: \(rise.date)")
            } else {
                print("dawn: ---")
            }
            if let set = twilights.set {
                print("dusk: \(set.date)")
            } else {
                print("dusk: ---")
            }

            date = gregorianCalendar.date(byAdding: .day, value: 1, to: date)!
        }
    }

Location coordinates (latitude: 48.866667, longitude: 2.350000) correspond to the following in the degrees: 48°52'00.0"N 2°21'00.0"E (that can be checked here)

This is the setup that I use for checking the data at http://aa.usno.navy.mil/:

Here is a peek into the results that I got from SwiftAA:

dawn: 2017-01-01 05:54:31 +0000
dusk: 2017-01-01 18:07:15 +0000
dawn: 2017-01-02 05:54:31 +0000
dusk: 2017-01-02 18:08:10 +0000
dawn: 2017-01-03 05:54:28 +0000
dusk: 2017-01-03 18:09:07 +0000
dawn: 2017-01-04 05:54:23 +0000
dusk: 2017-01-04 18:10:05 +0000
dawn: 2017-01-05 05:54:16 +0000
dusk: 2017-01-05 18:11:05 +0000
dawn: 2017-01-06 05:54:06 +0000
dusk: 2017-01-06 18:12:06 +0000
dawn: 2017-01-07 05:53:54 +0000
dusk: 2017-01-07 18:13:08 +0000
dawn: 2017-01-08 05:53:39 +0000
dusk: 2017-01-08 18:14:12 +0000
dawn: 2017-01-09 05:53:21 +0000
dusk: 2017-01-09 18:15:17 +0000
dawn: 2017-01-10 05:53:02 +0000
....

Here is what USNO gives: http://aa.usno.navy.mil/cgi-bin/aa_rstablew.pl?ID=AA&year=2017&task=4&place=&lon_sign=1&lon_deg=2&lon_min=21&lat_sign=1&lat_deg=48&lat_min=52&tz=0&tz_sign=-1

So it looks like we're a little bit off.

[onekiloparsec]

Thanks for reporting this! Indeed, we are a bit off. It will require some investigations.

[andris-zalitis]

Here I found a similar algorithm to the one used in twilights but there it is used with iterations for greater precision:

The orginal method only gives an approximate value of the Sun's rise/set times. The error rarely exceeds one or two minutes, but at high latitudes, when the Midnight Sun soon will start or just has ended, the errors may be much larger. If you want higher accuracy, you must then use the iteration feature. This feature is new as of version 0.7. Here is what I have tried to accomplish with this.

a) Compute sunrise or sunset as always, with one exception: to convert LHA from degrees to hours, divide by 15.04107 instead of 15.0 (this accounts for the difference between the solar day and the sidereal day.

b) Re-do the computation but compute the Sun's RA and Decl, and also GMST0, for the moment of sunrise or sunset last computed.

c) Iterate b) until the computed sunrise or sunset no longer changes significantly. Usually 2 iterations are enough, in rare cases 3 or 4 iterations may be needed.

[onekiloparsec]

Andris, thanks for checking this. If look closely, you'll see that the iterations hint is coming from the exact same source as mine (an algorithm from Paul Schlyter, see here.

I managed to make the iterations work, and fixed a by-one-day shift in the original algorithm. This is what I get now:

dawn: 2017-01-01 05:48:55 +0000
dusk: 2017-01-01 18:00:25 +0000
dawn: 2017-01-02 05:48:03 +0000
dusk: 2017-01-02 18:00:22 +0000
dawn: 2017-01-03 05:47:09 +0000
dusk: 2017-01-03 18:00:21 +0000
dawn: 2017-01-04 05:46:14 +0000
dusk: 2017-01-04 18:00:23 +0000
dawn: 2017-01-05 05:45:17 +0000
dusk: 2017-01-05 18:00:27 +0000
dawn: 2017-01-06 05:44:18 +0000
dusk: 2017-01-06 18:00:34 +0000
dawn: 2017-01-07 05:43:18 +0000
dusk: 2017-01-07 18:00:43 +0000
dawn: 2017-01-08 05:42:17 +0000
dusk: 2017-01-08 18:00:54 +0000
dawn: 2017-01-09 05:41:14 +0000
dusk: 2017-01-09 18:01:08 +0000
dawn: 2017-01-10 05:40:10 +0000
dusk: 2017-01-10 18:01:24 +0000

It looks like we are a lot less off compared to USNO, right? :-) In fact, accuracy is easily below the second. I need to cleanup a bit the code before it will appear here on GitHub, but I'm on my way.

[onekiloparsec]

I spoke a little too soon, business as usual... January first fits very well. However, along the year, it is not so fit. And my algorithm switch rise and set around April. More work to be done, but the code is pushed. You can try yourself.

[andris-zalitis]

Impressive work! I will try the new version. Interesting that it fits so well for January 1st, but then changes more rapidly for dusk and less rapidly for dawn than USNO values. Thank you for the link to Paul Schlyter's algorithm!

[andris-zalitis]

OK, I did spend some time looking into it.

I was able to come up with an updated version of Schlyter's algorithm implementation (see this branch in case you're curious). It did improve accuracy of the iteration implementation and also of the original non iteration algorithm. However, it was still not accurate enough (delta of 2-5 minutes).

I started looking into other variants of implementing this calculation because I suspected that Paul Schlyter's algorithm is not very accurate itself. Here's his quote:

These algortihms were developed by myself back in 1979, based on a preprint from T. van Flandern's and K. Pulkkinen's paper "Low precision formulae for planetary positions", published in the Astrophysical Journal Supplement Series, 1980. It's basically a simplification of these algorithms, while keeping a reasonable accuracy.

I even found this document which seems to include a detailed instruction of how to calculate sunrise and sunset times based on The Astronomical Algorithms.

But then... and here's the kicker - I found out that AA+ library already has this implemented! Moreover SwiftAA already has it wrapped! It's RiseTransitSetTimes and it works like magic 😃! The only thing I had to do was to add a possibility of providing custom altitude when doing the calculation (since we need -18° for Sun's astronomical dawn/dusk).

I also added a temporary test that just prints out the times. Here is the PR: #52. The results match USNO data precisely (too bad they round up to a minute), with the only exception being two dates when the period of no dawn/dusk starts.

[onekiloparsec]

Thanks Andris for that work! I knew there was something I overlooked. But I am glad it finally ended up with a fewer-lines-of-codes solution. I'll merge your PR, and refactor a bit the things to make the correct (AA+) twilights the default ones (that is, removing the Schylter code).

[onekiloparsec]

I just pushed some refactored code. It doesn't change much, but it is a lot better: we always get the transit time, not only rise and set. And I've added some better unit tests for above/below the artic/antarctic circles.