Update :
I have found a way to extract Chebyshev coefficients at a requested date for a particular body using CSPICE.
Here the program that allows to do this :
#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 );
}
And the meta-kernal file :
KPL/MK
\begindata
KERNELS_TO_LOAD = ( 'de430.bsp',
'naif0011.tls' )
\begintext
End of meta-kernel
