Index: samporbit.cpp
===================================================================
--- samporbit.cpp	(revision 4233)
+++ samporbit.cpp	(working copy)
@@ -30,17 +30,33 @@
 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 @@
 }
 
 
-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 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 @@
 
 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 @@
 
         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 @@
                 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 @@
                 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 @@
             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 @@
                 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 @@
 }
 
 
+// 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 @@
 }
 
 
+// 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);
+}
Index: samporbit.h
===================================================================
--- samporbit.h	(revision 4232)
+++ samporbit.h	(working copy)
@@ -27,5 +27,7 @@
 
 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_
Index: trajmanager.cpp
===================================================================
--- trajmanager.cpp	(revision 4232)
+++ trajmanager.cpp	(working copy)
@@ -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,22 +57,42 @@
     // 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;
 }