Added spice2xyzv, a tool for generating Celestia xyzv files from a set of SPICE SPK kernels.

ver1_6_1
Chris Laurel 2008-05-06 23:31:00 +00:00
parent 24967b4eda
commit dd31bae5d5
7 changed files with 815 additions and 0 deletions

View File

@ -0,0 +1,26 @@
StartDate "1997-Oct-15 09:27"
EndDate "2004-Jun-20 13:00"
Observer "SUN"
Target "CASSINI"
Frame "eclipJ2000"
Tolerance 20
MinStep 60
MaxStep 432000
KernelDirectory "c:/dev/cassini/new"
Kernels [
"de421.bsp"
"000331R_SK_LP0_V1P32.bsp"
"000331RB_SK_V1P32_V2P12.bsp"
"000331R_SK_V2P12_EP15.bsp"
"010420R_SCPSE_EP1_JP83.bsp"
"010423_SK_JP67_SP0.bsp"
"020425B_SK_SM812_T45.bsp"
"030201AP_OPK_SM546_T45.bsp"
"040622AP_OPK_04122_08222.bsp"
"041001AP_SK_04275_08222.bsp"
"041210AP_OPK_04329_08189.bsp"
"050720AP_SCPSE_05119_08222.bsp"
"060323AP_SCPSE_06082_08222.bsp"
"070209AP_SCPSE_06356_08222.bsp"
"070918AP_SCPSE_07261_10191.bsp"
]

View File

@ -0,0 +1,27 @@
StartDate "2004-Jun-20 12:00"
EndDate "2010-Jul-10 11:00"
Observer "SATURN"
Target "CASSINI"
Frame "eclipJ2000"
Tolerance 20
MinStep 60
MaxStep 432000
KernelDirectory "c:/dev/cassini/new"
Kernels [
"de421.bsp"
"sat252s.bsp"
"000331R_SK_LP0_V1P32.bsp"
"000331RB_SK_V1P32_V2P12.bsp"
"000331R_SK_V2P12_EP15.bsp"
"010420R_SCPSE_EP1_JP83.bsp"
"010423_SK_JP67_SP0.bsp"
"020425B_SK_SM812_T45.bsp"
"030201AP_OPK_SM546_T45.bsp"
"040622AP_OPK_04122_08222.bsp"
"041001AP_SK_04275_08222.bsp"
"041210AP_OPK_04329_08189.bsp"
"050720AP_SCPSE_05119_08222.bsp"
"060323AP_SCPSE_06082_08222.bsp"
"070209AP_SCPSE_06356_08222.bsp"
"070918AP_SCPSE_07261_10191.bsp"
]

View File

@ -0,0 +1,14 @@
StartDate "2004-Dec-25 02:01:04.183 TDB"
EndDate "2005-Jan-14 09:07:00.000 TDB"
Observer "SATURN"
Target "CASP"
Frame "eclipJ2000"
Tolerance 0.1
MinStep 60
MaxStep 86400
KernelDirectory "c:/dev/cassini/new"
Kernels [
"de421.bsp"
"sat252s.bsp"
"041210AP_OPK_04329_08189.bsp"
]

View File

@ -0,0 +1,11 @@
# Visual C++ makefile for spice2xyzv.exe
TARGET=spice2xyzv.exe
SPICE_INSTALL=c:\dev\spice\cspice
SOURCE_FILES=\
spice2xyzv.cpp
$(TARGET): $(SOURCE_FILES)
cl /EHsc /Ox /MT $(SOURCE_FILES) /Fe$(TARGET) /I $(SPICE_INSTALL)\include $(SPICE_INSTALL)\lib\cspice.lib

View File

