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."?
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$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$