L2 point ephemeris (celestial mechanics)

Question

I'm a master student and I am trying to get the L2 ephemeris for some calculations in my master's project. It was kinda hard to find a file with the L2 ephemeris, but once I found it, I tried the next syntax.

from jplephem.spk import SPK
kernel = SPK.open('L2_de431.bsp')
position = kernel[3,392].compute(2457061.5)

but I'm getting an exception: "Only SPK data types 2 and 3 are supported".

I tried the same syntax with a planetary ephemeris file (de432s.bsp) and it works fine.

Can someone help me with that or knows any other ephemeris file for L2?

L2 file is here: https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/lagrange_point/

Answer 1

Bad news, this type of SPK file has a different sort of interpolation that is not supported by the jplephem package (Hermite interpolation vs Chebyshev polynomials). You can find this out by doing:

In [1]: print(len(kernel.segments))                                                                                         
1
In [2]: print(kernel.segments[0].describe())                                                                                
2415020.50..2506696.50  Earth Barycenter (3) -> Unknown Target (392)
  frame=1 data_type=12 source=Sun-EarthMoon L2

which shows that the segment is of data type 12. Looking this up in the NAIF SPK documentation shows that this is Type 12: Hermite Interpolation --- Equal Time Steps; looking in the jplephem code in the spk.Segment._load() routine, shows the code that produces the error you are seeing and the supported types (2 and 3). According to the NAIF SPK docs these are Type 2: Chebyshev (position only) and Type 3: Chebyshev (position and velocity).

I think your options are either:

1. add support for this type of interpolation to jplephem (the math is at least documented in NAIF docs),
2. see if SpiceyPy has support for this more unusual type of SPK/BSP,
3. find an alternative source for an L2 ephemeris.

EDIT/UPDATE 2022-01-20:

I have found a way to do it with SpiceyPy and with JPL HORIZONS. In JPL HORIZONS, when you specify Target Body, if you set it to Search only major bodies (planets, satellites, spacecraft etc) and then put in SEMB-L2 (Sun, Earth-Moon Barycenter L2 Lagrange point) you can can use all of the normal HORIZONS functionality to calculate vectors and observables etc. (Some of the info on using SpiceyPy came from this tutorial series)

To achieve it in SpicePy you will need the following SPICE kernels (the main DE430 SPK ephemeris, the L2 SPK from the naif link in the question and the LSK leap second kernel) listed in this kernel meta file:

\begintext

    This meta file contains the relative paths to all needed SPICE kernels.

\begindata

KERNELS_TO_LOAD = (
    'spk/de430.bsp',
    'lsk/naif0012.tls',
    'pck/pck00010.tpc',
    'spk/L2_de431.bsp'
                  )

Some basic code to calculate vectors and RA, Dec position for the first time of a set of times; the code can be extended to make a numpy or Pandas Dataframe for multiple positions as needed:

from datetime import datetime
import spiceypy
import numpy as np
import astropy.units as u
from astropy.coordinates import SkyCoord

# Load the SPICE kernels via a meta file
spiceypy.furnsh('kernel_meta.txt')

# Create dictionary of body names and their SPICE ids
solsys_dict = {'SSB' : 0, 'SUN' : 10, 'EARTH' : 399, 'L2' : 392 }

# Create start and end datetime objects and convert to strings
start_time_UTC = datetime(year=2022, month=1, day=1, \
                                  hour=0, minute=0, second=0)
start_time_UTC_str = start_time_UTC.strftime('%Y-%m-%dT%H:%M:%S')

end_time_UTC = datetime(2027,1,1,0,0,0)
end_time_UTC_str = end_time_UTC.strftime('%Y-%m-%dT%H:%M:%S')

# Convert to Ephemeris Time (ET) using the SPICE function utc2et
start_time_ET = spiceypy.utc2et(start_time_UTC_str)
end_time_ET = spiceypy.utc2et(end_time_UTC_str)
# Create daily span of ET time intervals
delta_days = end_time_UTC-start_time_UTC
time_interval_et = np.linspace(start_time_ET,end_time_ET, delta_days.days)

# Calculate X,Y,Z position vector (and light travel time) in ecliptic frame
posvec_ecl, ltt = spiceypy.spkezp(targ=solsys_dict['L2'], et=time_interval_et[0], ref='ECLIPJ2000',abcorr='LT',obs=solsys_dict['EARTH'])
# Print vector (in km) and LTT (in seconds)
print(posvec_ecl)
print(ltt)
[-2.64321172e+05  1.45478338e+06 -1.66970236e+02]
4.932081466911206

# Calculate X,Y,Z position vector (and light travel time) in J2000 frame
posvec, ltt = spiceypy.spkezp(targ=solsys_dict['L2'], et=time_interval_et[0], ref='J2000',abcorr='LT',obs=solsys_dict['EARTH'])
# Convert position vector to RA, Dec and distance
dist,ra,dec = spiceypy.recrad(posvec)
# Convert to an AstroPy SkyCoordinate and print out
coord = SkyCoord(ra,dec,unit=u.rad)
print(coord.to_string("hmsdms"))
06h44m48.226497s +23d01m59.89402451s

I've not done extensive checking but the numbers seem to match the vector and observer table outputs from HORIZONS - usual caveats about getting correct time systems and coordinate frames apply...