previous_rising() skips sunrise times at the north pole

[ninneman]

Using the latest PyEphem pulled directly from Github, I'm attempting to compute sunrise times at Earth's north pole, which occur only once per year (in March). I've discovered that previous_rising() skips over the sunrise for one or two years, depending on which year I use as the starting point.

This program demonstrates the error:

#!/usr/bin/python3

import ephem
import datetime

def get_previous_rise(observer, body):
    while True:
        try:
            rise_time = observer.previous_rising(body)
            break
        except (ephem.AlwaysUpError, ephem.NeverUpError):
            observer.date = observer.previous_antitransit(body)
    return rise_time

np = ephem.Observer()
np.lon = 0 * ephem.degree
np.lat = 90 * ephem.degree
np.elevation = 0
sun = ephem.Sun()

for year in range(1900, 2100):
    start_time = ephem.Date(datetime.datetime(year, 6, 1))
    np.date = start_time
    rise_time = get_previous_rise(np, sun)
    difference = start_time.datetime() - rise_time.datetime()
    print("for {0}, sunrise was {1} previous".format(year, difference))

Which gives output like:

for 1900, sunrise was 804 days, 15:41:08.341447 previous
for 1901, sunrise was 73 days, 22:24:44.625040 previous
for 1902, sunrise was 73 days, 16:32:29.915412 previous
for 1903, sunrise was 438 days, 16:32:29.915412 previous
for 1904, sunrise was 804 days, 16:32:29.915412 previous
for 1905, sunrise was 73 days, 22:52:27.044984 previous
for 1906, sunrise was 73 days, 16:55:54.082062 previous
for 1907, sunrise was 438 days, 16:55:54.082062 previous
for 1908, sunrise was 804 days, 16:55:54.082062 previous
for 1909, sunrise was 73 days, 23:33:51.867867 previous
for 1910, sunrise was 73 days, 17:45:51.699219 previous
for 1911, sunrise was 73 days, 11:54:03.381875 previous
for 1912, sunrise was 439 days, 11:54:03.381875 previous
for 1913, sunrise was 804 days, 11:54:03.381875 previous
for 1914, sunrise was 73 days, 18:36:47.585368 previous
for 1915, sunrise was 73 days, 12:56:34.473033 previous
for 1916, sunrise was 439 days, 12:56:34.473033 previous

and so on. Note that "previous" means previous to June 1, a date I hard-coded into the program above. The sun is definitely up at that date. The values around 73 days are probably correct, since June 1 - 73 days is roughly the March equinox. The others obviously aren't right at all. Interestingly, changing the latitude to 89.9 makes this problem disappear. I'm continuing my investigation.

[ninneman]

Disabling refraction (pressure = 0) stops the skipping. That doesn't directly answer why it's happening in the first place. The investigation continues.

[brandon-rhodes]

Having thought about this issue for two weeks: I think the problem is that the rising logic expects a rising, if it occurs at all, to occur between the antitransit and the transit. But at the north pole, where the Sun's motion northward is simply a gradual uninterrupted climb into the sky, the sunrise is as likely to occur on the west side of the sky as in the east — and thus will, half the time, be missed by the algorithm?

[ninneman]

I think you're on to something. I whipped up this program to brute-force the rise times for 1900-1916, and then get the sun's azimuth at that time.

#!/usr/bin/python3

import sys
import os
import ephem
import datetime
import math

for year in range(1900, 1917):
    np.date = "{0}/4/1".format(str(year))
    old_sign = None
    for step in range(40000):
        sun = ephem.Sun(np)
        limb_alt = sun.alt + sun.radius - np.horizon
        new_sign = math.copysign(1, limb_alt)
        if (new_sign != old_sign) and old_sign != None:
            print(np.date, sun.az)
            break
        old_sign = new_sign
        np.date -= 0.5 * ephem.minute

The output is:

1900/3/18 19:49:30 295:19:45.7
1901/3/19 01:35:00 21:42:08.6
1902/3/19 07:27:30 109:49:40.6
1903/3/19 13:25:30 199:20:01.4
1904/3/18 19:11:00 285:42:16.1
1905/3/19 01:07:30 14:50:00.2
1906/3/19 07:04:00 103:57:46.6
1907/3/19 12:44:30 189:04:51.4
1908/3/18 18:36:30 277:04:59.8
1909/3/19 00:26:00 4:27:27.9
1910/3/19 06:14:00 91:27:28.9
1911/3/19 12:05:30 179:20:08.4
1912/3/18 17:42:30 263:34:42.8
1913/3/18 23:28:00 349:57:12.3
1914/3/19 05:23:00 78:42:31.0
1915/3/19 11:03:00 163:42:17.1
1916/3/18 16:57:30 252:20:02.8

