3
$\begingroup$

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?

$\endgroup$
4
$\begingroup$

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]
$\endgroup$
3
  • 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 '17 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 '17 at 22:39
  • $\begingroup$ Wow, there are lots of goodies there! OK I will take a look, thanks! $\endgroup$
    – uhoh
    Sep 4 '17 at 23:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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