NAIF y SpiceyPy arrojan resultados aparentemente inexactos

Soy otro aficionado entusiasta de la astrofísica que está creando un simulador del sistema solar. Empecé codificando cuerpos planetarios y sus características de órbita satelital a través de tablas de datos JPL como las que se encuentran en el sitio web jpl . Mi proyecto ahora se ha convertido en la creación de una API REST que interactúa con SpiceyPyguiones sobre la marcha. Desafortunadamente, me encontré con un obstáculo desconcertante con bastante rapidez. Al tratar de obtener los elementos orbitales de la Luna que gira alrededor de la Tierra (o el baricentro Tierra-Luna para el caso), veo una ligera discrepancia en algunos de los valores, a saber, la inclinación. Parece estar errado en 1/10 de grado, y no estoy seguro de por qué o cómo solucionarlo. Dado que la API simplemente ejecuta secuencias de comandos de Python, puedo probar las secuencias de comandos directamente en la línea de comandos. Aquí está el script en cuestión ( 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()

El elixir_formatmódulo es un módulo que simplemente formatea la salida en una representación de cadena de un mapa Elixir (el lenguaje del lado del servidor que se usa para la API, junto con el marco Phoenix).

El meta_kernelmódulo es un módulo personalizado que es más o menos un meta selector de kernel, ya que diferentes puntos finales de API requerirán diferentes combinaciones de kernels para cargar. En este caso, el meta kernel que se está cargando es:

\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

Y aquí está mi naifmódulo personalizado, o al menos las funciones relevantes, a las que se hace referencia en el archivo anterior:

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

¡Uf! Ahora que todo el código está ahí, aquí están las ejecuciones de prueba:

$ 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
}

Estoy usando el ECLIPJ2000marco (en lugar del predeterminado J2000) por conveniencia ya que la prueba es para la luna alrededor de la tierra.

Convirtiendo a grados, la inclinación anterior es de 5,29 grados, que es 0,13 grados por encima del valor dado en la tabla JPL (5,16 grados en el momento de esta pregunta). Wikipedia confirma el número en la tabla JPL. Experimenté un poco con la fecha, y retrocediendo 2 años se obtiene:

$ 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
}

Convirtiendo el valor en radianes a grados, la inclinación se devuelve como 5,05 grados, que es 0,11 grados menos que lo que se muestra en la tabla JPL.

No hay cambio si cambio el observador al baricentro Tierra-Luna:

$ 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
}

Me he estado tirando de los pelos durante días con esto. ¿Qué me estoy perdiendo? ¿La inclinación realmente varía con el tiempo con la nutación/precesión? ¿Debería tomar en serio el consejo del JPL cuando dice: " Estos parámetros orbitales medios no están destinados al cálculo de efemérides "?