¿Cuál es el formato exacto de los archivos de efemérides JPL?

Puede analizar y utilizar esos archivos de efemérides JPL ascii.

Advertencia : este es un ejercicio para personas a las que les gusta torturarse a sí mismas. Si no te gusta la autoflagelación, deberías hacer lo que escribió Mark Adler y usar el kit de herramientas SPICE. Pasé por esta autoflagelación hace más de una década, cuando SPICE estaba cerrado y no era lo que es ahora. Ahora la gente está pensando en acabar con mil líneas de código reemplazando mi trabajo con SPICE. Y eso es exactamente lo que deberían hacer.

Considere el archivo ascp2200.405 . (Elegí ese archivo porque es relativamente pequeño). La primera línea es

     1  1018

El "1" denota el número de bloque dentro del archivo. El "1018" denota el número de números en el bloque. Las siguientes 340 líneas contienen esos 1018 números (más dos más para relleno). Los dos primeros números representan las fechas julianas de inicio y finalización del bloque, en el propio tiempo de efemérides del JPL. Si no le importan los errores de milisegundos, puede usar TAI en lugar de JPL$T_{\text{eph}}$.

Los primeros 168 números de ese bloque pertenecen a Mercurio, los siguientes 60 a Venus, los siguientes 78 al baricentro Tierra-Luna, y así sucesivamente. Esto se especifica en la GROUP 1050sección del archivo de cabecera. Puede deducir la dimensionalidad de cada elemento de esa sección. Alternativamente, debe saber que la posición es tridimensional, mientras que las libraciones son bidimensionales.

Los 168 números que describen la posición de Mercurio están organizados en cuatro sub-bloques, cada uno compuesto por 3*14=42 números. (Con respecto a esos números mágicos: 14 es de los datos del GRUPO 1050 para Mercurio, y 3 es la dimensionalidad de los datos de posición). De esos 42 números, los primeros 14 son los coeficientes para los coeficientes polinómicos de Chebyshev de orden 13 que describen la ^x posición, los siguientes 14 describen la posición y, y los últimos 14 describen la posición z.

Esos 14 números son coeficientes polinómicos de Chebyshev con el tiempo (apropiadamente escalado) como variable independiente. Los polinomios de Chebyshev son funciones limitadas al dominio [-1,1]. Para encontrar la posición de Mercurio en algún momento$t$, primero deberá encontrar el bloque que pertenece a ese momento. Luego, deberá encontrar el subbloque (Mercurio tiene cuatro conjuntos de coeficientes) que pertenece a ese momento. A continuación, deberá compensar el tiempo desde la mitad de ese subbloque y escalar el resultado para que -1 represente el inicio del subbloque y 1 represente el final. Finalmente, deberás calcular$\sum_{k=0}^{13} C_n T_n(t)$para cada uno de los tres elementos del vector de posición, donde$C_n$son los coeficientes polinómicos de Chebyshev y el$T_n(t)$son los polinomios de Chebyshev:$T_0(t)=1$,$T_1(t)=t$, y$T_{n+1}(t) = 2tT_n(t)-T_{n-1}(t)$. Los resultados están en kilómetros, a excepción de los elementos relacionados con la rotación. Esto puede o no ser lo que quieres.

Para hacer las cosas aún más divertidas (la autoflagelación es divertida), ¿qué pasa si quieres saber la posición de Marte en relación con la Tierra? El archivo de efemérides del JPL no tiene coeficientes para la Tierra. Tiene coeficientes para el baricentro Tierra-Luna (relativo al baricentro del sistema solar) y para la Luna (relativo al centro de la Tierra). Deberá calcular la posición de la Tierra calculando la posición de Marte, la posición del baricentro Tierra-Luna y la posición de la Luna en relación con la Tierra. A continuación, deberá calcular la posición de la Tierra (en relación con el baricentro del sistema solar) utilizando la relación de masas Tierra-Luna (un número mágico que se encuentra en el archivo de encabezado). Con esto, finalmente puedes calcular la posición de Marte en relación con la Tierra. Esa es la posición relativa instantánea.

