Discrepancy in State objects depending on how JPL Horizons is queried

[jensimjenseits]

High-level problem description
Hi, I have been using sbpy to generate synchrones and syndynes and noticed some discrepancy in the creation of State objects that I could not find any explanation for.

What did you do?
Basically, I tried creating a State object for my comet in two different ways:

1. I used astroquery to obtain the cartesian heliocentric coordinates of the comet directly from JPL Horizons via Horizons(...).vectors().
2. I used the sbpy convenience functions sbpy.data.from_horizons() and created a State object from that using sbpy.dynamics.state.from_ephem()

Example:

from astroquery.jplhorizons import Horizons
from sbpy.data import Ephem
from sbpy.dynamics import State
from astropy.time import Time
import astropy.units as u

epoch=Time(2459946.0,format='jd')

# using astroquery
q = Horizons(id='311P', location='@sun', epochs=epoch, id_type='smallbody')
v = q.vectors(refplane='earth')
r = u.Quantity([v['x'][0],  v['y'][0],  v['z'][0]],  u.au)
rdot = u.Quantity([v['vx'][0], v['vy'][0], v['vz'][0]], u.au / u.day)
comet = State(r,rdot,epoch,frame='heliocentriceclipticiau76')
print(comet)

# using sbpy functions
eph = Ephem.from_horizons("311P",
                          epochs=epoch,
                          location='@sun',
                          id_type="smallbody")
comet = State.from_ephem(eph,frame='heliocentriceclipticiau76')[0]
print(comet)

Output:
<State (<HeliocentricEclipticIAU76 Frame (obstime=2459946.0)>):
r
[ 1.16927516 -1.85001645 -0.72280141] AU
v
[15.28357029 9.8210348 5.85998821] km / s
t
2459946.0>
<State (<HeliocentricEclipticIAU76 Frame (obstime=2459946.0)>):
r
[ 1.16916454 -1.85008735 -0.72284385] AU
v
[14.56214698 10.21163345 3.68954787] km / s
t
2459946.0>

What did you expect?
The same state vector for both methods.

What did really happen?
The resulting state vectors are quite similar but not identical in their positions (difference about 0.01%), but the velocity vectors differ quite drastically. These differences in the state vectors change the resulting Synchrones and Syndynes significantly. In my case, method 1 yielded the physically plausible projected synchrones and syndynes for my observations (unrelated to the example given here).

When checking the source code of from_ephem(), I saw that the cartesian coordinates are calculated from JPL Horizons' ephemerides using astropy functions (SphericalRepresentation, SphericalDifferential,...). It seems as if these astropy transformations return different cartesian vectors than Horizons gives out when using its vector mode. If I am missing something, maybe an additional input to the from_ephem() function or some frame transformation that I didn't know I have to apply, I would be happy about any hint. Otherwise, since sbpy's from_horizons() currently uses astroquery.jplhorizons.HorizonsClass.ephemerides(), it might also make sense to include another convenience function using astroquery.jplhorizons.HorizonsClass.vectors() for direct access to Horizons' state vectors.

Provide information on your environment:
operating system and version: Debian GNU/Linux 12
sbpy version: 0.6.0
astropy version: 7.2.0
numpy version: 1.26.4

Thank you for your work on sbpy!

[github-actions]

Welcome to sbpy 👋 and thank you for your first issue!

A project member will respond to you as soon as possible; in the meantime, please double-check the (Astropy) guidelines for submitting issues and make sure you've provided the requested details.

GitHub issues in the sbpy repository are used to track bug reports and feature requests; If your issue poses a question about how to use sbpy, please instead raise your question in sbpy Discussions and close this issue.

[mkelley]

Thanks for the report. I'm not sure what is happening yet, but I will continue to work on it.

