# Discrepancy between HORIZONS and SPICE

I am trying to reproduce the output of the HORIZONS system using the SPICE toolkit from NASA, but I am not getting exactly the same results. I am using this script to compute the position and velocity of the Venus barycenter in the J2000 frame, using the spkez function (Python docs).

import spiceypy as spice

SSB_ID = 0
VENUS_ID = 299

FRAME = 'J2000'  # Earth mean equator and equinox of J2000
#FRAME = 'ECLIPJ2000'  # Mean ecliptic and equinox of J2000

print(spice.tkvrsn("TOOLKIT"))

spice.furnsh("./cassMetaK.txt")

utc = '2017-09-01 12:05:50'
et = spice.str2et(utc)

state, lightTimes = spice.spkez(VENUS_ID, et, 'J2000', 'NONE', SSB_ID)

position = state[:3]
velocity = state[3:]

print(position)
print(velocity)


My kernels:

\$ cat cassMetaK.txt
\begindata
'naif0009.tls',
'de431.bsp',
)
\begintext


These are the values of position and velocity that I am getting:

CSPICE_N0066
[ 19217297.26117753  97944611.05090818  42842518.49871409]
[-34.60894367   4.64339943   4.27861046]


And notice the slight differences with the HORIZONS output:

*******************************************************************************
Ephemeris / WWW_USER Mon Sep  4 00:46:19 2017 Pasadena, USA      / Horizons
*******************************************************************************
Target body name: Venus (299)                     {source: DE431mx}
Center body name: Solar System Barycenter (0)     {source: DE431mx}
Center-site name: BODY CENTER
*******************************************************************************
Start time      : A.D. 2017-Sep-01 12:05:50.0000 TDB
Stop  time      : A.D. 2017-Sep-02 00:00:00.0000 TDB
Step-size       : 1440 minutes
*******************************************************************************
Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 0.00000000,0.00000000,0.0000000 {E-lon(deg),Dxy(km),Dz(km)}
Output units    : KM-S
Output type     : GEOMETRIC cartesian states
Output format   : 3 (position, velocity, LT, range, range-rate)
Reference frame : ICRF/J2000.0
Coordinate systm: Earth Mean Equator and Equinox of Reference Epoch
*******************************************************************************
JDTDB
X     Y     Z
VX    VY    VZ
LT    RG    RR
*******************************************************************************
$$SOE 2457998.004050926 = A.D. 2017-Sep-01 12:05:50.0000 TDB X = 1.921958776713952E+07 Y = 9.794430371601710E+07 Z = 4.284223531920820E+07 VX=-3.460881129595693E+01 VY= 4.644081825248464E+00 VZ= 4.278909130085439E+00 LT= 3.623116682955181E+02 RG= 1.086183056003940E+08 RR=-2.484635847985155E-01$$EOE


I know the errors are under 0.1 %, but it is very important to me to know how exactly is HORIZONS computing those values, if not with the code I am using. What can be the source of those discrepancies?

Based on this answer in spice-discussion, one must input the time explicitly as TDB, like this:

tdb = '2017-09-01 12:05:50 TDB'
et = spice.str2et(tdb)


Then the results match exactly:

CSPICE_N0066
[ 19219587.76713952  97944303.7160171   42842235.3192082 ]
[-34.6088113    4.64408183   4.27890913]

• Bravo! It's great to see python posted. You can also explore different timescales in python using Skyfield's Timescale class.
– uhoh
Sep 4 '17 at 17:43
• I usually use astropy.time for that. If you want to see more Python, make sure you check out this ;) docs.poliastro.space Sep 4 '17 at 22:39
• Wow, there are lots of goodies there! OK I will take a look, thanks!
– uhoh
Sep 4 '17 at 23:15