Split EllipticalOrbit and HyperbolicOrbit classes

- Store semimajor axis instead of pericenterDistance in EllipticalOrbit
- Store semiminor axis to avoid some sqrt calls during evaluation
- Support zero and negative mean anomaly in hyperbolic orbits
- Fix velocity calculation for hyperbolic orbits
- Consolidate state vector to orbital elements calculations
- Extend orbital elements calculation to support hyperbolic orbits
- Add tests for orbital elements calculations
- Fix hyperbolic orbital elements display in Qt info panel

Author: ajtribick

src/celengine/astro.cpp

@@ -8,6 +8,7 @@
 // as published by the Free Software Foundation; either version 2
 // of the License, or (at your option) any later version.
 
+#include <algorithm>
 #include <cmath>
 #include <ostream>
 #include <utility>
@@ -91,6 +92,13 @@ static const astro::LeapSecondRecord LeapSeconds[] =
 };
 
 
+static inline void negateIf(double& d, bool condition)
+{
+    if (condition)
+        d = -d;
+}
+
+
 static celestia::util::array_view<astro::LeapSecondRecord> g_leapSeconds = LeapSeconds;
 
 
@@ -765,3 +773,95 @@ void astro::setLeapSeconds(celestia::util::array_view<astro::LeapSecondRecord> l
 {
     g_leapSeconds = leapSeconds;
 }
+
+
+astro::KeplerElements astro::StateVectorToElements(const Eigen::Vector3d& r,
+                                                   const Eigen::Vector3d& v,
+                                                   double mu)
+{
+    constexpr double tolerance = 1e-9;
+
+    Eigen::Vector3d h = r.cross(v);
+    double rNorm = r.norm();
+
+    KeplerElements result;
+
+    // Compute eccentricity
+    Eigen::Vector3d evec = v.cross(h) / mu - r / rNorm;
+    result.eccentricity = evec.norm();
+
+    // Compute inclination
+    result.inclination = std::acos(std::clamp(h.y() / h.norm(), -1.0, 1.0));
+
+    // Normal vector (UnitY x h)
+    Eigen::Vector3d nvec(h[2], 0, -h[0]);
+    double nNorm = nvec.norm();
+
+    // compute longAscendingNode and argPericenter
+    if (result.inclination < tolerance)
+    {
+        // handle face-on orbit: by convention Omega = 0.0
+        if (result.eccentricity >= tolerance)
+        {
+            result.argPericenter = std::acos(evec.x() / result.eccentricity);
+            negateIf(result.argPericenter, evec.z() >= 0.0);
+        }
+    }
+    else
+    {
+        result.longAscendingNode = std::acos(nvec.x() / nNorm);
+        negateIf(result.longAscendingNode, nvec.z() >= 0.0);
+        if (result.eccentricity >= tolerance)
+        {
+            result.argPericenter = std::acos(std::clamp(nvec.dot(evec) / (nNorm * result.eccentricity), -1.0, 1.0));
+            negateIf(result.argPericenter, evec.y() < 0.0);
+        }
+    }
+
+    // compute true anomaly
+    double nu;
+    if (result.eccentricity >= tolerance)
+    {
+        nu = std::acos(std::clamp(evec.dot(r) / (result.eccentricity * rNorm), -1.0, 1.0));
+        negateIf(nu, r.dot(v) < 0.0);
+    }
+    else
+    {
+        if (result.inclination < tolerance)
+        {
+            // circular face-on orbit
+            nu = std::acos(r.x() / rNorm);
+            negateIf(nu, v.x() > 0.0);
+        }
+        else
+        {
+            nu = std::acos(std::clamp(nvec.dot(r) / (nNorm * rNorm), -1.0, 1.0));
+            negateIf(nu, nvec.dot(v) > 0.0);
+        }
+    }
+
+    double s_nu;
+    double c_nu;
+    celmath::sincos(nu, s_nu, c_nu);
+
+    // compute mean anomaly
+    double e2 = celmath::square(result.eccentricity);
+    if (result.eccentricity < 1.0)
+    {
+        double E = std::atan2(std::sqrt(1.0 - e2) * s_nu,
+                              result.eccentricity + c_nu);
+        result.meanAnomaly = E - result.eccentricity * std::sin(E);
+    }
+    else
+    {
+        double sinhE = std::sqrt(e2 - 1.0) * s_nu / (1.0 + result.eccentricity * c_nu);
+        double E = std::asinh(sinhE);
+        result.meanAnomaly = result.eccentricity * sinhE - E;
+    }
+
+    // compute semimajor axis
+    result.semimajorAxis = 1.0 / (2.0 / rNorm - v.squaredNorm() / mu);
+    result.period = 2.0 * celestia::numbers::pi * std::sqrt(celmath::cube(std::abs(result.semimajorAxis)) / mu);
+
+    return result;
+}

src/celengine/astro.h

@@ -320,4 +320,19 @@ namespace astro
         return astro::speedOfLight * n;
     }
     }
