SPK type 3 interface returns inconsistent and confusing velocity units, causes downstream bugs (example in astropy)

[vancky]

Problem

In the current design of jplephem's SPK interface, the compute and compute_and_differentiate methods behave inconsistently between type 2 and type 3 segments. This inconsistency leads to confusion and actual bugs in user code, such as in astropy's solar system ephemeris calculations.

Details

• For type 2 segments, both methods work as expected—returning 3D position and 3D velocity, respectively, with well-documented units.
• For type 3 segments, compute returns 6 components: [x, y, z, dx, dy, dz], where dx, dy, dz are velocity components. However, the velocity units are always in km/s, contrary to the user's assumption (often expecting km/day, as with type 2).
• Moreover, the compute_and_differentiate method returns the 6D components and their derivatives, but only the first 6D output is the correct position and velocity; the derivatives output is not suitable for physical velocity interpretation.

Real Bug Example in Astropy

In astropy/astropy/files/astropy/coordinates/solar_system.py, in the function:

def _get_body_barycentric_posvel(body, time, ephemeris=None, get_velocity=True)

the code for type 3 segments is:

if spk.data_type == 3:
    # Type 3 kernels contain both position and velocity.
    posvel = spk.compute(jd1, jd2)
    if get_velocity:
        body_posvel_bary += posvel.reshape(body_posvel_bary.shape)
    else:
        body_posvel_bary[0] += posvel[:3]

Here, the velocity (posvel[3:6]) is incorrectly assumed to be in km/day, resulting in erroneous output because the actual unit is km/s.

Why This Is Harmful

• The API exposes the internal type-specific format in the external interface, requiring users to understand SPK internals to correctly interpret output.
• Downstream consumers (like astropy) create bugs from incorrect unit assumptions.
• Documentation does not call out this subtlety loudly enough or offer unified, type-independent results.

Suggested Solution

• Ensure compute always returns just the position vector (3D), with consistent units, for both type 2 and type 3 segments.
• Ensure compute_and_differentiate always returns a (position, velocity) tuple, with velocity units explicitly documented (preferably standardized, e.g. always km/s or km/day, with correct conversion for type 3).
• Update documentation to highlight this change and warn about prior inconsistencies for downstream libraries.
• Consider providing automated unit conversion internally so users do not have to worry about SPK kernel types or units.

Impacted Users

• Any code that relies on jplephem/SPK, including major astronomy projects like astropy, is affected.
• Example code in astropy directly demonstrates the real-world impact and propagation of the bug.

[brandon-rhodes]

Since ephemeris files don't, as far as I know, specify any units, I don't see any way that jplephem could change its behavior — if you'll take a look at the Segment code in spk.py, you'll see that it's simply reading out the numbers in the ephemeris arrays and reporting them back; it doesn't know what they actually mean.

But let's double-check: which Type 3 ephemeris are you loading up? That will give us a specific case to look at. Thanks!

[vancky]

Thank you for the quick reply! I think there is a misunderstanding, so let me clarify the issue more precisely.

Clarification About Units in jplephem vs. SPK Kernels

According to jplephem’s own documentation, the velocity will by default be distance traveled per day, in whatever units for distance the ephemeris happens to use.
However, according to SPICE documentation, SPK kernel data always use second as its time unit.

You have also noted this distinction Velocities units · Issue #19 · brandon-rhodes/python-jplephem.

Where the Problem Actually Comes From

Inside the generate() method:

• When you yield components, you do not apply any unit conversion.
  For type-3 segments, this means the last 3 components (dx, dy, dz) are left in km/s, straight from the kernel.

• When you yield rates, you do apply time-unit conversion:

  rates /= intlen
  rates *= 2.0
  rates *= S_PER_DAY    # converts from per-second to per-day

  So the differential output is implicitly converted to per day, not per second.

This leads to an internal inconsistency inside jplephem:

• components[3:6] → km/s

• rates[...] → km/day

Yet the documentation states that results use day-based units.

Why This Becomes a Real Bug

Downstream code (e.g., in Astropy) assumes that:

• If jplephem returns velocity, it is in km/day, because that is the documented behavior for type-2 segments.

So in Astropy’s type-3 code, they simply take:

posvel = spk.compute(jd1, jd2)
velocity = posvel[3:6]   # assumed to be km/day

But with a Type 3 kernel, this is actually km/s, leading to wrong barycentric velocities by a factor of 86,400.

Summary

The issue is not that SPK kernels specify units inconsistently — they don't.
The issue is that jplephem mixes units internally:

• components → seconds

• rates → days

while documentation indicates day-based units.

For Type3, https://ssd.jpl.nasa.gov/ftp/eph/satellites/bsp/ura116.bsp for a ref.
I hope this clarifies the problem.