For every year where previous_rising() overshot sunrise, the azimuth was greater than 180, and less than 180 for all other years.

[ninneman]

It occurred to me that if this hypothesis is correct, then changing the longitude of the observer would change the overshooting behavior. Sure enough, with lon = 180:

for 1900, sunrise was 74 days, 4:10:24.125156 after
for 1901, sunrise was 439 days, 4:10:24.125156 after
for 1902, sunrise was 804 days, 4:10:24.125156 after
for 1903, sunrise was 73 days, 10:34:16.844250 after
for 1904, sunrise was 74 days, 4:48:45.140556 after
for 1905, sunrise was 439 days, 4:48:45.140556 after
for 1906, sunrise was 804 days, 4:48:45.140556 after
for 1907, sunrise was 73 days, 11:15:04.181975 after
for 1908, sunrise was 74 days, 5:23:13.585147 after
for 1909, sunrise was 439 days, 5:23:13.585147 after
for 1910, sunrise was 804 days, 5:23:13.585147 after
for 1911, sunrise was 1169 days, 5:23:13.585147 after
for 1912, sunrise was 74 days, 6:17:29.478559 after
for 1913, sunrise was 74 days, 0:31:54.131598 after
for 1914, sunrise was 439 days, 0:31:54.131598 after
for 1915, sunrise was 804 days, 0:31:54.131598 after
for 1916, sunrise was 74 days, 7:02:17.576049 after

That's the inverse of the behavior when lon = 0.

[ghost]

Possibly related:

http://astronomy.stackexchange.com/questions/937/can-moon-set-after-being-up-when-due-north-pseudo-circumpolarity

On Tue, 9 Dec 2014, Brandon Rhodes wrote:

