0

I'm refactoring some old code that used PyEphem to use Skyfield, and I'm getting a slight difference in the results of the GHA/dec of a body.

def sf(year):
    from skyfield.api import N, W, wgs84
    from skyfield.api import load
    from skyfield.units import Angle

    ts = load.timescale()
    eph = load('de421.bsp')
    greenwich = wgs84.latlon(58.47722 * N, 0.0 * W)
    t = ts.utc(year)
    e = eph['earth'] + greenwich
    v = eph['venus']
    ha = e.at(t).observe(v).hadec()
    gha = Angle(degrees=360.0 + ha[0]._degrees)

    print(f'{gha}, {ha[1]}')


def pe(year):
    import ephem

    greenwich = ephem.Observer()
    greenwich.lon = 0.0
    greenwich.lat = ephem.degrees('51:28:38')
    greenwich.pressure = 0.0
    greenwich.horizon = '-0:34'

    t = ephem.date(str(year))
    v = ephem.Venus()
    greenwich.date = t
    greenwich.epoch = t
    v.compute(greenwich)
    gha = ephem.twopi - v.g_ra + greenwich.sidereal_time()

    print(f'{ephem.degrees(gha)}, {ephem.degrees(v.g_dec)}')


if __name__ == '__main__':
    sf(2016)
    pe(2016)

Yields the results:

219deg 41' 57.5", -18deg 37' 03.9"
219:42:15.7, -18:36:56.0

The old PyEphem code agrees with the Nautical Almanac I have in front of me.

This is definitely a PEBKAC, but I'm scratching my head to find what temporal or spatial transformation I've missed.

MerseyViking
  • 389
  • 3
  • 19

1 Answers1

1

It looks like you are asking Skyfield for astrometric positions, but PyEphem for apparent positions. According to the PyEphem documentation:

"g_ra and ra — Apparent geocentric right ascension for the epoch-of-date"

https://rhodesmill.org/pyephem/quick.html#body-compute-date

Whereas with Skyfield, you have to call the .apparent() method on a position to learn the corresponding apparent position; it does not happen automatically:

https://rhodesmill.org/skyfield/positions.html#barycentric-astrometric-apparent

See if that change eliminates most of the difference between coordinates.

Brandon Rhodes
  • 83,755
  • 16
  • 106
  • 147
  • I can't believe I missed that! I had been reading about `.apparent()` as well, and it didn't click for some reason. The HA is close enough for my purposes, being about 1 arcsec different in the example, but the declination is still off by 5 arcsecs which is just a tad too much. Any other thoughts? – MerseyViking Apr 12 '22 at 08:17
  • You might be asking Skyfield for J2000 (ICRS) coordinates, while the quick reference suggests that PyEphem `g_ra` and `g_dec` are always in the coordinates of the specific date you're asking about. Maybe try a date argument to `.radec()`? – Brandon Rhodes Apr 12 '22 at 13:28
  • My bad, I had already tried that but didn't update the question. With current epoch I get -18:37:01.1 and with J2000 I get -18:34:30.7 PyEphem and the Nautical Almanac give me -18:36:56.0 – MerseyViking Apr 13 '22 at 09:08
  • PyEphem doesn't know about leap seconds. Try the year 2000, and see if the results are closer. I suspect the two libraries might agree about that year, then maybe diverge by 1 second each time since 2000 that there's been a leap second? – Brandon Rhodes Apr 14 '22 at 05:18