Was ist das genaue Format der JPL-Ephemeridendateien?

Sie können diese ASCII-JPL-Ephemeridendateien analysieren und verwenden.

Warnung : Dies ist eine Übung für Leute, die sich gerne selbst quälen. Wenn Sie sich nicht für Selbstgeißelung interessieren, sollten Sie tun, was Mark Adler geschrieben hat, und das SPICE-Toolkit verwenden. Ich habe diese Selbstgeißelung vor über einem Jahrzehnt durchgemacht, als SPICE geschlossen wurde und nicht das war, was es jetzt ist. Jetzt denken die Leute daran, tausend Codezeilen oder so zu töten, indem sie meine Arbeit durch SPICE ersetzen. Und genau das sollten sie tun.

Betrachten Sie die Datei ascp2200.405 . (Ich habe diese Datei gewählt, weil sie relativ klein ist.) Die erste Zeile ist

     1  1018

Die "1" bezeichnet die Blocknummer innerhalb der Datei. Die "1018" bezeichnet die Anzahl der Zahlen im Block. Die nächsten 340 Zeilen enthalten diese 1018 Zahlen (plus zwei weitere zum Auffüllen). Die ersten beiden Zahlen stellen die julianischen Start- und Enddaten des Blocks in JPLs eigener Ephemeridenzeit dar. Wenn Ihnen Fehler im Millisekundenbereich egal sind, können Sie TAI anstelle von JPLs verwenden$T_{\text{eph}}$.

Die ersten 168 Zahlen in diesem Block beziehen sich auf Merkur, die nächsten 60 auf die Venus, die nächsten 78 auf das Erd-Mond-Schwerzentrum und so weiter. Dies wird im GROUP 1050Abschnitt der Header-Datei angegeben. Sie können die Dimensionalität jedes Elements aus diesem Abschnitt ableiten. Alternativ sollten Sie wissen, dass die Position dreidimensional ist, während Librationen zweidimensional sind.

Die 168 Zahlen, die die Position des Merkur beschreiben, sind in vier Unterblöcken organisiert, die jeweils 3*14=42 Zahlen umfassen. (In Bezug auf diese magischen Zahlen: 14 stammt aus den GRUPPE 1050-Daten für Merkur, und 3 ist die Dimensionalität der Positionsdaten.) Von diesen 42 Zahlen sind die ersten 14 die Koeffizienten für die Tschebyscheff-Polynomkoeffizienten 13. ^Ordnung , die das x beschreiben Position, die nächsten 14 beschreiben die y-Position und die letzten 14 beschreiben die z-Position.

Diese 14 Zahlen sind Tschebyscheff-Polynomkoeffizienten mit Zeit (korrekt skaliert) als unabhängige Variable. Tschebyscheff-Polynome sind Funktionen, die auf den Bereich [-1,1] beschränkt sind. Um irgendwann die Position des Merkur zu finden$t$, müssen Sie zuerst den Block finden, der sich auf diese Zeit bezieht. Dann müssen Sie den Unterblock (Mercury hat vier Sätze von Koeffizienten) finden, der sich auf diese Zeit bezieht. Als nächstes müssen Sie die Zeit von der Mitte dieses Unterblocks verschieben und das Ergebnis so skalieren, dass -1 den Beginn des Unterblocks und 1 das Ende darstellt. Schließlich müssen Sie berechnen$\sum_{k=0}^{13} C_n T_n(t)$für jedes der drei Elemente des Positionsvektors, wobei$C_n$sind die Tschebyscheff-Polynomkoeffizienten und die$T_n(t)$sind die Tschebyscheff-Polynome:$T_0(t)=1$,$T_1(t)=t$, und$T_{n+1}(t) = 2tT_n(t)-T_{n-1}(t)$. Die Ergebnisse sind in Kilometern angegeben, mit Ausnahme von Elementen, die sich auf die Rotation beziehen. Dies kann oder kann nicht das sein, was Sie wollen.

Um die Sache noch lustiger zu machen (also Selbstgeißelung-Spaß), was ist, wenn Sie die Position des Mars relativ zur Erde wissen möchten? Die JPL-Ephemeridendatei enthält keine Koeffizienten für die Erde. Es hat Koeffizienten für den Schwerpunkt Erde-Mond (relativ zum Schwerpunkt des Sonnensystems) und für den Mond (relativ zum Erdmittelpunkt). Sie müssen die Position der Erde berechnen, indem Sie die Position des Mars, die Position des Erde-Mond-Schwerpunkts und die Position des Mondes relativ zur Erde berechnen. Als nächstes müssen Sie die Position der Erde (relativ zum Baryzentrum des Sonnensystems) mithilfe des Erd-Mond-Massenverhältnisses berechnen (eine magische Zahl, die sich in der Header-Datei befindet). Damit können Sie endlich die Position des Mars relativ zur Erde berechnen. Das ist die momentane relative Position.

