Implement orbital precession

src/celastro/astro.h

@@ -218,6 +218,8 @@ struct KeplerElements
     double argPericenter{ 0.0 };
     double meanAnomaly{ 0.0 };
     double period{ 0.0 };
+    double nodalPeriod{ 0.0 };
+    double apsidalPeriod{ 0.0 };
 };
 
 

src/celengine/parseobject.cpp

@@ -133,6 +133,11 @@ GetDefaultUnits(bool usePlanetUnits, double& distanceScale)
  *     # One or none of the following:
  *     MeanAnomaly <degrees>     (default: 0.0)
  *     MeanLongitude <degrees>   (default: 0.0)
+ * 
+ *     # For precessing orbits (default unit is Julian years):
+ *     # Default values result in no precession
+ *     NodalPeriod <number>   (default: 0.0)
+ *     ApsidalPeriod <number> (default: 0.0)
  * } \endcode
  *
  * If usePlanetUnits is true:
@@ -198,6 +203,10 @@ CreateKeplerianOrbit(const AssociativeArray* orbitData,
         return nullptr;
     }
 
+    elements.nodalPeriod = orbitData->getTime<double>("NodalPeriod", 1.0, astro::DAYS_PER_YEAR).value_or(0.0);
+
+    elements.apsidalPeriod = orbitData->getTime<double>("ApsidalPeriod", 1.0, astro::DAYS_PER_YEAR).value_or(0.0);
+
     elements.inclination = orbitData->getAngle<double>("Inclination").value_or(0.0);
 
     elements.longAscendingNode = orbitData->getAngle<double>("AscendingNode").value_or(0.0);
@@ -228,7 +237,12 @@ CreateKeplerianOrbit(const AssociativeArray* orbitData,
 
     if (elements.eccentricity < 1.0)
     {
-        return std::make_shared<ephem::EllipticalOrbit>(elements, epoch);
+        if (elements.nodalPeriod == 0.0 && elements.apsidalPeriod == 0.0)
+        {
+            return std::make_shared<ephem::EllipticalOrbit>(elements, epoch);
+        }
+
+        return std::make_shared<ephem::PrecessingOrbit>(elements, epoch);
     }
 
     return std::make_shared<ephem::HyperbolicOrbit>(elements, epoch);

src/celephem/orbit.cpp

@@ -127,7 +127,12 @@ std::unique_ptr<Orbit>
 CreateKeplerOrbit(const astro::KeplerElements& elements, double epoch)
 {
     if (elements.eccentricity < 1.0)
-        return std::make_unique<EllipticalOrbit>(elements, epoch);
+    {
+        if (elements.nodalPeriod == 0.0 && elements.apsidalPeriod == 0.0)
+            return std::make_unique<EllipticalOrbit>(elements, epoch);
+
+        return std::make_unique<PrecessingOrbit>(elements, epoch);
+    }
 
     return std::make_unique<HyperbolicOrbit>(elements, epoch);
 }
@@ -383,6 +388,138 @@ double EllipticalOrbit::getBoundingRadius() const
 }
 
 