-}
+
+    struct KeplerElements
+    {
+        double semimajorAxis{ 0.0 };
+        double eccentricity{ 0.0 };
+        double inclination{ 0.0 };
+        double longAscendingNode{ 0.0 };
+        double argPericenter{ 0.0 };
+        double meanAnomaly{ 0.0 };
+        double period{ 0.0 };
+    };
+
+    KeplerElements StateVectorToElements(const Eigen::Vector3d&,
+                                         const Eigen::Vector3d&,
+                                         double);
+} // end namespace astro

src/celengine/parseobject.cpp

@@ -143,8 +143,8 @@ ParseDate(const Hash* hash, const string& name, double& jd)
  *     SemiMajorAxis or PericenterDistance is in kilometers.
  */
 static std::unique_ptr<celestia::ephem::Orbit>
-CreateEllipticalOrbit(const Hash* orbitData,
-                      bool usePlanetUnits)
+CreateKeplerianOrbit(const Hash* orbitData,
+                     bool usePlanetUnits)
 {
 
     // default units for planets are AU and years, otherwise km and days
@@ -153,78 +153,85 @@ CreateEllipticalOrbit(const Hash* orbitData,
     double timeScale;
     GetDefaultUnits(usePlanetUnits, distanceScale, timeScale);
 
+    astro::KeplerElements elements;
+
+    elements.eccentricity = orbitData->getNumber<double>("Eccentricity").value_or(0.0);
+    if (elements.eccentricity < 0.0)
+    {
+        GetLogger()->error("Negative eccentricity is invalid.\n");
+        return nullptr;
+    }
+    else if (elements.eccentricity == 1.0)
+    {
+        GetLogger()->error("Parabolic orbits are not supported.\n");
+        return nullptr;
+    }
+
     // SemiMajorAxis and Period are absolutely required; everything
     // else has a reasonable default.
-    double pericenterDistance = 0.0;
-    std::optional<double> semiMajorAxis = orbitData->getLength<double>("SemiMajorAxis", 1.0, distanceScale);
-    if (!semiMajorAxis.has_value())
+    if (auto semiMajorAxisValue = orbitData->getLength<double>("SemiMajorAxis", 1.0, distanceScale); semiMajorAxisValue.has_value())
     {
-        if (auto pericenter = orbitData->getLength<double>("PericenterDistance", 1.0, distanceScale); pericenter.has_value())
-        {
-            pericenterDistance = *pericenter;
-        }
-        else
-        {
-            GetLogger()->error("SemiMajorAxis/PericenterDistance missing!  Skipping planet . . .\n");
-            return nullptr;
-        }
+        elements.semimajorAxis = *semiMajorAxisValue;
+    }
+    else if (auto pericenter = orbitData->getLength<double>("PericenterDistance", 1.0, distanceScale); pericenter.has_value())
+    {
+        elements.semimajorAxis = *pericenter / (1.0 - elements.eccentricity);
+    }
+    else
+    {
+        GetLogger()->error("SemiMajorAxis/PericenterDistance missing from orbit definition.\n");
+        return nullptr;
     }
 
-    double period = 0.0;
     if (auto periodValue = orbitData->getTime<double>("Period", 1.0, timeScale); periodValue.has_value())
     {
-        period = *periodValue;
+        elements.period = *periodValue;
+        if (elements.period == 0.0)
+        {
+            GetLogger()->error("Period cannot be zero.\n");
+            return nullptr;
+        }
     }
     else
     {
-        GetLogger()->error("Period missing!  Skipping planet . . .\n");
+        GetLogger()->error("Period must be specified in EllipticalOrbit.\n");
         return nullptr;
     }
 
-    auto eccentricity = orbitData->getNumber<double>("Eccentricity").value_or(0.0);
+    elements.inclination = orbitData->getAngle<double>("Inclination").value_or(0.0);
 
-    auto inclination = orbitData->getAngle<double>("Inclination").value_or(0.0);
+    elements.longAscendingNode = orbitData->getAngle<double>("AscendingNode").value_or(0.0);
 
-    double ascendingNode = orbitData->getAngle<double>("AscendingNode").value_or(0.0);
-
-    double argOfPericenter = 0.0;
     if (auto argPeri = orbitData->getAngle<double>("ArgOfPericenter"); argPeri.has_value())
     {
-        argOfPericenter = *argPeri;
+        elements.argPericenter = *argPeri;
     }
     else if (auto longPeri = orbitData->getAngle<double>("LongOfPericenter"); longPeri.has_value())
     {
-        argOfPericenter = *longPeri - ascendingNode;
+        elements.argPericenter = *longPeri - elements.longAscendingNode;
     }
 
     double epoch = astro::J2000;
     ParseDate(orbitData, "Epoch", epoch);
 
     // Accept either the mean anomaly or mean longitude--use mean anomaly
     // if both are specified.
-    double anomalyAtEpoch = 0.0;
     if (auto meanAnomaly = orbitData->getAngle<double>("MeanAnomaly"); meanAnomaly.has_value())
-    {
-        anomalyAtEpoch = *meanAnomaly;
-    }
+        elements.meanAnomaly = *meanAnomaly;
     else if (auto meanLongitude = orbitData->getAngle<double>("MeanLongitude"); meanLongitude.has_value())
+        elements.meanAnomaly = *meanLongitude - (elements.argPericenter + elements.longAscendingNode);
+
+    elements.inclination = celmath::degToRad(elements.inclination);
+    elements.longAscendingNode = celmath::degToRad(elements.longAscendingNode);
+    elements.argPericenter = celmath::degToRad(elements.argPericenter);
+    elements.meanAnomaly = celmath::degToRad(elements.meanAnomaly);
+
+    if (elements.eccentricity < 1.0)
     {
-        anomalyAtEpoch = *meanLongitude - (argOfPericenter + ascendingNode);
+        return std::make_unique<celestia::ephem::EllipticalOrbit>(elements, epoch);
     }
 
-    // If we read the semi-major axis, use it to compute the pericenter
-    // distance.
-    if (semiMajorAxis.has_value())
-        pericenterDistance = *semiMajorAxis * (1.0 - eccentricity);
-
-    return std::make_unique<celestia::ephem::EllipticalOrbit>(pericenterDistance,
-                                                              eccentricity,
-                                                              degToRad(inclination),
-                                                              degToRad(ascendingNode),
-                                                              degToRad(argOfPericenter),
-                                                              degToRad(anomalyAtEpoch),
-                                                              period,
-                                                              epoch);
+    return std::make_unique<celestia::ephem::HyperbolicOrbit>(elements, epoch);
 }
 
 
@@ -310,7 +317,7 @@ CreateFixedPosition(const Hash* trajData, const Selection& centralObject, bool u
     {
         if (centralObject.getType() != SelectionType::Body)
         {
-            GetLogger()->error("FixedPosition planetographic coordinates aren't valid for stars.\n");
+            GetLogger()->error("FixedPosition planetographic coordinates are not valid for stars.\n");
             return nullptr;
         }
 
@@ -770,7 +777,7 @@ CreateOrbit(const Selection& centralObject,
             return nullptr;
         }
 
-        return CreateEllipticalOrbit(orbitData, usePlanetUnits).release();
+        return CreateKeplerianOrbit(orbitData, usePlanetUnits).release();
     }
 
     // Create an 'orbit' that places the object at a fixed point in its

src/celengine/render.cpp

@@ -1090,42 +1090,24 @@ void Renderer::renderOrbit(const OrbitPathListEntry& orbitPath,
     if (cachedOrbit == nullptr)
     {
         double startTime = t;
-#if 0
-        int nSamples = detailOptions.orbitPathSamplePoints;
-#endif
 
         // Adjust the number of samples used for aperiodic orbits--these aren't
         // true orbits, but are sampled trajectories, generally of spacecraft.
         // Better control is really needed--some sort of adaptive sampling would
         // be ideal.
-        if (!orbit->isPeriodic())
+        if (orbit->isPeriodic())
+        {
+            startTime = t - orbit->getPeriod();
+        }
+        else
         {
             double begin = 0.0, end = 0.0;
             orbit->getValidRange(begin, end);
 
             if (begin != end)
             {
                 startTime = begin;
-#if 0
-                nSamples = (int) (orbit->getPeriod() * 100.0);
-                nSamples = std::clamp(nSamples, 100, 1000);
-#endif
             }
-#if 0
-            else
-            {
-                // If the orbit is aperiodic and doesn't have a
-                // finite duration, we don't render it. A compromise
-                // would be to pick some time window centered at the
-                // current time, but we'd have to pick some arbitrary
-                // duration.
-                nSamples = 0;
-            }
-#endif
-        }
-        else
-        {
-            startTime = t - orbit->getPeriod();
         }
 
         cachedOrbit = new CurvePlot(*this);