How can I get Skyfield to agree with the Nautical Almanac for the Sun's Declination?

Question

Here's a sample script that calculates the sun's declination on 2016/7/23 at 00:00:00 GMT using both PyEphem and Skyfield:

import ephem

sun1 = ephem.Sun('2016/7/23 00:00:00')
dec1 = sun1.dec

print 'PyEphem Declination:', dec1

#-------------------------------------
from skyfield.api import load

planets = load('de421.bsp')

earth = planets['earth']
sun2 = planets['sun']

ts = load.timescale()
time2 = ts.utc(2016, 7, 23)

dec2 = earth.at(time2).observe(sun2).apparent().radec()[1]

print 'Skyfield Declination:', dec2

When I run this I get:

PyEphem Declination: 20:01:24.0
Skyfield Declination: +20deg 04' 30.0"

The Nautical Almanac gives 20deg 01.4' for that time. I'm not sure what I'm doing incorrectly to cause this discrepancy. Thanks!

P.S. I'm using Python 2.7, and Skyfield 0.8.

Answer 1

The answer that PyEphem is giving you is in exact agreement with the almanac, but expressed in traditional hours-minutes-seconds of arc instead of in hours and decimal minutes. The fraction .4 that is part of the minutes-of-arc quantity 1.4 becomes, if expressed as arcseconds, 60 × 0.4 = 24 arcseconds. So:

20°1.4′ = 20°1′24″

Skyfield is giving you, by default, coordinates in the permanent GCRS coordinate system that is the updated replacement for J2000. But the Almanac does not use the equator-and-equinox coordinate system of the year 2000, but of the current year — and in fact the exact moment ­— of each datum that it reports. To ask Skyfield to express the declination in year-2016 coordinates, give it the epoch value "date":

from skyfield.api import load

ts = load.timescale()
planets = load('de421.bsp')

earth = planets['earth']
sun = planets['sun']

t = ts.utc(2016, 7, 23)
ra, dec, distance = earth.at(t).observe(sun).apparent().radec(epoch='date')
print(dec)

The result:

+20deg 01' 24.0"