Horizons System vs Skyfield - Difference in Geocentric Vector Table Values

[ghost]

Hi there!

I was hoping someone could help me please. I'm trying to obtain ephemeris information, specifically the position of Mars with a Geocentric Coordinate Centre. In Horizons, this looks as follows:

For Ecliptic J2000 Reference Frame, this gives:

and for ICRF this gives:

I'm trying to replicate the [X, Y, Z] results for both of these in Skyfield but I'm getting slightly different values. My Skyfield code is as follows:

from skyfield.api import Loader
from skyfield.framelib import ecliptic_J2000_frame
load = Loader('.')
ts = load.timescale()
eph = load('de421.bsp')

earth, mars = [eph[x] for x in ('earth', 'mars')]
time = ts.tdb(2023, 10, 4, 0, 0, 0)

geocentric = mars.at(time) - earth.at(time)
print(geocentric)
barycentric = mars.at(time)
print(barycentric)

print("--- Geocentric Frame ---")
x, y, z = geocentric.position.km
print(x)
print(-3.471720294584179E+08)
print(y)
print(-1.439381140200791E+08)
print(z)
print(-6.001575746852054E+07)

print("--- Geocentric Ecliptic J2000 Frame ---")
x, y, z = geocentric.frame_xyz(ecliptic_J2000_frame).km
print(x)
print(-3.471720294584179E+08)
print(y)
print(-1.559335349784188E+08)
print(z)
print(2.191912706248462E+06)

which gives the output

<Geocentric ICRS position and velocity at date t center=399 target=499>
<Barycentric BCRS position and velocity at date t center=0 target=499>
--- Geocentric Frame ---
-347172029.72241116
-347172029.4584179
-143938113.43708077
-143938114.0200791
-60015757.70761283
-60015757.46852054
--- Geocentric Ecliptic J2000 Frame ---
-347172029.72241116
-347172029.4584179
-155933534.53863373
-155933534.9784188
2191912.2549821567
2191912.706248462

I just wanted to check if I'd set up the problem correctly and if you had any solution to this? I've done the same for Barycentric and get different results to Horizons too.

Thank you!

[brandon-rhodes]

I just wanted to check if I'd set up the problem correctly and if you had any solution to this? I've done the same for Barycentric and get different results to Horizons too.

You have set up the problem correctly!

The differences you see — of 264 meters on the x axis, of 440 meters on the y axis, and 451 meters on the z axis — are the differences in the Mars position estimated by the DE421 ephemeris you are loading, and a more recent ephemeris like DE440.

We should note that these differences fit almost inside of the error bars of the DE421 ephemeris — the report on the ephemeris says that ‘The error in the plane-of-sky position of Mars relative to Earth through the end of 2008 is about 300 m.’ So accuracy of better than 300 meters is not to be had from DE421.

The documentation of DE440 cites our several ongoing spacecraft missions to Mars for tightening our estimate of its position — which is apparently now within a meter! One diagram has the caption: ‘Residuals of the Mars orbiter range data against DE440. The rms residual of the MEX ranges is about 2 m and the rms residuals of the MGS, ODY, and MRO ranges are about 0.7 m.’

The more accurate ephemeris is a pretty big download; here is the adjustment you would make to your code:

eph = load('de440.bsp')
earth, mars = [eph[x] for x in ('Earth', 'Mars Barycenter')]

Note that a 440 meter difference amounts to around 0.013% of the Mars radius, so it’s not likely that the difference would affect an observation of Mars made with a telescope, but only radio-telescope interactions with the transmitters of the various Mars probes — and, I suppose, the difference could also be rather important to anyone trying to land something on Mars.

[ghost]

Makes so much sense, thank you very much!