However, I do see two problems with your example. We have a relevant update to the documentation that is currently in review (#434):

Ephem.from_horizons returns equatorial coordinates in the ICRF reference frame, which has its origin at the Solar System barycenter. For State to correctly convert the ephemeris object to vectors, we need to set the Horizons observer to the Solar System barycenter ("@ssb"). However, dynamical integrations are done in a heliocentric reference frame, so we transform the result to "heliocentriceclipticiau76":

   >>> eph = Ephem.from_horizons(
   ...     "48P",
   ...     id_type="designation",
   ...     closest_apparition=True,
   ...     epochs=Time("2004-10-13T21:08:23.894"),
   ...     location="@ssb",
   ... )
   >>> comet = State.from_ephem(eph, frame="icrs")
   >>> comet = comet.transform_to("heliocentriceclipticiau76")

And these comments about the origin and reference frame also apply to your astroquery example. The location should be "@ssb" and the state object's reference frame should initially be "icrs" (https://astroquery.readthedocs.io/en/latest/jplhorizons/jplhorizons.html#reference-frames):

Here is your example edited to this effect:

from astroquery.jplhorizons import Horizons
from sbpy.data import Ephem
from sbpy.dynamics import State
from astropy.time import Time
import astropy.units as u

epoch=Time(2459946.0,format='jd')

# using astroquery
q = Horizons(id='311P', location='@ssb', epochs=epoch, id_type='smallbody')
v = q.vectors(refplane='earth')
r = u.Quantity([v['x'][0],  v['y'][0],  v['z'][0]],  u.au)
rdot = u.Quantity([v['vx'][0], v['vy'][0], v['vz'][0]], u.au / u.day)
comet = State(r,rdot,epoch,frame='icrs').transform_to("heliocentriceclipticiau76")
print(comet)

# using sbpy functions
eph = Ephem.from_horizons("311P",
                          epochs=epoch,
                          location='@ssb',
                          id_type="smallbody")
comet = State.from_ephem(eph,frame='icrs').transform_to("heliocentriceclipticiau76")[0]
print(comet)

Output:

<State (<HeliocentricEclipticIAU76 Frame (obstime=2459946.0)>):
 r
  [ 1.16927445 -1.9848711   0.07273703] AU
 v
  [15.28357033 11.34159258  1.46985071] km / s
 t
  2459946.0>
<State (<HeliocentricEclipticIAU76 Frame (obstime=2459946.0)>):
 r
  [ 1.16916419 -1.98495278  0.07272649] AU
 v
  [14.505259   10.89243103  1.66418022] km / s
 t
  2459946.0>

The position vectors are still in agreement, so that's good, but the velocity vectors are not. I've tested astropy's conversion with a different solar system ephemeris, but the result was effectively the same.

And I agree on that a path to using Horizons.vectors would be useful. Thanks for the suggestion.

[mkelley]

I've made some progress. I can confirm that the results from Horizons in Cartesian and spherical representations are in agreement, but that passing the spherical representation through State.from_ephem causes the issue:

import numpy as np
import astropy.units as u
from astropy.time import Time
from astroquery.jplhorizons import Horizons
from spiceypy import furnsh, spkez

import sbpy.time
from sbpy.data import Ephem
from sbpy.dynamics import State

epoch = Time(2459946.0, format='jd')

# cartesian
q = Horizons(id='311P', location='@0', epochs=epoch, id_type='smallbody')
tab1 = q.vectors(refplane='earth')
tab1["k1"].unit = u.dimensionless_unscaled
tab1["date"] = Time(tab1["datetime_jd"], format="jd")
eph = Ephem.from_table(tab1)
a = State.from_ephem(eph, frame="icrs")
av = np.linalg.norm(a.v.unmasked.value[0])

# spherical
eph = Ephem.from_horizons("311P",
                          epochs=epoch,
                          location='@ssb',
                          id_type="smallbody")
b = State.from_ephem(eph, frame='icrs')
bv = np.linalg.norm(b.v.unmasked.value[0])

delta = eph["delta"].unmasked[0]
ephv = np.linalg.norm([
    (delta * np.sin(eph["RA*cos(Dec)_rate"].unmasked * u.hr) / u.hr).to_value("km/s"),
    (delta * np.sin(eph["DEC_rate"].unmasked * u.hr) / u.hr).to_value("km/s"),
    eph["delta_rate"].unmasked.value
])

print("Speed from Cartesian vectors:", av)
print("Speed from Spherical vectors via State:", bv)
print("Speed from Spherical vectors:", ephv)

Output:

Speed from Cartesian vectors: 19.08080607836211
Speed from Spherical vectors via State: 18.20785571585951
Speed from Spherical vectors: 19.080684597624735

[jensimjenseits]

Thank you for your corrections about observer location and reference frames.

Yes, I had also confirmed that the disagreement occurs when State.from_ephem is applied. The astropy functions used in State.from_ephem are SphericalRepresentation.to_cartesian and SphericalDifferential.to_cartesian (and also implicitly SphericalRepresentation.represent_as when putting a base in the argument of to_cartesian for the conversion of the velocities), which itself refer to the ERFA library using ufunc.s2p for the conversion from spherical to cartesian. Maybe somewhere along this line of dependencies there is a difference in convention to the way that Horizons is transforming between representations?

[mkelley]

I've got it now. The issue is that sbpy was using the RA*cos(Dec) values and SphericalDifferential. It should instead be using SphericalCosLatDifferential. Continuing from the last example:

import astropy.coordinates.representation as cr

c = cr.SphericalRepresentation(eph["ra"], eph["dec"], eph["Delta"])
d = cr.SphericalDifferential(
    eph["RA*cos(Dec)_rate"],
    eph["Dec_rate"],
    eph["deltadot"],
).to_cartesian(base=c)[0]
print("Speed from SphericalDifferential:", d.norm().to("km/s"))

d = cr.SphericalCosLatDifferential(
    eph["RA*cos(Dec)_rate"],
    eph["Dec_rate"],
    eph["deltadot"],
).to_cartesian(base=c)[0]
print("Speed from SphericalCosLatDifferential:", d.norm().to("km/s"))

Speed from SphericalDifferential: 18.20785571585951 km / s
Speed from SphericalCosLatDifferential: 19.08068470247059 km / s

I'll make a fix. Thanks again for the report.