@ -0,0 +1,136 @@
KPL/LSK
LEAPSECONDS KERNEL FILE
===========================================================================
Modifications:
--------------
2005, Aug. 3 NJB Modified file to account for the leapsecond that
will occur on December 31, 2005.
1998, Jul 17 WLT Modified file to account for the leapsecond that
will occur on December 31, 1998.
1997, Feb 22 WLT Modified file to account for the leapsecond that
will occur on June 30, 1997.
1995, Dec 14 KSZ Corrected date of last leapsecond from 1-1-95
to 1-1-96.
1995, Oct 25 WLT Modified file to account for the leapsecond that
will occur on Dec 31, 1995.
1994, Jun 16 WLT Modified file to account for the leapsecond on
June 30, 1994.
1993, Feb. 22 CHA Modified file to account for the leapsecond on
June 30, 1993.
1992, Mar. 6 HAN Modified file to account for the leapsecond on
June 30, 1992.
1990, Oct. 8 HAN Modified file to account for the leapsecond on
Dec. 31, 1990.
Explanation:
------------
The contents of this file are used by the routine DELTET to compute the
time difference
[1] DELTA_ET = ET - UTC
the increment to be applied to UTC to give ET.
The difference between UTC and TAI,
[2] DELTA_AT = TAI - UTC
is always an integral number of seconds. The value of DELTA_AT was 10
seconds in January 1972, and increases by one each time a leap second
is declared. Combining [1] and [2] gives
[3] DELTA_ET = ET - (TAI - DELTA_AT)
= (ET - TAI) + DELTA_AT
The difference (ET - TAI) is periodic, and is given by
[4] ET - TAI = DELTA_T_A + K sin E
where DELTA_T_A and K are constant, and E is the eccentric anomaly of the
heliocentric orbit of the Earth-Moon barycenter. Equation [4], which ignores
small-period fluctuations, is accurate to about 0.000030 seconds.
The eccentric anomaly E is given by
[5] E = M + EB sin M
where M is the mean anomaly, which in turn is given by
[6] M = M + M t
0 1
where t is the number of ephemeris seconds past J2000.
Thus, in order to compute DELTA_ET, the following items are necessary.
DELTA_TA
K
EB
M0
M1
DELTA_AT after each leap second.
The numbers, and the formulation, are taken from the following sources.
1) Moyer, T.D., Transformation from Proper Time on Earth to
Coordinate Time in Solar System Barycentric Space-Time Frame
of Reference, Parts 1 and 2, Celestial Mechanics 23 (1981),
33-56 and 57-68.
2) Moyer, T.D., Effects of Conversion to the J2000 Astronomical
Reference System on Algorithms for Computing Time Differences
and Clock Rates, JPL IOM 314.5--942, 1 October 1985.
The variable names used above are consistent with those used in the
Astronomical Almanac.
\begindata
DELTET/DELTA_T_A = 32.184
DELTET/K = 1.657D-3
DELTET/EB = 1.671D-2
DELTET/M = ( 6.239996D0 1.99096871D-7 )
DELTET/DELTA_AT = ( 10, @1972-JAN-1
11, @1972-JUL-1
12, @1973-JAN-1
13, @1974-JAN-1
14, @1975-JAN-1
15, @1976-JAN-1
16, @1977-JAN-1
17, @1978-JAN-1
18, @1979-JAN-1
19, @1980-JAN-1
20, @1981-JUL-1
21, @1982-JUL-1
22, @1983-JUL-1
23, @1985-JUL-1
24, @1988-JAN-1
25, @1990-JAN-1
26, @1991-JAN-1
27, @1992-JUL-1
28, @1993-JUL-1
29, @1994-JUL-1
30, @1996-JAN-1
31, @1997-JUL-1
32, @1999-JAN-1
33, @2006-JAN-1 )
\begintext

View File

