Каков точный формат файлов эфемерид JPL?

Вы можете анализировать и использовать эти эфемеридные файлы ascii JPL.

Предупреждение : это упражнение для тех, кто любит мучить себя. Если вы не занимаетесь самобичеванием, вам следует сделать то, что написал Марк Адлер, и использовать инструментарий SPICE. Я прошел через это самобичевание более десяти лет назад, когда SPICE была закрыта и не была тем, чем она является сейчас. Теперь люди думают убить тысячу строк кода или около того, заменив мою работу на SPICE. И это именно то, что они должны делать.

Рассмотрим файл ascp2200.405 . (Я выбрал этот файл, потому что он относительно небольшой.) Первая строка

     1  1018

«1» обозначает номер блока в файле. «1018» обозначает количество чисел в блоке. Следующие 340 строк содержат эти 1018 чисел (плюс еще два для заполнения). Первые два числа представляют начальную и конечную юлианские даты блока в собственном эфемеридном времени JPL. Если вас не волнуют миллисекундные ошибки, вы можете использовать TAI вместо JPL.$T_{\text{eph}}$.

Первые 168 чисел в этом блоке относятся к Меркурию, следующие 60 — к Венере, следующие 78 — к барицентру Земля-Луна и так далее. Это указывается в GROUP 1050разделе заголовочного файла. Вы можете вывести размерность каждого элемента из этого раздела. В качестве альтернативы вы должны знать, что позиция трехмерна, а либрации двумерны.

168 чисел, описывающих положение Меркурия, организованы в четыре подблока, каждый из которых состоит из 3*14=42 чисел. (Относительно этих магических чисел: 14 — из данных GROUP 1050 для Меркурия, а 3 — размерность данных о положении.) Из этих 42 чисел первые 14 — это коэффициенты полинома Чебышева 13 ^-го порядка, описывающие x положение, следующие 14 описывают положение y, а последние 14 описывают положение z.

Эти 14 чисел являются полиномиальными коэффициентами Чебышева со временем (правильно масштабированным) в качестве независимой переменной. Многочлены Чебышева — это функции, ограниченные областью определения [-1,1]. Найти положение Меркурия в какой-то момент$t$, вам сначала нужно найти блок, относящийся к этому времени. Затем вам нужно будет найти подблок (у Меркурия есть четыре набора коэффициентов), относящийся к этому времени. Затем вам нужно будет сместить время от середины этого подблока и масштабировать результат так, чтобы -1 представляло начало подблока, а 1 представляло конец. Наконец, вам нужно вычислить$\sum_{k=0}^{13} C_n T_n(t)$для каждого из трех элементов вектора положения, где$C_n$– коэффициенты полинома Чебышева и$T_n(t)$полиномы Чебышева:$T_0(t)=1$,$T_1(t)=t$, а также$T_{n+1}(t) = 2tT_n(t)-T_{n-1}(t)$. Результаты указаны в километрах, за исключением пунктов, относящихся к вращению. Это может быть или не быть тем, что вы хотите.

Чтобы сделать ситуацию еще более забавной (самобичевание, то есть), что, если вы хотите узнать положение Марса относительно Земли? В файле эфемерид JPL нет коэффициентов для Земли. Он имеет коэффициенты для барицентра Земля-Луна (относительно барицентра Солнечной системы) и для Луны (относительно центра Земли). Вам нужно будет рассчитать положение Земли, рассчитав положение Марса, положение барицентра Земля-Луна и положение Луны относительно Земли. Далее вам нужно рассчитать положение Земли (относительно барицентра Солнечной системы), используя отношение масс Земли и Луны (магическое число, которое находится в заголовочном файле). Благодаря этому вы, наконец, можете вычислить положение Марса относительно Земли. Это мгновенное относительное положение.

В качестве альтернативы вы можете последовать совету Марка Адлера: используйте инструментарий SPICE.

Нет, первым шагом является загрузка SPICE Toolkit . Он будет считывать эфемериды и генерировать нужные вам данные, а также делать многое, многое другое . В файлах ядра SPICE используется много форматов, а не только тот, на который вы смотрите, который, вероятно, представляет собой массив Чебышева только с позициями (тип 2). Есть 16 видов. Значения могут храниться с прямым или обратным порядком байтов. В файле могут быть большие разделы комментариев. В каждом файле несколько тел. Тела будут иметь разную степень разрешения во времени. Вам потребуется некоторая работа, чтобы расшифровать это самостоятельно, а это пустая трата времени, поскольку инструментарий сделает все это за вас.

Загрузите набор инструментов.

Чтобы ответить на ваш вопрос, если вы успешно нашли одну запись массива Чебышева типа 2, состоящую только из позиций, и пришли к полиномиальным коэффициентам (прочитав MID и RADIUS для этой записи), это: последовательность двойников начинается с константы член многочлена для$x$, и продолжается линейным коэффициентом$t$, пока этот полином не будет завершен, за которым следует то же самое для$y$а также$z$. Двойники — это 64-битные числа IEEE с плавающей запятой в порядке младшего или старшего порядка, в зависимости от того, что написано в заголовке.

$t$это эфемеридное время в секундах, а результаты в км.

Формат задокументирован здесь .

Однако вы должны скачать инструментарий.

Точный формат SPK-файлов SPICE Toolkit (также называемых BSP-файлами) называется DAF ( файл массива двойной точности ) . Это двоичный формат, состоящий из последовательных блоков по 1024 байта. Первой записью всегда является запись файла , затем несколько блоков, служащих областью комментариев , затем записи массива и, наконец, набор записей данных .

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