Alternativamente, puede seguir el consejo de Mark Adler: use el kit de herramientas SPICE.

No, el primer paso es descargar SPICE Toolkit . Leerá las efemérides y generará los datos que necesita para usted, además de hacer mucho, mucho más . Se utilizan muchos formatos en los archivos del núcleo SPICE, no solo el que está viendo en particular, que probablemente sea una matriz de solo posición de Chebyshev (tipo 2). Hay 16 tipos. Los valores se pueden almacenar en little o big endian. Puede haber grandes secciones de comentarios en el archivo. Hay varios cuerpos en cada archivo. Los cuerpos tendrán diferentes grados de resolución en el tiempo. Le costará algo de trabajo decodificarlo usted mismo, lo cual es una pérdida de tiempo, ya que el conjunto de herramientas hará todo eso por usted.

Descarga el kit de herramientas.

Para responder a su pregunta, si ha encontrado con éxito un solo registro de matriz de solo posición de Chebyshev de tipo 2 y ha llegado a los coeficientes polinómicos (después de haber leído MID y RADIUS para ese registro), es que: la secuencia de dobles comienza con la constante término del polinomio para$x$, y continúa con el coeficiente lineal de$t$, hasta que ese polinomio esté completo, seguido de lo mismo para$y$y$z$. Los dobles son números de coma flotante de 64 bits IEEE en orden pequeño o grande, según lo que diga en el encabezado.

$t$es el tiempo de las efemérides en segundos, y los resultados están en km.

El formato está documentado aquí .

Sin embargo, debe descargar el kit de herramientas.

El formato exacto de los archivos SPK de SPICE Toolkit (también llamados archivos BSP) se denomina DAF ( archivo de matriz de doble precisión ) . Es un formato binario formado por bloques consecutivos de 1.024 bytes. El primer registro es siempre el registro del archivo , luego tiene algunos bloques que sirven como área de comentarios , luego tiene los registros de matriz y, finalmente, una colección de registros de datos .

[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)

Todos los datos SPICE binarios se almacenan en formato DAF (no solo las efemérides). El tipo específico de datos (es decir, la semántica de los registros) depende de la información almacenada en los registros de resumen y los registros de elementos. Necesita saber qué tipo de datos está leyendo para interpretar correctamente los registros.

Hay muchos tipos de efemérides . Los planetas son a menudo del Tipo II ; los satélites son a menudo del tipo III ; las naves espaciales son a menudo Tipo I. Los tipos II y III representan polinomios de Chebyshev, y su diferencia es si usan la derivada del polinomio de Chebyshev para calcular la velocidad (Tipo II) o si se usa un polinomio separado para calcular la velocidad (Tipo III).

En otras palabras: un archivo SPK Tipo II le proporcionará un polinomio de cierto orden cuyo valor representa la posición y cuya derivada representa la velocidad. Un archivo SPK Tipo III le proporcionará un polinomio para la posición y un polinomio diferente para la velocidad. Ambos tipos son polinomios de Chebyshev.

Para los archivos SPK, los registros de resumen le indicarán dónde comienza y termina cada registro, el tipo de efemérides, el intervalo de tiempo, el marco, el objetivo y el observador. Cada registro específico (que muy probablemente contendrá más de un bloque) tiene el siguiente aspecto:

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

y cada registro parece

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

para Tipo II y similares

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

para el Tipo III.

El grado del polinomio es ( RSIZE - 2 ) / 3 - 1. El MIDle indica el punto medio del intervalo (tiempo de efemérides) y RADIUSle indica la mitad del tamaño del intervalo. Entonces su lapso de tiempo es [MID - RADIUS, MID + RADIUS].