@ -0,0 +1,129 @@
spice2xyzv
Author: Chris Laurel (claurel@gmail.com)
Spice2xyzv is a command line tool that is used to produce Celestia xyzv files
from a set of SPICE SPK kernels. An xyzv file is simply an ASCII file with
position and velocity states at unequal time steps. Each record in the file
consists of seven double precision values, arranged as such:
<Julian date (TDB)> <x> <y> <z> <vx> <vy> <vz>
X, y, and z are components of the position vector; vx, vy, and vz are the
velocity vector components.
Why?
----
Celestia supports using SPICE SPK kernels directly, and when high accuracy is
required SPICE kernels are the preferred ephemeris source. But the large
amount of data required to represent an entire mission with SPK kernels can
make this approach less attractive. This is the case for the spacecraft
trajectories in the default Celestia package. Spice2xyzv offers a tolerance
parameter that lets you pick the optimal balance between size and accuracy for
the situation. Spice2xyzv can thus be thought of as a sort of lossy
compression tool for SPK files.
Usage
-----
The spice2xyzv command line is extremely simple:
spice2xyzv <config file>
The xyzv file is written to standard output, thus you'll generally use the
tool with output redirection, e.g:
spice2xyzv cassini-cruise.cfg > cruise.xyzv
The configuration file is a text file with a list of named parameters. These
parameters have either string, numeric, or string list values. Some of the
parameters have defaults and can be omitted from the file. The order in which
the parameters appear is not important. The following example is used to
produce an xyzv file for the Huygens probe:
StartDate "2004-Dec-25 02:01:04.183 TDB"
EndDate "2005-Jan-14 09:07:00.000 TDB"
Observer "SATURN"
Target "CASP"
Frame "eclipJ2000"
Tolerance 0.1
MinStep 60
MaxStep 86400
KernelDirectory "c:/dev/cassini/new"
Kernels [
"de421.bsp"
"sat252s.bsp"
"041210AP_OPK_04329_08189.bsp"
]
Configuration Parameters
------------------------
StartDate and EndDate (string)
These parameters specify the time range for the xyzv file. The set of SPK
kernels used must provide coverage for the target object (and any necessary
frame centers) for any time between start and end. The date strings are
parsed by SPICE, which accepts a wide range of different formats. By default,
all times are UTC, but the TDB or TDT suffix can be appended (as in the
example) to specify a different time system.
Observer (string)
The name or NAIF integer ID of the coordinate center. The choice of
coordinate center depends on how the resulting xyzv file will be used: it
must match the frame center chosen for the SSC file.
Target (string)
The name or NAIF integer ID of the object that you're generating a trajectory
for.
Frame (string)
The name of the reference frame for the xyzv coordinates. Two common frames
are:
"J2000" - Earth mean equator and equinox of J2000
"eclipJ2000" - Earth ecliptic and equinox of J2000
If the frame is not one of SPICE's built-in frames, it is necessary to
add a frame kernel to the Kernels list in the configuration file.
Tolerance (number)
Maximum permitted error in units of kilometers. Making this value smaller will
result in a larger xyzv file.
MinStep and MaxStep (number)
The minimum and maximum time interval (in seconds) between states. The
default values are 60 for MinStep and 432000 (5 days) for MaxStep.
KernelDirectory (string)
Directory containing the SPICE kernel files. The default value is ".",
indicating that the kernel files reside in the current directory.
Kernels (string list)
List of kernel files that need to be loaded to generate the xyzv file. This
list may include more than just SPK files. For instance, a frame kernel file
may be necessary to get the position of an object in the desired frame.
Method
------
Spice2xyzv uses a very simple method to choose sample times for xyzv states.
It calls SPICE to generate a state at a base time t0. It then generates two
more states: one at t0+dt/2 and one at t0+dt. Next, the position at t0+dt/2
is compared to the result of cubic Hermite interpolation of the SPICE
computed positions at t0 and t0+dt. If the distance is within the tolerance
specified in the configuration file, the test is repeated with dt*2. This
continues until either MaxStep is reached or the interpolated and SPICE
calculated positions are further than Tolerance kilometers apart. The last
value of dt for which the interpolated position was close enough the the
SPICE calculated position is used as the time step, t0 is incremented by
dt, and the process is repeated over the entire time span. This adaptive
sampling results in a low number of samples in slowly varying parts of the
trajectory and more samples at times when the trajectory changes more
dramatically.

View File

