how accurate are PyEphem rise and set calculations over several hundred years?

Question

Calculating daily hours of sunlight (sunset-sunrise) and finding the maximum for each year (generally on the solstice but not always), an interesting pattern emerges. About 5 seconds of sunlight are lost each century.

Is this an error factor within PyEphem? Is this accurate and PyEphem taking into account variations in Earth's orbit? some other reason?

import pandas as pd
import ephem

sun = ephem.Sun()
raleigh = ephem.Observer()
raleigh.lon, raleigh.lat = "-78.6382", '35.7796'
raleigh.horizon = '-0:34'  # USNO standard atmospheric diffraction
raleigh.pressure = 0       # atmospheric refraction parameters

def riseset(date, f):
    # compute passed function (sunrise or sunset)
    raleigh.date = date
    sun.compute(raleigh)
    sr = ephem.localtime(f(sun))
    return sr

def createdataframe(start, end):
    # create a dataframe index by daily dates, add columns for the
    # sunrise, sunset, and their delta
    df = pd.DataFrame(index=pd.date_range(start=start, end=end,  freq='D'))
    df['date'] = df.index.map(lambda d: d.strftime("%b %d"))
    df['sunrise'] = df.index.map(lambda d: riseset(d, raleigh.next_rising))
    df['sunset'] = df.index.map(lambda d: riseset(d, raleigh.next_setting))
    df['daylightdelta'] = df['sunset'] - df['sunrise']
    return df

def outputmax(df, year):
    i = df['daylightdelta'].idxmax()  # index of the day where the sun is visible above the horizon for the most time
    return "solstice: %s longest day sunrise: %s sunset: %s daylight: %s" % (
        ephem.localtime(ephem.next_solstice(str(year))).strftime("%Y %b %d %X"),
        df.loc[i]['sunrise'].strftime("%b %d %X"),
        df.loc[i]['sunset'].strftime("%T"),
        df.loc[i]['daylightdelta'])

if __name__ == "__main__":
    for year in range(1900,2201):
        # looping through 1900-2200, find the date with the most hours of sunlight
        start = '%d-01-01 04:00:00' % year # compensating for UTC which can throw off pandas columnar math
        end = '%d-12-31 23:59:00' % year
        print outputmax(createdataframe(start, end), year)

Could you supply a small Python script demonstrating this result over several centuries? Otherwise other programmers reading this question will not know where to start to replicate your result. — Brandon Rhodes, Jun 22 '17 at 00:29
Of course, I should have done that from the start. Edited above. Uses PyEphem (of course) and a pandas dataframe for simplified calculation of the longest day each year. — rtphokie, Jun 23 '17 at 01:54

score 2 · Accepted Answer · answered Jul 19 '17 at 09:03

My guess is that PyEphem is showing you a real phenomenon. While I am not enough of an expert to list all of the many swirling factors that affect a number like the length of the longest day, one that stands out to me is that the tilt of the Earth’s pole varies over the ages and is currently decreasing:

https://en.wikipedia.org/wiki/Milankovitch_cycles#Axial_tilt_.28obliquity.29

Let's attempt a very rough back-of-an-envelope guess as to the size this effect might have. If over 41,000 years the tilt goes from maximum to minimum and back, then the current half-cycle from maximum tilt back to minimum must take around 20,500 years. Though of course the real adjustment is sinusoidal, with slow change near the maxima and then more rapid change in the middle between the extremes, what if it were simply linear, as a first approximation? Then the rate of change over 20,500 years = 205 centuries would be roughly:

(24.5 - 22.1) degrees / 205 centuries ≅ 0.01 degrees

So, the tilt of the axis might be expected to change by about a hundredth of a degree per year. How many seconds of daylight does the longest day change in Raleigh if you vary the axis tilt by 0.01 degrees? PyEphem doesn't let us change the axial tilt arbitrarily, so let's adjust the position of Raleigh instead. Change the bottom clause of your program to:

if __name__ == "__main__":
    year = 2000
    start = '%d-01-01 04:00:00' % year # compensating for UTC which can throw off pandas columnar math
    end = '%d-12-31 23:59:00' % year

    raleigh.lat = '35.76'
    print outputmax(createdataframe(start, end), year)

    raleigh.lat = '35.77'
    print outputmax(createdataframe(start, end), year)

    raleigh.lat = '35.78'
    print outputmax(createdataframe(start, end), year)

And as output you should get:

solstice: 2000 Jun 20 21:47:51 longest day sunrise: Dec 21 07:20:52 sunset: 17:04:27 daylight: 0 days 14:16:24.970989
solstice: 2000 Jun 20 21:47:51 longest day sunrise: Dec 21 07:20:54 sunset: 17:04:26 daylight: 0 days 14:16:28.206467
solstice: 2000 Jun 20 21:47:51 longest day sunrise: Dec 21 07:20:55 sunset: 17:04:24 daylight: 0 days 14:16:31.442902

Around three or four seconds per century, which is roughly the magnitude of the effect you are seeing. My guess, then, is that your program is uncovering the gradual decrease in the planet's axial tilt that is gradually making the seasons less extreme and the longest day and longest night less extreme.

how accurate are PyEphem rise and set calculations over several hundred years?

1 Answers1