Numerische Ephemeridenrechnung

Zusammenfassung

Der beschriebene Rechenkern berechnet scheinbare topozentrische Ephemeriden von Asteroiden, Kometen und interstellaren Objekten aus oskulierenden Bahnelementen. Elliptische, parabolische und hyperbolische Anfangszustände werden ohne externen Orbitdienst behandelt. Die anschließende numerische Integration berücksichtigt solare Gravitation, die erste relativistische Korrektur, Planeten, Erde und Mond, Erd-J2, optional 16 massereiche Asteroiden sowie Marsden-A1/A2/A3. Planeten- und Asteroidenzustände werden mit kubischer Hermite-Interpolation aus kompakten Tabellen gewonnen.

Die beobachtete Richtung entsteht nach UTC-nach-TDB-Näherung, Lichtlaufzeit, Erdrotation, Nutation und Präzession. Sechs Extremobjekte mit 20 eingefrorenen JPL-Horizons-Zeilen bilden den gemeinsamen Abnahmetest aller 49 eigenständigen Zielprogramme.

1. Fachbegriffe und Symbole

Alle Fachbegriffe und Abkürzungen, die für das Verständnis der folgenden Rechnung erforderlich sind, werden hier zusammenhängend erklärt. Bei der ersten Verwendung besonders wichtiger Begriffe folgt zusätzlich eine kurze Erklärung im Fließtext.

Glossar der astronomischen und numerischen Fachbegriffe.

Die Symbole μ⊙ und μE werden im gesamten Artikel eindeutig unterschieden.

2. Konventionen, μ und Grundgrößen

Der dynamische Zustand y=[x,y,z,vx,vy,vz] ist heliozentrisch, rechtshändig und auf das äquatoriale J2000-System bezogen. Positionen sind AE, Geschwindigkeiten AE/Tag, Beschleunigungen AE/Tag². Interne Winkel sind Radiant. Bahnelemente, geographische Eingaben und Ergebniswinkel sind Grad.

Das Symbol μ bezeichnet grundsätzlich einen Gravitationsparameter μ=G·M, also das Produkt aus Gravitationskonstante G und Masse M. Ohne Körperindex wurde es in älteren Fassungen missverständlich für die Sonne verwendet. Dieser Artikel schreibt deshalb immer μ⊙=G·M⊙ für die Sonne und μE=G·ME für die Erde. Für einen Störkörper b gilt ηb=μb/μ⊙=Mb/M⊙.

3. Anfangszustand aus oskulierenden Elementen

Eingaben sind Epoche t0, Exzentrizität e, Periheldistanz q, Perihelzeit Tp, Knoten Ω, Perihelargument ω und Inklination i. Oskulierend bedeutet, dass diese Elemente an der Epoche dieselbe Position und Geschwindigkeit wie die gestörte Bahn darstellen. Mit μ⊙=G·M⊙ und Δ=t0−Tp wird zuerst die wahre Anomalie ν bestimmt.

3.1 Ellipse

Gültigkeit: e < 1 - 1e-8
a = q / (1 - e)
n = sqrt(μ⊙ / a^3)
M = normalize(n Δ, -π, +π)
E0 = M                    für e < 0.8
E0 = sign(M) π            sonst
Ek+1 = Ek - (Ek - e sin Ek - M) / (1 - e cos Ek)
Abbruch: |Korrektur| < 2e-14 oder nach 80 Iterationen
ν = 2 atan2(sqrt(1+e) sin(E/2), sqrt(1-e) cos(E/2))

3.2 Hyperbel

Gültigkeit: e > 1 + 1e-8
a = q / (e - 1)                    (positive Betragskonvention)
M = sqrt(μ⊙ / a^3) Δ
H0 = asinh(M / e)
Hk+1 = Hk - (e sinh Hk - Hk - M) / (e cosh Hk - 1)
Abbruch: |Korrektur| < 2e-14 oder nach 80 Iterationen
ν = 2 atan2(sqrt(e+1) sinh(H/2), sqrt(e-1) cosh(H/2))

3.3 Parabel und nahezu parabolischer Grenzbereich

Gültigkeit: |e - 1| <= 1e-8
s = sqrt(2 q^3 / μ⊙)
B = Δ / s
D = 2 sinh(asinh(1.5 B) / 3)
ν = 2 atan D

