From 912cf9af06d731ef86dc82e971a77b1e1ccf03d8 Mon Sep 17 00:00:00 2001 From: "Pedro J. Garcia" Date: Wed, 16 Jul 2025 23:27:49 -0300 Subject: [PATCH] Implement orbital precession --- src/celastro/astro.h | 2 + src/celengine/parseobject.cpp | 16 +++- src/celephem/orbit.cpp | 139 +++++++++++++++++++++++++++++++++- src/celephem/orbit.h | 31 ++++++++ test/unit/kepler_test.cpp | 56 ++++++++++++++ 5 files changed, 242 insertions(+), 2 deletions(-) diff --git a/src/celastro/astro.h b/src/celastro/astro.h index e504b473ba..d236cd0e35 100644 --- a/src/celastro/astro.h +++ b/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 }; }; diff --git a/src/celengine/parseobject.cpp b/src/celengine/parseobject.cpp index 4b1b79db9c..17ebb01101 100644 --- a/src/celengine/parseobject.cpp +++ b/src/celengine/parseobject.cpp @@ -133,6 +133,11 @@ GetDefaultUnits(bool usePlanetUnits, double& distanceScale) * # One or none of the following: * MeanAnomaly (default: 0.0) * MeanLongitude (default: 0.0) + * + * # For precessing orbits (default unit is Julian years): + * # Default values result in no precession + * NodalPeriod (default: 0.0) + * ApsidalPeriod (default: 0.0) * } \endcode * * If usePlanetUnits is true: @@ -198,6 +203,10 @@ CreateKeplerianOrbit(const AssociativeArray* orbitData, return nullptr; } + elements.nodalPeriod = orbitData->getTime("NodalPeriod", 1.0, astro::DAYS_PER_YEAR).value_or(0.0); + + elements.apsidalPeriod = orbitData->getTime("ApsidalPeriod", 1.0, astro::DAYS_PER_YEAR).value_or(0.0); + elements.inclination = orbitData->getAngle("Inclination").value_or(0.0); elements.longAscendingNode = orbitData->getAngle("AscendingNode").value_or(0.0); @@ -228,7 +237,12 @@ CreateKeplerianOrbit(const AssociativeArray* orbitData, if (elements.eccentricity < 1.0) { - return std::make_shared(elements, epoch); + if (elements.nodalPeriod == 0.0 && elements.apsidalPeriod == 0.0) + { + return std::make_shared(elements, epoch); + } + + return std::make_shared(elements, epoch); } return std::make_shared(elements, epoch); diff --git a/src/celephem/orbit.cpp b/src/celephem/orbit.cpp index 1e2e2ecbed..8a1246a2ba 100644 --- a/src/celephem/orbit.cpp +++ b/src/celephem/orbit.cpp @@ -127,7 +127,12 @@ std::unique_ptr CreateKeplerOrbit(const astro::KeplerElements& elements, double epoch) { if (elements.eccentricity < 1.0) - return std::make_unique(elements, epoch); + { + if (elements.nodalPeriod == 0.0 && elements.apsidalPeriod == 0.0) + return std::make_unique(elements, epoch); + + return std::make_unique(elements, epoch); + } return std::make_unique(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), diff --git a/src/celephem/orbit.h b/src/celephem/orbit.h index b4a4b60c70..8c68907e86 100644 --- a/src/celephem/orbit.h +++ b/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: diff --git a/test/unit/kepler_test.cpp b/test/unit/kepler_test.cpp index d16c45f900..4c005d5156 100644 --- a/test/unit/kepler_test.cpp +++ b/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(fixedZero) + : celestia::util::array_view(testAngles); + for (double nodeDeg : testNodes) + { + auto testNodalPeriods = nodeDeg == 0.0 + ? celestia::util::array_view(fixedZero) + : celestia::util::array_view(testPeriods); + for (double nodalPeriod : testNodalPeriods) + for (double eccentricity : testEccentricities) + { + auto testPericenters = eccentricity == 0.0 + ? celestia::util::array_view(fixedZero) + : celestia::util::array_view(testAngles); + for (double pericenterDeg : testPericenters) + { + auto testApsidalPeriods = pericenterDeg == 0.0 + ? celestia::util::array_view(fixedZero) + : celestia::util::array_view(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 };