2
$\begingroup$

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

from meta_kernel import load as load_mk, unload as unload_mk


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',
    'state. Units are km, rad, rad/sec. The same elements',
    '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
)
parser.add_argument('date', metavar='date',
                    help='a utc date')
parser.add_argument('obs', metavar='observer',
                    help='name of primary (observing) body/barycenter')
parser.add_argument('targ', metavar='target',
                    help='name of orbiting (target) body/barycenter')
parser.add_argument('--frame', default='J2000',
                    help='frame of reference')
parser.add_argument('--abcorr', default='LT+S',
                    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():
    load_mk( meta_kernel_name )

    # 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 )

    unload_mk( meta_kernel_name )


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_TO_LOAD = (
                  '$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."?

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – barrycarter Jan 15 at 16:53
  • $\begingroup$ @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. =] $\endgroup$ – smola Jan 16 at 1:05
2
$\begingroup$

Should I actually be taking the JPL's advice to heart when it says, "These mean orbital parameters are not intended for ephemeris computation."?

When you request orbital elements from Horizons, it does not return those mean orbital elements. It instead returns osculating orbital elements. The conversion between Cartesian states (position and velocity) and osculating orbital elements is rather mechanical. The expressions used assume a Keplerian orbit. What if the object isn't in a Keplerian orbit? For example, one could apply those equations to calculate Neptune's orbital elements as it orbits the Earth. The results will of course be pure garbage because Neptune isn't in a Keplerian orbit about the Earth.

Here's the problem: While the Moon does indeed orbit the Earth, that orbit is markedly non-Keplerian due to perturbations by the Sun.

$\endgroup$
  • $\begingroup$ Your example involving Neptune seems rather extreme, but what you are saying implies that none of the satellite orbits for any planet in any solar system can have accurate ephemeris (or Keplerian orbits) determined from a single state vector because the star(s) and other bodies in the system create perturbations in the orbits? I suppose this makes sense as the positions/velocities of all bodies in a system are in constant flux and probably never (or very seldom) repeat. Maybe I should have asked how the sample size is determined for the mean elements in the JPL tables... $\endgroup$ – smola Jan 16 at 1:00
  • $\begingroup$ @smola - An accurate ephemeris and Keplerian orbits are two different things. None of the planets have Keplerian orbits. The solar system is an example of the N-body problem. (It is the quintessential N-body problem.) The planets interact with one another gravitationally as well as with the Sun. That the planets' orbits are not Keplerian does not mean all hope is lost. Space agencies around the world have sent probes to other bodies. They wouldn't be able to do that if those bodies' orbits couldn't be accurately predicted. $\endgroup$ – David Hammen Jan 16 at 4:39
  • $\begingroup$ Ah yes, that's what I was deducing there. And, in reading the example given in the documentation of oscltx_c, the example given mentions that the result "...produces a set of osculating elements describing the nominal orbit that the spacecraft would follow in the absence of all other bodies in the solar system." I am now convinced that trying to match the numbers in the JPL exactly was a futile exercise. Thanks! $\endgroup$ – smola Jan 16 at 17:37

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.