Alternativ können Sie dem Rat von Mark Adler folgen: Verwenden Sie das SPICE-Toolkit.

Nein, der erste Schritt besteht darin, das SPICE Toolkit herunterzuladen . Es wird die Ephemeriden lesen und die Daten generieren, die Sie für Sie benötigen, und noch viel, viel mehr tun . Es gibt viele Formate, die in SPICE-Kernel-Dateien verwendet werden, nicht nur dasjenige, das Sie gerade betrachten, das wahrscheinlich ein Nur-Chebyshev-Positions-Array (Typ 2) ist. Es gibt 16 Typen. Die Werte können Little oder Big Endian gespeichert werden. Die Datei kann große Kommentarabschnitte enthalten. Jede Datei enthält mehrere Körper. Die Körper haben zeitlich unterschiedliche Auflösungsgrade. Es wird einige Arbeit für Sie erfordern, das Ding selbst zu entschlüsseln, was nur verschwendete Zeit ist, da das Toolkit all das für Sie erledigt.

Laden Sie das Toolkit herunter.

Um auf Ihre Frage zu antworten, wenn Sie erfolgreich einen einzelnen Chebyshev-Array-Datensatz vom Typ 2 gefunden haben und zu den Polynomkoeffizienten gekommen sind (nachdem Sie MID und RADIUS für diesen Datensatz gelesen haben), lautet Folgendes: Die Sequenz der Doppel beginnt mit der Konstante Term des Polynoms für$x$, und fährt fort mit dem linearen Koeffizienten von$t$, bis dieses Polynom vollständig ist, gefolgt von der gleichen Sache für$y$und$z$. Die Doubles sind IEEE-64-Bit-Gleitkommazahlen in Little- oder Big-Endian-Reihenfolge, je nachdem, was im Header steht.

$t$ist die Ephemeridenzeit in Sekunden und die Ergebnisse in km.

Das Format ist hier dokumentiert .

Sie sollten jedoch das Toolkit herunterladen.

Das genaue Format der SPICE Toolkit SPK-Dateien (auch BSP-Dateien genannt) heißt DAF ( Double-precision Array File ) . Es ist ein Binärformat, das aus aufeinanderfolgenden Blöcken von 1.024 Bytes besteht. Der erste Datensatz ist immer der Dateidatensatz , dann hat er einige Blöcke, die als Kommentarbereich dienen , dann hat er die Array-Datensätze und schließlich eine Sammlung von Datensätzen .

[FILE RECORD] (1024 bytes) [FIRST COMMENT BLOCK] (1024 bytes) ... [LAST COMMENT BLOCK] (1024 bytes) [FIRST BLOCK OF SUMMARY RECORDS] (1024 bytes) ... [LAST BLOCK OF SUMMARY RECORDS] (1024 bytes) [FIRST BLOCK OF NAME RECORDS] (1024 bytes) ... [LAST BLOCK OF NAME RECORDS] (1024 bytes) [FIRST BLOCK OF ELEMENT RECORDS] (1024 bytes) ... [LAST BLOCK OF ELEMENT RECORDS] (1024 bytes)

Alle binären SPICE-Daten werden im DAF-Format gespeichert (nicht nur Ephemeriden). Die spezifische Art von Daten (d. h. die Semantik der Datensätze) hängt von den Informationen ab, die in den Übersichtsdatensätzen und den Elementdatensätzen gespeichert sind. Sie müssen wissen, welche Art von Daten Sie lesen, um die Aufzeichnungen richtig interpretieren zu können.

Es gibt viele Arten von Ephemeriden . Planeten sind oft Typ II ; Satelliten sind oft vom Typ III ; Raumfahrzeuge sind oft Typ I. Die Typen II und III stellen Chebyshev-Polynome dar, und ihr Unterschied besteht darin, ob sie die Ableitung des Chebyshev-Polynoms zur Berechnung der Geschwindigkeit verwenden (Typ II) oder ob ein separates Polynom zur Berechnung der Geschwindigkeit verwendet wird (Typ III).

Mit anderen Worten: Eine SPK-Datei vom Typ II liefert Ihnen ein Polynom bestimmter Ordnung, dessen Wert die Position und dessen Ableitung die Geschwindigkeit darstellt. Eine SPK-Datei vom Typ III liefert Ihnen ein Polynom für die Position und ein anderes Polynom für die Geschwindigkeit. Beide Arten sind Tschebyscheff-Polynome.