Die Barker-Form verhindert die numerisch schlechte Division durch 1−e oder e−1 im unmittelbaren Parabelbereich. Anschließend gelten für alle drei Fälle:

p = q (1 + e)
r = p / (1 + e cos ν)
v0 = sqrt(μ⊙ / p)
r_orb = [r cos ν, r sin ν, 0]
v_orb = [-v0 sin ν, v0 (e + cos ν), 0]

Die Drehung aus der Bahnebene verwendet Rz(Ω) Rx(i) Rz(ω). Der ausgeschriebene zweidimensionale Eingang [x,y] wird zu:

X = (cosΩ cosω - sinΩ sinω cosi) x + (-cosΩ sinω - sinΩ cosω cosi) y
Y = (sinΩ cosω + cosΩ sinω cosi) x + (-sinΩ sinω + cosΩ cosω cosi) y
Z = (sinω sini) x + (cosω sini) y

Die feste Ekliptik-nach-Äquator-Rotation ist eine x-Achsenrotation um 84381.448 Bogensekunden und wird identisch auf Position und Geschwindigkeit angewandt.

4. Störkörperdaten und Hermite-Interpolation

Jeder Tabellenknoten enthält Position r0 und Geschwindigkeit v0. Für zwei Knoten im Abstand h und u=(t−t0)/h interpoliert der Kern mit den folgenden Polynomen. Eine kubische Hermite-Interpolation übernimmt an beiden Enden Position und Geschwindigkeit exakt; dies sind die vollständigen Reihen und ihre exakten numerischen Koeffizienten:

h00(u) =  2 u^3 - 3 u^2 + 1
h10(u) =    u^3 - 2 u^2 + u
h01(u) = -2 u^3 + 3 u^2
h11(u) =    u^3 -   u^2

r(u) = h00 r0 + h10 h v0 + h01 r1 + h11 h v1

dh00/dt = ( 6 u^2 - 6 u) / h
dh10/dt =   3 u^2 - 4 u + 1
dh01/dt = (-6 u^2 + 6 u) / h
dh11/dt =   3 u^2 - 2 u
v(u) = (dh00/dt) r0 + dh10 v0 + (dh01/dt) r1 + dh11 v1

5. Bewegungsgleichung

d/dt [r, v] = [v, a(r,v,t)]

Die solare Grundbeschleunigung ist aS=−μ⊙ r/|r|³. Für einen Störkörper mit heliozentrischer Position rb und Massenverhältnis ηb wird der direkte und indirekte Term gemeinsam addiert:

ab = μ⊙ ηb [(rb - r)/|rb - r|^3 - rb/|rb|^3]

Erde und Mond werden innerhalb 0,2 AE getrennt verwendet, außerhalb als Erde-Mond-Baryzentrum. Die acht übrigen Körper Merkur, Venus, Mars, Jupiter, Saturn, Uranus, Neptun und Pluto wirken stets. N16 wird nur bei gesetztem Schalter und vollständigen Daten ergänzt.

5.1 Relativistische Sonnenkorrektur

arel = μ⊙ / (c^2 |r|^3) * [ (4 μ⊙/|r| - |v|^2) r + 4 (r·v) v ]

5.2 Erd-J2

Innerhalb 0,02 AE von der Erde wird die übliche J2-Beschleunigung im erdfesten Näherungsrahmen des Ports als direkter minus indirekter Term addiert:

μE = G·ME = 8.887692445125148326e-10 AE^3/Tag^2
k = 1.5 J2 μE RE^2 / ρ^5
aJ2,x = k x (5 z^2/ρ^2 - 1)
aJ2,y = k y (5 z^2/ρ^2 - 1)
aJ2,z = k z (5 z^2/ρ^2 - 3)

5.3 Nichtgravitative Beschleunigung

g(r) = α (r/r0)^(-m) [1 + (r/r0)^n]^(-k)
Standard: r0=2.808, α=0.1112620426, m=2.15, n=5.093, k=4.6142
ang = g(r) [A1 er + A2 et + A3 en]

er ist radial, en liegt normal zur Bahnebene in Richtung r×v, und et=en×er ist transversal. Sind A1=A2=A3=0, wird exakt der Nullvektor geliefert.

6. Dormand-Prince 5(4)