Para obtener la posición, simplemente evalúe el Chebyshev en el momento correspondiente. Para la velocidad, evalúas la derivada o un polinomio diferente. Recuerde normalizar el tiempo antes de evaluar y desnormalizar el tiempo después de evaluar (esto me ha afectado muchas veces).

Tratar directamente con DAF puede ser una experiencia gratificante (por ejemplo: tengo un evaluador de efemérides de alto rendimiento seguro para subprocesos), pero es fácil equivocarse. Y después de mucho trabajo, solo tiene una posición o velocidad en un marco determinado que aún necesita rotar/traducir/etc. (para lo cual necesita reinventar SPICE o confiar en él).

A menos que tenga una muy buena razón para no hacerlo, simplemente apéguese a SPICE.

Un simple script de Python

Acabo de hacer un script de Python que podría ayudar a alguien a comenzar con el formato binario DAF. Simplemente volcará la información de registro del archivo y los contenidos (resúmenes de matriz) contenidos en un DAF.

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)

Por ejemplo: ejecutarlo para jup310.bsp da

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)

Eso significa que jup310 contiene registros que abarcan ET (-3155716758.8160305, 3155716867.183885) para algunas lunas de Júpiter (501 a 519) en relación con el Baricentro de Júpiter (5), en formato J2000 (1), y son de tipo 3. Otros registros son para los Baricentros de la Tierra, Júpiter y el Sol (3, 5 y 10) en relación con el Baricentro del Sistema Solar (0) en J200 (1), y son de tipo 2.

Cada registro le dice en qué parte del archivo comienza y termina cada matriz. Puede extender la secuencia de comandos para saltar a cada ubicación y leer el registro completo para una interpolación posterior.

Software disponible

Rápidamente creé una implementación simple en Python que puede leer y trabajar directamente con los Tipos II y III. Se compara bien con CSPICE.

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

No estoy de acuerdo con todas las respuestas anteriores. Yo diría que implementar su propio programa Ephemeris es casi trivial una vez que conoce el formato de los archivos y tiene las ecuaciones relativamente simples para usar los polinomios para calcular las posiciones. La tarea está dentro del ámbito de un programador principiante. El formato del archivo es la parte más difícil, solo porque hay muchas piezas para armar, pero no exactamente difíciles de hacer.

El artículo Format of the JPL Ephemeris Files tiene un recorrido bastante completo, con ejemplos y una implementación de JavaScript como ejemplo. Tenga en cuenta que toda la clase que realiza las búsquedas de coeficientes y los cálculos tiene menos de 100 líneas de código. Además, hay implementaciones de ejemplo en el repositorio de Github gmiller123456/jpl-development-ephemeris que actualmente tiene implementaciones en Python, Perl, C# y JavaScript.

Como ya se dio cuenta, la sección "GRUPO 1050" del encabezado contiene muchas de las claves para el formato del archivo. Esta es la sección para 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

El orden y el número de componentes a los que corresponden se dan en el documento en formato Ascii de JPL. Cada columna representa los valores para un planeta dado (u otra propiedad):

#    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)

Cada uno de los archivos ASCII se divide en bloques, válidos para el número de días indicado en el archivo de cabecera. Los dos primeros números del bloque son las fechas julianas para las que ese bloque es válido, por lo que encontrar el bloque correspondiente a una fecha juliana dada es trivial.

Del GRUPO 1050 de la cabecera, los números de cada columna corresponden a:

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

Entonces, para calcular el desplazamiento en el bloque para un planeta dado:

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

Entonces, para los planetas, tendrá todos los coeficientes X, luego todos los coeficientes Y y finalmente todos los coeficientes Z. Por ejemplo, los coeficientes de Mercurio para JD=2458850,5 se dan en el siguiente código.

Tomando prestado el ejemplo del artículo Format of the JPL Ephemeris Files , esto en realidad calcula las posiciones:

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));
}

Llamar a la función computeExamplePolynomials() anterior, produce los valores a continuación:

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