# NAIF and SpiceyPy returning seemingly inaccurate results

I am yet another amateur astrophysics enthusiast creating a solar system simulator. I began by hard-coding planetary bodies and their satellite orbiting characteristics via JPL data tables such as those found at the jpl website. My project has now pivoted into creating a REST API that interacts with SpiceyPy scripts on the fly. Unfortunately, I ran into a perplexing hurdle fairly quickly. In trying to attain the orbital elements for the moon revolving around Earth (or the Earth-Moon barycenter for that matter) I am seeing a slight discrepancy for some of the values, namely the inclination. It appears to be off by 1/10 of a degree, and I am not sure why or how to fix it. Given that the api simply executes python scripts, I am able to test the scripts directly on the command line. Here is the script in question (orbital_elements.py):

import argparse

import naif
import elixir_format as fmt

epi = "\n".join([
'Outputs an Elixir map with the following keys:',
'  pa   Perifocal distance. (periapsis)',
'  e    Eccentricity.',
'  i    Inclination.',
'  O    Longitude of the ascending node.',
'  w    Argument of periapsis.',
'  M    Mean anomaly at epoch.',
'  t0   Epoch.',
'  mu   Gravitational parameter.',
'  nu   True anomaly at epoch.',
'  a    Semi-major axis. A is set to zero if',
'       it is not computable.',
'  T    Orbital period. Applicable only for',
'       elliptical orbits. Set to zero otherwise.',
'',
'The epoch of the elements is the epoch of the input',
'are used to describe all three types (elliptic,',
'hyperbolic, and parabolic) of conic orbits.',
])

parser = argparse.ArgumentParser(
formatter_class=argparse.RawDescriptionHelpFormatter,
description='Get orbital elements for given observer and target bodies.',
epilog=epi
)
help='a utc date')
help='name of primary (observing) body/barycenter')
help='name of orbiting (target) body/barycenter')
help='frame of reference')
choices=['NONE', 'LT', 'LT+S', 'CN', 'CN+S', 'XLT', 'XLT+S', 'XCN', 'XCN+S'],
help='aberrational correction method')

args = parser.parse_args()

meta_kernel_name = 'meta_kernel'

def orbital_elements():

# get elements
elements = naif.orbital_elements( args.date, args.obs, args.targ,
args.frame, args.abcorr         )

# grab output
elements_map = fmt.orbital_elements_map( elements )

#
# Display the results.
#
print( elements_map )

if __name__ == '__main__':
orbital_elements()

The elixir_format module is a module that simply formats output into a string representation of an Elixir map (the server-side language being used for the api, along with the Phoenix framework).

The meta_kernel module is a custom module that is more or less a meta kernel selector, as different api endpoints will require different combinations of kernels to be loaded. In this case the meta kernel being loaded is:

\begindata
PATH_VALUES     = ( '/path/to/kernels' )
PATH_SYMBOLS    = ( 'KERNELS' )
'$$KERNELS/lsk/naif0012.tls', '$$KERNELS/pck/gm_de431.tpc',
'$$KERNELS/pck/pck00010.tpc', '$$KERNELS/spk/planets/de432s.bsp',
)
\begintext

And here is my custom naif module--or at least the relevant functions--referenced in the above file:

import math
import spiceypy
from spiceypy.utils.support_types import SpiceyError

def get_state(date, observer, target, frame='J2000', abcorr='LT+S'):
et = spiceypy.str2et( date )

#
# Compute the apparent state of target as seen from
# observer in the J2000 frame.
#
# targ (str) – Target body name.
# et (Union[float,Iterable[float]]) – Observer epoch.
# ref (str) – Reference frame of output state vector.
# abcorr (str) – Aberration correction flag.
# obs (str) – Observing body name.
#
[state, ltime] = spiceypy.spkezr( target, et, frame, abcorr, observer )

return state

def orbital_elements(date, observer, target, frame='J2000', abcorr='LT+S'):
state = get_state( date, observer, target, frame, abcorr )

et = spiceypy.str2et( date )

mu = spiceypy.bodvrd(observer, 'GM', 1)[1][0]

#
# Compute the orbital elements
#
elements = spiceypy.oscltx(state, et, mu)

return elements

Whew! Now that the code is all there, here are the test runs:

$python3 priv/scripts/orbital_elements.py 2019-01-06T00:00:00 3 Moon --frame=ECLIPJ2000 %{ pa: 342071.326649, e: 0.076903, i: 0.092276, O: 2.032973, w: 0.088103, M: 2.794102, t0: 600004869.184060, mu: 403503.235502, nu: 2.842107, a: 370569.235442, T: 2231312.507324 } I'm using the ECLIPJ2000 frame (instead of the default J2000) for convenience since the test is for the moon around earth. Converting to degrees, the inclination above is 5.29deg which is 0.13deg above the value given in the JPL table (5.16deg at the time of this question). Wikipedia confirms the number in the JPL table. I experimented with the date a bit, and moving 2 years back yields:$ python3 priv/scripts/orbital_elements.py 2016-01-06T00:00:00 Earth Moon --frame=ECLIPJ2000
%{
pa:       367491.005309,
e:            0.051208,
i:            0.088237,
O:            3.044581,
w:            3.195189,
M:            4.256709,
t0:    505310468.184053,
mu:       398600.435436,
nu:            4.167384,
a:       387325.259045,
T:      2398969.430945
}

Converting the radian value to degrees, the inclination is returned as 5.05deg, which is 0.11deg less than what is shown in the JPL table.

There is no change if I change the observer to the Earth-Moon barycenter:

\$ python3 priv/scripts/orbital_elements.py 2016-01-06T00:00:00 3 Moon --frame=ECLIPJ2000
%{
pa:       335220.207195,
e:            0.084074,
i:            0.088237,
O:            3.044581,
w:            3.702010,
M:            3.748739,
t0:    505310468.184053,
mu:       403503.235502,
nu:            3.660563,
a:       365990.590986,
T:      2190086.350198
}

I've been pulling my hair out for days on this. What am I missing? Does the inclination actually vary over time with nutation/precession? Should I actually be taking the JPL's advice to heart when it says, "These mean orbital parameters are not intended for ephemeris computation."?

• Your code includes elements = spiceypy.oscltx(state, et, mu). The caveats for naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/oscltx_c.html explain why this function doesn't work well for non-Keplerian orbits as @david-hammen notes. astronomy.stackexchange.com/questions/28691/… is a similar example of why osculating elements aren't always ideal. – barrycarter Jan 15 at 16:53
• @barrycarter I am aware of special circumstances as I began this project by learning about orbital mechanics from an old college textbook and completed calculations by hand for given exercises. In my case, however, the eccentricity of the Moon's orbit is nowhere near 1 nor is it near Earth's equator. Thank you for the second link, though. The author points out that NASA uses oscelt for argument of perifocus, so I may explore the use of that function instead. This project is at its early stages so complete accuracy is only a distant goal at this point. =] – smola Jan 16 at 1:05