Все бинарные данные SPICE хранятся в формате DAF (не только эфемериды). Конкретный вид данных (то есть семантика записей) зависит от информации, хранящейся в сводных записях и записях элементов. Вам нужно знать, какие данные вы читаете, чтобы правильно интерпретировать записи.

Существует множество типов эфемерид . Планеты часто относятся к типу II ; спутники часто относятся к типу III ; космические корабли часто относятся к типу I. Типы II и III представляют полиномы Чебышева, и их различие заключается в том, используют ли они производную полинома Чебышева для расчета скорости (тип II) или используют ли отдельный полином для расчета скорости (тип III).

Другими словами: файл SPK типа II предоставит вам один полином определенного порядка, значение которого представляет позицию, а производная — скорость. Файл SPK типа III предоставит вам один полином для положения и другой полином для скорости. Оба типа являются полиномами Чебышева.

Для файлов SPK сводные записи сообщают вам, где начинается и заканчивается каждая запись, тип эфемерид, временной интервал, кадр, цель и наблюдатель. Каждая конкретная запись (которая, скорее всего, будет содержать более одного блока) выглядит следующим образом:

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

и каждая запись выглядит как

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

для типа II и т.п.

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

для типа III.

Степень полинома равна ( RSIZE - 2 ) / 3 - 1. Сообщает MIDвам среднюю точку интервала (время эфемерид) и RADIUSсообщает вам половину размера интервала. Таким образом, его временной интервал равен [MID - RADIUS, MID + RADIUS].

Для получения позиции вы просто оцениваете Чебышева в соответствующее время. Для скорости вы либо оцениваете производную, либо другой полином. Не забудьте нормализовать время перед оценкой и денормализировать время после оценки (мне приходилось это делать много раз).

Работа напрямую с DAF может быть полезным опытом (например: у меня есть многопоточный высокопроизводительный оценщик эфемерид), но легко ошибиться. И после большой работы у вас есть только положение или скорость в данном кадре, которые вам все еще нужно вращать/перемещать/и т. д. (для чего вам нужно либо заново изобрести SPICE, либо положиться на него).

Если у вас нет очень веской причины не делать этого, просто придерживайтесь SPICE.

Простой скрипт Python

Я только что сделал скрипт Python, который может помочь кому-то начать работу с двоичным форматом DAF. Он просто выгружает информацию о записи файла и содержимое (сводки массива), содержащиеся в 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)

Например: запуск его для jup310.bsp дает

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)

Это означает, что jup310 содержит записи, охватывающие ET (-3155716758.8160305, 3155716867.183885) для некоторых спутников Юпитера (от 501 до 519) относительно барицентра Юпитера (5), в формате J2000 (1) и имеют тип 3. Другие записи относятся к барицентрам. Земли, Юпитера и Солнца (3, 5 и 10) относительно барицентра Солнечной системы (0) в J200 (1) и относятся к типу 2.

Каждая запись сообщает вам, где в файле начинается и заканчивается каждый массив. Вы можете расширить сценарий, чтобы пропустить каждое место и прочитать всю запись для последующей интерполяции.

Доступное программное обеспечение

Я быстро создал простую реализацию на Python, которая может читать и работать напрямую с типами II и III. Это хорошо сравнивается с CSPICE.

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

Я собираюсь не согласиться со всеми предыдущими ответами. Я бы сказал, что реализация вашей собственной программы Ephemeris является тривиальной задачей, если вы знаете формат файлов и имеете относительно простые уравнения для использования полиномов для вычисления позиций. Задача вполне по силам начинающему программисту. Формат файла - самая сложная часть, только потому, что нужно собрать много частей, но это не так сложно.

В статье « Формат файлов эфемерид JPL» есть довольно подробное описание с примерами и реализацией JavaScript в качестве примера. Обратите внимание, что весь класс, выполняющий поиск и вычисление коэффициентов, занимает менее 100 строк кода. Кроме того, в репозитории Github есть примеры реализации gmiller123456/jpl-development-ephemeris , которые в настоящее время имеют реализации на Python, Perl, C# и JavaScript.

Как вы уже поняли, раздел «ГРУППА 1050» заголовка содержит множество ключей к формату файла. Это раздел для 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

Порядок и количество компонентов, которым они соответствуют, приведены в документе формата Ascii от JPL. Каждый столбец представляет значения для данной планеты (или другого свойства):

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

Каждый из файлов ASCII разделен на блоки, действительные в течение количества дней, указанного в заголовочном файле. Первые два числа в блоке — это даты по юлианскому календарю, для которых этот блок действителен, поэтому найти блок, соответствующий данной дате по юлианскому календарю, несложно.

Начиная с ГРУППЫ 1050 в заголовке, числа в каждом столбце соответствуют:

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

Итак, чтобы вычислить смещение в блоке для данной планеты:

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

Итак, для планет у вас будут все коэффициенты X, затем все коэффициенты Y и, наконец, все коэффициенты Z. Например, коэффициенты для Меркурия для JD=2458850,5 приведены в приведенном ниже коде.

Чтобы позаимствовать пример из статьи « Формат файлов эфемерид JPL» , это фактически вычисляет позиции:

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

Вызов функции calculateExamplePolynomials(), приведенной выше, дает следующие значения:

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