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.