Ein eingebettetes Runge-Kutta-Paar erzeugt gleichzeitig eine Lösung fünfter und vierter Ordnung. Die Differenz schätzt den lokalen Fehler. Eine Butcher-Tafel ist die tabellarische Sammlung der Stufenzeitpunkte ci, der Zwischengewichte aij und der Endgewichte bi eines Runge-Kutta-Verfahrens. Sie beschreibt damit den Integrationsschritt vollständig. Sämtliche im Code verwendeten Butcher-Koeffizienten sind:

c2=1/5;       a21=1/5
c3=3/10;      a31=3/40; a32=9/40
c4=4/5;       a41=44/45; a42=-56/15; a43=32/9
c5=8/9;       a51=19372/6561; a52=-25360/2187;
              a53=64448/6561; a54=-212/729
c6=1;         a61=9017/3168; a62=-355/33; a63=46732/5247;
              a64=49/176; a65=-5103/18656
c7=1;         a71=35/384; a72=0; a73=500/1113;
              a74=125/192; a75=-2187/6784; a76=11/84

Ordnung 5: b = [35/384, 0, 500/1113, 125/192,
                -2187/6784, 11/84, 0]
Ordnung 4: b*= [5179/57600, 0, 7571/16695, 393/640,
                -92097/339200, 187/2100, 1/40]

Für jede der sechs Komponenten gilt sc=2e−15+2e−13·max(|yalt|,|yneu|). RMS bedeutet quadratischer Mittelwert; die normierte RMS-Fehlerzahl ist sqrt(sum((Fehler/sc)²)/6). Ein Schritt wird bei Fehlerzahl ≤1 oder |h|≤1e−8 Tag angenommen. Der neue Faktor ist clamp(0.9·Fehlerzahl^(−0.2),0.1,5). clamp begrenzt den Faktor auf das geschlossene Intervall [0,1;5]; bei Fehlerzahl null wird 5 verwendet. Die Schrittweite bleibt zwischen 1e−8 und 0,25 Tag. Eine Integration darf höchstens 1.000.000 Schritte ausführen.

7. Zeitskalen und Beobachter

Die öffentliche Zeit ist UTC. Für die Dynamik wird näherungsweise TDB=UTC+(TAI−UTC+32,184 s)/86400 verwendet. Die implementierten Schaltsekundenstufen lauten:

JD >= 2457754.5 : TAI-UTC = 37 s
JD >= 2457204.5 : TAI-UTC = 36 s
JD >= 2456109.5 : TAI-UTC = 35 s
JD >= 2454832.5 : TAI-UTC = 34 s
JD >= 2453736.5 : TAI-UTC = 33 s
sonst          : TAI-UTC = 32 s

Mit T=(JDUTC−2451545)/36525 wird GMST in Grad berechnet. Die folgenden Polynome und kurzen trigonometrischen Reihen sind vollständig und enthalten alle im Port verwendeten Zahlen:

GMST = 280.46061837 + 360.98564736629 (JD-2451545)
       + 0.000387933 T^2 - T^3/38710000

L_sun  = 280.4665  + 36000.7698 T              [Grad]
L_moon = 218.3165  + 481267.8813 T             [Grad]
Ω_moon = 125.04452 - 1934.136261 T             [Grad]

Δψ = -17.20 sin Ω_moon -1.32 sin(2 L_sun)
     -0.23 sin(2 L_moon) +0.21 sin(2 Ω_moon)   [Bogensekunden]
Δε =  9.20 cos Ω_moon +0.57 cos(2 L_sun)
     +0.10 cos(2 L_moon) -0.09 cos(2 Ω_moon)   [Bogensekunden]

ε0 = 23 + 26/60 + 21.448/3600 Grad
     -(46.8150 T + 0.00059 T^2 - 0.001813 T^3) Bogensekunden

ζ = 2306.2181 T + 0.30188 T^2 + 0.017998 T^3  [Bogensekunden]
z = 2306.2181 T + 1.09468 T^2 + 0.018203 T^3  [Bogensekunden]
θ = 2004.3109 T - 0.42665 T^2 - 0.041833 T^3 [Bogensekunden]

Der geozentrische Standortvektor wird aus MPC-ρ-Faktoren oder geodätischer Länge, Breite und Höhe gebildet, mit scheinbarer Sternzeit rotiert und über Nutation sowie Präzession nach J2000 zurückgeführt.

8. Lichtlaufzeit und scheinbare Richtung

Der Erdzustand und Beobachter werden zur Empfangszeit bestimmt. Für die Emissionszeit wird höchstens fünfmal iteriert:

