PyEphem: dates of solstices and equinoxes and their validity on longer time scales

Question

My question is what the time span is over which PyEphem provides accurate results for the dates of the solstices and equinoxes and the solar geometry.

So far, I found the limits B.C. 9998-03-20 to A.D. 9999-12-31 in this very informative post on GitHub https://github.com/brandon-rhodes/pyephem/issues/61 and an overall indication that results become erratic when moving beyond +/- 20,000 years on either site of present.

I'd like to have this clarified for I am trying to obtain the position of the sun over longer periods -typically going back to 10,000 years BP- in order to compute the incoming solar radiation from the altitude and azimuth of the sun at a given location. PyEphem seems to offer a good alternative to the functions provided by for example Berger, 1978 (J. Atmosph. Sc., 35: 2362-2367). To me, an essential advantage of PyEphem over these algorithms would be that it also keeps track of time whereas the earth orbit is usually fixed at one particular moment in the aforementioned algorithms (e.g., vernal equinox at March 21).

Typically, the variations in the earth orbit with the algorithms of Berger and others are evaluated on the scale of the last glacial (up to 126,000 years BP). While evaluating PyEphem over that range, I came across some strange behaviour for the dates of the solstices when the dates lie well before the present day:

import ephem
date= ephem.date((-59000,1,1))

orbitPoints= ['vernal_equinox_start','summer_solstice',\
  'autumnal_equinox','winter_solstice','vernal_equinox_end']

dates= {}
dates['vernal_equinox_start']= ephem.next_vernal_equinox(date)
dates['summer_solstice']= ephem.next_summer_solstice(dates['vernal_equinox_start'])
dates['autumnal_equinox']= ephem.next_autumnal_equinox(dates['vernal_equinox_start'])
dates['winter_solstice']= ephem.next_winter_solstice(dates['vernal_equinox_start'])
dates['vernal_equinox_end']= ephem.next_vernal_equinox(dates['winter_solstice'])

for orbitPoint in orbitPoints:
  date= dates[orbitPoint]
  distance= body_distance(sun, date)
  hlon, hlat= body_hpos(sun, date)
  print '%-20s %30s %8.4f %12s %12s' %\
    (orbitPoint, date, distance, hlon, hlat)

print 'days between between equinoxes: %.1f, year length (vernal equinox): %.1f' %\
  (dates['autumnal_equinox']-dates['vernal_equinox_start'],\
   dates['vernal_equinox_end']-dates['vernal_equinox_start'])

Gives:

vernal_equinox_start          -59001/12/20 01:33:21   1.6255  180:01:53.6    0:00:01.7
summer_solstice                -59000/2/19 04:16:18   1.3380  180:03:56.7    0:00:09.6
autumnal_equinox                -59000/7/1 09:20:58   0.3704    0:00:42.2    0:00:09.6
winter_solstice                -59000/10/9 07:35:39   1.1335  179:56:59.4   -0:00:12.4
vernal_equinox_end            -58999/10/10 07:40:47   1.1348  180:00:52.7    0:00:16.2
days between between equinoxes: 193.3, year length (vernal equinox): 659.3

I still get a reasonable (?) value for the year length when I set the date to date= ephem.date((-25000,1,1)).

If PyEphem indeed gives accurate results over my period of interest (10,000 years BP) it would be serve my purpose. However, I'd like to have this clarified and am open for suggestions to extend this range, even if it were just for validation. I was looking at SkyField as an alternative but it does not appear to offer an extended range.

A clear description of the range of validity of PyEphem and any suggestions would be greatly appreciated.

Answer 1

The results you are getting are reasonable, though I cannot speak to their specific accuracy. It is only briefly that the Earth’s rate of rotation is similar to that of the present day. Just a few hundred years ago days were shorter, and in the future they will be ever longer (as the error accumulates we have leap seconds), and the accumulated difference each and every year of history makes our current Gregorian calendar a worse and worse approximation to the seasons each thousand years you jump backwards. See:

https://eclipse.gsfc.nasa.gov/SEcat5/deltat.html

I have not sat down and calculated whether the accumulated difference really accounts for the difference you see above, or whether another effect is at work as well, but I would certainly not be surprised to find our current Gregorian calendar of days exactly the length of ours, when projected back to 59,000 BC, lands somewhere else in the seasons entirely than where it does today.