Interoperable BC timescale for tandem use of pyEphem / skyField

[tostt]

This is not really an issue, rather the results of my explorations with pyEphem, skyfield and Stellarium that I'd like to share. Any comments or corrections about my experiments are welcome.

I use pyEphem for its speed, skyField for its accuracy and "grammar" and Stellarium to get a more hands-on feeling.

When exploring an event in the distant past (BC), finding a common timescale to set the event on the 3 sofware package was not immediate. As Brandon Rhodes writes about libastro, the library that underlies pyEphem :"It provides no facility for consistently handling dates, but leaves it up to the user to understand the underlying models and call everything correctly."

For planetary events, I found that using the following method yield a timescale that can be used in all 3 programs for BC years.

1) Negative years

First, a reminder of how years are counted.

AD/BC   |    pyEphem    | skyField/Stellarium
--------+---------------+--------------------
N BC    |   -N  |   -N+1
... |   ... |   ...
2 BC    |   -2  |   -1
1 BC    |    -1 and 0   |    0
1 AD    |    1  |    1
... |   ... |   ...

If phEpyem reports an event with a negative year, you will have to add 1 to it to use it in the other packages and conversely.

2) Setting BC dates and times

After adjusting the year per the previous paragraph, use the following calls.

with pyEphem:

• Use ephem.Date()

with skyField:

• Set ts.julian_calendar_cutoff = GREGORIAN_START
• Use timescale.UT1()

with Stellarium:

• In location panel:
  • Set Custom Time Zone: checked
  • Time zone: UTC
• In Date and time panel:
  • Set time in UTC

3) Python example

In the following python example :

• test 1 shows that ephem.Date() is the same thing as Timescale.tt() for julian day computations
• test 2 shows that ephem.Date() is the same thing as Timescale.ut1() for planetary computations

# Timescale test on pyEphem and skyField

import ephem
date1ephem = ephem.Date('-100/01/01 13:00:00.0')

from skyfield.api import load, GREGORIAN_START
ts = load.timescale()
ts.julian_calendar_cutoff = GREGORIAN_START
date1skyfield = ts.tt(-99, 1, 1, 13, 0, 0.0)

# Test 1: Julian day
print('Julian day test (BC)')
print('--------------------')
print(f'pyEphem  julian day date 1 : {ephem.julian_date(date1ephem)}')
print(f'skyField julian day date 1 : {date1skyfield.tt}')
print()

# Test 2 : Moon position
print('Moon position test (BC)')
print('-----------------------')

## pyEphem 
moon = ephem.Moon()
moon.compute(date1ephem)
print('pyEphem  moon at date1', moon.a_ra, moon.a_dec)

## skyField
from skyfield.api import N,S,E,W, wgs84
ts = load.timescale()
ts.julian_calendar_cutoff = GREGORIAN_START # repeated for clarity
planets = load('de406.bsp')
earth, moon = planets['earth'], planets['moon']
observer = earth

def showSkyfieldPosition(typeDate, date):
    moonAstrometric = observer.at(date).observe(moon)
    moonRa, moonDec, _ = moonAstrometric.radec()
    print(f'skyField moon at date 1 {typeDate}', moonAstrometric.radec())

showSkyfieldPosition('UT1', ts.ut1(-99, 1, 1, 13, 0, 0.0))
showSkyfieldPosition('TT',  ts.tt (-99, 1, 1, 13, 0, 0.0))
showSkyfieldPosition('TAI', ts.tai(-99, 1, 1, 13, 0, 0.0))
showSkyfieldPosition('TDB', ts.tdb(-99, 1, 1, 13, 0, 0.0))

[brandon-rhodes]

This is not really an issue, rather the results of my explorations with pyEphem, skyfield and Stellarium that I'd like to share. Any comments or corrections about my experiments are welcome.

I don't see any corrections to offer, so unless there are mistakes that you've found Skyfield to be making, my guess is that we can go ahead and close this issue, and maybe other folks in the future who are curious will find the issue thanks to web searches and arrive here.

test 1 shows that ephem.Date() is the same thing as Timescale.tt() for julian day computations

This does not strike me as unique to tt — I think it's true that either a JD or a calendar date can be used as input or output for a Skyfield Time for the scales TAI, TT, and TDB. For example:

date1skyfield = ts.tai(-99, 1, 1, 13, 0, 0.0)
print(f'skyField julian day date 1 : {date1skyfield.tai}')

Doesn’t that have the same power to convert calendar dates to JD's as your example code that uses .tt instead?

[tostt]

Indeed, no correction is needed.

Regarding your final point, skyField's TT timescale appears to be the closest match to pyEphem. I wrote a more complete version of the first test to demonstrate.

print(f'pyEphem  julian day date 1     : {pyEphemDate}')
print(f'skyField julian day date 1 TT  : {date1skyfield.tt} \t {pyEphemDate-date1skyfield.tt}')
print(f'skyField julian day date 1 TDB : {date1skyfield.tdb} \t {pyEphemDate-date1skyfield.tdb}')
print(f'skyField julian day date 1 TAI : {date1skyfield.tai} \t {pyEphemDate-date1skyfield.tai}')
print(f'skyField julian day date 1 UT1 : {date1skyfield.ut1} \t {pyEphemDate-date1skyfield.ut1}')

Running this snippet, I get the following values and sorted time differences, in days. It appears that the skyField TT timescale is the closest match to pyEphem : the difference is of the order of 10^-10 day.

pyEphem  julian day date 1     : 1684899.0416666665
skyField julian day date 1 TT  : 1684899.0416666667      -2.3283064365386963e-10
skyField julian day date 1 TDB : 1684899.0416666768      -1.0244548320770264e-08
skyField julian day date 1 TAI : 1684899.0412941666       0.00037249992601573467
skyField julian day date 1 UT1 : 1684898.90803716         0.13362950645387173

[brandon-rhodes]

Running this snippet, I get the following values and sorted time differences, in days. It appears that the skyField TT timescale is the closest match to pyEphem : the difference is of the order of 10^-10 day.

That might be because you created your date1skyfield using ts.tt(). What if you use ts.tai() instead?

date1skyfield = ts.tai(-99, 1, 1, 13, 0, 0.0)

Then wouldn’t it be the .tai JD that's a close match?

[tostt]

You are absolutely right. When doing what you suggest, the match is excellent.

date1skyfield = ts.tai(-99, 1, 1, 13, 0, 0.0)
print(f'skyField julian day date 1 : {date1skyfield.tai}')

gives skyField julian day date 1 : 1684899.0416666667 which is extremely close to the target: 1684899.0416666665

And the same goes with all timescales. Thank you for taking the time to help me.

[brandon-rhodes]

And the same goes with all timescales. Thank you for taking the time to help me.

I'm glad you were able to work out the date conversions you needed! I hope your project goes well.