+PrecessingOrbit::PrecessingOrbit(const astro::KeplerElements& _elements, double _epoch) :
+    semiMajorAxis(_elements.semimajorAxis),
+    eccentricity(_elements.eccentricity),
+    inclination(_elements.inclination),
+    longAscendingNodeAtEpoch(_elements.longAscendingNode),
+    argPericenterAtEpoch(_elements.argPericenter),
+    meanAnomalyAtEpoch(_elements.meanAnomaly),
+    period(_elements.period),
+    nodalPeriod(_elements.nodalPeriod),
+    apsidalPeriod(_elements.apsidalPeriod),
+    epoch(_epoch)
+{
+    assert(eccentricity >= 0.0 && eccentricity < 1.0);
+    assert(semiMajorAxis >= 0.0);
+    assert(period != 0.0);
+    semiMinorAxis = semiMajorAxis * std::sqrt(1.0 - math::square(eccentricity));
+}
+
+
+double PrecessingOrbit::eccentricAnomaly(double M) const
+{
+    if (eccentricity == 0.0)
+    {
+        // Circular orbit
+        return M;
+    }
+    if (eccentricity < 0.2)
+    {
+        // Low eccentricity, so use the standard iteration technique
+        return math::solve_iteration_fixed(SolveKeplerFunc1(eccentricity, M), M, 5).first;
+    }
+    if (eccentricity < 0.9)
+    {
+        // Higher eccentricity elliptical orbit; use a more complex but
+        // much faster converging iteration.
+        return math::solve_iteration_fixed(SolveKeplerFunc2(eccentricity, M), M, 6).first;
+    }
+
+    // Extremely stable Laguerre-Conway method for solving Kepler's
+    // equation.  Only use this for high-eccentricity orbits, as it
+    // requires more calcuation.
+    double E = M + 0.85 * eccentricity * math::sign(std::sin(M));
+    return math::solve_iteration_fixed(SolveKeplerLaguerreConway(eccentricity, M), E, 8).first;
+}
+
+
+// Compute the position at the specified eccentric
+// anomaly E.
+Eigen::Vector3d PrecessingOrbit::positionAtE(double E, double longAscendingNode, double argPericenter) const
+{
+    double x = semiMajorAxis * (std::cos(E) - eccentricity);
+    double y = semiMinorAxis * std::sin(E);
+
+    Eigen::Matrix3d orbitPlaneRotation = (math::ZRotation(longAscendingNode) *
+                                          math::XRotation(inclination) *
+                                          math::ZRotation(argPericenter)).toRotationMatrix();
+    Eigen::Vector3d p = orbitPlaneRotation * Eigen::Vector3d(x, y, 0);
+
+    // Convert to Celestia's internal coordinate system
+    return Eigen::Vector3d(p.x(), p.z(), -p.y());
+}
+
+
+// Compute the velocity at the specified eccentric
+// anomaly E.
+Eigen::Vector3d PrecessingOrbit::velocityAtE(double E, double longAscendingNode, double argPericenter, double meanMotion) const
+{
+    double sinE;
+    double cosE;
+    math::sincos(E, sinE, cosE);
+
+    double edot = meanMotion / (1.0 - eccentricity * cosE);
+
+    double x = -semiMajorAxis * sinE * edot;
+    double y =  semiMinorAxis * cosE * edot;
+
+    Eigen::Matrix3d orbitPlaneRotation = (math::ZRotation(longAscendingNode) *
+                                          math::XRotation(inclination) *
+                                          math::ZRotation(argPericenter)).toRotationMatrix();
+    Eigen::Vector3d v = orbitPlaneRotation* Eigen::Vector3d(x, y, 0);
+
+    // Convert to Celestia's coordinate system
+    return Eigen::Vector3d(v.x(), v.z(), -v.y());
+}
+
+
+// Return the offset from the center
+Eigen::Vector3d PrecessingOrbit::positionAtTime(double t) const
+{
+    t = t - epoch;
+    double meanMotion = 2.0 * celestia::numbers::pi / period;
+    double nodalPrecessionRate = nodalPeriod != 0.0 ? 2.0 * celestia::numbers::pi / nodalPeriod : 0.0;
+    double apsidalPrecessionRate = apsidalPeriod != 0.0 ? 2.0 * celestia::numbers::pi / apsidalPeriod : 0.0;
+
+    double meanAnomaly = meanAnomalyAtEpoch + t * (meanMotion - apsidalPrecessionRate - nodalPrecessionRate);
+    double E = eccentricAnomaly(meanAnomaly);
+    double longAscendingNode = longAscendingNodeAtEpoch + t * nodalPrecessionRate;
+    double argPericenter = argPericenterAtEpoch + t * (apsidalPrecessionRate - nodalPrecessionRate);
+
+    return positionAtE(E, longAscendingNode, argPericenter);
+}
+
+
+Eigen::Vector3d PrecessingOrbit::velocityAtTime(double t) const
+{
+    t = t - epoch;
+    double meanMotion = 2.0 * celestia::numbers::pi / period;
+    double nodalPrecessionRate = nodalPeriod != 0.0 ? 2.0 * celestia::numbers::pi / nodalPeriod : 0.0;
+    double apsidalPrecessionRate = apsidalPeriod != 0.0 ? 2.0 * celestia::numbers::pi / apsidalPeriod : 0.0;
+
+    double meanAnomaly = meanAnomalyAtEpoch + t * (meanMotion - apsidalPrecessionRate - nodalPrecessionRate);
+    double E = eccentricAnomaly(meanAnomaly);
+    double longAscendingNode = longAscendingNodeAtEpoch + t * nodalPrecessionRate;
+    double argPericenter = argPericenterAtEpoch + t * (apsidalPrecessionRate - nodalPrecessionRate);
+
+    return velocityAtE(E, longAscendingNode, argPericenter, meanMotion);
+}
+
+
+double PrecessingOrbit::getPeriod() const
+{
+    return period;
+}
+
+
+double PrecessingOrbit::getBoundingRadius() const
+{
+    // TODO: watch out for unbounded parabolic and hyperbolic orbits
+    return semiMajorAxis * (1.0 + eccentricity);
+}
+
+
 HyperbolicOrbit::HyperbolicOrbit(const astro::KeplerElements& _elements, double _epoch) :
     semiMajorAxis(_elements.semimajorAxis),
     eccentricity(_elements.eccentricity),

