Index: orbit.cpp
===================================================================
--- orbit.cpp	(revision 4190)
+++ orbit.cpp	(working copy)
@@ -42,6 +42,9 @@
     period(_period),
     epoch(_epoch)
 {
+    orbitPlaneRotation = (Mat3d::zrotation(_ascendingNode) *
+                          Mat3d::xrotation(_inclination) *
+                          Mat3d::zrotation(_argOfPeriapsis));
 }
 
 
@@ -178,35 +181,64 @@
 
 Point3d EllipticalOrbit::positionAtE(double E) const
 {
-    double x, z;
+    double x, y;
 
     if (eccentricity < 1.0)
     {
         double a = pericenterDistance / (1.0 - eccentricity);
         x = a * (cos(E) - eccentricity);
-        z = a * sqrt(1 - square(eccentricity)) * -sin(E);
+        y = a * sqrt(1 - square(eccentricity)) * sin(E);
     }
     else if (eccentricity > 1.0)
     {
         double a = pericenterDistance / (1.0 - eccentricity);
         x = -a * (eccentricity - cosh(E));
-        z = -a * sqrt(square(eccentricity) - 1) * -sinh(E);
+        y = -a * sqrt(square(eccentricity) - 1) * sinh(E);
     }
     else
     {
         // TODO: Handle parabolic orbits
         x = 0.0;
-        z = 0.0;
+        y = 0.0;
     }
 
-    Mat3d R = (Mat3d::yrotation(ascendingNode) *
-               Mat3d::xrotation(inclination) *
-               Mat3d::yrotation(argOfPeriapsis));
+    Point3d p = orbitPlaneRotation * Point3d(x, y, 0);
 
-    return R * Point3d(x, 0, z);
+    // Convert to Celestia's internal coordinate system
+    return Point3d(p.x, p.z, -p.y);
 }
 
 
+Vec3d EllipticalOrbit::velocityAtE(double E) const
+{
+    double x, y;
+
+    if (eccentricity < 1.0)
+    {
+        double a = pericenterDistance / (1.0 - eccentricity);
+        x = a * (-sin(E) - eccentricity);
+        y = a * sqrt(1 - square(eccentricity)) * cos(E);
+    }
+    else if (eccentricity > 1.0)
+    {
+        double a = pericenterDistance / (1.0 - eccentricity);
+        x = -a * (eccentricity - cosh(E));
+        y = -a * sqrt(square(eccentricity) - 1) * sinh(E);
+    }
+    else
+    {
+        // TODO: Handle parabolic orbits
+        x = 0.0;
+        y = 0.0;
+    }
+
+    Vec3d v = orbitPlaneRotation * Vec3d(x, y, 0);
+
+    // Convert to Celestia's coordinate system
+    return Vec3d(v.x, v.z, -v.y);
+}
+
+
 // Return the offset from the center
 Point3d EllipticalOrbit::positionAtTime(double t) const
 {
@@ -219,6 +251,17 @@
 }
 
 
+Vec3d EllipticalOrbit::velocityAtTime(double t) const
+{
+    t = t - epoch;
+    double meanMotion = 2.0 * PI / period;
+    double meanAnomaly = meanAnomalyAtEpoch + t * meanMotion;
+    double E = eccentricAnomaly(meanAnomaly);
+
+    return velocityAtE(E);
+}
+
+
 double EllipticalOrbit::getPeriod() const
 {
     return period;
@@ -405,16 +448,6 @@
     double T = 2 * PI * sqrt(cube(a) / GM);
     T = T / 86400.0; // Convert from seconds to days
 
-#if 0
-    cout << "a: " << astro::kilometersToAU(a) << '\n';
-    cout << "e: " << e << '\n';
-    cout << "i: " << radToDeg(i) << '\n';
-    cout << "Om: " << radToDeg(Om) << '\n';
-    cout << "om: " << radToDeg(om) << '\n';
-    cout << "M: " << radToDeg(M) << '\n';
-    cout << "T: " << T << '\n';
-#endif
-
     return new EllipticalOrbit(a * (1 - e), e, i, Om, om, M, T, t);
 }
 
Index: orbit.h
===================================================================
--- orbit.h	(revision 4190)
+++ orbit.h	(working copy)
@@ -12,7 +12,6 @@
 
 #include <celmath/vecmath.h>
 
-
 class OrbitSampleProc;
 
 class Orbit
@@ -54,7 +53,7 @@
 
     // Compute the orbit for a specified Julian date
     virtual Point3d positionAtTime(double) const;
-    //virtual Vec3d velocityAtTime(double) const;
+    virtual Vec3d velocityAtTime(double) const;
     double getPeriod() const;
     double getBoundingRadius() const;
     virtual void sample(double, double, int, OrbitSampleProc&) const;
@@ -62,6 +61,7 @@
  private:
     double eccentricAnomaly(double) const;
     Point3d positionAtE(double) const;
+    Vec3d velocityAtE(double) const;
 
     double pericenterDistance;
     double eccentricity;
@@ -71,6 +71,8 @@
     double meanAnomalyAtEpoch;
     double period;
     double epoch;
+
+    Mat3d orbitPlaneRotation;
 };