t_emit,0 = t_TDB
ρk = r_object(t_emit,k) - r_earth(t_TDB) - r_observer(t_UTC)
t_emit,k+1 = t_TDB - |ρk| * 0.0057755183
Abbruch: |t_emit,k+1 - t_emit,k| < 1e-12 Tag

RA  = normalize(atan2(ρy,ρx), 0, 2π)
Dec = atan2(ρz, sqrt(ρx^2+ρy^2))

9. Binärtabellen

Planetenblöcke sind von 2000 bis 2040 in neun überlappende Zeitabschnitte geteilt und verwenden 0,5 Tag Schrittweite. N16 verwendet ebenfalls neun Blöcke und exakt 1,9999817474948436 Tage. Die Überlappung erlaubt sichere Hermite-Interpolation an Abschnittsgrenzen. Alle Werte sind Little-Endian float64.

10. Extrembeispiele

Die nachfolgende Tabelle nennt UTC, eingefrorene JPL-Horizons-Referenz, die C-Port-Rechnung und den sphärisch genäherten Winkelabstand. Station A68 und dieselben nominalen Elemente wurden verwendet. Werte sind Grad.

Bahnform und Periheldistanz der sechs Prüffälle.

Perihelzeit und Orientierung der Bahnebene.

Nichtgravitative Parameter; null bedeutet, dass der jeweilige Term entfällt.

Eingefrorene JPL-Horizons-Referenzpositionen.

Berechnete Positionen und sphärischer Winkelabstand zur Referenz.

433 Eros prüft die gewöhnliche Ellipse. C/2023 A3 liegt fast auf der Parabel und enthält A1/A2. 2P/Encke ist ein periodischer Komet mit Nichtgravitation. 2I/Borisov erzwingt e=3,356475782676596 und A1/A2/A3. Apophis und Duende prüfen die rasche Richtungsänderung bei dichter Erdpassage.

11. Reproduzierbarkeit und Grenzen

• Die 1-Bogensekunden-Grenze gilt gegen dieselben eingefrorenen nominalen Elemente, nicht gegen eine später aktualisierte Bahnlösung.
• Die kurze Nutations-, Präzessions- und TDB-Näherung ist auf den dokumentierten Zeitraum und das angestrebte Genauigkeitsniveau zugeschnitten.
• Schaltsekunden nach der letzten eingebauten Stufe müssen vor Nutzung neuerer Epochen ergänzt werden.
• Die Tabellen decken 2000 bis 2040 ab; außerhalb dieses Bereichs ist ein neuer Datensatz erforderlich.
• Unsicherheitsellipse und Kovarianzfortpflanzung sind nicht Bestandteil dieses reinen Nominalkerns.
• Find_Orb ist kein Laufzeitbestandteil und wird nur privat als zusätzlicher lokaler Vergleich verwendet.

Anhang A. Vergleich der Programmiersprachen

Alle erfolgreichen Fassungen bestanden dieselben 20 Referenzzeilen. Die Genauigkeitstabelle ist alphabetisch nach Sprache sortiert. Die getrennte Geschwindigkeitstabelle ist nach der Medianzeit aufsteigend sortiert, also mit der schnellsten Fassung zuerst. Die Zeiten beziehen sich auf sechs Fälle, die je Sprache siebenmal in sprachübergreifend wechselnder Reihenfolge gerechnet wurden; sie zeigen auf diesem Rechner nur die relative Laufzeit.

A.1 Genauigkeit, alphabetisch

Maximaler sphärischer Winkelabstand über alle 20 Referenzzeilen.

A.2 Geschwindigkeit, schnellste zuerst

Der Faktor 1,000 bezeichnet die schnellste Medianzeit des Prüflaufs.

Anhang B. Lizenzabgrenzung

Der projektgeschriebene Rechenkern und die zugehörigen Generatoren und Hilfen stehen, soweit LICENSE-SCOPE.md sie erfasst, unter Apache License 2.0. Datentabellen, externe Ephemeriden, JPL-Kernel und Drittprogramme sind getrennt zu beurteilen. Eine Weitergabe muss LICENSE und NOTICE beilegen, Copyright- und Lizenzvermerke erhalten und geänderte Dateien kenntlich machen. Patentlizenz und Kündigungsregel ergeben sich aus Abschnitt 3 der Apache-Lizenz.

Projektwebseite: https://www.jostjahn.de/amrum-software/ephemeris/