Я хочу создать файл эфемерид, содержащий коэффициенты Чебышева, аналогично тому, как мы можем создать файл эфемерид позиций с помощью инструмента Horizon ( https://ssd.jpl.nasa.gov/horizons/app.html#/ ). Я хотел бы сделать то же самое, то есть, скажем, иметь возможность выбирать целевое тело, центр координат, даты и генерировать файл эфемерид, содержащий связанные коэффициенты Чебышева.
Знаете ли вы инструмент или программу, которая может это сделать? Я провел много исследований, но не нашел ни одного.
Я пытался использовать SPICE, но не нашел метода, позволяющего извлекать коэффициенты Чебышева таким образом.
Ваш вопрос немного неоднозначен, файл, содержащий коэффициенты Чебышева, не будет эфемеридой, но его можно использовать для создания эфемерид. Несмотря на это, статья « Формат файлов эфемерид JPL» ответит как на то, как получить коэффициенты, так и на то, как их использовать для создания эфемерид.
В статье содержится более подробная информация, но краткий ответ: посмотрите на заголовочный файл для JPL Development Ephemeris , который вы хотите использовать. В разделе «Группа 1050» содержится самая важная информация. Первый столбец — Меркурий, затем Венера, барицентр Земля-Луна… 9-й — Плутон, Луна и Солнце.
Коэффициенты сгруппированы в 32-дневные блоки, отмеченные юлианскими днями, для которых они действительны. Приведенная выше информация о группе 1050 показывает, как разбивается каждый блок. Первая строка показывает смещение в блоке, где начинаются коэффициенты для этой планеты, вторая — количество коэффициентов для каждого свойства (например, X, Y и Z), а последняя строка — количество подинтервалов каждые 32 дневной блок разбит на.
Например, ряд Меркурия равен 3, 14, 4. Таким образом, он начинается со смещения 3, имеет 14 коэффициентов на свойство (x, y, z) = 3 * 14 = 42 и делится на 4 подинтервала, что в сумме из 42 * 3 = 168 коэффициентов. Обратите внимание, что столбец для Венеры имеет смещение 171, что является смещением Меркурия плюс общие коэффициенты.
Когда у вас есть 168 коэффициентов, вам нужно определить, какой подинтервал вам нужен, поскольку Меркурий делится на 4 подинтервала, или 4/32 = 8-дневные интервалы. Первые две записи каждого блока предоставляют действительный диапазон юлианских дней, поэтому просто определите, для какого интервала вы хотите выполнить расчет, и выберите соответствующие 42 коэффициента для этого диапазона.
Из этих 42 коэффициентов первые 14 относятся к координате x, следующие 14 — к координате y, а последние — к координате z. Это используемые коэффициенты Чебышева, первая ссылка выше содержит пример извлечения коэффициентов и выполнения вычислений, а также исходный код на JavaScript.
Этот репозиторий github содержит исходный код на нескольких языках, JavaScript, Python, Java, C#, Perl и, возможно, на других.
Лучший способ, который я нашел, — это использовать «jplephem» Python Брэндона Роудса отсюда , в частности, используя _load()функцию сегмента SPK. Вы получаете список сегментов SPK при загрузке файла BSP.
В качестве быстрого плагина я являюсь частью команды из четырех человек, которые только начали работать над ANISE несколько дней назад: https://github.com/anise-toolkit/ . План состоит в том, чтобы создать версию SPICE с открытым исходным кодом (Mozilla Public License 2.0), готовую к полету. Я выполнял подобную работу раньше для частной компании (и, следовательно, эта работа является собственностью, поэтому у меня нет к ней доступа), потому что CSPICE не был жизнеспособной альтернативой в то время. Что касается ANISE, мы ищем новых людей, которые присоединятся к обсуждению и команде разработчиков, поэтому, если вы заинтересованы, обращайтесь в нашу группу Element/Matrix ( https://matrix.to/#/#anise:matrix. орг ).
Я нашел способ извлечь коэффициенты Чебышева на запрошенную дату для конкретного тела с помощью CSPICE.
Вот программа, которая позволяет это сделать:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "SpiceUsr.h"
/*
Author : G. Juif
spk_to_asc.c
Computation : gcc spk_to_asc.c -Iinclude lib/cspice.a -lm -o spk_to_asc
Input : meta-kernel file (name must be modified in the code)
This script extract Chebyshev coefficients from a SPICE bsp file.
The body for which we want the coefficients must be defined in the code. It is by default MARS BARYCENTER.
List of name bodies can be find here : https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/naif_ids.html#Planets%20and%20Satellites
Care to have the associated bsp file.
The coefficients are written in 3 differents files for X, Y and Z components.
At execution several parameters must be given :
- Start date in calendar format AAAA/MM/JJ. It is date from which the coefficients will be extract.
By default this date is in TDB time scale but it can be modified in the code.
- Duration in days (integer) : duration from the start date
- Time step in days (integer) : Chebyshev coefficients will be extract at each of theses days step
Output 3 files : coeffs_x ; ceffs_y ; coeffs_z
In these files, at each time step, it gives :
# Step_number Date Order Date_start Date_end
Coefficients at each order line per line
Dates are written in CJD format (number of julian days since January 1st, 1950 0h).
It can be modified to have J2000 julian days by deleting 18262.5 additions in this code.
- Step_number : gives the step number
- Date : Date (TDB by default) corresponding of Chebyshev coefficients given below
- Order : Order of the Chebyshev polynom
- Date_start Date_end : dates where Chebyshev coefficients are defined in JPL file in CJD days.
It is usefull in particular to compute scaled time for the Chebyshev polynom.
*/
int main()
{
/*
Local parameters
*/
#define ND 2
#define NI 6
#define DSCSIZ 5
#define SIDLEN1 41
#define MAXLEN 50
#define utclen 35
/*
Local variables
*/
SpiceBoolean found;
SpiceChar segid [ SIDLEN1 ];
SpiceChar utcstr[ utclen ];
SpiceChar calendar[ utclen ];
SpiceChar frname [MAXLEN];
SpiceChar cname [MAXLEN];
SpiceChar bname [MAXLEN];
SpiceDouble dc [ ND ];
SpiceDouble descr [ DSCSIZ ];
SpiceDouble et;
SpiceDouble record [99];
SpiceDouble date [99];
SpiceDouble dateCJD;
SpiceDouble Cheby_order;
SpiceDouble Interv_length;
SpiceDouble Start_interv_cheby;
SpiceDouble End_interv_cheby;
SpiceInt handle;
SpiceInt ic [ NI ];
SpiceInt idcode;
SpiceInt temp;
SpiceInt i = 0;
SpiceInt j = 0;
SpiceInt nbJours = 0;
SpiceInt nbInterv = 0;
SpiceInt nbCoeffs = 0;
FILE * eph_file_x;
FILE * eph_file_y;
FILE * eph_file_z;
eph_file_x = fopen("coeffs_x", "w");
eph_file_y = fopen("coeffs_y", "w");
eph_file_z = fopen("coeffs_z", "w");
/*
Load a meta-kernel that specifies a planetary SPK file
and leapseconds kernel. The contents of this meta-kernel
are displayed above.
*/
furnsh_c ( "spksfs_ex1.tm" );
printf("Retrieve Chebyshev coefficients at a given date with duration and time step in days\n");
/*
Get the NAIF ID code for the Pluto system barycenter.
This is a built-in ID code, so something's seriously
wrong if we can't find the code.
*/
bodn2c_c ( "MARS BARYCENTER", &idcode, &found );
if ( !found )
{
sigerr_c( "SPICE(BUG)" );
}
/*
Pick a request time; convert to seconds past J2000 TDB.
*/
printf("Enter start date aaaa/mm/jj (TDB time scale):\n");
scanf("%12s", calendar);
strcat(calendar," TDB");
str2et_c ( calendar, &et );
et2utc_c ( et, "J", 7, utclen, utcstr );
printf ( "Date : %s \n",calendar);
printf("Enter duration in days :\n");
scanf("%d", &nbInterv);
printf("Enter time step in days :\n");
scanf("%d", &nbJours);
/* Loop on et */
nbInterv /= nbJours;
date[0] = et;
for (i = 0 ; i < nbInterv ; i++)
{
/*
Find a loaded segment for the specified body and time.
*/
spksfs_c ( idcode, date[i], SIDLEN1, &handle, descr, segid, &found );
if ( !found )
{
printf ( "No descriptor was found for ID %d at "
"TDB %24.17e\n",
(int) idcode,
et );
}
else
{
/* Convert date in CJD CNES date */
dateCJD = (date[i]/86400 ) + 18262.5;
/*
Display the segment ID.
Unpack the descriptor. Display the contents.
*/
dafus_c ( descr, ND, NI, dc, ic );
temp = spkr02_(&handle,descr,&date[i],record);
/* Chebyshev polynom order (minus 2 because length of record doesn't consider first element, see fortran spice doc)*/
Cheby_order = (record[0] - 2)/3;
/* Interval length of chebyshev coefficients in days */
Interv_length = (record[2]/86400)*2;
/* Start and end interval dates where Chebyshev coefficients are defined in JPL file in CJD days*/
Start_interv_cheby = (record[1]/86400) + 18262.5 - Interv_length/2 ;
End_interv_cheby = (record[1]/86400) + 18262.5 + Interv_length/2 ;
/* Print informations in files */
fprintf(eph_file_x, "# %ld %lf %lf %lf %lf\n", i+1, dateCJD, Cheby_order, Start_interv_cheby, End_interv_cheby);
fprintf(eph_file_y, "# %ld %lf %lf %lf %lf\n", i+1, dateCJD, Cheby_order, Start_interv_cheby, End_interv_cheby);
fprintf(eph_file_z, "# %ld %lf %lf %lf %lf\n", i+1, dateCJD, Cheby_order, Start_interv_cheby, End_interv_cheby);
nbCoeffs = (int) Cheby_order ;
/* Coeffs for X,Y,Z component */
for (j = 0 ; j < nbCoeffs ; j++)
{
fprintf(eph_file_x, "%24.17e\n", record[3+j]);
fprintf(eph_file_y, "%24.17e\n", record[3+nbCoeffs+j]);
fprintf(eph_file_z, "%24.17e\n", record[3+2*nbCoeffs+j]);
}
}
/* Compute next date in seconds past J2000 TDB */
date[i+1] = date[i] + 86400*nbJours;
}
/* Translate SPICE codes into common names */
frmnam_c ( (int) ic[2], MAXLEN, frname );
bodc2n_c ( (int) ic[1], MAXLEN, cname, &found );
bodc2n_c ( (int) ic[0], MAXLEN, bname, &found );
/* Print configuration*/
printf ( "Segment ID: %s\n"
"\n--------Configuration-------\n"
"Body ID code: %s\n"
"Center ID code: %s\n"
"Frame ID code: %s\n"
"SPK data type: %d\n"
"Start ephemeris file time (TDB): %24.17e\n"
"Stop ephemeris file time (TDB): %24.17e\n"
"\n--------Chenbyshev polynom informations-------\n"
"Chebyshev polynom order: %lf\n"
"Time step in days where Chebyshev coefficients are defined: %lf\n",
segid,
bname,
cname,
frname,
(int) ic[3],
dc[0],
dc[1],
Cheby_order,
Interv_length );
fclose(eph_file_x);
fclose(eph_file_y);
fclose(eph_file_z);
return ( 0 );
}
И файл мета-ядра:
KPL/MK
\begindata
KERNELS_TO_LOAD = ( 'de430.bsp',
'naif0011.tls' )
\begintext
End of meta-kernel
Может быть, это помогает?
https://www.orekit.org/static/xref/org/orekit/bodies/PosVelChebyshev.html
Это из orekit, java-библиотеки с открытым исходным кодом.
Альфонсо Гонсалес
Дэвид Хаммен
Гийом Ж
Гийом Ж
ооо