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



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:]


My kernels:

$ cat cassMetaK.txt 

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

[ 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)}
Center radii    : (undefined)                                                  
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            
   X     Y     Z
   VX    VY    VZ
   LT    RG    RR
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

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?


1 Answer 1


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:

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.