spacecruft
/
celestia
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Implemented xyz+velocity sampled trajectories (xyzv files) for better accuracy.

ver1_6_1
Chris Laurel 2008-04-03 01:26:36 +00:00
parent f86ba2b6dd
commit 9017964c44
3 changed files with 425 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

409
src/celengine/samporbit.cpp
View File

@ -30,17 +30,33 @@ static const double MinSampleInterval = 1.0 / 1440.0; // one minute
static const double MaxSampleInterval = 50.0;
static const double SampleThresholdAngle = 2.0;
// Position-only sample
template <typename T> struct Sample
{
double t;
T x, y, z;
};
// Position + velocity sample
template <typename T> struct SampleXYZV
{
double t;
Vector3<T> position;
Vector3<T> velocity;
};
template <typename T> bool operator<(const Sample<T>& a, const Sample<T>& b)
{
return a.t < b.t;
}
template <typename T> bool operator<(const SampleXYZV<T>& a, const SampleXYZV<T>& b)
{
return a.t < b.t;
}
template <typename T> class SampledOrbit : public CachingOrbit
{
public:
@ -123,9 +139,9 @@ template <typename T> double SampledOrbit<T>::getBoundingRadius() const
}
static Point3d cubicInterpolate(const Point3d& p0, const Vec3d& v0,
const Point3d& p1, const Vec3d& v1,
double t)
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)) +
@ -133,8 +149,8 @@ static Point3d cubicInterpolate(const Point3d& p0, const Vec3d& v0,
}
static Vec3d cubicInterpolateVelocity(const Point3d& p0, const Vec3d& v0,
const Point3d& p1, const Vec3d& v1,
static Vec3d cubicInterpolateVelocity(const Vec3d& p0, const Vec3d& v0,
const Vec3d& p1, const Vec3d& v1,
double t)
{
return ((2.0 * (p0 - p1) + v1 + v0) * (3.0 * t * t)) +
@ -145,14 +161,14 @@ static Vec3d cubicInterpolateVelocity(const Point3d& p0, const Vec3d& v0,
template <typename T> Point3d SampledOrbit<T>::computePosition(double jd) const
{
Point3d pos;
Vec3d pos;
if (samples.size() == 0)
{
pos = Point3d(0.0, 0.0, 0.0);
pos = Vec3d(0.0, 0.0, 0.0);
}
else if (samples.size() == 1)
{
pos = Point3d(samples[0].x, samples[1].y, samples[2].z);
pos = Vec3d(samples[0].x, samples[1].y, samples[2].z);
}
else
{
@ -175,7 +191,7 @@ template <typename T> Point3d SampledOrbit<T>::computePosition(double jd) const
if (n == 0)
{
pos = Point3d(samples[n].x, samples[n].y, samples[n].z);
pos = Vec3d(samples[n].x, samples[n].y, samples[n].z);
}
else if (n < (int) samples.size())
{
@ -185,9 +201,9 @@ template <typename T> Point3d SampledOrbit<T>::computePosition(double jd) const
Sample<T> s1 = samples[n];
double t = (jd - s0.t) / (s1.t - s0.t);
pos = Point3d(Mathd::lerp(t, (double) s0.x, (double) s1.x),
Mathd::lerp(t, (double) s0.y, (double) s1.y),
Mathd::lerp(t, (double) s0.z, (double) s1.z));
pos = Vec3d(Mathd::lerp(t, (double) s0.x, (double) s1.x),
Mathd::lerp(t, (double) s0.y, (double) s1.y),
Mathd::lerp(t, (double) s0.z, (double) s1.z));
}
else if (interpolation == TrajectoryInterpolationCubic)
{
@ -206,8 +222,8 @@ template <typename T> Point3d SampledOrbit<T>::computePosition(double jd) const
double h = s2.t - s1.t;
double ih = 1.0 / h;
double t = (jd - s1.t) * ih;
Point3d p0(s1.x, s1.y, s1.z);
Point3d p1(s2.x, s2.y, s2.z);
Vec3d p0(s1.x, s1.y, s1.z);
Vec3d p1(s2.x, s2.y, s2.z);
Vec3d v10((double) s1.x - (double) s0.x,
(double) s1.y - (double) s0.y,
@ -228,12 +244,12 @@ template <typename T> Point3d SampledOrbit<T>::computePosition(double jd) const
else
{
// Unknown interpolation type
pos = Point3d(0.0, 0.0, 0.0);
pos = Vec3d(0.0, 0.0, 0.0);
}
}
else
{
pos = Point3d(samples[n - 1].x, samples[n - 1].y, samples[n - 1].z);
pos = Vec3d(samples[n - 1].x, samples[n - 1].y, samples[n - 1].z);
}
}
@ -298,8 +314,8 @@ template <typename T> Vec3d SampledOrbit<T>::computeVelocity(double jd) const
double h = s2.t - s1.t;
double ih = 1.0 / h;
double t = (jd - s1.t) * ih;
Point3d p0(s1.x, s1.y, s1.z);
Point3d p1(s2.x, s2.y, s2.z);
Vec3d p0(s1.x, s1.y, s1.z);
Vec3d p1(s2.x, s2.y, s2.z);
Vec3d v10((double) s1.x - (double) s0.x,
(double) s1.y - (double) s0.y,
@ -385,6 +401,292 @@ template <typename T> void SampledOrbit<T>::sample(double start, double t, int,
}
// Sampled orbit with positions and velocities
template <typename T> class SampledOrbitXYZV : public CachingOrbit
{
public:
SampledOrbitXYZV(TrajectoryInterpolation);
virtual ~SampledOrbitXYZV();
void addSample(double t, const Vec3d& position, const Vec3d& velocity);
void setPeriod();
double getPeriod() const;
double getBoundingRadius() const;
Point3d computePosition(double jd) const;
Vec3d computeVelocity(double jd) const;
bool isPeriodic() const;
void getValidRange(double& begin, double& end) const;
virtual void sample(double, double, int, OrbitSampleProc& proc) const;
private:
vector<SampleXYZV<T> > samples;
double boundingRadius;
double period;
mutable int lastSample;
TrajectoryInterpolation interpolation;
};
template <typename T> SampledOrbitXYZV<T>::SampledOrbitXYZV(TrajectoryInterpolation _interpolation) :
boundingRadius(0.0),
period(1.0),
lastSample(0),
interpolation(_interpolation)
{
}
template <typename T> SampledOrbitXYZV<T>::~SampledOrbitXYZV()
{
}
// Add a new sample to the trajectory:
// Position in km
// Velocity in km/Julian day
template <typename T> void SampledOrbitXYZV<T>::addSample(double t, const Vec3d& position, const Vec3d& velocity)
{
double r = sqrt(position.length());
if (r > boundingRadius)
boundingRadius = r;
SampleXYZV<T> samp;
samp.t = t;
samp.position = Vector3<T>((T) position.x, (T) position.y, (T) position.z);
samp.velocity = Vector3<T>((T) velocity.x, (T) velocity.y, (T) velocity.z);
samples.push_back(samp);
}
template <typename T> double SampledOrbitXYZV<T>::getPeriod() const
{
return samples[samples.size() - 1].t - samples[0].t;
}
template <typename T> bool SampledOrbitXYZV<T>::isPeriodic() const
{
return false;
}
template <typename T> void SampledOrbitXYZV<T>::getValidRange(double& begin, double& end) const
{
begin = samples[0].t;
end = samples[samples.size() - 1].t;
}
template <typename T> double SampledOrbitXYZV<T>::getBoundingRadius() const
{
return boundingRadius;
}
template <typename T> Point3d SampledOrbitXYZV<T>::computePosition(double jd) const
{
Vec3d pos;
if (samples.size() == 0)
{
pos = Vec3d(0.0, 0.0, 0.0);
}
else if (samples.size() == 1)
{
pos = Vec3d(samples[0].position.x, samples[1].position.y, samples[2].position.z);
}
else
{
SampleXYZV<T> samp;
samp.t = jd;
int n = lastSample;
if (n < 1 || n >= (int) samples.size() || jd < samples[n - 1].t || jd > samples[n].t)
{
typename vector<SampleXYZV<T> >::const_iterator iter = lower_bound(samples.begin(),
samples.end(),
samp);
if (iter == samples.end())
n = samples.size();
else
n = iter - samples.begin();
lastSample = n;
}
if (n == 0)
{
pos = Vec3d(samples[n].position.x, samples[n].position.y, samples[n].position.z);
}
else if (n < (int) samples.size())
{
SampleXYZV<T> s0 = samples[n - 1];
SampleXYZV<T> s1 = samples[n];
if (interpolation == TrajectoryInterpolationLinear)
{
double t = (jd - s0.t) / (s1.t - s0.t);
pos = Vec3d(Mathd::lerp(t, (double) s0.position.x, (double) s1.position.x),
Mathd::lerp(t, (double) s0.position.y, (double) s1.position.y),
Mathd::lerp(t, (double) s0.position.z, (double) s1.position.z));
}
else if (interpolation == TrajectoryInterpolationCubic)
{
double h = s1.t - s0.t;
double ih = 1.0 / h;
double t = (jd - s0.t) * ih;
Vec3d p0(s0.position.x, s0.position.y, s0.position.z);
Vec3d v0(s0.velocity.x, s0.velocity.y, s0.velocity.z);
Vec3d p1(s1.position.x, s1.position.y, s1.position.z);
Vec3d v1(s1.velocity.x, s1.velocity.y, s1.velocity.z);
pos = cubicInterpolate(p0, v0 * h, p1, v1 * h, t);
}
else
{
// Unknown interpolation type
pos = Vec3d(0.0, 0.0, 0.0);
}
}
else
{
pos = Vec3d(samples[n - 1].position.x, samples[n - 1].position.y, samples[n - 1].position.z);
}
}
// Add correction for Celestia's coordinate system
return Point3d(pos.x, pos.z, -pos.y);
}
// Velocity is computed as the derivative of the interpolating function
// for position.
template <typename T> Vec3d SampledOrbitXYZV<T>::computeVelocity(double jd) const
{
Vec3d vel(0.0, 0.0, 0.0);
if (samples.size() >= 2)
{
SampleXYZV<T> samp;
samp.t = jd;
int n = lastSample;
if (n < 1 || n >= (int) samples.size() || jd < samples[n - 1].t || jd > samples[n].t)
{
typename vector<SampleXYZV<T> >::const_iterator iter = lower_bound(samples.begin(),
samples.end(),
samp);
if (iter == samples.end())
n = samples.size();
else
n = iter - samples.begin();
lastSample = n;
}
if (n > 0 && n < (int) samples.size())
{
SampleXYZV<T> s0 = samples[n - 1];
SampleXYZV<T> s1 = samples[n];
if (interpolation == TrajectoryInterpolationLinear)
{
double h = s1.t - s0.t;
vel = Vec3d(s1.position.x - s0.position.x,
s1.position.y - s0.position.y,
s1.position.z - s0.position.z) * (1.0 / h) * astro::daysToSecs(1.0);
}
else if (interpolation == TrajectoryInterpolationCubic)
{
double h = s1.t - s0.t;
double ih = 1.0 / h;
double t = (jd - s0.t) * ih;
Vec3d p0(s0.position.x, s0.position.y, s0.position.z);
Vec3d p1(s1.position.x, s1.position.y, s1.position.z);
Vec3d v0(s0.velocity.x, s0.velocity.y, s0.velocity.z);
Vec3d v1(s1.velocity.x, s1.velocity.y, s1.velocity.z);
vel = cubicInterpolate(p0, v0 * h, p1, v1 * h, t) * ih;
}
else
{
// Unknown interpolation type
vel = Vec3d(0.0, 0.0, 0.0);
}
}
}
// Add correction for Celestia's coordinate system
return Vec3d(vel.x, vel.z, -vel.y);
}
template <typename T> void SampledOrbitXYZV<T>::sample(double start, double t, int,
OrbitSampleProc& proc) const
{
double cosThresholdAngle = cos(degToRad(SampleThresholdAngle));
double dt = MinSampleInterval;
double end = start + t;
double current = start;
proc.sample(current, positionAtTime(current));
while (current < end)
{
double dt2 = dt;
Point3d goodpt;
double gooddt = dt;
Point3d pos0 = positionAtTime(current);
goodpt = positionAtTime(current + dt2);
while (1)
{
Point3d pos1 = positionAtTime(current + dt2);
Point3d pos2 = positionAtTime(current + dt2 * 2.0);
Vec3d vec1 = pos1 - pos0;
Vec3d vec2 = pos2 - pos0;
vec1.normalize();
vec2.normalize();
double dot = vec1 * vec2;
if (dot > 1.0)
dot = 1.0;
else if (dot < -1.0)
dot = -1.0;
if (dot > cosThresholdAngle && dt2 < MaxSampleInterval)
{
gooddt = dt2;
goodpt = pos1;
dt2 *= 2.0;
}
else
{
proc.sample(current + gooddt, goodpt);
break;
}
}
current += gooddt;
}
}
// Load an ASCII xyz trajectory file. The file contains records with 4 double
// precision values each:
//
// 1: TDB time
// 2: Position x
// 3: Position y
// 4: Position z
//
// Positions are in kilometers.
template <typename T> SampledOrbit<T>* LoadSampledOrbit(const string& filename, TrajectoryInterpolation interpolation, T)
{
ifstream in(filename.c_str());
@ -428,13 +730,84 @@ template <typename T> SampledOrbit<T>* LoadSampledOrbit(const string& filename,
}
// Load an xyzv sampled trajectory file. The file contains records with 7 double
// precision values:
//
// 1: TDB time
// 2: Position x
// 3: Position y
// 4: Position z
// 5: Velocity x
// 6: Velocity y
// 7: Velocity z
//
// Positions are in kilometers, velocities are kilometers per second.
template <typename T> SampledOrbitXYZV<T>* LoadSampledOrbitXYZV(const string& filename, TrajectoryInterpolation interpolation, T)
{
ifstream in(filename.c_str());
if (!in.good())
return NULL;
SampledOrbitXYZV<T>* orbit = NULL;
orbit = new SampledOrbitXYZV<T>(interpolation);
int nSamples = 0;
while (in.good())
{
double tdb = 0.0;
Vec3d position;
Vec3d velocity;
in >> tdb;
in >> position.x;
in >> position.y;
in >> position.z;
in >> velocity.x;
in >> velocity.y;
in >> velocity.z;
// Convert velocities from km/sec to km/Julian day
velocity = velocity * astro::daysToSecs(1.0);
if (in.good())
{
orbit->addSample(tdb, position, velocity);
nSamples++;
}
}
return orbit;
}
/*! Load a trajectory file containing single precision positions.
*/
Orbit* LoadSampledTrajectorySinglePrec(const string& filename, TrajectoryInterpolation interpolation)
{
return LoadSampledOrbit(filename, interpolation, 0.0f);
}
/*! Load a trajectory file containing double precision positions.
*/
Orbit* LoadSampledTrajectoryDoublePrec(const string& filename, TrajectoryInterpolation interpolation)
{
return LoadSampledOrbit(filename, interpolation, 0.0);
}
/*! Load a trajectory file with single precision positions and velocities.
*/
Orbit* LoadXYZVTrajectorySinglePrec(const string& filename, TrajectoryInterpolation interpolation)
{
return LoadSampledOrbitXYZV(filename, interpolation, 0.0f);
}
/*! Load a trajectory file with double precision positions and velocities.
*/
Orbit* LoadXYZVTrajectoryDoublePrec(const string& filename, TrajectoryInterpolation interpolation)
{
return LoadSampledOrbitXYZV(filename, interpolation, 0.0);
}

2
src/celengine/samporbit.h
View File

@ -27,5 +27,7 @@ enum TrajectoryPrecision
extern Orbit* LoadSampledTrajectoryDoublePrec(const std::string& name, TrajectoryInterpolation interpolation);
extern Orbit* LoadSampledTrajectorySinglePrec(const std::string& name, TrajectoryInterpolation interpolation);
extern Orbit* LoadXYZVTrajectoryDoublePrec(const std::string& name, TrajectoryInterpolation interpolation);
extern Orbit* LoadXYZVTrajectorySinglePrec(const std::string& name, TrajectoryInterpolation interpolation);
#endif // _CELENGINE_SAMPORBIT_H_

43
src/celengine/trajmanager.cpp
View File

@ -1,6 +1,6 @@
// trajmanager.cpp
//
// Copyright (C) 2001-2007 Chris Laurel <claurel@gmail.com>
// Copyright (C) 2001-2008 Chris Laurel <claurel@gmail.com>
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
@ -13,6 +13,7 @@
#include "celestia.h"
#include <celutil/debug.h>
#include <celutil/filetype.h>
#include "samporbit.h"
#include "trajmanager.h"
@ -56,21 +57,41 @@ Orbit* TrajectoryInfo::load(const string& filename)
// strip off the uniquifying suffix
string::size_type uniquifyingSuffixStart = filename.rfind(UniqueSuffixChar);
string strippedFilename(filename, 0, uniquifyingSuffixStart);
ContentType filetype = DetermineFileType(strippedFilename);
DPRINTF(1, "Loading trajectory: %s\n", strippedFilename.c_str());
Orbit* sampTrajectory = NULL;
switch (precision)
if (filetype == Content_CelestiaXYZVTrajectory)
{
case TrajectoryPrecisionSingle:
sampTrajectory = LoadSampledTrajectorySinglePrec(strippedFilename, interpolation);
break;
case TrajectoryPrecisionDouble:
sampTrajectory = LoadSampledTrajectoryDoublePrec(strippedFilename, interpolation);
break;
default:
assert(0);
break;
switch (precision)
{
case TrajectoryPrecisionSingle:
sampTrajectory = LoadXYZVTrajectorySinglePrec(strippedFilename, interpolation);
break;
case TrajectoryPrecisionDouble:
sampTrajectory = LoadXYZVTrajectoryDoublePrec(strippedFilename, interpolation);
break;
default:
assert(0);
break;
}
}
else
{
switch (precision)
{
case TrajectoryPrecisionSingle:
sampTrajectory = LoadSampledTrajectorySinglePrec(strippedFilename, interpolation);
break;
case TrajectoryPrecisionDouble:
sampTrajectory = LoadSampledTrajectoryDoublePrec(strippedFilename, interpolation);
break;
default:
assert(0);
break;
}
}
return sampTrajectory;