3

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
KERNELS_TO_LOAD=(
   '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)}
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            
*******************************************************************************
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?

astrojuanlu
  • 1,052
  • 7
  • 19

1 Answers1

5

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