NAIF und SpiceyPy liefern scheinbar ungenaue Ergebnisse

Ich bin ein weiterer Amateur-Astrophysik-Enthusiast, der einen Sonnensystem-Simulator entwickelt. Ich begann damit, planetare Körper und ihre Satellitenumlaufeigenschaften über JPL-Datentabellen, wie sie auf der jpl-Website zu finden sind, fest zu codieren . Mein Projekt hat sich nun darauf konzentriert, eine REST-API zu erstellen, die mit SpiceyPy interagiertSkripte im laufenden Betrieb. Leider stieß ich ziemlich schnell auf eine verwirrende Hürde. Beim Versuch, die Orbitalelemente für den Mond zu ermitteln, der sich um die Erde dreht (oder das Erde-Mond-Schwerzentrum für diese Angelegenheit), sehe ich eine leichte Diskrepanz für einige der Werte, nämlich die Neigung. Es scheint um 1/10 Grad ausgeschaltet zu sein, und ich bin mir nicht sicher, warum oder wie ich es beheben soll. Da die API einfach Python-Skripte ausführt, kann ich die Skripte direkt auf der Befehlszeile testen. Hier ist das betreffende Skript ( 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()

Das elixir_formatModul ist ein Modul, das die Ausgabe einfach in eine Zeichenfolgendarstellung einer Elixir-Karte formatiert (die serverseitige Sprache, die zusammen mit dem Phoenix-Framework für die API verwendet wird).

Das meta_kernelModul ist ein benutzerdefiniertes Modul, das mehr oder weniger ein Meta-Kernel-Selektor ist, da für verschiedene API-Endpunkte unterschiedliche Kombinationen von Kerneln geladen werden müssen. In diesem Fall ist der geladene Meta-Kernel:

\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

Und hier ist mein benutzerdefiniertes naifModul – oder zumindest die relevanten Funktionen – auf die in der obigen Datei verwiesen wird:

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

Wütend! Nun, da der Code vollständig ist, hier sind die Testläufe:

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

Ich verwende den ECLIPJ2000Rahmen (anstelle des Standardwerts J2000) der Einfachheit halber, da der Test für den Mond um die Erde ist.

In Grad umgerechnet beträgt die obige Neigung 5,29 Grad, was 0,13 Grad über dem in der JPL-Tabelle angegebenen Wert liegt (5,16 Grad zum Zeitpunkt dieser Frage). Wikipedia bestätigt die Zahl in der JPL-Tabelle. Ich habe ein wenig mit dem Datum experimentiert, und das Verschieben um 2 Jahre zurück ergibt:

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

Wenn man den Bogenmaßwert in Grad umwandelt, wird die Neigung als 5,05 Grad zurückgegeben, was 0,11 Grad weniger ist als in der JPL-Tabelle angezeigt.

Es gibt keine Änderung, wenn ich den Beobachter auf das Baryzentrum Erde-Mond ändere:

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

Ich reiße mir seit Tagen damit die Haare aus. Was vermisse ich? Ändert sich die Neigung tatsächlich im Laufe der Zeit mit Nutation/Präzession? Sollte ich mir eigentlich den Rat des JPL zu Herzen nehmen, wenn es heißt: " Diese mittleren Orbitalparameter sind nicht für die Berechnung von Ephemeriden gedacht "?