<p>Bad news, this type of SPK file has a different sort of interpolation that is not supported by  the <code>jplephem</code> package (Hermite interpolation vs Chebyshev polynomials). You can find this out by doing:</p>
<pre><code>In [1]: print(len(kernel.segments))                                                                                         
1
In [2]: print(kernel.segments[0].describe())                                                                                
2415020.50..2506696.50  Earth Barycenter (3) -&gt; Unknown Target (392)
  frame=1 data_type=12 source=Sun-EarthMoon L2
</code></pre>
<p>which shows that the segment is of data type <code>12</code>. Looking this up in the <a href="https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/FORTRAN/req/spk.html#Supported%20Data%20Types" rel="nofollow noreferrer">NAIF SPK documentation</a> shows that this is <code>Type 12: Hermite Interpolation --- Equal Time Steps</code>; looking in the <code>jplephem</code> code in the <code>spk.Segment._load()</code> routine, shows the code that produces the error you are seeing and the supported types (<code>2</code> and <code>3</code>). According to the NAIF SPK docs these are <code>Type 2: Chebyshev (position only)</code> and <code>Type 3: Chebyshev (position and velocity)</code>.</p>
<p>I think your options are either:</p>
<ol>
<li>add support for this type of interpolation to <code>jplephem</code> (the math is at least documented in NAIF docs),</li>
<li>see if <a href="https://github.com/AndrewAnnex/SpiceyPy" rel="nofollow noreferrer">SpiceyPy</a> has support for this more unusual type of SPK/BSP,</li>
<li>find an alternative source for an L2 ephemeris.</li>
</ol>
<p>EDIT/UPDATE 2022-01-20:</p>
<p>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 <code>Search only major bodies (planets, satellites, spacecraft etc)</code> and then put in <code>SEMB-L2</code> (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 <a href="https://towardsdatascience.com/tagged/space-science-with-python" rel="nofollow noreferrer">this tutorial series</a>)</p>
<p>To achieve it in <a href="https://github.com/AndrewAnnex/SpiceyPy" rel="nofollow noreferrer">SpicePy</a> 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:</p>
<pre><code>\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'
                  )
</code></pre>
<p>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:</p>
<pre><code>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(&quot;hmsdms&quot;))
06h44m48.226497s +23d01m59.89402451s
</code></pre>
<p>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...</p>

CC BY-SA 4.0
2019-10-29T08:39:59.023Z
Yeah... those are bad news indeed. I'm trying Spice toolkit, hope I will get same results and get used to the new syntax. Thank you!
Luis García Iglesias
CC BY-SA 4.0
2022-01-03T17:58:58.687Z
@LuisGarcíaIglesiasThanks for raising this question. And astrosnapper, thanks for your insights. Did either of you end up finding a good Python package for working with Hermite interpolations in SPICE?
nealmcb
CC BY-SA 4.0
2022-01-20T16:56:14.207Z
@nealmcb I have managed to find a way to get this working via HORIZONS and SpiceyPy if you are still interested.
astrosnapper
CC BY-SA 4.0
2022-01-22T20:31:12.147Z
Good work, @astrosnapper - thank you for updating the answer!
nealmcb
CC BY-SA 4.0
2022-02-18T12:23:43.197Z
If you know how the JPL DE actually *defines* the Lagrange points, please take a look at https://space.stackexchange.com/q/58359/38535
PM 2Ring
CC BY-SA 4.0
2019-10-29T03:20:25.783Z
2022-01-20T16:55:10.333Z
2022-01-20T16:55:10.333Z
astrosnapper
<p>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.</p>

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

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

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

<p>Can someone help me with that or knows any other ephemeris file for L2?</p>

<p>L2 file is here: <a href="https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/lagrange_point/" rel="nofollow noreferrer">https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/lagrange_point/</a></p>

<ephemeris><lagrange-point><astropy>
L2 point ephemeris (celestial mechanics)
Luis García Iglesias
# L2 point ephemeris (celestial mechanics)

(Asked by: username_0 on 2019-10-28)

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/>

**Tags:** ephemeris, lagrange-point, astropy


---

## 1 answers

### Answer by username_1 on 2019-10-29. Score: 5 | Accepted answer

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](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/FORTRAN/req/spk.html#Supported%20Data%20Types) 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](https://github.com/AndrewAnnex/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](https://towardsdatascience.com/tagged/space-science-with-python))

To achieve it in [SpicePy](https://github.com/AndrewAnnex/SpiceyPy) 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...

> **Comments:**
> **username_0**: Yeah... those are bad news indeed. I'm trying Spice toolkit, hope I will get same results and get used to the new syntax. Thank you!  
> **username_2**: @LuisGarcíaIglesiasThanks for raising this question. And username_1, thanks for your insights. Did either of you end up finding a good Python package for working with Hermite interpolations in SPICE?  
> **username_1**: @username_2 I have managed to find a way to get this working via HORIZONS and SpiceyPy if you are still interested.  
> **username_2**: Good work, @username_1 - thank you for updating the answer!  
> **username_3**: If you know how the JPL DE actually *defines* the Lagrange points, please take a look at https://space.stackexchange.com/q/58359/38535  

---