Bei SPK-Dateien erfahren Sie in den zusammenfassenden Datensätzen, wo jeder Datensatz beginnt und endet, den Ephemeridentyp, die Zeitspanne, den Frame, das Ziel und den Beobachter. Jeder spezifische Datensatz (der höchstwahrscheinlich mehr als einen Block enthalten wird) sieht wie folgt aus:

+---------------+ | Record 1 | +---------------+ | Record 2 | +---------------+ . . . +---------------+ | Record N | +---------------+ | INIT | +---------------+ | INTLEN | +---------------+ | RSIZE | +---------------+ | N | +---------------+

und jeder Datensatz aussieht

+------------------+ | MID | +------------------+ | RADIUS | +------------------+ | X coefficients | +------------------+ | Y coefficients | +------------------+ | Z coefficients | +------------------+

für Typ II und dergleichen

+------------------+ | MID | +------------------+ | RADIUS | +------------------+ | X coefficients | +------------------+ | Y coefficients | +------------------+ | Z coefficients | +------------------+ | X' coefficients | +------------------+ | Y' coefficients | +------------------+ | Z' coefficients | +------------------+

für Typ III.

Der Grad des Polynoms ist ( RSIZE - 2 ) / 3 - 1. Der MIDgibt Ihnen den Intervallmittelpunkt (Ephemeridenzeit) und RADIUSdie halbe Größe des Intervalls an. Seine Zeitspanne ist also [MID - RADIUS, MID + RADIUS].

Zur Erlangung der Position wertet man einfach Chebyshev zum entsprechenden Zeitpunkt aus. Für die Geschwindigkeit werten Sie entweder die Ableitung oder ein anderes Polynom aus. Denken Sie daran, die Zeit vor der Bewertung zu normalisieren und die Zeit nach der Bewertung zu denormalisieren (das hat mich oft erwischt).

Der direkte Umgang mit DAF kann eine lohnende Erfahrung sein (z. B.: Ich habe einen Thread-sicheren Hochleistungs-Ephemeriden-Evaluator), aber es ist leicht, Fehler zu machen. Und nach viel Arbeit haben Sie nur eine Position oder Geschwindigkeit in einem bestimmten Frame, die Sie noch drehen/übersetzen müssen usw. (für die Sie entweder SPICE neu erfinden oder sich darauf verlassen müssen).

Wenn Sie keinen triftigen Grund dagegen haben, bleiben Sie einfach bei SPICE.

Ein einfaches Python-Skript

Ich habe gerade ein Python-Skript erstellt, das jemandem beim Einstieg in das DAF-Binärformat helfen könnte. Es wird einfach die Dateiaufzeichnungsinformationen und den Inhalt (Array-Zusammenfassungen) ausgeben, die in einem DAF enthalten sind.

import mmap
import struct

RECLEN = 1024

def peek_spk(mem):
    # file record is always the first one

    # string data
    locidw = mem[0:8]
    locifn = mem[16:16 + 60]
    locfmt = mem[88:88+8]
    ftpstr = mem[699:699+28]

    # endianness
    fmt     = "<" if locfmt == "LTL-IEEE" else ">"
    int_fmt = fmt + "I"

    nd,    = struct.unpack_from(int_fmt, mem, offset=8)
    ni,    = struct.unpack_from(int_fmt, mem, offset=12)
    fward, = struct.unpack_from(int_fmt, mem, offset=76)
    bward, = struct.unpack_from(int_fmt, mem, offset=80)
    free,  = struct.unpack_from(int_fmt, mem, offset=84)

    print "locidw {0}".format(locidw)
    print "nd     {0}".format(nd)
    print "ni     {0}".format(ni)
    print "locifn {0}".format(locifn)
    print "fward  {0}".format(fward)
    print "bward  {0}".format(bward)
    print "free   {0}".format(free)
    print "locfmt {0}".format(locfmt)
    print "ftpstr {0}".format(repr(ftpstr))

    # let's do the first summary record only... we need to read nd
    # doubles and ni integers starting at offset fward

    offset = (fward - 1) * RECLEN
    sum_fmt = fmt + nd * "d" + ni * "I"
    size = nd + (ni + 1) / 2 # integer division
    while True:
        nxt, prv, nsum = struct.unpack_from(fmt + "ddd", mem, offset=offset)
        offset += 24 # skip three doubles
        for n in range(int(nsum)):
            print struct.unpack_from(sum_fmt, mem, offset = offset)
            offset += size * 8
        if nxt == 0:
            break


