Quiero generar un archivo de efemérides que contenga los coeficientes de Chebyshev, similar a cómo podemos generar un archivo de efemérides de posiciones con la herramienta Horizon ( https://ssd.jpl.nasa.gov/horizons/app.html#/ ). Me gustaría hacer lo mismo, es decir, poder elegir el cuerpo objetivo, el centro de coordenadas, las fechas y generar el archivo de efemérides que contiene los coeficientes de Chebyshev asociados.
¿Conoces una herramienta o un programa que sea capaz de hacer esto? He investigado mucho, pero no he encontrado ninguno.
Probé con SPICE pero no encontré un método que permitiera extraer coeficientes de Chebyshev como este.
Su pregunta es un poco ambigua, un archivo que contiene los coeficientes de Chebyshev no sería una efeméride, pero es algo que podría usar para generar una efeméride. Independientemente, el artículo Formato de los archivos de efemérides JPL responderá tanto cómo obtener los coeficientes como cómo usarlos para generar una efemérides.
El artículo entra en más detalles, pero la respuesta corta: mire el archivo de encabezado para las efemérides de desarrollo JPL que desea usar. La sección denominada "Grupo 1050" tiene la información más importante. La primera columna es para Mercurio, luego Venus, el Baricentro Tierra-Luna... la novena es Plutón, la Luna y el Sol.
Los coeficientes se agrupan en bloques de 32 días, marcados con los Días Julianos para los que son válidos. La información anterior del Grupo 1050 muestra cómo se desglosa cada bloque. La primera fila muestra el desplazamiento en el bloque donde comienzan los coeficientes para ese planeta, la segunda es la cantidad de coeficientes para cada propiedad (por ejemplo, X, Y y Z), y la última fila es la cantidad de subintervalos cada 32 el bloque de días se desglosa en.
Por ejemplo, la fila de Mercurio es 3, 14, 4. Entonces comienza en el desplazamiento 3, tiene 14 coeficientes por propiedad (x, y, z) = 3 * 14 = 42 y se divide en 4 subintervalos, para un total de 42 * 3 = 168 coeficientes. Observe que la columna de Venus tiene una compensación de 171, que es la compensación de Mercurio más los coeficientes totales.
Una vez que tenga los 168 coeficientes, debe determinar qué subintervalo necesita, ya que el mercurio se divide en 4 subintervalos, o 4/32 = intervalos de 8 días. Las dos primeras entradas de cada bloque proporcionan el rango de días julianos válido, así que simplemente determine para qué intervalo desea calcular y elija los 42 coeficientes correspondientes para ese rango.
Con esos 42 coeficientes, los primeros 14 son para la coordenada x, los siguientes 14 para la coordenada y y los últimos son para la coordenada z. Estos son los coeficientes de Chebyshev a usar, el primer enlace de arriba proporciona un ejemplo de cómo extraer los coeficientes y realizar el cálculo, así como el código fuente en JavaScript.
Este repositorio de github contiene código fuente en varios lenguajes, JavaScript, Python, Java, C#, Perl y quizás otros.
El mejor método que he encontrado es usar Python "jplephem" de Brandon Rhodes desde aquí , específicamente usando la _load()función de un segmento SPK. Obtiene una lista de segmentos SPK al cargar un archivo BSP.
Como complemento rápido, formo parte de un equipo de cuatro personas que comenzaron a trabajar en ANISE hace unos días: https://github.com/anise-toolkit/ . El plan es crear una versión de código abierto (Mozilla Public License 2.0) de SPICE que esté lista para el software de vuelo. He realizado un trabajo similar antes para una empresa privada (y, por lo tanto, ese trabajo es propietario, por lo que no tengo acceso a él) porque CSPICE no era una alternativa viable en ese momento. Para ANISE, estamos buscando a más personas para que se unan al equipo de conversación y desarrollo, así que si está interesado, comuníquese con nuestro espacio Element/Matrix ( https://matrix.to/#/#anise:matrix. org ).
Encontré una manera de extraer los coeficientes de Chebyshev en una fecha solicitada para un cuerpo en particular usando CSPICE.
Aquí el programa que permite hacer esto:
#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 );
}
Y el archivo meta-kernal:
KPL/MK
\begindata
KERNELS_TO_LOAD = ( 'de430.bsp',
'naif0011.tls' )
\begintext
End of meta-kernel
¿Quizás esto ayude?
https://www.orekit.org/static/xref/org/orekit/bodies/PosVelChebyshev.html
Esto es de orekit, una biblioteca java de código abierto
alfonso gonzalez
david hamen
GuillaumeJ
GuillaumeJ
UH oh