Different results for the same JPL Horizons query

Question

On JPL Horizons website, I'm using the following query to get vector state for Earth:

MAKE_EPHEM=YES
COMMAND=3
EPHEM_TYPE=VECTORS
CENTER='500@0'
START_TIME='2023-02-06'
STOP_TIME='2023-02-07'
STEP_SIZE='1 DAYS'
VEC_TABLE='3'
REF_SYSTEM='ICRF'
REF_PLANE='ECLIPTIC'
VEC_CORR='NONE'
OUT_UNITS='KM-S'
VEC_LABELS='YES'
VEC_DELTA_T='NO'
CSV_FORMAT='NO'
OBJ_DATA='YES'

I'm getting two values for 'x', which I guess make perfect sense, being in km:

X =-1.085163786984627E+08
X =-1.103097877035814E+08

I'm trying to have the same query in Python, using jplhorizons from astroquery library, like this:

from astroquery.jplhorizons import Horizons
from datetime import datetime, timedelta
today = datetime.now()
yesterday = today - timedelta(days=1)
obj_earth = Horizons(id='3', location='0', epochs={'start':str(yesterday), 'stop':str(today),'step':'1d'})
print(obj_earth.vectors().columns['x'])

Running this code, I get the following output:

          x           
          AU          
----------------------
-3.582503526151746e-05
-3.948114036370441e-05

The values are apparently in AU, but honestly, they make no sense to me. Shouldn't I get the exact same values running this through web Horizons' interface and through the script --even if the values are in different units?

Answer 1

The issue is with differences in how the defaults behave between the code and the website. In the code version of the API id=3 gives the Earth-Moon Barycenter (as opposed to id=399 which is the Earth itself; not sure which is wanted but it's not the most significant problem). The major problem is that location=0 in the Python code produces the position of Greenwich Observatory relative to the geocenter, which is why the numbers are so small (~2850 km).

To get the desired position of the Earth relative to the Solar System Barycenter requires using location='500@0' in the code (I also used astropy.units to do the conversion to km):

from astroquery.jplhorizons import Horizons
import astropy.units as u
from datetime import datetime, timedelta
today = datetime(2022,2,7)
yesterday = today - timedelta(days=1)
obj_earth = Horizons(id='399', location='500@0', epochs={'start':str(yesterday), 'stop':str(today),'step':'1d'})
print(obj_earth.vectors().columns['x'].to(u.km))

This produces the following output, matching the web query:

[-1.08939344e+08 -1.10724448e+08] km

(As an aside to help diagnose this, you can do obj_earth.uri which shows the query URL string (urlencoded) in the same format as originally posted. This showed that the original Python code was submitting CENTER=0 not the desired CENTER=500@0)