src/celephem/orbit.h

@@ -95,6 +95,37 @@ class EllipticalOrbit final : public Orbit
 };
 
 
+class PrecessingOrbit final : public Orbit
+{
+public:
+    PrecessingOrbit(const astro::KeplerElements&, double _epoch = 2451545.0);
+    ~PrecessingOrbit() override = default;
+
+    // Compute the orbit for a specified Julian date
+    Eigen::Vector3d positionAtTime(double) const override;
+    Eigen::Vector3d velocityAtTime(double) const override;
+    double getPeriod() const override;
+    double getBoundingRadius() const override;
+
+private:
+    double eccentricAnomaly(double) const;
+    Eigen::Vector3d positionAtE(double, double, double) const;
+    Eigen::Vector3d velocityAtE(double, double, double, double) const;
+
+    double semiMajorAxis;
+    double semiMinorAxis;
+    double eccentricity;
+    double inclination;
+    double longAscendingNodeAtEpoch;
+    double argPericenterAtEpoch;
+    double meanAnomalyAtEpoch;
+    double period;
+    double nodalPeriod;
+    double apsidalPeriod;
+    double epoch;
+};
+
+
 class HyperbolicOrbit final : public Orbit
 {
 public:

test/unit/kepler_test.cpp

@@ -99,6 +99,62 @@ TEST_CASE("Elliptical orbits")
     }
 }
 
+TEST_CASE("Precessing orbits")
+{
+    constexpr std::array testEccentricities{ 0.0, 0.2, 0.6 };
+
+    for (double period : testPeriods)
+    {
+        double semimajor = std::cbrt(GMsun * math::square(period) / fourpi2);
+        for (double meanAnomalyDeg : testAngles)
+        for (double inclinationDeg : testInclinations)
+        {
+            auto testNodes = inclinationDeg == 0.0 || inclinationDeg == 180.0
+                ? celestia::util::array_view<double>(fixedZero)
+                : celestia::util::array_view<double>(testAngles);
+            for (double nodeDeg : testNodes)
+            {
+                auto testNodalPeriods = nodeDeg == 0.0
+                    ? celestia::util::array_view<double>(fixedZero)
+                    : celestia::util::array_view<double>(testPeriods);
+                for (double nodalPeriod : testNodalPeriods)
+                for (double eccentricity : testEccentricities)
+                {
+                    auto testPericenters = eccentricity == 0.0
+                        ? celestia::util::array_view<double>(fixedZero)
+                        : celestia::util::array_view<double>(testAngles);
+                    for (double pericenterDeg : testPericenters)
+                    {
+                        auto testApsidalPeriods = pericenterDeg == 0.0
+                            ? celestia::util::array_view<double>(fixedZero)
+                            : celestia::util::array_view<double>(testPeriods);
+                        for (double apsidalPeriod : testApsidalPeriods)
+                        {
+                            astro::KeplerElements expected;
+                            expected.period = period;
+                            expected.semimajorAxis = semimajor;
+                            expected.eccentricity = eccentricity;
+                            expected.inclination = math::degToRad(inclinationDeg);
+                            expected.longAscendingNode = math::degToRad(nodeDeg);
+                            expected.argPericenter = math::degToRad(pericenterDeg);
+                            expected.meanAnomaly = math::degToRad(meanAnomalyDeg);
+                            expected.nodalPeriod = nodalPeriod;
+                            expected.apsidalPeriod = apsidalPeriod;
+                            auto orbit = celestia::ephem::PrecessingOrbit(expected, 0.0);
+                            auto position = orbit.positionAtTime(0.0);
+                            auto velocity = orbit.velocityAtTime(0.0);
+
+                            auto actual = astro::StateVectorToElements(position, velocity, GMsun);
+
+                            TestElements(expected, actual);
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
 TEST_CASE("Hyperbolic orbits")
 {
     constexpr std::array testEccentricities{ 1.5, 2.4 };