Date: Tue, 09 Dec 2014 20:22:13 -0800 From: Brandon Rhodes notifications@github.com Reply-To: brandon-rhodes/pyephem <reply+0004c41d51a95e76917f9f3488a365ad05ba63a805e8bfe192cf00000001109f8c7 592a169ce02fa6c79@reply.github.com> To: brandon-rhodes/pyephem pyephem@noreply.github.com Subject: Re: [pyephem] previous_rising() skips sunrise times at the north pole (#65)

Having thought about this issue for two weeks: I think the problem is that the rising logic expects a rising, if it occurs at all, to occur between the antitransit and the transit. But at the north pole, where the Sun's motion northward is simply a gradual uninterrupted climb into the sky, the sunrise is as likely to occur on the west side of the sky as in the east ? and thus will, half the time, be missed by the algorithm?

---

Reply to this email directly or view it on GitHub: https://github.com/brandon-rhodes/pyephem/issues/65#issuecomment-66402835

[brandon-rhodes]

Excellent work, @ninneman! The two tests you put together have conclusively demonstrated the failure mode. I will try to figure out how to make the search routines smarter about the situation.

[brandon-rhodes]

And, @barrycarter, that Stack Exchange post does have a nice picture of what it looks like when something sets in an unusual direction from near the Pole.

[ninneman]

I've also been trying to think up a "global" fix that can be used for any observer. I'm not coming up with one, but using a transit-antitransit window doesn't seem to be salvageable. Incidentally, I'm assuming that "transit-antitransit" is used in the search algorithm as "upper culmination-lower culmination", but as a link at the bottom of that Stack Exchange post makes clear, those aren't generally the same in a way that's even more subtle than the error discovered here. There may actually be other missed rises/sets in highly particular body-observer pairs.

[ninneman]

It occurred to me that my "incidental" comment was, in fact, a global solution. I cooked up a draft fix in the PR listed above. I tested it using my test code above, and it works correctly! Maybe what I've done here can inspire something better.

[brandon-rhodes]

Interesting! That is an interesting alternative way to code up a search, and I am reading over it. How did you choose 10 hours as the particular step size when doing this search?

[ninneman]

That 10-hour thing was supposed to be a placeholder. My reasoning was that 10 hours is less than half of 24 hours, so I'm guaranteed not to skip over a zero by walking that interval.

[ninneman]

My latest PR only fails two of the regression tests. I made a dumb mistake in the last one that gave wrong results everywhere except the north and south poles.

[brandon-rhodes]

I have just tried out your new algorithm, and I think that it might still miss certain sunrises. For example, the sunrise in 1905 if viewed from close to, but not exactly at, the pole:

import ephem

polar = ephem.Observer()
polar.lon = 0 * ephem.degree
polar.lat = 89.9 * ephem.degree
polar.elevation = 0

sun = ephem.Sun()

for date in [
        '1905/3/18',
        '1905/3/19',
        '1905/3/20',
        '1905/3/21',
        ]:
    polar.date = date
    try:
        polar.previous_rising(sun)
    except Exception as e:
        print e

The output that I get from this script, when on your branch (at least as shown in the pull request), is:

'Sun' culminates below the horizon at 1905/3/17 14:40:00
'Sun' culminates below the horizon at 1905/3/18 14:41:54
'Sun' is still above the horizon at 1905/3/19 21:31:46
'Sun' is still above the horizon at 1905/3/20 21:30:00

Do you see the same thing, or do you get a different result?

[ghost]

Thought: the sun's changing declination (which causes all this odd behavior) can never increase by more than .12 degrees per day.

Perhaps change the algorithm to assume an up-to-.01 degrees/hour (ie, twice the max value) change in declination to find the check interval, instead of a fixed 6 hour period?

Sort of like a Newton's method.

On Tue, 13 Jan 2015, Brandon Rhodes wrote:

Date: Tue, 13 Jan 2015 21:07:00 -0800 From: Brandon Rhodes notifications@github.com Reply-To: brandon-rhodes/pyephem <reply+0004c41dced237ef7332a3a4ba1f9560cae0982070ce6cac92cf0000000110cdbb7 492a169ce02fa6c79@reply.github.com> To: brandon-rhodes/pyephem pyephem@noreply.github.com Cc: barrycarter github@barrycarter.info Subject: Re: [pyephem] previous_rising() skips sunrise times at the north pole (#65)

I have just tried out your new algorithm, and I think that it might still miss certain sunrises. For example, the sunrise in 1905 if viewed from close to, but not exactly at, the pole:

import ephem

polar = ephem.Observer()
polar.lon = 0 * ephem.degree
polar.lat = 89.9 * ephem.degree
polar.elevation = 0

sun = ephem.Sun()

for date in [
       '1905/3/18',
       '1905/3/19',
       '1905/3/20',
       '1905/3/21',
       ]:
   polar.date = date
   try:
       polar.previous_rising(sun)
   except Exception as e:
       print e

The output that I get from this script, when on your branch (at least as shown in the pull request), is:

'Sun' culminates below the horizon at 1905/3/17 14:40:00
'Sun' culminates below the horizon at 1905/3/18 14:41:54
'Sun' is still above the horizon at 1905/3/19 21:31:46
'Sun' is still above the horizon at 1905/3/20 21:30:00

Do you see the same thing, or do you get a different result?

---

Reply to this email directly or view it on GitHub: https://github.com/brandon-rhodes/pyephem/issues/65#issuecomment-69869862

[brandon-rhodes]

I suspect that for high latitudes we are supposed to adjust the algorithm to find the minimum of the second derivative each day, to find the local maximum that the sun creates as it dips and then starts upward again.

I wonder if it is time for me to invest in the book behind the Astronomical Almanac, in the hopes that they already have a perfect algorithm for doing sunrises from even high latitudes?

http://www.amazon.com/Explanatory-Supplement-Astronomical-Almanac-Urban/dp/1891389858/

[ninneman]

Do you see the same thing, or do you get a different result?

Yep, I see the same behavior. However, try this:

import ephem

polar = ephem.Observer()
polar.lon = 0 * ephem.degree
polar.lat = 89.9 * ephem.degree
polar.elevation = 0

sun = ephem.Sun()

for date in [
        '1905/3/18 06:00',
        '1905/3/19 06:00',
        '1905/3/20 06:00',
        '1905/3/21 06:00',
        ]:
    polar.date = date
    try:
        risetime = polar.previous_rising(sun)
        print("rise detected at {0}".format(risetime))
    except Exception as e:
        print(e)

Which for me produces:

'Sun' culminates below the horizon at 1905/3/17 14:40:00
rise detected at 1905/3/19 04:08:50
'Sun' is still above the horizon at 1905/3/19 21:31:46
'Sun' is still above the horizon at 1905/3/20 21:30:00

What this tells me is that the root-finding is failing or succeeding depending on the initial guesses.

[ninneman]

Nope, it wasn't the root-finding algorithm. It's much simpler. Since sunrise is at 1905/3/19 at 04:08, it would only be caught by searching backward from 1905/3/20 in your test code. But since the sun upper-culminates above the horizon at 1905/3/19 21:31:46 and errors out accordingly, the search never sees 04:08.

So the problem--if there is one--is in the CircumpolarError decision code. However, there is similar behavior in the current transit-based rise and set calculator. The canonical way to work around it is to re-search starting from the datetime where the error left off. That could work here too.

[brandon-rhodes]

I am going to close this issue for now. While I do plan to keep PyEphem updated so that it can compile under new versions of Python as they come out, I am not planning to invest my time in its old solvers, particularly for edge cases like the Poles; my attention instead needs to be focused on the more robust routines in Skyfield.

If anyone reading this later does develop a solution along with tests, though, I'm not averse to having them comment here and submitting the code for consideration!