@ -0,0 +1,472 @@
// spice2xyzv.cpp
//
// Copyright (C) 2008, Chris Laurel <claurel@shatters.net>
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 2
// of the License, or (at your option) any later version.
//
// Create a Celestia xyzv file from a pool of SPICE SPK files
#include "SpiceUsr.h"
#include <string>
#include <vector>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
using namespace std;
const double J2000 = 2451545.0;
// Default values
// Units are seconds
const double MIN_STEP_SIZE = 60.0;
const double MAX_STEP_SIZE = 5 * 86400.0;
// Units are kilometers
const double TOLERANCE = 20.0;
class Configuration
{
public:
Configuration() :
kernelDirectory("."),
frameName("eclipJ2000"),
minStepSize(MIN_STEP_SIZE),
maxStepSize(MAX_STEP_SIZE),
tolerance(TOLERANCE)
{
}
string kernelDirectory;
vector<string> kernelList;
string startDate;
string endDate;
string observerName;
string targetName;
string frameName;
double minStepSize;
double maxStepSize;
double tolerance;
};
// Very basic 3-vector class
class Vec3d
{
public:
Vec3d() : x(0.0), y(0.0), z(0.0) {}
Vec3d(double _x, double _y, double _z) : x(_x), y(_y), z(_z) {}
Vec3d(const double v[]) : x(v[0]), y(v[1]), z(v[2]) {}
double length() const { return std::sqrt(x * x + y * y + z * z); }
double x, y, z;
};
// Vector add
Vec3d operator+(const Vec3d& v0, const Vec3d& v1)
{
return Vec3d(v0.x + v1.x, v0.y + v1.y, v0.z + v1.z);
}
// Vector subtract
Vec3d operator-(const Vec3d& v0, const Vec3d& v1)
{
return Vec3d(v0.x - v1.x, v0.y - v1.y, v0.z - v1.z);
}
// Additive inverse
Vec3d operator-(const Vec3d& v)
{
return Vec3d(-v.x, -v.y, -v.z);
}
// Scalar multiply
Vec3d operator*(const Vec3d& v, double d)
{
return Vec3d(v.x * d, v.y * d, v.z * d);
}
// Scalar multiply
Vec3d operator*(double d, const Vec3d& v)
{
return Vec3d(v.x * d, v.y * d, v.z * d);
}
ostream& operator<<(ostream& o, const Vec3d& v)
{
return o << v.x << ' ' << v.y << ' ' << v.z;
}
// The StateVector object just contains the position and velocity
// in 3D.
class StateVector
{
public:
// Construct a new StateVector from an array of 6 doubles
// (as used by SPICE.)
StateVector(const double v[]) :
position(v), velocity(v + 3) {};
Vec3d position;
Vec3d velocity;
};
// QuotedString is used read a double quoted string from a C++ input
// stream.
class QuotedString
{
public:
string value;
};
istream& operator>>(istream& in, QuotedString& qs)
{
char c = '\0';
in >> c;
while (in && isspace(c))
{
in >> c;
}
if (c != '"')
{
in.setstate(ios::failbit);
return in;
}
string s;
in >> c;
while (in && c != '"')
{
s += c;
in.get(c);
}
if (in)
qs.value = s;
return in;
}
// QuoteStringList is used to read a list of double quoted strings from
// a C++ input stream. The string list must be enclosed by square brackets.
class QuotedStringList
{
public:
vector<string> value;
};
istream& operator>>(istream& in, QuotedStringList& qsl)
{
qsl.value.clear();
char c = '\0';
in >> c;
if (c != '[')
{
in.setstate(ios::failbit);
return in;
}
in >> c;
while (in && c == '"')
{
in.unget();
QuotedString qs;
if (in >> qs)
{
qsl.value.push_back(qs.value);
in >> c;
}
}
if (c != ']')
in.setstate(ios::failbit);
return in;
}
static Vec3d cubicInterpolate(const Vec3d& p0, const Vec3d& v0,
const Vec3d& p1, const Vec3d& v1,
double t)
{
return p0 + (((2.0 * (p0 - p1) + v1 + v0) * (t * t * t)) +
((3.0 * (p1 - p0) - 2.0 * v0 - v1) * (t * t)) +
(v0 * t));
}
double et2jd(double et)
{
return J2000 + et / 86400.0;
}
string etToString(double et)
{
char buf[200];
et2utc_c(et, "C", 3, sizeof(buf), buf);
return string(buf);
}
void printRecord(ostream& out, double et, const StateVector& state)
{
// < 1 second error around J2000
out << setprecision(12) << et2jd(et) << " ";
// < 1 meter error at 1 billion km
out << setprecision(12) << state.position << " ";
// < 0.1 mm/s error at 10 km/s
out << setprecision(8) << state.velocity << endl;
}
StateVector getStateVector(SpiceInt targetID,
double et,
const string& frameName,
SpiceInt observerID)
{
double stateVector[6];
double lightTime = 0.0;
spkgeo_c(targetID, et, frameName.c_str(), observerID,
stateVector, &lightTime);
return StateVector(stateVector);
}
bool convertSpkToXyzv(const Configuration& config,
ostream& out)
{
// Load the required SPICE kernels
for (vector<string>::const_iterator iter = config.kernelList.begin();
iter != config.kernelList.end(); iter++)
{
string pathname = config.kernelDirectory + "/" + *iter;
furnsh_c(pathname.c_str());
}
double startET = 0.0;
double endET = 0.0;
str2et_c(config.startDate.c_str(), &startET);
str2et_c(config.endDate.c_str(), &endET);
SpiceBoolean found = SPICEFALSE;
SpiceInt observerID = 0;
SpiceInt targetID = 0;
bodn2c_c(config.observerName.c_str(), &observerID, &found);
if (!found)
{
cerr << "Observer object " << config.observerName << " not found. Aborting.\n";
return false;
}
bodn2c_c(config.targetName.c_str(), &targetID, &found);
if (!found)
{
cerr << "Target object " << config.targetName << " not found. Aborting.\n";
return false;
}
StateVector lastState = getStateVector(targetID, startET, config.frameName, observerID);
double et = startET;
printRecord(out, et, lastState);
while (et + config.minStepSize < endET)
{
double dt = config.minStepSize;
StateVector s0 = getStateVector(targetID, et + dt, config.frameName, observerID);
double et0 = et + dt;
while (dt < config.maxStepSize && et + dt * 2.0 < endET)
{
dt *= 2.0;
StateVector s1 = getStateVector(targetID, et + dt, config.frameName, observerID);
double et1 = et + dt;
Vec3d pInterp = cubicInterpolate(lastState.position,
lastState.velocity * dt,
s1.position,
s1.velocity * dt,
0.5);
double positionError = (pInterp - s0.position).length();
if (positionError > config.tolerance || dt > config.maxStepSize)
break;
s0 = s1;
et0 = et1;
}
lastState = s0;
et = et0;
printRecord(out, et0, lastState);
}
lastState = getStateVector(targetID, endET, config.frameName, observerID);
printRecord(out, endET, lastState);
return true;
}
bool readConfig(istream& in, Configuration& config)
{
QuotedString qs;
while (in && !in.eof())
{
string key;
in >> key;
if (in.eof())
return true;
if (!in.eof())
{
if (key == "StartDate")
{
if (in >> qs)
config.startDate = qs.value;
}
else if (key == "EndDate")
{
if (in >> qs)
config.endDate = qs.value;
}
else if (key == "Observer")
{
if (in >> qs)
config.observerName = qs.value;
}
else if (key == "Target")
{
if (in >> qs)
config.targetName = qs.value;
}
else if (key == "Frame")
{
if (in >> qs)
config.frameName = qs.value;
}
else if (key == "MinStep")
{
in >> config.minStepSize;
}
else if (key == "MaxStep")
{
in >> config.maxStepSize;
}
else if (key == "Tolerance")
{
in >> config.tolerance;
}
else if (key == "KernelDirectory")
{
if (in >> qs)
config.kernelDirectory = qs.value;
}
else if (key == "Kernels")
{
QuotedStringList qsl;
if (in >> qsl)
{
config.kernelList = qsl.value;
}
}
}
}
return in.good();
}
int main(int argc, char* argv[])
{
// Load the leap second kernel
furnsh_c("naif0008.tls");
if (argc < 2)
{
cerr << "Usage: spice2xyzv <config filename> [output filename]\n";
return 1;
}
ifstream configFile(argv[1]);
if (!configFile)
{
cerr << "Error opening configuration file.\n";
return 1;
}
Configuration config;
if (!readConfig(configFile, config))
{
cerr << "Error in configuration file.\n";
return 1;
}
// Check that all required parameters are present.
if (config.startDate.empty())
{
cerr << "StartDate missing from configuration file.\n";
return 1;
}
if (config.endDate.empty())
{
cerr << "EndDate missing from configuration file.\n";
return 1;
}
if (config.targetName.empty())
{
cerr << "Target missing from configuration file.\n";
return 1;
}
if (config.observerName.empty())
{
cerr << "Observer missing from configuration file.\n";
return 1;
}
if (config.kernelList.empty())
{
cerr << "Kernels missing from configuration file.\n";
return 1;
}
convertSpkToXyzv(config, cout);
return 0;
}