%pylab inline

import ephem

start1 = ephem.date('2013/11/26 18:37')
start2 = ephem.date('2013/11/26 18:38')

arctic = ephem.Observer()
arctic.pressure = 0 
arctic.horizon = '-0:34' 
arctic.lat, arctic.lon = '89:30', '0' 

arctic.date = start1
try:
    arctic.next_setting(ephem.Moon())
except Exception as e:
    print e

m = ephem.Moon()
arctic.date = start2
setting = arctic.next_setting(m)
print setting, 'gives altitude', m.alt

m = ephem.Moon()
m.compute(arctic)
m.alt

arctic.date = '2013/11/26 00:00'
dates = []
altitudes = []
azimuths = []
for n in range(24 * 60 * 3):
    arctic.date += ephem.minute
    m.compute(arctic)
    dates.append(arctic.date)
    altitudes.append(m.alt)
    azimuths.append(m.az)

plot(dates, altitudes)
vlines([start1, setting], [-0.199], [0.099])
hlines([arctic.horizon - m.radius], [4.1602e4 + 0.5], [4.1602e4 + 3.99])

class MyMoon(ephem.Moon):
    def compute(self, observer):
        ephem.Moon.compute(self, observer)
        print '*', observer.date, self.alt

# From start1, the algorithm barely explores at all

mm = MyMoon()
arctic.date = start1
try:
    arctic.next_setting(mm)
except Exception as e:
    print e

# From start2, it moves boldly forward to the setting time

arctic.date = start2
setting = arctic.next_setting(mm)

# But why is the first search ending so soon?
# Let's watch.

import os, trace
t = trace.Trace(count=0, trace=1, ignoredirs=[os.path.dirname(trace.__file__)])
arctic.date = start1
try:
    t.runfunc(arctic.next_setting, mm)
except Exception as e:
    print e

# Wait! Sunrise and sunset are searched for in the current
# day in LOCAL time, not Universal Time, where the day
# changes when the sun passes zero azimuth. Is that the
# crucial difference between these two dates?

plt.subplot(2, 1, 1)
plot(dates, altitudes)
hlines([arctic.horizon - m.radius], [4.1602e4 + 0.5], [4.1602e4 + 3.99])
vlines([start1, setting], [-0.199], [0.099])
plt.subplot(2, 1, 2)
plot(dates, azimuths)

# Yes! That is the solution: start1 and start2
# are on different days, if we count days in local
# solar time.

arctic.date = start1
m.compute(arctic)
print m.az

arctic.date = start2
m.compute(arctic)
print m.az