def peek(path):

    print "peeking into {0}".format(path)

    with open(path, "rb") as fp:
        mem = mmap.mmap(fp.fileno(), 0, access=mmap.PROT_READ)

    peek_spk(mem)


if __name__ == "__main__":
    import sys
    for path in sys.argv[1:]:
        peek(path)

Zum Beispiel: läuft es für jup310.bsp gibt

peeking into /data/spice/jup310.bsp locidw DAF/SPK nd 2 ni 6 locifn NIO2SPK
fward 6 bward 6 free 122077514 locfmt LTL-IEEE ftpstr 'FTPSTR:\r:\n:\r\n:\r\x00:\x81:\x10\xce:ENDFTP' (-3155716758.8160305, 3155716867.183885, 501, 5, 1, 3, 897, 7208500) (-3155716758.8160305, 3155716867.183885, 502, 5, 1, 3, 7208501, 14367404) (-3155716758.8160305, 3155716867.183885, 503, 5, 1, 3, 14367405, 19140106) (-3155716758.8160305, 3155716867.183885, 504, 5, 1, 3, 19140107, 22451778) (-3155716758.8160305, 3155716867.183885, 505, 5, 1, 3, 22451779, 44074360) (-3155716758.8160305, 3155716867.183885, 514, 5, 1, 3, 44074361, 65696942) (-3155716758.8160305, 3155716867.183885, 515, 5, 1, 3, 65696943, 90825888) (-3155716758.8160305, 3155716867.183885, 516, 5, 1, 3, 90825889, 115954834) (-3155716758.8160305, 3155716867.183885, 599, 5, 1, 3, 115954835, 120922238) (-3155716758.8160305, 3155716867.183885, 3, 0, 1, 2, 120922239, 121109489) (-3155716758.8160305, 3155716867.183885, 5, 0, 1, 2, 121109490, 121168877) (-3155716758.8160305, 3155716867.183885, 10, 0, 1, 2, 121168878, 121328726) (-3155716758.8160305, 3155716867.183885, 399, 3, 1, 2, 121328727, 122077513)

Das bedeutet, dass jup310 Aufzeichnungen über ET (-3155716758,8160305, 3155716867,183885) für einige Jupitermonde (501 bis 519) relativ zu Jupiters Baryzentrum (5) im J2000-Format (1) enthält und vom Typ 3 ist. Andere Aufzeichnungen sind für die Baryzentren von Erde, Jupiter und Sonne (3, 5 und 10) relativ zum Baryzentrum des Sonnensystems (0) in J200 (1) und sind Typ 2.

Jeder Datensatz sagt Ihnen, wo in der Datei jedes Array beginnt und endet. Sie können das Skript erweitern, um zu jeder Position zu springen und den gesamten Datensatz für eine spätere Interpolation zu lesen.

Verfügbare Software

Ich habe schnell eine einfache Implementierung in Python erstellt, die die Typen II und III lesen und direkt damit arbeiten kann. Es lässt sich gut mit CSPICE vergleichen.

https://gist.github.com/arrieta/c2b56f1e2277a6fede6d1afbc85095fb

Ich werde allen vorherigen Antworten widersprechen. Ich würde argumentieren, dass die Implementierung Ihres eigenen Ephemeridenprogramms grenzwertig trivial ist, sobald Sie das Format der Dateien kennen und die relativ einfachen Gleichungen für die Verwendung der Polynome zur Berechnung der Positionen haben. Die Aufgabe liegt durchaus im Bereich eines Programmieranfängers. Das Format der Datei ist der schwierigste Teil, nur weil es viele Teile zusammenzusetzen gibt, aber nicht gerade schwierig zu machen.

Der Artikel Format of the JPL Ephemeris Files enthält eine ziemlich gründliche Anleitung mit Beispielen und einer JavaScript-Implementierung als Beispiel. Beachten Sie, dass die gesamte Klasse, die die Koeffizientensuchen und Berechnungen durchführt, weniger als 100 Codezeilen umfasst. Zusätzlich gibt es Beispielimplementierungen im Github-Repository gmiller123456/jpl-development-ephemeris , das derzeit Implementierungen in Python, Perl, C# und JavaScript enthält.

Wie Sie bereits herausgefunden haben, enthält der Abschnitt "GROUP 1050" des Headers viele Schlüssel zum Format der Datei. Dies ist der Abschnitt für DE405

GROUP   1050

 3   171   231   309   342   366   387   405   423   441   753   819   899
