From 6956a4e18fabf2a52c0b3aafc7efc7dbc5ef3eb2 Mon Sep 17 00:00:00 2001 From: Alex Liang <35860220+alex-liang3@users.noreply.github.com> Date: Tue, 17 Dec 2024 17:50:49 +0100 Subject: [PATCH] feat: frozen orbit conditions (#473) * feat: add method to get coe for frozen orbit * feat: add frozen Orbit constructor * test: add coe gtest * chore: improve error handling * refactor: make re, j2, and j3 args * test: add tests for orbit * chore: bump ostk-physics dependency version * test: add python tests * perf: loosen criticality requirement * style: use camelCase * style: add units and long for to var names * style: make lambda captures explicit * perf: use ostk Arrays * test: add more tests * doc: add more documentation * test: improve tests * doc: add docs for Orbit constructor * feat: add overloaded ctor for COE * test: add test for ctor with Celestial * test: add python test * perf: use reference * chore: fix typos * doc: add doc to orbit.hpp * test: add negative tests for celestial ctor --- bindings/python/requirements.txt | 2 +- .../Trajectory/Orbit.cpp | 48 +++++ .../Trajectory/Orbit/Model/Kepler/COE.cpp | 113 ++++++++++ .../orbit/models/kepler/test_coe.py | 21 ++ bindings/python/test/trajectory/test_orbit.py | 9 + .../Astrodynamics/Trajectory/Orbit.hpp | 43 ++++ .../Trajectory/Orbit/Model/Kepler/COE.hpp | 83 ++++++++ .../Astrodynamics/Trajectory/Orbit.cpp | 36 ++++ .../Trajectory/Orbit/Model/Kepler/COE.cpp | 154 ++++++++++++++ .../Astrodynamics/Trajectory/Orbit.test.cpp | 74 +++++++ .../Orbit/Model/Kepler/COE.test.cpp | 200 ++++++++++++++++++ 11 files changed, 782 insertions(+), 1 deletion(-) diff --git a/bindings/python/requirements.txt b/bindings/python/requirements.txt index 512443264..758ac7233 100644 --- a/bindings/python/requirements.txt +++ b/bindings/python/requirements.txt @@ -3,4 +3,4 @@ open-space-toolkit-core~=4.1 open-space-toolkit-io~=4.0 open-space-toolkit-mathematics~=4.3 -open-space-toolkit-physics~=11.1 +open-space-toolkit-physics~=11.2 diff --git a/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit.cpp b/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit.cpp index 8db8dc3fd..72ef875f3 100644 --- a/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit.cpp +++ b/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit.cpp @@ -375,6 +375,54 @@ inline void OpenSpaceToolkitAstrodynamicsPy_Trajectory_Orbit(pybind11::module& a )doc" ) + .def_static( + "frozen", + &Orbit::Frozen, + arg("epoch"), + arg("altitude"), + arg("celestial_object"), + arg_v("eccentricity", Real::Undefined(), "Real.undefined()"), + arg_v("inclination", Angle::Undefined(), "Angle.undefined()"), + arg_v("raan", Angle::Degrees(0.0), "Angle.degrees(0.0)"), + arg_v("aop", Angle::Undefined(), "Angle.undefined()"), + arg_v("true_anomaly", Angle::Degrees(0.0), "Angle.degrees(0.0)"), + R"doc( + Create a frozen `Orbit` object. + + The critical angles for inclination are 63.4349 degrees and 116.5651 degrees. + The critical angles for AoP are 90.0 degrees and 270.0 degrees. + + At a minimum, an epoch, altitude, and celestial body with a defined J2 and J3 must be provided. + In this case, the inclination and AoP are set to critical angles, and the eccentricity is derived + from inclination. RAAN and true anomaly default to zero degrees. + + Additionally, the following combinations of inputs are supported: + - AoP (inclination set to critical value, eccentricity derived) + - AoP and eccentricity (inclination derived) + - AoP and inclination, but at least one of them must be a critical value (eccentricity derived) + - Inclination (AoP set to critical value, eccentricity derived) + - Eccentricity (AoP set to critical value, inclination derived) + + Note that inclination and eccentricity cannot both be provided. + + RAAN and True Anomaly may be provided alongside any of these arguments, and will be passed through + to the resulting Orbit as they do not impact the frozen orbit condition. + + Args: + epoch (Instant): The epoch. + altitude (Length): The altitude. + celestial_object (Celestial): The celestial object. + eccentricity (float): The eccentricity. + inclination (Angle): The inclination. + raan (Angle): The right ascension of the ascending node. + aop (Angle): The argument of periapsis. + true_anomaly (Angle): The true anomaly. + + Returns: + Orbit: The frozen `Orbit` object. + )doc" + ) + .def_static( "compute_passes", &Orbit::ComputePasses, diff --git a/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit/Model/Kepler/COE.cpp b/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit/Model/Kepler/COE.cpp index 7b8b6eaeb..7c08d6455 100644 --- a/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit/Model/Kepler/COE.cpp +++ b/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit/Model/Kepler/COE.cpp @@ -8,6 +8,7 @@ inline void OpenSpaceToolkitAstrodynamicsPy_Trajectory_Orbit_Model_Kepler_COE(py using ostk::core::type::Real; + using ostk::physics::environment::object::Celestial; using ostk::physics::unit::Angle; using ostk::physics::unit::Length; @@ -387,6 +388,118 @@ inline void OpenSpaceToolkitAstrodynamicsPy_Trajectory_Orbit_Model_Kepler_COE(py arg("anomaly_type") ) + .def_static( + "frozen_orbit", + overload_cast< + const Length&, + const Shared&, + const Real&, + const Angle&, + const Angle&, + const Angle&, + const Angle&>(&COE::FrozenOrbit), + R"doc( + Build a `COE` model of a frozen orbit. + + The critical angles for inclination are 63.4349 degrees and 116.5651 degrees. + The critical angles for AoP are 90.0 degrees and 270.0 degrees. + + At a minimum, a semi-major axis and shared pointer to a central celestial body with a defined J2 and J3 + must be provided. In this case, the inclination and AoP are set to critical angles, and the eccentricity + is derived from inclination. RAAN and true anomaly default to zero degrees. + + Additionally, the following combinations of inputs are supported: + - AoP (inclination set to critical value, eccentricity derived) + - AoP and eccentricity (inclination derived) + - AoP and inclination, but at least one of them must be a critical value (eccentricity derived) + - Inclination (AoP set to critical value, eccentricity derived) + - Eccentricity (AoP set to critical value, inclination derived) + + Note that inclination and eccentricity cannot both be provided. + + RAAN and True Anomaly may be provided alongside any of these arguments, and will be passed through + to the resulting COE as they do not impact the frozen orbit condition. + + Args: + semi_major_axis (Length): The semi-major axis. + celestial_object (Celestial): The celestial object. + eccentricity (float): The eccentricity. + inclination (Angle): The inclination. + raan (Angle): The right ascension of the ascending node. + aop (Angle): The argument of periapsis. + true_anomaly (Angle): The true anomaly. + + Returns: + COE: The `COE` model. + )doc", + arg("semi_major_axis"), + arg("celestial_object"), + arg_v("eccentricity", Real::Undefined(), "Real.undefined()"), + arg_v("inclination", Angle::Undefined(), "Angle.undefined()"), + arg_v("raan", Angle::Degrees(0.0), "Angle.degrees(0.0)"), + arg_v("aop", Angle::Undefined(), "Angle.undefined()"), + arg_v("true_anomaly", Angle::Degrees(0.0), "Angle.degrees(0.0)") + ) + + .def_static( + "frozen_orbit", + overload_cast< + const Length&, + const Length&, + const Real&, + const Real&, + const Real&, + const Angle&, + const Angle&, + const Angle&, + const Angle&>(&COE::FrozenOrbit), + R"doc( + Build a `COE` model of a frozen orbit. + + The critical angles for inclination are 63.4349 degrees and 116.5651 degrees. + The critical angles for AoP are 90.0 degrees and 270.0 degrees. + + At a minimum, a semi-major axis, equatorial radius, J2, and J3 must be provided. In this case, + the inclination and AoP are set to critical angles, and the eccentricity is derived from inclination. + RAAN and true anomaly default to zero degrees. + + Additionally, the following combinations of inputs are supported: + - AoP (inclination set to critical value, eccentricity derived) + - AoP and eccentricity (inclination derived) + - AoP and inclination, but at least one of them must be a critical value (eccentricity derived) + - Inclination (AoP set to critical value, eccentricity derived) + - Eccentricity (AoP set to critical value, inclination derived) + + Note that inclination and eccentricity cannot both be provided. + + RAAN and True Anomaly may be provided alongside any of these arguments, and will be passed through + to the resulting COE as they do not impact the frozen orbit condition. + + Args: + semi_major_axis (Length): The semi-major axis. + equatorial_radius (Length): The equatorial radius. + j2 (float): The second zonal harmonic coefficient. + j3 (float): The third zonal harmonic coefficient. + eccentricity (float): The eccentricity. + inclination (Angle): The inclination. + raan (Angle): The right ascension of the ascending node. + aop (Angle): The argument of periapsis. + true_anomaly (Angle): The true anomaly. + + Returns: + COE: The `COE` model. + )doc", + arg("semi_major_axis"), + arg("equatorial_radius"), + arg("j2"), + arg("j3"), + arg_v("eccentricity", Real::Undefined(), "Real.undefined()"), + arg_v("inclination", Angle::Undefined(), "Angle.undefined()"), + arg_v("raan", Angle::Degrees(0.0), "Angle.degrees(0.0)"), + arg_v("aop", Angle::Undefined(), "Angle.undefined()"), + arg_v("true_anomaly", Angle::Degrees(0.0), "Angle.degrees(0.0)") + ) + .def_static( "eccentric_anomaly_from_true_anomaly", &COE::EccentricAnomalyFromTrueAnomaly, diff --git a/bindings/python/test/trajectory/orbit/models/kepler/test_coe.py b/bindings/python/test/trajectory/orbit/models/kepler/test_coe.py index a8f023d9f..b1da95df8 100644 --- a/bindings/python/test/trajectory/orbit/models/kepler/test_coe.py +++ b/bindings/python/test/trajectory/orbit/models/kepler/test_coe.py @@ -7,6 +7,7 @@ from ostk.physics.time import Instant from ostk.physics.environment.gravitational import Earth from ostk.physics.environment.object.celestial import Sun +from ostk.physics import Environment from ostk.astrodynamics.trajectory.orbit.model.kepler import COE @@ -74,6 +75,26 @@ def test_constructors(self): assert isinstance(coe, COE) assert coe.is_defined() is False + coe: COE = COE.frozen_orbit( + a, + Earth.EGM2008.equatorial_radius, + Earth.EGM2008.J2, + Earth.EGM2008.J3, + inclination=i, + ) + + assert coe is not None + assert isinstance(coe, COE) + assert coe.is_defined() + + coe: COE = COE.frozen_orbit( + a, Environment.default().access_celestial_object_with_name("Earth") + ) + + assert coe is not None + assert isinstance(coe, COE) + assert coe.is_defined() + def test_comparators( self, coe: COE, diff --git a/bindings/python/test/trajectory/test_orbit.py b/bindings/python/test/trajectory/test_orbit.py index 2fd8c189a..796cbeb52 100644 --- a/bindings/python/test/trajectory/test_orbit.py +++ b/bindings/python/test/trajectory/test_orbit.py @@ -181,6 +181,15 @@ def test_sun_synchronous(self, earth): argument_of_latitude=Angle.degrees(50.0), ).is_defined() + def test_frozen(self, earth, epoch): + altitude = Length.kilometers(500.0) + + orbit: Orbit = Orbit.frozen(epoch, altitude, earth) + + assert orbit is not None + assert isinstance(orbit, Orbit) + assert orbit.is_defined() + def test_compute_passes(self, states: list[State]): passes: list[tuple[int, Pass]] = Orbit.compute_passes(states, 1) assert passes is not None diff --git a/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.hpp b/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.hpp index 0f616d937..764c5098f 100644 --- a/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.hpp +++ b/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.hpp @@ -250,6 +250,49 @@ class Orbit : public Trajectory const Angle& anArgumentOfLatitude = Angle::Zero() ); + /// @brief Constructs a Frozen orbit + /// + /// Model: Kepler (J2 Perturbation) + /// + /// The critical angles for inclination are 63.4349 degrees and 116.5651 degrees. + /// The critical angles for AoP are 90.0 degrees and 270.0 degrees. + /// + /// At a minimum, an epoch, altitude, and celestial body with a defined J2 and J3 must be provided. + /// In this case, the inclination and AoP are set to critical angles, and the eccentricity is derived + /// from inclination. RAAN and true anomaly default to zero degrees. + /// + /// Additionally, the following combinations of inputs are supported: + /// - AoP (inclination set to critical value, eccentricity derived) + /// - AoP and eccentricity (inclination derived) + /// - AoP and inclination, but at least one of them must be a critical value (eccentricity derived) + /// - Inclination (AoP set to critical value, eccentricity derived) + /// - Eccentricity (AoP set to critical value, inclination derived) + /// + /// Note that inclination and eccentricity cannot both be provided. + /// + /// RAAN and True Anomaly may be provided alongside any of these arguments, and will be passed through + /// to the resulting Orbit as they do not impact the frozen orbit condition. + /// + /// @param anEpoch An orbit epoch + /// @param anAltitude An orbit altitude (wrt. equatorial radius) + /// @param aCelestialObjectSPtr A shared pointer to a central celestial body + /// @param anEccentricity An eccentricity + /// @param anInclination An inclination + /// @param aRaan A raan + /// @param anAop An aop + /// @param aTrueAnomaly A true anomaly + /// @return Frozen orbit + static Orbit Frozen( + const Instant& anEpoch, + const Length& anAltitude, + const Shared& aCelestialObjectSPtr, + const Real& anEccentricity = Real::Undefined(), + const Angle& anInclination = Angle::Undefined(), + const Angle& aRaan = Angle::Degrees(0.0), + const Angle& anAop = Angle::Undefined(), + const Angle& aTrueAnomaly = Angle::Degrees(0.0) + ); + /// @brief Get the string representation of a frame type /// /// @param aFrameType Type of the frame to get the string representation of diff --git a/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Kepler/COE.hpp b/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Kepler/COE.hpp index 81af6a3e7..6d07e9516 100644 --- a/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Kepler/COE.hpp +++ b/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Kepler/COE.hpp @@ -42,6 +42,7 @@ using ostk::mathematics::object::Vector6d; using ostk::physics::coordinate::Frame; using ostk::physics::coordinate::Position; using ostk::physics::coordinate::Velocity; +using ostk::physics::environment::object::Celestial; using ostk::physics::environment::object::celestial::Sun; using ostk::physics::time::Duration; using ostk::physics::time::Instant; @@ -252,6 +253,88 @@ class COE /// @return COE static COE FromSIVector(const Vector6d& aCOEVector, const AnomalyType& anAnomalyType); + /// @brief Construct a frozen orbit from an incomplete set of COEs + /// + /// The critical angles for inclination are 63.4349 degrees and 116.5651 degrees. + /// The critical angles for AoP are 90.0 degrees and 270.0 degrees. + /// + /// At a minimum, a semi-major axis and shared pointer to a central celestial body with a defined J2 and J3 + /// must be provided. In this case, the inclination and AoP are set to critical angles, and the eccentricity + /// is derived from inclination. RAAN and true anomaly default to zero degrees. + /// + /// Additionally, the following combinations of inputs are supported: + /// - AoP (inclination set to critical value, eccentricity derived) + /// - AoP and eccentricity (inclination derived) + /// - AoP and inclination, but at least one of them must be a critical value (eccentricity derived) + /// - Inclination (AoP set to critical value, eccentricity derived) + /// - Eccentricity (AoP set to critical value, inclination derived) + /// + /// Note that inclination and eccentricity cannot both be provided. + /// + /// RAAN and True Anomaly may be provided alongside any of these arguments, and will be passed through + /// to the resulting COE as they do not impact the frozen orbit condition. + /// + /// @param aSemiMajorAxis A semi-major axis + /// @param aCelestialObjectSPtr A shared pointer to a central celestial body + /// @param anEccentricity An eccentricity + /// @param anInclination An inclination + /// @param aRaan A raan + /// @param anAop An aop + /// @param aTrueAnomaly A true anomaly + /// @return COE + static COE FrozenOrbit( + const Length& aSemiMajorAxis, + const Shared& aCelestialObjectSPtr, + const Real& anEccentricity = Real::Undefined(), + const Angle& anInclination = Angle::Undefined(), + const Angle& aRaan = Angle::Degrees(0.0), + const Angle& anAop = Angle::Undefined(), + const Angle& aTrueAnomaly = Angle::Degrees(0.0) + ); + + /// @brief Construct a frozen orbit from an incomplete set of COEs + /// + /// The critical angles for inclination are 63.4349 degrees and 116.5651 degrees. + /// The critical angles for AoP are 90.0 degrees and 270.0 degrees. + /// + /// At a minimum, a semi-major axis, equatorial radius, J2, and J3 must be provided. In this case, + /// the inclination and AoP are set to critical angles, and the eccentricity is derived from inclination. + /// RAAN and true anomaly default to zero degrees. + /// + /// Additionally, the following combinations of inputs are supported: + /// - AoP (inclination set to critical value, eccentricity derived) + /// - AoP and eccentricity (inclination derived) + /// - AoP and inclination, but at least one of them must be a critical value (eccentricity derived) + /// - Inclination (AoP set to critical value, eccentricity derived) + /// - Eccentricity (AoP set to critical value, inclination derived) + /// + /// Note that inclination and eccentricity cannot both be provided. + /// + /// RAAN and True Anomaly may be provided alongside any of these arguments, and will be passed through + /// to the resulting COE as they do not impact the frozen orbit condition. + /// + /// @param aSemiMajorAxis A semi-major axis + /// @param anEquatorialRadius An equatorial radius + /// @param aJ2 A J2 + /// @param aJ3 A J3 + /// @param anEccentricity An eccentricity + /// @param anInclination An inclination + /// @param aRaan A raan + /// @param anAop An aop + /// @param aTrueAnomaly A true anomaly + /// @return COE + static COE FrozenOrbit( + const Length& aSemiMajorAxis, + const Length& anEquatorialRadius, + const Real& aJ2, + const Real& aJ3, + const Real& anEccentricity = Real::Undefined(), + const Angle& anInclination = Angle::Undefined(), + const Angle& aRaan = Angle::Degrees(0.0), + const Angle& anAop = Angle::Undefined(), + const Angle& aTrueAnomaly = Angle::Degrees(0.0) + ); + /// @brief Convert True anomaly to Eccentric anomaly /// /// @param aTrueAnomaly A true anomaly diff --git a/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.cpp b/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.cpp index 9168265bb..ccaec3029 100644 --- a/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.cpp +++ b/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.cpp @@ -1078,6 +1078,42 @@ Orbit Orbit::SunSynchronous( return {orbitalModel, aCelestialObjectSPtr}; } +Orbit Orbit::Frozen( + const Instant& anEpoch, + const Length& anAltitude, + const Shared& aCelestialObjectSPtr, + const Real& anEccentricity, + const Angle& anInclination, + const Angle& aRaan, + const Angle& anAop, + const Angle& aTrueAnomaly +) +{ + if (!anEpoch.isDefined()) + { + throw ostk::core::error::runtime::Undefined("Epoch"); + } + if (!anAltitude.isDefined()) + { + throw ostk::core::error::runtime::Undefined("Altitude"); + } + if ((aCelestialObjectSPtr == nullptr) || (!aCelestialObjectSPtr->isDefined())) + { + throw ostk::core::error::runtime::Undefined("Celestial object"); + } + + const Length equatorialRadius = aCelestialObjectSPtr->getEquatorialRadius(); + const Length semiMajorAxis = equatorialRadius + anAltitude; + + const COE coe = COE::FrozenOrbit( + semiMajorAxis, aCelestialObjectSPtr, anEccentricity, anInclination, aRaan, anAop, aTrueAnomaly + ); + + const Kepler orbitalModel = {coe, anEpoch, {*aCelestialObjectSPtr}, Kepler::PerturbationType::J2, false}; + + return {orbitalModel, aCelestialObjectSPtr}; +} + String Orbit::StringFromFrameType(const Orbit::FrameType& aFrameType) { switch (aFrameType) diff --git a/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Kepler/COE.cpp b/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Kepler/COE.cpp index 452c37d51..cc2c98931 100644 --- a/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Kepler/COE.cpp +++ b/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Kepler/COE.cpp @@ -24,6 +24,7 @@ namespace model namespace kepler { +using ostk::core::type::Array; using ostk::core::type::Integer; using ostk::core::type::Real; using ostk::core::type::Shared; @@ -641,6 +642,159 @@ COE COE::FromSIVector(const Vector6d& aCOEVector, const AnomalyType& anAnomalyTy }; } +COE COE::FrozenOrbit( + const Length& aSemiMajorAxis, + const Shared& aCelestialObjectSPtr, + const Real& anEccentricity, + const Angle& anInclination, + const Angle& aRaan, + const Angle& anAop, + const Angle& aTrueAnomaly +) +{ + if ((aCelestialObjectSPtr == nullptr) || (!aCelestialObjectSPtr->isDefined())) + { + throw ostk::core::error::runtime::Undefined("Celestial object"); + } + + const Length equatorialRadius = aCelestialObjectSPtr->getEquatorialRadius(); + const Real j2 = aCelestialObjectSPtr->accessGravitationalModel()->getParameters().J2_; + const Real j3 = aCelestialObjectSPtr->accessGravitationalModel()->getParameters().J3_; + + return COE::FrozenOrbit( + aSemiMajorAxis, equatorialRadius, j2, j3, anEccentricity, anInclination, aRaan, anAop, aTrueAnomaly + ); +} + +COE COE::FrozenOrbit( + const Length& aSemiMajorAxis, + const Length& anEquatorialRadius, + const Real& aJ2, + const Real& aJ3, + const Real& anEccentricity, + const Angle& anInclination, + const Angle& aRaan, + const Angle& anAop, + const Angle& aTrueAnomaly +) +{ + if (!aSemiMajorAxis.isDefined()) + { + throw ostk::core::error::runtime::Undefined("Semi-major Axis"); + } + if (!anEquatorialRadius.isDefined()) + { + throw ostk::core::error::runtime::Undefined("Equatorial Radius"); + } + if (!aJ2.isDefined() || aJ2 == 0.0) + { + throw ostk::core::error::runtime::Undefined("J2"); + } + if (!aJ3.isDefined() || aJ3 == 0.0) + { + throw ostk::core::error::runtime::Undefined("J3"); + } + if (!aRaan.isDefined()) + { + throw ostk::core::error::runtime::Undefined("Raan"); + } + if (!aTrueAnomaly.isDefined()) + { + throw ostk::core::error::runtime::Undefined("True anomaly"); + } + if (anEccentricity.isDefined() && anInclination.isDefined()) + { + throw ostk::core::error::RuntimeError("Cannot define both eccentricity and inclination"); + } + + const Array criticalInclinations = {Angle::Degrees(63.4349), Angle::Degrees(116.5651)}; + const Array criticalAops = {Angle::Degrees(90.0), Angle::Degrees(270.0)}; + + const Real equatorialRadius_meters = anEquatorialRadius.inMeters(); + const Real j2 = aJ2; + const Real j3 = aJ3; + + // ecc =~ eccCoefficient * sin(incl) + const Real eccCoefficient = -j3 * equatorialRadius_meters / 2.0 / j2 / aSemiMajorAxis.inMeters(); + + const auto eccentricityFromInclination = [&eccCoefficient](const Angle& inclination) -> Real + { + return eccCoefficient * std::sin(inclination.inRadians()); + }; + + const auto inclinationFromEccentricity = [&eccCoefficient](const Real& eccentricity) -> Angle + { + return Angle::Radians(asin(eccentricity / eccCoefficient)); + }; + + // TBI: this only works because neither set of critical angles are 0/360 degrees + const auto isCritical = [](const Angle& angle, const Array& criticalAngles) -> bool + { + const Real epsilon = Angle::Arcseconds(1.0).inRadians(); // TBI: make configurable? + + return (angle.inRadians().isNear(criticalAngles[0].inRadians(), epsilon)) || + (angle.inRadians().isNear(criticalAngles[1].inRadians(), epsilon)); + }; + + // Use the provided AoP, or default to a critical value + const Angle aop = anAop.isDefined() ? anAop : criticalAops[0]; + + // If AoP matches a critical value + if (isCritical(aop, criticalAops)) + { + Angle inclination = Angle::Undefined(); + Real eccentricity = Real::Undefined(); + + // If ecc not defined, use the given inclination or default to critical inclination, and calculate ecc + if (!anEccentricity.isDefined()) + { + inclination = anInclination.isDefined() ? anInclination : criticalInclinations[0]; + eccentricity = eccentricityFromInclination(inclination); + } + // If ecc defined, use it to calculate inc + else + { + // If the eccentricity is larger than this value, the approximation isn't valid + if (anEccentricity > eccCoefficient) + { + throw ostk::core::error::RuntimeError( + "Provided eccentricity [" + anEccentricity.toString() + "] cannot be greater than " + + eccCoefficient.toString() + ); + } + inclination = inclinationFromEccentricity(anEccentricity); + eccentricity = anEccentricity; + } + + return { + aSemiMajorAxis, + eccentricity, + inclination, + aRaan, + aop, + aTrueAnomaly, + }; + } + + // If AoP is (given) and not a critical angle, then the inclination must be critical + const Angle inclination = anInclination.isDefined() ? anInclination : criticalInclinations[0]; + if (!isCritical(inclination, criticalInclinations)) + { + throw ostk::core::error::RuntimeError("Provided inclination must be a critical value"); + } + + const Real eccentricity = eccentricityFromInclination(inclination); + + return { + aSemiMajorAxis, + eccentricity, + inclination, + aRaan, + aop, + aTrueAnomaly, + }; +} + Angle COE::EccentricAnomalyFromTrueAnomaly(const Angle& aTrueAnomaly, const Real& anEccentricity) { if (!aTrueAnomaly.isDefined()) diff --git a/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.test.cpp b/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.test.cpp index f7d547643..e8e9afa32 100644 --- a/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.test.cpp +++ b/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.test.cpp @@ -2127,3 +2127,77 @@ TEST(OpenSpaceToolkit_Astrodynamics_Trajectory_Orbit, SunSynchronous) // } // } } + +TEST(OpenSpaceToolkit_Astrodynamics_Trajectory_Orbit, Frozen) +{ + const Instant epoch = Instant::DateTime(DateTime::Parse("2018-01-01 00:00:00"), Scale::UTC); + const Length altitude = Length::Kilometers(500.0); + { + { + const Orbit orbit = + Orbit::Frozen(epoch, altitude, Environment::Default().accessCelestialObjectWithName("Earth")); + + EXPECT_TRUE(orbit.isDefined()); + } + + // Provide target inclination + { + const Orbit orbit = Orbit::Frozen( + epoch, + altitude, + Environment::Default().accessCelestialObjectWithName("Earth"), + Real::Undefined(), + Angle::Degrees(60.0) + ); + + EXPECT_TRUE(orbit.isDefined()); + EXPECT_EQ( + orbit.accessModel().as().getClassicalOrbitalElements().getInclination(), Angle::Degrees(60.0) + ); + } + } + + { + // Undefined epoch + { + EXPECT_THROW( + Orbit::Frozen( + Instant::Undefined(), altitude, Environment::Default().accessCelestialObjectWithName("Earth") + ), + ostk::core::error::runtime::Undefined + ); + } + + // Undefined altitude + { + EXPECT_THROW( + Orbit::Frozen( + epoch, Length::Undefined(), Environment::Default().accessCelestialObjectWithName("Earth") + ), + ostk::core::error::runtime::Undefined + ); + } + + // Undefined celestial object + { + { + EXPECT_THROW(Orbit::Frozen(epoch, altitude, nullptr), ostk::core::error::runtime::Undefined); + } + + { + EXPECT_THROW( + Orbit::Frozen(epoch, altitude, std::make_shared(Celestial::Undefined())), + ostk::core::error::runtime::Undefined + ); + } + } + + // Celestial object model with no J3 term + { + EXPECT_THROW( + Orbit::Frozen(epoch, altitude, std::make_shared(Earth::WGS84())), + ostk::core::error::runtime::Undefined + ); + } + } +} diff --git a/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Kepler/COE.test.cpp b/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Kepler/COE.test.cpp index d2bb01c89..3c9556311 100644 --- a/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Kepler/COE.test.cpp +++ b/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Kepler/COE.test.cpp @@ -5,6 +5,7 @@ #include +#include #include #include #include @@ -28,7 +29,9 @@ using ostk::mathematics::object::Vector6d; using ostk::physics::coordinate::Frame; using ostk::physics::coordinate::Position; using ostk::physics::coordinate::Velocity; +using ostk::physics::Environment; using ostk::physics::environment::gravitational::Earth; +using ostk::physics::environment::object::celestial::Celestial; using ostk::physics::environment::object::celestial::Sun; using ostk::physics::time::DateTime; using ostk::physics::time::Duration; @@ -675,6 +678,203 @@ TEST_F(OpenSpaceToolkit_Astrodynamics_Trajectory_Orbit_Model_Kepler_COE, FromSIV } } +TEST_F(OpenSpaceToolkit_Astrodynamics_Trajectory_Orbit_Model_Kepler_COE, FrozenOrbit) +{ + const Length semiMajorAxis = Length::Kilometers(8000); + + const Length re = Earth::EGM2008.equatorialRadius_; + const Real j2 = Earth::EGM2008.J2_; + const Real j3 = Earth::EGM2008.J3_; + + // Only provide SMA + { + const COE coe = COE::FrozenOrbit(semiMajorAxis, re, j2, j3); + + EXPECT_NEAR(coe.getSemiMajorAxis().inMeters(), semiMajorAxis.inMeters(), 1e-10); + EXPECT_NEAR(coe.getEccentricity(), 8.340085671757e-4, 1e-10); + EXPECT_NEAR(coe.getInclination().inDegrees(), 63.4349, 1e-10); + EXPECT_NEAR(coe.getRaan().inDegrees(), 0.0, 1e-10); + EXPECT_NEAR(coe.getAop().inDegrees(), 90.0, 1e-10); + EXPECT_NEAR(coe.getTrueAnomaly().inDegrees(), 0.0, 1e-10); + } + + // Provide inclination + // @ref: https://www.mathworks.com/matlabcentral/fileexchange/39119-frozen-orbit-design?s_tid=FX_rc3_behav + { + const COE coe = COE::FrozenOrbit(semiMajorAxis, re, j2, j3, Real::Undefined(), Angle::Degrees(45.0)); + + EXPECT_NEAR(coe.getEccentricity(), 6.594e-4, 1e-6); + EXPECT_NEAR(coe.getAop().inDegrees(), 90.0, 1e-10); + } + + // Provide critical AoP and eccentricity + // @ref: https://www.mathworks.com/matlabcentral/fileexchange/39119-frozen-orbit-design?s_tid=FX_rc3_behav + { + const COE coe = COE::FrozenOrbit( + semiMajorAxis, re, j2, j3, 6.5941377e-4, Angle::Undefined(), Angle::Degrees(0.0), Angle::Degrees(90.0) + ); + + EXPECT_NEAR(coe.getInclination().inDegrees(), 45.0, 1e-2); + } + + // Provide AoP + { + const COE coe = COE::FrozenOrbit( + semiMajorAxis, + re, + j2, + j3, + Real::Undefined(), + Angle::Undefined(), + Angle::Degrees(45.0), + Angle::Degrees(45.0), + Angle::Degrees(0.0) + ); + + EXPECT_NEAR(coe.getInclination().inDegrees(), 63.4349, 1e-10); + } + + // Non-exact critical angle supplied + { + const Angle non_critical_aop = Angle::Degrees(80.0); + const COE coe = COE::FrozenOrbit( + semiMajorAxis, + re, + j2, + j3, + Real::Undefined(), + Angle::Degrees(63.43490001), + Angle::Degrees(0.0), + non_critical_aop + ); + + EXPECT_NEAR(coe.getAop().inDegrees(), non_critical_aop.inDegrees(), 1e-10); + } + + // Construct with Celestial + { + COE coe = COE::FrozenOrbit(semiMajorAxis, Environment::Default().accessCelestialObjectWithName("Earth")); + + EXPECT_TRUE(coe.isDefined()); + } + + { + { + // No SMA provided + { + EXPECT_THROW(COE::FrozenOrbit(Length::Undefined(), re, j2, j3), ostk::core::error::runtime::Undefined); + } + + // No equatorial radius provided + { + EXPECT_THROW( + COE::FrozenOrbit(semiMajorAxis, Length::Undefined(), j2, j3), ostk::core::error::runtime::Undefined + ); + } + + // No J2 provided + { + { + EXPECT_THROW( + COE::FrozenOrbit(semiMajorAxis, re, Real::Undefined(), j3), + ostk::core::error::runtime::Undefined + ); + } + + { + EXPECT_THROW(COE::FrozenOrbit(semiMajorAxis, re, 0.0, j3), ostk::core::error::runtime::Undefined); + } + } + + // No J3 provided + { + { + EXPECT_THROW( + COE::FrozenOrbit(semiMajorAxis, re, j2, Real::Undefined()), + ostk::core::error::runtime::Undefined + ); + } + + { + EXPECT_THROW(COE::FrozenOrbit(semiMajorAxis, re, j2, 0.0), ostk::core::error::runtime::Undefined); + } + } + + // No RAAN provided + { + EXPECT_THROW( + COE::FrozenOrbit( + semiMajorAxis, re, j2, j3, Real::Undefined(), Angle::Undefined(), Angle::Undefined() + ), + ostk::core::error::runtime::Undefined + ); + } + + // No true anomaly provided + { + EXPECT_THROW( + COE::FrozenOrbit( + semiMajorAxis, + re, + j2, + j3, + Real::Undefined(), + Angle::Undefined(), + Angle::Degrees(0.0), + Angle::Undefined(), + Angle::Undefined() + ), + ostk::core::error::runtime::Undefined + ); + } + + // Null Celestial provided + { + { + EXPECT_THROW( + COE::FrozenOrbit(semiMajorAxis, std::make_shared(Celestial::Undefined())), + ostk::core::error::runtime::Undefined + ); + } + + { + EXPECT_THROW(COE::FrozenOrbit(semiMajorAxis, nullptr), ostk::core::error::runtime::Undefined); + } + } + } + + // Excessively large eccentricity + { + EXPECT_THROW(COE::FrozenOrbit(semiMajorAxis, re, j2, j3, 0.1), ostk::core::error::RuntimeError); + } + + // Inclination and eccentricity both provided + { + EXPECT_THROW( + COE::FrozenOrbit(semiMajorAxis, re, j2, j3, 0.0001, Angle::Degrees(45.0)), + ostk::core::error::RuntimeError + ); + } + + // AoP and inclination both provided and both non-critical + { + EXPECT_THROW( + COE::FrozenOrbit( + semiMajorAxis, + re, + j2, + j3, + Real::Undefined(), + Angle::Degrees(45.0), + Angle::Degrees(0.0), + Angle::Degrees(45.0) + ), + ostk::core::error::RuntimeError + ); + } + } +} + // TEST (OpenSpaceToolkit_Astrodynamics_Trajectory_Orbit_Model_Kepler_COE, EccentricAnomalyFromTrueAnomaly) // {