Astropy's get_body gives slightly different RA and Dec values than those retrieved from JPL horizons #252
Description
Maybe this issue is really for Astropy, but I want people to cross-check with me before we report this to Astropy. Astropy has a nice integration to get solar system ephemeris quickly in an offline manner (after downloading a ~10-MB database). See https://docs.astropy.org/en/stable/coordinates/solarsystem.html. If this can replace our current read_horizons() function, which invokes an online query, it would be great. Once could use astropy's build-in ephemeris or use the more accurate JPL version after installing jplephem (just do “pip install jplephem”). While it works nicely and fast, I am puzzled that I always get ~20" difference in RA and ~6" difference in Dec. To reproduce my results, here are the steps. [I am using suncasa’s Python 3.8 version. My astropy version is 5.1. ]
# First, use suncasa's helioimage2fits to query JPL Horizons. For simplicity, a geocentric observer is used.
>>> import numpy as np
>>> from astropy.time import Time
>>> t0 = Time('2023-05-02T17:42:55')
>>> from suncasa.utils import helioimages2fits as hf
>>> out = hf.read_horizons(t0, dur=1./60./24., observatory='500')
# out is a dictionary that has ra, dec, etc. parse from JPL horizons query results
>>> print(out)
{'time': [60066.73813657416, 60066.73883101903], 'ra': [0.6861736017226758, 0.6861852028389474], 'dec': [0.2680396049676472, 0.26804322245922796], 'p0': [336.1374, 336.1376], 'delta': [1.00774195240406, 1.00774212515192]}
>>> print(Time(out['time'], format='mjd').isot)
['2023-05-02T17:42:55.000' '2023-05-02T17:43:55.000']
# Now use Astropy to try to get the solar ephemeris
>>> from astropy.coordinates import solar_system_ephemeris, EarthLocation
>>> from astropy.coordinates import get_body
>>> sun_loc = get_body('sun', t0, ephemeris='jpl')
>>> print(sun_loc)
<SkyCoord (GCRS: obstime=2023-05-02T17:42:55.000, obsgeoloc=(0., 0., 0.) m, obsgeovel=(0., 0., 0.) m / s): (ra, dec, distance) in (deg, deg, km)
(39.30927983, 15.35580065, 1.5075605e+08)>
# Compare the two results
# difference in RA and DEC, in arcsec
>>> print((sun_loc.ra.value - np.degrees(out['ra'][0]))*3600)
-20.057608507445934
>>> print((sun_loc.dec.value - np.degrees(out['dec'][0]))*3600)
-6.254842632380075
I verified that helioimages2fits gives exactly the same results as the web API of JPL Horizons (https://ssd.jpl.nasa.gov/horizons/app.html#/). Changing the ephemeris to "buildin" or the latest used by JPL Horizons ("de440s") gives similar results. One difference is that Astropy reports coordinates in the Geocentric Celestial Reference System (GCRS) system (see https://docs.astropy.org/en/stable/api/astropy.coordinates.GCRS.html#astropy.coordinates.GCRS), while JPL Horizons uses the ICRF (International Celestial Reference Frame, equatorial-aligned radio frame) coordinate system (i.e., J2000), which uses the solar system's barycenter as the origin - see https://en.wikipedia.org/wiki/International_Celestial_Reference_System_and_its_realizations for more information. The axes of GCRS and ICRS should share the same orientation and are both non-rotating w.r.t. extragalactic sources. The only difference is where the origin is. [Although there may be very subtle differences in dynamically non-rotating rather than kinematically non-rotating, which, if it is present, can only change by 1.9" per century]. See this document for official definitions and detailed discussions (https://arxiv.org/abs/astro-ph/0602086). Therefore, the differences we see here can not be understood from the coordinate systems used but should have other origins. Given that JPL Horizons have been used widely in the studies of solar system bodies (including us!), I am inclined to believe that the Astropy version has errors...