14    10    13    11     8     7     6     6     6    13    11    10    10
 4     2     2     1     1     1     1     1     1     8     2     4     4

Die Reihenfolge und Anzahl der Komponenten, denen diese entsprechen, ist im Dokument im Ascii-Format von JPL angegeben. Jede Spalte stellt die Werte für einen bestimmten Planeten (oder eine andere Eigenschaft) dar:

#    Properties Units        Center   Name
1    3          km           SSB      Mercury
2    3          km           SSB      Venus
3    3          km           SSB      Earth-Moon barycenter
4    3          km           SSB      Mars
5    3          km           SSB      Jupiter
6    3          km           SSB      Saturn
7    3          km           SSB      Uranus
8    3          km           SSB      Neptune
9    3          km           SSB      Pluto
10   3          km           Earth    Moon (geocentric)
11   3          km           SSB      Sun
12   2          radians               Earth Nutations in longitude and obliquity (IAU 1980 model)
13   3          radians               Lunar mantle libration
14   3          radians/day           Lunar mantle angular velocity
15   1          Seconds               TT-TDB (at geocenter)

Jede der ASCII-Dateien ist in Blöcke unterteilt, die für die in der Header-Datei angegebene Anzahl von Tagen gültig sind. Die ersten beiden Zahlen im Block sind die julianischen Daten, für die dieser Block gültig ist, daher ist es trivial, den entsprechenden Block zu einem bestimmten julianischen Datum zu finden.

Ab GRUPPE 1050 in der Kopfzeile entsprechen die Zahlen in jeder Spalte:

1. The start offset in the block (starting at 1)
2. Number of coefficients for each property
3. The number of subintervals

Um also den Offset in den Block für einen bestimmten Planeten zu berechnen:

lengthOfSubinterval = daysPerBlock / numberOfSubintervals
subinterval = floor( (JD-blockStartDate)/lengthOfSubinterval )
offset=subinterval*numberOfCoefficients*numberOfProperties+seriesStartOffset;

Für Planeten haben Sie also alle X-Koeffizienten, dann alle Y-Koeffizienten und schließlich alle Z-Koeffizienten. Beispielsweise sind die Koeffizienten für Merkur für JD=2458850,5 im folgenden Code angegeben.

Um das Beispiel aus dem Artikel Format des JPL-Ephemeris-Dateien auszuleihen , werden die Positionen tatsächlich berechnet:

function computePolynomial(x,coefficients){
  let T=new Array();

  T[0]=1;
  T[1]=x;
  for(let n=2;n<14;n++)  {
    T[n]=2*x*T[n-1] - T[n-2];
  }

  let v=0;
  for(let i=coefficients.length-1;i>=0;i--){
    v+=T[i]*coefficients[i];
  }
  return v;
}    

function computeExamplePolynomials(){   
    let X=[0.230446411715880504E+04,  0.133726736662702635E+08, -0.782187090879053358E+04, -0.267678745522568279E+05,  
           -0.227070698075548364E+03, -0.142012340261296774E+02, -0.924872006275108544E-01,  0.431659104815666252E-02,
            0.356917634561652571E-03,  0.302564651657819373E-04, 0.980701702776103911E-06,  0.505819702568259545E-07,
            0.113034198242195379E-08,  0.323800745882515925E-10];

    let Y=[-0.593914454531169161E+08,  0.138391312173493067E+07, 0.725419090211108793E+06,  0.139471465250126903E+04,
           -0.290917422263861397E+03, -0.635064566332839320E+01, -0.646844700926034299E+00, -0.120797394835047579E-01, 
           -0.681164244772722110E-03, -0.783160742259704191E-05, -0.953933699143903451E-07,  0.170514411319974421E-07,
            0.132846579503924915E-08,  0.625629348278007546E-10];

    let Z=[-0.318846685234071501E+08, -0.647159726192409638E+06,  0.388325111500594590E+06,  0.351975238047553557E+04,
           -0.131868705094903135E+03, -0.192040059987689204E+01, -0.335952534459033059E+00, -0.690035617434751804E-02,
           -0.400870301836372738E-03, -0.731991272299537233E-05, -0.152615685994738755E-06,  0.386553675297770635E-08,  
            0.592487943320094233E-09,  0.300642066655273442E-10];

    let x=-0.5;
    console.log(computePolynomial(x,X));
    console.log(computePolynomial(x,Y));
    console.log(computePolynomial(x,Z));
}

Der Aufruf der obigen Funktion computeExamplePolynomials() erzeugt die folgenden Werte:

X = -6706768.766943997 km
Y = -60444568.85087551 km
Z = -31751664.901437085 km