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

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,36 +153,51 @@ CreateEllipticalOrbit(const Hash* orbitData,
     double timeScale;
     GetDefaultUnits(usePlanetUnits, distanceScale, timeScale);
 
+    auto eccentricity = orbitData->getNumber<double>("Eccentricity").value_or(0.0);
+    if (eccentricity < 0.0)
+    {
+        GetLogger()->error("Negative eccentricity is invalid.\n");
+        return nullptr;
+    }
+    else if (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())
+    double semiMajorAxis = 0.0;
+    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;
-        }
+        semiMajorAxis = *semiMajorAxisValue;
+    }
+    else if (auto pericenter = orbitData->getLength<double>("PericenterDistance", 1.0, distanceScale); pericenter.has_value())
+    {
+        semiMajorAxis = *pericenter / (1.0 - 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;
+        if (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);
-
     auto inclination = orbitData->getAngle<double>("Inclination").value_or(0.0);
 
     double ascendingNode = orbitData->getAngle<double>("AscendingNode").value_or(0.0);
@@ -204,20 +219,23 @@ CreateEllipticalOrbit(const Hash* orbitData,
     // if both are specified.
     double anomalyAtEpoch = 0.0;
     if (auto meanAnomaly = orbitData->getAngle<double>("MeanAnomaly"); meanAnomaly.has_value())
-    {
         anomalyAtEpoch = *meanAnomaly;
-    }
     else if (auto meanLongitude = orbitData->getAngle<double>("MeanLongitude"); meanLongitude.has_value())
-    {
         anomalyAtEpoch = *meanLongitude - (argOfPericenter + ascendingNode);
-    }
 
-    // If we read the semi-major axis, use it to compute the pericenter
-    // distance.
-    if (semiMajorAxis.has_value())
-        pericenterDistance = *semiMajorAxis * (1.0 - eccentricity);
+    if (eccentricity < 1.0)
+    {
+        return std::make_unique<celestia::ephem::EllipticalOrbit>(semiMajorAxis,
+                                                                  eccentricity,
+                                                                  degToRad(inclination),
+                                                                  degToRad(ascendingNode),
+                                                                  degToRad(argOfPericenter),
+                                                                  degToRad(anomalyAtEpoch),
+                                                                  period,
+                                                                  epoch);
+    }
 
-    return std::make_unique<celestia::ephem::EllipticalOrbit>(pericenterDistance,
+    return std::make_unique<celestia::ephem::HyperbolicOrbit>(semiMajorAxis,
                                                               eccentricity,
                                                               degToRad(inclination),
                                                               degToRad(ascendingNode),
@@ -310,7 +328,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 +788,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);