diff --git a/bindings/python/requirements.txt b/bindings/python/requirements.txt index b6952ce95..29a33a5c0 100644 --- a/bindings/python/requirements.txt +++ b/bindings/python/requirements.txt @@ -4,3 +4,4 @@ open-space-toolkit-core~=3.0 open-space-toolkit-io~=3.0 open-space-toolkit-mathematics~=3.0 open-space-toolkit-physics~=7.0 +numpy~=1.26 \ No newline at end of file diff --git a/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit.cpp b/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit.cpp index 5e4c9dbe4..810c05fcf 100644 --- a/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit.cpp +++ b/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit.cpp @@ -229,6 +229,20 @@ inline void OpenSpaceToolkitAstrodynamicsPy_Trajectory_Orbit(pybind11::module& a arg("revolution_number"), arg_v("step_duration", Duration::Minutes(10.0), "Duration.minutes(10.0)") ) + .def( + "get_passes_within_interval", + &Orbit::getPassesWithinInterval, + R"doc( + Get the passes within a given interval. + + Args: + interval (Interval): The interval. + + Returns: + list[Pass]: The passes. + )doc", + arg("interval") + ) .def( "get_orbital_frame", &Orbit::getOrbitalFrame, @@ -374,7 +388,28 @@ inline void OpenSpaceToolkitAstrodynamicsPy_Trajectory_Orbit(pybind11::module& a initial_revolution_number (Integer): The initial revolution number. Returns: - list[tuple[int, Pass]]: List of index-pass pairs + list[tuple[int, Pass]]: List of index-pass pairs. + )doc" + ) + + .def_static( + "compute_passes_with_model", + &Orbit::ComputePassesWithModel, + arg("model"), + arg("start_instant"), + arg("end_instant"), + arg("initial_revolution_number"), + R"doc( + Compute passes with the given model for the provided interval. + + Args: + model (orbit.Model): The model. + start_instant (Instant): The start instant. + end_instant (Instant): The end instant. + initial_revolution_number (int): The initial revolution number. + + Returns: + list[Pass]: List of passes. )doc" ) diff --git a/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit/Model/SGP4.cpp b/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit/Model/SGP4.cpp index 874d90462..d7e7a44ab 100644 --- a/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit/Model/SGP4.cpp +++ b/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit/Model/SGP4.cpp @@ -99,22 +99,6 @@ inline void OpenSpaceToolkitAstrodynamicsPy_Trajectory_Orbit_Model_SGP4(pybind11 )doc" ) - .def( - "calculate_revolution_number_at", - &SGP4::calculateRevolutionNumberAt, - arg("instant"), - R"doc( - Calculate the revolution number of the `SGP4` model at a given instant. - - Args: - instant (Instant): The instant. - - Returns: - int: The revolution number. - - )doc" - ) - ; } diff --git a/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit/Model/Tabulated.cpp b/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit/Model/Tabulated.cpp index eb55da7a3..0dbc06042 100644 --- a/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit/Model/Tabulated.cpp +++ b/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit/Model/Tabulated.cpp @@ -139,21 +139,5 @@ inline void OpenSpaceToolkitAstrodynamicsPy_Trajectory_Orbit_Model_Tabulated(pyb arg("instants") ) - .def( - "calculate_revolution_number_at", - &Tabulated::calculateRevolutionNumberAt, - R"doc( - Calculate the revolution number of the `Tabulated` model at a given instant. - - Args: - instant (Instant): The instant. - - Returns: - int: The revolution number. - - )doc", - arg("instant") - ) - ; } diff --git a/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit/Pass.cpp b/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit/Pass.cpp index e8628d127..afb510dce 100644 --- a/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit/Pass.cpp +++ b/bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Orbit/Pass.cpp @@ -157,6 +157,16 @@ inline void OpenSpaceToolkitAstrodynamicsPy_Trajectory_Orbit_Pass(pybind11::modu )doc" ) + .def( + "get_interval", + &Pass::getInterval, + R"doc( + Get the interval of the pass. Undefined if the pass is not complete. + + Returns: + Interval: The interval of the pass. + )doc" + ) .def( "get_instant_at_ascending_node", &Pass::accessInstantAtAscendingNode, diff --git a/bindings/python/test/trajectory/orbit/test_pass.py b/bindings/python/test/trajectory/orbit/test_pass.py index 99f80b5b7..935b65c42 100644 --- a/bindings/python/test/trajectory/orbit/test_pass.py +++ b/bindings/python/test/trajectory/orbit/test_pass.py @@ -47,6 +47,9 @@ def test_get_start_instant(self, pass_: Pass): def test_get_end_instant(self, pass_: Pass): assert pass_.get_end_instant() is not None + def test_get_interval(self, pass_: Pass): + assert pass_.get_interval() is not None + def test_get_instant_at_ascending_node(self, pass_: Pass): assert pass_.get_instant_at_ascending_node() is not None diff --git a/bindings/python/test/trajectory/test_orbit.py b/bindings/python/test/trajectory/test_orbit.py index f94ef3b5b..2fd8c189a 100644 --- a/bindings/python/test/trajectory/test_orbit.py +++ b/bindings/python/test/trajectory/test_orbit.py @@ -98,7 +98,19 @@ def test_get_pass_with_revolution_number(self, orbit: Orbit): assert isinstance(pass_, Pass) assert pass_.is_defined() - assert orbit.get_pass_with_revolution_number(2, Duration.minutes(10.0)) is not None + assert ( + orbit.get_pass_with_revolution_number(2, Duration.minutes(10.0)) is not None + ) + + def test_get_passes_within_interval(self, orbit: Orbit): + passes: list[Pass] = orbit.get_passes_within_interval( + Interval.closed( + Instant.date_time(DateTime(2018, 1, 1, 0, 0, 0), Scale.UTC), + Instant.date_time(DateTime(2018, 1, 1, 0, 10, 0), Scale.UTC), + ) + ) + + assert len(passes) > 0 def test_undefined(self): assert Orbit.undefined().is_defined() is False @@ -169,6 +181,16 @@ def test_sun_synchronous(self, earth): argument_of_latitude=Angle.degrees(50.0), ).is_defined() - def test_compute_passes(self, orbit: Orbit, states: list[State]): - passes: list[tuple[int, Pass]] = orbit.compute_passes(states, 1) + def test_compute_passes(self, states: list[State]): + passes: list[tuple[int, Pass]] = Orbit.compute_passes(states, 1) assert passes is not None + + def test_compute_passes_with_model(self, orbit: Orbit): + passes: list[tuple[int, Pass]] = Orbit.compute_passes_with_model( + model=orbit.access_kepler_model(), + start_instant=Instant.date_time(DateTime(2018, 1, 1, 0, 0, 0), Scale.UTC), + end_instant=Instant.date_time(DateTime(2018, 1, 1, 0, 10, 0), Scale.UTC), + initial_revolution_number=1, + ) + + assert len(passes) > 0 diff --git a/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.hpp b/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.hpp index 7b2c8a4d5..b87b91a30 100644 --- a/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.hpp +++ b/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.hpp @@ -7,6 +7,7 @@ #include #include +#include #include #include #include @@ -35,6 +36,7 @@ namespace trajectory using ostk::core::container::Array; using ostk::core::container::Map; using ostk::core::container::Pair; +using ostk::core::container::Tuple; using ostk::core::type::Index; using ostk::core::type::Integer; using ostk::core::type::Real; @@ -55,7 +57,7 @@ using ostk::astrodynamics::trajectory::State; /// @brief Gravitationally curved trajectory of an object /// -/// @ref https://en.wikipedia.org/wiki/Orbit +/// @ref https://en.wikipedia.org/wiki/Orbit class Orbit : public Trajectory { public: @@ -78,46 +80,101 @@ class Orbit : public Trajectory typedef Array::ConstIterator ConstPassIterator; + /// @brief Constructor. + /// + /// @param aModel An orbit model. + /// @param aCelestialObjectSPtr A shared pointer to a central celestial body. Orbit(const orbit::Model& aModel, const Shared& aCelestialObjectSPtr); + /// @brief Constructor. + /// + /// @param aStateArray An array of states. + /// @param anInitialRevolutionNumber An initial revolution number. + /// @param aCelestialObjectSPtr A shared pointer to a central celestial body. Orbit( const Array& aStateArray, const Integer& anInitialRevolutionNumber, const Shared& aCelestialObjectSPtr ); + /// @brief Copy constructor. + /// + /// @param anOrbit Orbit to copy from. Orbit(const Orbit& anOrbit); + /// @brief Destructor. ~Orbit(); + /// @brief Copy assignment operator. + /// + /// @param anOrbit Orbit to copy from. + /// @return Reference to the updated Orbit. Orbit& operator=(const Orbit& anOrbit); + /// @brief Equality operator. + /// + /// @param anOrbit Orbit to compare to. + /// @return True if Orbits are equal, false otherwise. bool operator==(const Orbit& anOrbit) const; + /// @brief Inequality operator. + /// + /// @param anOrbit Orbit to compare to. + /// @return True if Orbits are not equal, false otherwise. bool operator!=(const Orbit& anOrbit) const; + /// @brief Check if the Orbit is defined. + /// + /// @return True if the Orbit is defined, false otherwise. bool isDefined() const; + /// @brief Get the revolution number at a given instant. + /// + /// @param anInstant Instant to get the revolution number at. + /// @return Revolution number at the given instant. Integer getRevolutionNumberAt(const Instant& anInstant) const; + /// @brief Get the pass at a given instant. + /// + /// @param anInstant Instant to get the pass at. + /// @return Pass at the given instant. Pass getPassAt(const Instant& anInstant) const; + /// @brief Get the pass with a given revolution number. + /// + /// @param aRevolutionNumber Revolution number to get the pass for. + /// @param aStepDuration Step duration to use. Defaults to Duration::Minutes(10.0). + /// @return Pass with the given revolution number. Pass getPassWithRevolutionNumber( const Integer& aRevolutionNumber, const Duration& aStepDuration = Duration::Minutes(10.0) ) const; + /// @brief Get all passes within a given interval. + /// + /// @param anInterval Interval to get the passes within. + /// @return Array of passes within the given interval. + Array getPassesWithinInterval(const Interval& anInterval) const; + + /// @brief Get the orbital frame of a given type. + /// + /// @param aFrameType Type of the frame to get. + /// @return Shared pointer to the orbital frame of the given type. Shared getOrbitalFrame(const Orbit::FrameType& aFrameType) const; + /// @brief Print the Orbit to an output stream. + /// + /// @param anOutputStream Output stream to print to. + /// @param displayDecorator Whether to display the decorator. virtual void print(std::ostream& anOutputStream, bool displayDecorator = true) const override; /// @brief Constructs an undefined orbit /// - /// Undefined orbit + /// @return Undefined orbit static Orbit Undefined(); /// @brief Constructs a circular orbit /// - /// Model: Kepler (No Perturbation). + /// Model: Kepler (No Perturbation). /// /// @param anEpoch An orbit epoch /// @param anAltitude An orbit altitude (wrt. equatorial radius) @@ -133,7 +190,7 @@ class Orbit : public Trajectory /// @brief Constructs an equatorial orbit /// - /// Model: Kepler (No Perturbation). + /// Model: Kepler (No Perturbation). /// /// @param anEpoch An orbit epoch /// @param anApoapsisAltitude An orbit apoapsis altitude (wrt. equatorial radius) @@ -149,7 +206,7 @@ class Orbit : public Trajectory /// @brief Constructs a circular-equatorial orbit /// - /// Model: Kepler (No Perturbation). + /// Model: Kepler (No Perturbation). /// /// @param anEpoch An orbit epoch /// @param anAltitude An orbit altitude (wrt. equatorial radius) @@ -161,7 +218,7 @@ class Orbit : public Trajectory /// @brief Constructs a geosynchronous orbit /// - /// Model: Kepler (J2 Perturbation). + /// Model: Kepler (J2 Perturbation). /// /// @param anEpoch An orbit epoch /// @param anInclination An orbit inclination @@ -177,7 +234,7 @@ class Orbit : public Trajectory /// @brief Constructs a Sun-synchronous orbit /// - /// Model: Kepler (J2 Perturbation). + /// Model: Kepler (J2 Perturbation). /// /// @param anEpoch An orbit epoch /// @param anAltitude An orbit altitude (wrt. equatorial radius) @@ -193,12 +250,36 @@ class Orbit : public Trajectory const Angle& anArgumentOfLatitude = Angle::Zero() ); + /// @brief Get the string representation of a frame type + /// + /// @param aFrameType Type of the frame to get the string representation of + /// @return String representation of the frame type static String StringFromFrameType(const Orbit::FrameType& aFrameType); + /// @brief Compute passes for a given state array and initial revolution number + /// + /// @param aStateArray Array of states + /// @param anInitialRevolutionNumber Initial revolution number + /// @return Array of Pair of Index + Passes static Array> ComputePasses( const Array& aStateArray, const Integer& anInitialRevolutionNumber ); + /// @brief Compute passes for a given model and interval. + /// Note: An Interval object is not used, as the start and end Instant can be forward or backward intervals. + /// + /// @param aModel An Orbit model. + /// @param aStartInstant Start instant of the interval. + /// @param anEndInstant End instant of the interval. + /// @param anInitialRevolutionNumber Initial revolution number at the start of the interval. + /// @return Array of passes + static Array ComputePassesWithModel( + const orbit::Model& aModel, + const Instant& aStartInstant, + const Instant& anEndInstant, + const Integer& anInitialRevolutionNumber + ); + private: const orbit::Model* modelPtr_; @@ -213,17 +294,35 @@ class Orbit : public Trajectory /// Use a bisection search to find the Instant between `previousInstant` and `nextInstant` at which the return value /// of `getValue` crosses zero. /// - /// @param anEpoch An orbit epoch from which to measure - /// @param previousInstant lower bound of search - /// @param currentInstant upper bound of search - /// @param getValue function that receives a duration [s] from `anEpoch` and returns a real value - /// @return Instant at which the crossing occurs + /// @param anEpoch An orbit epoch from which to measure. + /// @param previousInstant lower bound of search. + /// @param currentInstant upper bound of search. + /// @param getValue function that receives a duration [s] from `anEpoch` and returns a real value. + /// @return Instant at which the crossing occurs. static Instant GetCrossingInstant( const Instant& anEpoch, const Instant& previousInstant, const Instant& currentInstant, const std::function& getValue ); + + /// @brief Compute the instants at which the orbit crosses the north point, descending node, south point and pass + /// break + /// + /// @param aModel An orbit model. + /// @param previousInstant The previous instant. + /// @param stepDuration The step duration. + /// @param isForwardPropagated Whether the orbit is forward propagated. + /// @param anEndInstant The end instant. + /// @return Tuple of instants at which the orbit crosses the north point, descending node, south point and pass + /// break. + static Tuple ComputeCrossings( + const orbit::Model& aModel, + Instant previousInstant, + const Duration& stepDuration, + const bool& isForwardPropagated, + const Instant& anEndInstant = Instant::Undefined() + ); }; } // namespace trajectory diff --git a/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model.hpp b/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model.hpp index cbbe71ad6..6daa07238 100644 --- a/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model.hpp +++ b/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model.hpp @@ -42,7 +42,7 @@ class Model : public virtual ostk::astrodynamics::trajectory::Model virtual State calculateStateAt(const Instant& anInstant) const = 0; - virtual Integer calculateRevolutionNumberAt(const Instant& anInstant) const = 0; // [TBR] + virtual Integer calculateRevolutionNumberAt(const Instant& anInstant) const; virtual void print(std::ostream& anOutputStream, bool displayDecorator = true) const = 0; }; diff --git a/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/SGP4.hpp b/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/SGP4.hpp index 04156ebe4..bb07d0547 100644 --- a/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/SGP4.hpp +++ b/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/SGP4.hpp @@ -70,8 +70,6 @@ class SGP4 : public ostk::astrodynamics::trajectory::orbit::Model virtual State calculateStateAt(const Instant& anInstant) const override; - virtual Integer calculateRevolutionNumberAt(const Instant& anInstant) const override; // [TBR] ? - virtual void print(std::ostream& anOutputStream, bool displayDecorator = true) const override; protected: diff --git a/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Tabulated.hpp b/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Tabulated.hpp index 608db6f4b..e3bdf66ba 100644 --- a/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Tabulated.hpp +++ b/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Tabulated.hpp @@ -58,8 +58,6 @@ class Tabulated : public virtual trajectory::orbit::Model, public trajectory::mo virtual State calculateStateAt(const Instant& anInstant) const override; - virtual Integer calculateRevolutionNumberAt(const Instant& anInstant) const override; - virtual void print(std::ostream& anOutputStream, bool displayDecorator) const override; protected: diff --git a/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Pass.hpp b/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Pass.hpp index 81f38d5fd..f74c9caf1 100644 --- a/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Pass.hpp +++ b/include/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Pass.hpp @@ -8,6 +8,7 @@ #include #include +#include namespace ostk { @@ -23,6 +24,7 @@ using ostk::core::type::String; using ostk::physics::time::Duration; using ostk::physics::time::Instant; +using ostk::physics::time::Interval; /// @class Pass /// @brief A revolution of an orbiting object. @@ -123,6 +125,11 @@ class Pass /// @return The end instant of the pass. Instant getEndInstant() const; + /// @brief Gets the interval of the pass. + /// + /// @return The interval of the pass. + Interval getInterval() const; + /// @brief Accesses the instant at the ascending node of the pass. /// /// @return The instant at the ascending node of the pass. diff --git a/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.cpp b/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.cpp index 741909d95..5d1c8c3fe 100644 --- a/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.cpp +++ b/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.cpp @@ -43,8 +43,6 @@ using ostk::mathematics::geometry::d3::transformation::rotation::RotationMatrix; using ostk::mathematics::geometry::d3::transformation::rotation::RotationVector; using ostk::mathematics::object::Vector3d; -using orbit::model::Kepler; -using orbit::model::kepler::COE; using ostk::physics::environment::object::celestial::Sun; using ostk::physics::time::Duration; using ostk::physics::time::Interval; @@ -54,6 +52,8 @@ using ostk::physics::unit::Length; using ostk::physics::unit::Mass; using ostk::astrodynamics::RootSolver; +using ostk::astrodynamics::trajectory::orbit::model::Kepler; +using ostk::astrodynamics::trajectory::orbit::model::kepler::COE; static const Derived::Unit GravitationalParameterSIUnit = Derived::Unit::GravitationalParameter(Length::Unit::Meter, ostk::physics::unit::Time::Unit::Second); @@ -181,7 +181,7 @@ Pass Orbit::getPassWithRevolutionNumber(const Integer& aRevolutionNumber, const const std::lock_guard lock {this->mutex_}; - const auto getClosestPass = [this](const Integer& aRevolutionNumber) -> Pass + const auto getClosestPass = [this](const Integer& revolutionNumber) -> Pass { if (this->passMap_.empty()) { @@ -190,15 +190,15 @@ Pass Orbit::getPassWithRevolutionNumber(const Integer& aRevolutionNumber, const // exact revolution number exists - if (this->passMap_.count(aRevolutionNumber)) + if (this->passMap_.count(revolutionNumber)) { - return this->passMap_.at(aRevolutionNumber); + return this->passMap_.at(revolutionNumber); } - const auto lowerBoundMapIt = this->passMap_.lower_bound(aRevolutionNumber); + const auto lowerBoundMapIt = this->passMap_.lower_bound(revolutionNumber); // Revolution number is greater than any existing revolution number in map - // {5, 6, 9, 10} -> aRevolutionNumber=12 -> return 10 + // {5, 6, 9, 10} -> revolutionNumber=12 -> return 10 if (lowerBoundMapIt == this->passMap_.end()) { @@ -206,7 +206,7 @@ Pass Orbit::getPassWithRevolutionNumber(const Integer& aRevolutionNumber, const } // Revolution number is lesser than any existing revolution number in map - // {5, 6, 9, 10} -> aRevolutionNumber=4 -> return 5 + // {5, 6, 9, 10} -> revolutionNumber=4 -> return 5 if (lowerBoundMapIt == this->passMap_.begin()) { @@ -217,15 +217,15 @@ Pass Orbit::getPassWithRevolutionNumber(const Integer& aRevolutionNumber, const const auto closestPassMapIt = std::prev(lowerBoundMapIt); - // {5, 6, 9, 10} -> aRevolutionNumber=7 -> return 6 + // {5, 6, 9, 10} -> revolutionNumber=7 -> return 6 // lowerBoundMapIt = 9, closestPassMapIt = 6 - if ((aRevolutionNumber - closestPassMapIt->first) < (lowerBoundMapIt->first - aRevolutionNumber)) + if ((revolutionNumber - closestPassMapIt->first) < (lowerBoundMapIt->first - revolutionNumber)) { return closestPassMapIt->second; } - // {5, 6, 9, 10} -> aRevolutionNumber=8 -> return 9 + // {5, 6, 9, 10} -> revolutionNumber=8 -> return 9 // lowerBoundMapIt = 9, closestPassMapIt = 6 return lowerBoundMapIt->second; @@ -239,114 +239,19 @@ Pass Orbit::getPassWithRevolutionNumber(const Integer& aRevolutionNumber, const const bool isForwardPropagated = currentRevolutionNumber <= aRevolutionNumber; const Integer propagationSign = isForwardPropagated ? 1 : -1; - const auto computeCrossings = [this, &isForwardPropagated]( - Instant previousInstant, const Duration& stepDuration - ) -> Tuple - { - Instant northPointCrossing = Instant::Undefined(); - Instant southPointCrossing = Instant::Undefined(); - Instant descendingNodeCrossing = Instant::Undefined(); - Instant passBreakCrossing = Instant::Undefined(); - - const Instant epoch = this->modelPtr_->getEpoch(); - - const State previousState = this->modelPtr_->calculateStateAt(previousInstant); - Real previousStateCoordinates_ECI_z = previousState.getPosition().accessCoordinates().z(); - Real previousStateCoordinates_ECI_zdot = previousState.getVelocity().accessCoordinates().z(); - - const auto getZ = [this, &epoch](const double& aDurationInSeconds) -> Real - { - return this->modelPtr_->calculateStateAt(epoch + Duration::Seconds(aDurationInSeconds)) - .getPosition() - .accessCoordinates() - .z(); - }; - - const auto getZDot = [this, &epoch](const double& aDurationInSeconds) -> Real - { - return this->modelPtr_->calculateStateAt(epoch + Duration::Seconds(aDurationInSeconds)) - .getVelocity() - .accessCoordinates() - .z(); - }; - - while (!passBreakCrossing.isDefined()) - { - const Instant currentInstant = previousInstant + stepDuration; - - const State currentState = this->modelPtr_->calculateStateAt(currentInstant); - const Real currentStateCoordinates_ECI_z = currentState.getPosition().accessCoordinates().z(); - const Real currentStateCoordinates_ECI_zdot = currentState.getVelocity().accessCoordinates().z(); - - if ((previousStateCoordinates_ECI_z == 0.0) && (currentStateCoordinates_ECI_z == 0.0)) - { - throw ostk::core::error::runtime::ToBeImplemented("Equatorial orbit support."); - } - - if (((previousStateCoordinates_ECI_zdot > 0.0) && (currentStateCoordinates_ECI_zdot <= 0.0)) || - ((previousStateCoordinates_ECI_zdot < 0.0) && (currentStateCoordinates_ECI_zdot >= 0.0))) - { - if (currentStateCoordinates_ECI_z > 0.0) - { - northPointCrossing = Orbit::GetCrossingInstant(epoch, previousInstant, currentInstant, getZDot); - } - else - { - southPointCrossing = Orbit::GetCrossingInstant(epoch, previousInstant, currentInstant, getZDot); - } - } - - if ((previousStateCoordinates_ECI_z > 0.0) && (currentStateCoordinates_ECI_z <= 0.0)) - { - const Instant crossingInstant = Orbit::GetCrossingInstant(epoch, previousInstant, currentInstant, getZ); - - if (isForwardPropagated) - { - descendingNodeCrossing = crossingInstant; - } - else - { - passBreakCrossing = crossingInstant; - } - } - - if ((previousStateCoordinates_ECI_z < 0.0) && (currentStateCoordinates_ECI_z >= 0.0)) - { - const Instant crossingInstant = Orbit::GetCrossingInstant(epoch, previousInstant, currentInstant, getZ); - if (isForwardPropagated) - { - passBreakCrossing = crossingInstant; - } - else - { - descendingNodeCrossing = crossingInstant; - } - } - - previousStateCoordinates_ECI_z = currentStateCoordinates_ECI_z; - previousStateCoordinates_ECI_zdot = currentStateCoordinates_ECI_zdot; - previousInstant = currentInstant; - } - - return { - northPointCrossing, - descendingNodeCrossing, - southPointCrossing, - passBreakCrossing, - }; - }; - Instant northPointCrossing = Instant::Undefined(); Instant southPointCrossing = Instant::Undefined(); Instant descendingNodeCrossing = Instant::Undefined(); Instant passBreakCrossing = Instant::Undefined(); + Duration stepDuration = aStepDuration; + if (!currentPass.isDefined()) { const Instant previousInstant = this->modelPtr_->getEpoch(); std::tie(northPointCrossing, descendingNodeCrossing, southPointCrossing, passBreakCrossing) = - computeCrossings(previousInstant, aStepDuration * Real::Integer(propagationSign)); + Orbit::ComputeCrossings(*this->modelPtr_, previousInstant, stepDuration, isForwardPropagated); const Instant crossingInstant = Real(this->modelPtr_->calculateStateAt(this->modelPtr_->getEpoch()).getPosition().accessCoordinates().z()) @@ -374,8 +279,6 @@ Pass Orbit::getPassWithRevolutionNumber(const Integer& aRevolutionNumber, const while ((currentPass.getRevolutionNumber() != aRevolutionNumber)) { - Duration stepDuration = aStepDuration; - if (currentPass.isComplete()) { Array durations = { @@ -392,7 +295,7 @@ Pass Orbit::getPassWithRevolutionNumber(const Integer& aRevolutionNumber, const : currentPass.getStartInstant() - Duration::Microseconds(1.0); std::tie(northPointCrossing, descendingNodeCrossing, southPointCrossing, passBreakCrossing) = - computeCrossings(previousInstant, stepDuration * Real::Integer(propagationSign)); + Orbit::ComputeCrossings(*this->modelPtr_, previousInstant, stepDuration, isForwardPropagated); currentRevolutionNumber += propagationSign; @@ -415,6 +318,55 @@ Pass Orbit::getPassWithRevolutionNumber(const Integer& aRevolutionNumber, const return currentPass; } +Array Orbit::getPassesWithinInterval(const Interval& anInterval) const +{ + if (!this->isDefined()) + { + throw ostk::core::error::runtime::Undefined("Orbit"); + } + + if (!anInterval.isDefined()) + { + throw ostk::core::error::runtime::Undefined("Interval"); + } + + Instant currentInstant = anInterval.getStart(); + + Integer revolutionNumber = Integer::Undefined(); + + // Check if the Pass is already within Pass map + for (const auto& passPair : this->passMap_) + { + const Pass& pass = passPair.second; + + if (pass.getInterval().contains(currentInstant)) + { + revolutionNumber = pass.getRevolutionNumber(); + break; + } + } + + // If not, then calculate using the model method + if (!revolutionNumber.isDefined()) + { + revolutionNumber = this->modelPtr_->calculateRevolutionNumberAt(currentInstant); + } + + Array passes = {}; + + while (currentInstant <= anInterval.accessEnd()) + { + const Pass pass = this->getPassWithRevolutionNumber(revolutionNumber); + + passes.add(pass); + + currentInstant = pass.accessInstantAtPassBreak(); + revolutionNumber++; + } + + return passes; +} + Shared Orbit::getOrbitalFrame(const Orbit::FrameType& aFrameType) const { using ostk::mathematics::geometry::d3::transformation::rotation::Quaternion; @@ -1177,10 +1129,6 @@ Array> Orbit::ComputePasses(const Array& aStateArray, c throw ostk::core::error::runtime::Undefined("Initial revolution number"); } - // [TBI] Deal w/ z_ECI always equal to 0.0 case (equatorial orbit) - - Array> passMap; - if (aStateArray.getSize() < 2) { throw ostk::core::error::RuntimeError( @@ -1196,121 +1144,128 @@ Array> Orbit::ComputePasses(const Array& aStateArray, c } } - const model::Tabulated tabulated = model::Tabulated(aStateArray, Interpolator::Type::BarycentricRational); + const orbit::model::Tabulated tabulated = + orbit::model::Tabulated(aStateArray, anInitialRevolutionNumber, Interpolator::Type::BarycentricRational); + + const Array passes = Orbit::ComputePassesWithModel( + tabulated, + aStateArray.accessFirst().accessInstant(), + aStateArray.accessLast().accessInstant(), + anInitialRevolutionNumber + ); + + Array> passMap; + passMap.reserve(passes.getSize()); - const Instant& epoch = aStateArray.accessFirst().accessInstant(); + Index stateIndex = 0; + Index passIndex = 0; - const auto getZ = [&tabulated, &epoch](const double& aDurationInSeconds) -> double + while (stateIndex < aStateArray.getSize()) { - return tabulated.calculateStateAt(epoch + Duration::Seconds(aDurationInSeconds)) - .getPosition() - .accessCoordinates() - .z(); - }; + const Instant& stateInstant = aStateArray[stateIndex].accessInstant(); + + if (passes[passIndex].accessInstantAtPassBreak().isDefined() && + stateInstant >= passes[passIndex].accessInstantAtPassBreak()) + { + passMap.add({stateIndex, passes[passIndex]}); + ++passIndex; + } + + stateIndex++; + } - const auto getZDot = [&tabulated, &epoch](const double& aDurationInSeconds) -> double + passMap.add({stateIndex, passes.accessLast()}); + + return passMap; +} + +Array Orbit::ComputePassesWithModel( + const orbit::Model& aModel, + const Instant& aStartInstant, + const Instant& anEndInstant, + const Integer& aRevolutionNumber +) +{ + // [TBI] Deal w/ z_ECI always equal to 0.0 case (equatorial orbit) + + if (!aModel.isDefined()) { - return tabulated.calculateStateAt(epoch + Duration::Seconds(aDurationInSeconds)) - .getVelocity() - .accessCoordinates() - .z(); - }; + throw ostk::core::error::runtime::Undefined("Model"); + } + + if (!aRevolutionNumber.isDefined()) + { + throw ostk::core::error::runtime::Undefined("Revolution number"); + } + + Array passes = Array::Empty(); - Index currentIndex = 0; - State const* previousStatePtr = nullptr; + Integer revolutionNumber = aRevolutionNumber; + Instant previousInstant = aStartInstant; + State previousState = aModel.calculateStateAt(previousInstant); - Integer revolutionNumber = anInitialRevolutionNumber; + Instant previousPassEndInstant = (Real(previousState.getPosition().accessCoordinates().z()).isNear(0.0, epsilon)) + ? previousState.accessInstant() + : Instant::Undefined(); - Instant previousPassEndInstant = - (Real(aStateArray.accessFirst().getPosition().accessCoordinates().z()).isNear(0.0, epsilon)) - ? aStateArray.accessFirst().accessInstant() - : Instant::Undefined(); Instant northPointCrossing = Instant::Undefined(); Instant descendingNodeCrossing = Instant::Undefined(); Instant southPointCrossing = Instant::Undefined(); - Instant passBreak = Instant::Undefined(); + Instant passBreakCrossing = Instant::Undefined(); - for (const auto& state : aStateArray) + const bool isForwardPropagated = aStartInstant < anEndInstant; + Duration stepDuration = Duration::Minutes(5.0); + + while (true) { - if (previousStatePtr != nullptr) + if ((isForwardPropagated && (previousInstant >= anEndInstant)) || + (!isForwardPropagated && (previousInstant <= anEndInstant))) { - const Vector3d previousPositionCoordinates_ECI = previousStatePtr->getPosition().accessCoordinates(); - const Vector3d previousVelocityCoordinates_ECI = previousStatePtr->getVelocity().accessCoordinates(); - const Vector3d currentPositionCoordinates_ECI = state.getPosition().accessCoordinates(); - const Vector3d currentVelocityCoordinates_ECI = state.getVelocity().accessCoordinates(); - - // North & South point crossings - if (((previousVelocityCoordinates_ECI.z() > 0.0) && (currentVelocityCoordinates_ECI.z() <= 0.0)) || - ((previousVelocityCoordinates_ECI.z() < 0.0) && (currentVelocityCoordinates_ECI.z() >= 0.0))) - { - if (currentPositionCoordinates_ECI.z() > 0.0) - { - northPointCrossing = Orbit::GetCrossingInstant( - epoch, previousStatePtr->accessInstant(), state.accessInstant(), getZDot - ); - } - else - { - southPointCrossing = Orbit::GetCrossingInstant( - epoch, previousStatePtr->accessInstant(), state.accessInstant(), getZDot - ); - } - } + break; + } - // Descending node - if ((previousPositionCoordinates_ECI.z() > 0.0) && (currentPositionCoordinates_ECI.z() <= 0.0)) - { - descendingNodeCrossing = - Orbit::GetCrossingInstant(epoch, previousStatePtr->accessInstant(), state.accessInstant(), getZ); - } + std::tie(northPointCrossing, descendingNodeCrossing, southPointCrossing, passBreakCrossing) = + Orbit::ComputeCrossings(aModel, previousInstant, stepDuration, isForwardPropagated, anEndInstant); - // Pass break - if ((previousPositionCoordinates_ECI.z() < 0.0) && (currentPositionCoordinates_ECI.z() >= 0.0)) - { - passBreak = - Orbit::GetCrossingInstant(epoch, previousStatePtr->accessInstant(), state.accessInstant(), getZ); - } + Pass pass = { + revolutionNumber, + isForwardPropagated ? previousPassEndInstant : passBreakCrossing, + northPointCrossing, + descendingNodeCrossing, + southPointCrossing, + isForwardPropagated ? passBreakCrossing : previousPassEndInstant, + }; - if (passBreak.isDefined()) - { - const Pass pass = { - revolutionNumber, - previousPassEndInstant, - northPointCrossing, - descendingNodeCrossing, - southPointCrossing, - passBreak, - }; + passes.add(pass); - passMap.add({currentIndex, pass}); + // if forward propagating and the pass has an undefined PassBreak, then we have crossed the end of the interval + // if backward propagating and the pass has an undefined AscendingNode, then we have crossed the end of the + // interval + if ((isForwardPropagated && !pass.accessInstantAtPassBreak().isDefined()) || + (!isForwardPropagated && !pass.accessInstantAtAscendingNode().isDefined())) + { + break; + } - revolutionNumber++; - previousPassEndInstant = passBreak; + revolutionNumber += isForwardPropagated ? 1 : -1; + previousPassEndInstant = passBreakCrossing; - southPointCrossing = Instant::Undefined(); - northPointCrossing = Instant::Undefined(); - descendingNodeCrossing = Instant::Undefined(); - passBreak = Instant::Undefined(); - } + if (pass.isComplete()) + { + Array durations = { + (pass.accessInstantAtNorthPoint() - pass.accessInstantAtAscendingNode()), + (pass.accessInstantAtDescendingNode() - pass.accessInstantAtNorthPoint()), + (pass.accessInstantAtSouthPoint() - pass.accessInstantAtDescendingNode()), + (pass.accessInstantAtPassBreak() - pass.accessInstantAtSouthPoint()), + }; + stepDuration = *std::min_element(durations.begin(), durations.end()) / 2.0; } - previousStatePtr = &state; - currentIndex++; + previousInstant = isForwardPropagated ? pass.accessInstantAtPassBreak() + Duration::Microseconds(1.0) + : pass.accessInstantAtAscendingNode() - Duration::Microseconds(1.0); } - // Add last partial pass - const Pass pass = { - revolutionNumber, - previousPassEndInstant, - northPointCrossing, - descendingNodeCrossing, - southPointCrossing, - passBreak, - }; - - passMap.add({currentIndex, pass}); - - return passMap; + return passes; } Instant Orbit::GetCrossingInstant( @@ -1331,6 +1286,130 @@ Instant Orbit::GetCrossingInstant( return anEpoch + Duration::Seconds(solution.root); } +Tuple Orbit::ComputeCrossings( + const orbit::Model& aModel, + Instant previousInstant, + const Duration& aStepDuration, + const bool& isForwardPropagated, + const Instant& anEndInstant +) +{ + Instant northPointCrossing = Instant::Undefined(); + Instant southPointCrossing = Instant::Undefined(); + Instant descendingNodeCrossing = Instant::Undefined(); + Instant passBreakCrossing = Instant::Undefined(); + + if (anEndInstant.isDefined()) + { + if ((isForwardPropagated && (previousInstant >= anEndInstant)) || + (!isForwardPropagated && (previousInstant <= anEndInstant))) + { + return {northPointCrossing, descendingNodeCrossing, southPointCrossing, passBreakCrossing}; + } + } + + const Instant epoch = aModel.getEpoch(); + + const State previousState = aModel.calculateStateAt(previousInstant); + Real previousStateCoordinates_ECI_z = previousState.getPosition().accessCoordinates().z(); + Real previousStateCoordinates_ECI_zdot = previousState.getVelocity().accessCoordinates().z(); + + const auto getZ = [&aModel, &epoch](const double& aDurationInSeconds) -> Real + { + return aModel.calculateStateAt(epoch + Duration::Seconds(aDurationInSeconds)) + .getPosition() + .accessCoordinates() + .z(); + }; + + const auto getZDot = [&aModel, &epoch](const double& aDurationInSeconds) -> Real + { + return aModel.calculateStateAt(epoch + Duration::Seconds(aDurationInSeconds)) + .getVelocity() + .accessCoordinates() + .z(); + }; + + const Duration stepDuration = isForwardPropagated ? aStepDuration : -aStepDuration; + + bool lastStep = false; + + while (!passBreakCrossing.isDefined() && !lastStep) + { + Instant currentInstant = previousInstant + stepDuration; + + if (anEndInstant.isDefined()) + { + if ((isForwardPropagated && (currentInstant >= anEndInstant)) || + (!isForwardPropagated && (currentInstant <= anEndInstant))) + { + currentInstant = anEndInstant; + lastStep = true; + } + } + + const State currentState = aModel.calculateStateAt(currentInstant); + const Real currentStateCoordinates_ECI_z = currentState.getPosition().accessCoordinates().z(); + const Real currentStateCoordinates_ECI_zdot = currentState.getVelocity().accessCoordinates().z(); + + if ((previousStateCoordinates_ECI_z == 0.0) && (currentStateCoordinates_ECI_z == 0.0)) + { + throw ostk::core::error::runtime::ToBeImplemented("Equatorial orbit support."); + } + + if (((previousStateCoordinates_ECI_zdot > 0.0) && (currentStateCoordinates_ECI_zdot <= 0.0)) || + ((previousStateCoordinates_ECI_zdot < 0.0) && (currentStateCoordinates_ECI_zdot >= 0.0))) + { + if (currentStateCoordinates_ECI_z > 0.0) + { + northPointCrossing = Orbit::GetCrossingInstant(epoch, previousInstant, currentInstant, getZDot); + } + else + { + southPointCrossing = Orbit::GetCrossingInstant(epoch, previousInstant, currentInstant, getZDot); + } + } + + if ((previousStateCoordinates_ECI_z > 0.0) && (currentStateCoordinates_ECI_z <= 0.0)) + { + const Instant crossingInstant = Orbit::GetCrossingInstant(epoch, previousInstant, currentInstant, getZ); + + if (isForwardPropagated) + { + descendingNodeCrossing = crossingInstant; + } + else + { + passBreakCrossing = crossingInstant; + } + } + + if ((previousStateCoordinates_ECI_z < 0.0) && (currentStateCoordinates_ECI_z >= 0.0)) + { + const Instant crossingInstant = Orbit::GetCrossingInstant(epoch, previousInstant, currentInstant, getZ); + if (isForwardPropagated) + { + passBreakCrossing = crossingInstant; + } + else + { + descendingNodeCrossing = crossingInstant; + } + } + + previousStateCoordinates_ECI_z = currentStateCoordinates_ECI_z; + previousStateCoordinates_ECI_zdot = currentStateCoordinates_ECI_zdot; + previousInstant = currentInstant; + } + + return { + northPointCrossing, + descendingNodeCrossing, + southPointCrossing, + passBreakCrossing, + }; +} + } // namespace trajectory } // namespace astrodynamics } // namespace ostk diff --git a/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model.cpp b/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model.cpp index 38f848237..17a73ce96 100644 --- a/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model.cpp +++ b/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model.cpp @@ -1,6 +1,15 @@ /// Apache License 2.0 +#include +#include +#include + +#include +#include + +#include #include +#include namespace ostk { @@ -11,10 +20,46 @@ namespace trajectory namespace orbit { +using ostk::core::container::Array; + +using ostk::physics::time::Instant; +using ostk::physics::time::Interval; + +using ostk::astrodynamics::trajectory::Orbit; +using ostk::astrodynamics::trajectory::orbit::Pass; + Model::Model() {} Model::~Model() {} +Integer Model::calculateRevolutionNumberAt(const Instant& anInstant) const +{ + if (!anInstant.isDefined()) + { + throw ostk::core::error::runtime::Undefined("Instant"); + } + + if (!this->isDefined()) + { + throw ostk::core::error::runtime::Undefined("Model"); + } + + if (anInstant == this->getEpoch()) + { + return this->getRevolutionNumberAtEpoch(); + } + + const Array passes = + Orbit::ComputePassesWithModel(*this, this->getEpoch(), anInstant, this->getRevolutionNumberAtEpoch()); + + if (passes.isEmpty()) + { + throw ostk::core::error::RuntimeError("No pass found"); + } + + return passes.accessLast().getRevolutionNumber(); +} + } // namespace orbit } // namespace trajectory } // namespace astrodynamics diff --git a/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/SGP4.cpp b/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/SGP4.cpp index ba82ad2d5..4d551af09 100644 --- a/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/SGP4.cpp +++ b/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/SGP4.cpp @@ -184,26 +184,6 @@ State SGP4::calculateStateAt(const Instant& anInstant) const return this->implUPtr_->calculateStateAt(anInstant); } -Integer SGP4::calculateRevolutionNumberAt(const Instant& anInstant) const -{ - if (!anInstant.isDefined()) - { - throw ostk::core::error::runtime::Undefined("Instant"); - } - - if (!this->isDefined()) - { - throw ostk::core::error::runtime::Undefined("SGP4"); - } - - if (anInstant == this->tle_.getEpoch()) - { - return this->getRevolutionNumberAtEpoch(); - } - - return Integer::Undefined(); -} - void SGP4::print(std::ostream& anOutputStream, bool displayDecorator) const { displayDecorator ? ostk::core::utils::Print::Header(anOutputStream, "SGP4") : void(); diff --git a/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Tabulated.cpp b/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Tabulated.cpp index 880631772..097c51e8d 100644 --- a/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Tabulated.cpp +++ b/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Tabulated.cpp @@ -78,25 +78,6 @@ State Tabulated::calculateStateAt(const Instant& anInstant) const return trajectory::model::Tabulated::calculateStateAt(anInstant); } -Integer Tabulated::calculateRevolutionNumberAt(const Instant& anInstant) const -{ - if (!anInstant.isDefined()) - { - throw ostk::core::error::runtime::Undefined("Instant"); - } - - if (!this->isDefined()) - { - throw ostk::core::error::runtime::Undefined("Tabulated"); - } - - throw ostk::core::error::runtime::ToBeImplemented("Kepler::calculateRevolutionNumberAt"); - - // [TBI] - - return Integer::Undefined(); -} - void Tabulated::print(std::ostream& anOutputStream, bool displayDecorator) const { trajectory::model::Tabulated::print(anOutputStream, displayDecorator); diff --git a/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Pass.cpp b/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Pass.cpp index 51e6e189c..232af17df 100644 --- a/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Pass.cpp +++ b/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Pass.cpp @@ -247,6 +247,21 @@ Instant Pass::getEndInstant() const return instantAtAscendingNode_; } +Interval Pass::getInterval() const +{ + if (!this->isDefined()) + { + throw ostk::core::error::runtime::Undefined("Pass"); + } + + if (type_ == Pass::Type::Complete) + { + return Interval::Closed(instantAtAscendingNode_, instantAtPassBreak_); + } + + return Interval::Undefined(); +} + const Instant& Pass::accessInstantAtAscendingNode() const { if (!this->isDefined()) diff --git a/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.test.cpp b/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.test.cpp index 7c06a126f..f7d547643 100644 --- a/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.test.cpp +++ b/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.test.cpp @@ -25,9 +25,12 @@ #include #include #include +#include #include #include #include +#include +#include #include @@ -67,10 +70,13 @@ using EarthGravitationalModel = ostk::physics::environment::gravitational::Earth using ostk::astrodynamics::trajectory::Orbit; using ostk::astrodynamics::trajectory::orbit::model::Kepler; using ostk::astrodynamics::trajectory::orbit::model::kepler::COE; +using ostk::astrodynamics::trajectory::orbit::model::Propagated; using ostk::astrodynamics::trajectory::orbit::model::SGP4; using ostk::astrodynamics::trajectory::orbit::model::sgp4::TLE; using ostk::astrodynamics::trajectory::orbit::Pass; +using ostk::astrodynamics::trajectory::Propagator; using ostk::astrodynamics::trajectory::State; +using ostk::astrodynamics::trajectory::state::NumericalSolver; TEST(OpenSpaceToolkit_Astrodynamics_Trajectory_Orbit, Constructor) { @@ -570,6 +576,263 @@ TEST(OpenSpaceToolkit_Astrodynamics_Trajectory_Orbit, ComputePasses) } } +TEST(OpenSpaceToolkit_Astrodynamics_Trajectory_Orbit, ComputePassesWithModel) +{ + // undefined model + { + const Kepler kepler = { + COE::Undefined(), + Instant::Undefined(), + Earth::Default(), + Kepler::PerturbationType::None, + }; + const Instant startInstant = Instant::J2000(); + const Instant endInstant = startInstant + Duration::Days(1.0); + EXPECT_THROW( + Orbit::ComputePassesWithModel(kepler, startInstant, endInstant, 1), ostk::core::error::runtime::Undefined + ); + } + + // undefined revolution number + { + const TLE tle = + TLE("1 43890U 18111Q 20195.55622117 +.00000241 +00000-0 +28512-4 0 9991", + "2 43890 097.6899 099.9703 0003551 181.1072 179.0140 14.92932318084236"); + const SGP4 sgp4 = SGP4(tle); + const Instant startInstant = Instant::J2000(); + const Instant endInstant = startInstant + Duration::Days(1.0); + EXPECT_THROW( + Orbit::ComputePassesWithModel(sgp4, startInstant, endInstant, Integer::Undefined()), + ostk::core::error::runtime::Undefined + ); + } + + // Kepler model + { + // Orbit setup + + const Length semiMajorAxis = Length::Kilometers(7000.0); + const Real eccentricity = 0.0; + const Angle inclination = Angle::Degrees(45.0); + const Angle raan = Angle::Degrees(0.0); + const Angle aop = Angle::Degrees(0.0); + const Angle trueAnomaly = Angle::Degrees(0.0); + + const COE coe = {semiMajorAxis, eccentricity, inclination, raan, aop, trueAnomaly}; + + const Instant epoch = Instant::DateTime(DateTime(2018, 1, 1, 0, 0, 0), Scale::UTC); + const Derived gravitationalParameter = EarthGravitationalModel::EGM2008.gravitationalParameter_; + const Length equatorialRadius = EarthGravitationalModel::EGM2008.equatorialRadius_; + const Real J2 = EarthGravitationalModel::EGM2008.J2_; + const Real J4 = EarthGravitationalModel::EGM2008.J4_; + + const Kepler keplerianModel = { + coe, epoch, gravitationalParameter, equatorialRadius, J2, J4, Kepler::PerturbationType::None + }; + + // partial + { + // forward + { + const Instant startInstant = Instant::DateTime(DateTime::Parse("2018-01-01 00:00:00"), Scale::UTC); + const Instant endInstant = Instant::DateTime(DateTime::Parse("2018-01-01 02:30:00"), Scale::UTC); + + const Array passes = Orbit::ComputePassesWithModel(keplerianModel, startInstant, endInstant, 1); + + for (const auto &pass : passes) + { + EXPECT_TRUE(pass.isDefined()); + } + EXPECT_EQ(passes.accessLast().getType(), Pass::Type::Partial); + } + + // backward + { + const Instant startInstant = Instant::DateTime(DateTime::Parse("2018-01-01 02:30:00"), Scale::UTC); + const Instant endInstant = Instant::DateTime(DateTime::Parse("2018-01-01 00:00:00"), Scale::UTC); + + const Array passes = Orbit::ComputePassesWithModel(keplerianModel, startInstant, endInstant, 1); + + for (const auto &pass : passes) + { + EXPECT_TRUE(pass.isDefined()); + } + + EXPECT_EQ(passes.accessLast().getType(), Pass::Type::Complete); + } + } + + { + // Reference data setup + + const Table referenceData = Table::Load( + File::Path(Path::Parse("/app/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Kepler/Test_1/" + "Satellite Passes.csv")), + Table::Format::CSV, + true + ); + + // Pass test + const Instant startInstant = Instant::DateTime(DateTime::Parse("2018-01-01 00:00:00"), Scale::UTC); + const Instant endInstant = Instant::DateTime(DateTime::Parse("2018-01-01 23:00:00"), Scale::UTC); + + const Array passes = Orbit::ComputePassesWithModel(keplerianModel, startInstant, endInstant, 1); + + Array reversePasses = Orbit::ComputePassesWithModel(keplerianModel, endInstant, startInstant, 15); + std::reverse(reversePasses.begin(), reversePasses.end()); + + EXPECT_EQ(passes.getSize(), reversePasses.getSize()); + + EXPECT_EQ( + referenceData.getRowCount(), passes.size() - 1 + ); // We're generating 1 pass over the reference data + + Index i = 0; + for (const auto &pass : passes) + { + if (pass.accessInstantAtAscendingNode().isDefined() && + reversePasses[i].accessInstantAtAscendingNode().isDefined()) + { + EXPECT_LT( + (pass.accessInstantAtAscendingNode() - reversePasses[i].accessInstantAtAscendingNode()) + .getAbsolute(), + Duration::Microseconds(1.0) + ); + } + if (pass.accessInstantAtNorthPoint().isDefined() && + reversePasses[i].accessInstantAtNorthPoint().isDefined()) + { + EXPECT_LT( + (pass.accessInstantAtNorthPoint() - reversePasses[i].accessInstantAtNorthPoint()).getAbsolute(), + Duration::Microseconds(1.0) + ); + } + if (pass.accessInstantAtDescendingNode().isDefined() && + reversePasses[i].accessInstantAtDescendingNode().isDefined()) + { + EXPECT_LT( + (pass.accessInstantAtDescendingNode() - reversePasses[i].accessInstantAtDescendingNode()) + .getAbsolute(), + Duration::Microseconds(1.0) + ); + } + if (pass.accessInstantAtSouthPoint().isDefined() && + reversePasses[i].accessInstantAtSouthPoint().isDefined()) + { + EXPECT_LT( + (pass.accessInstantAtSouthPoint() - reversePasses[i].accessInstantAtSouthPoint()).getAbsolute(), + Duration::Microseconds(1.0) + ); + } + if (pass.accessInstantAtPassBreak().isDefined() && + reversePasses[i].accessInstantAtPassBreak().isDefined()) + { + EXPECT_LT( + (pass.accessInstantAtPassBreak() - reversePasses[i].accessInstantAtPassBreak()).getAbsolute(), + Duration::Microseconds(1.0) + ); + } + + // Ignore the lass pass, as it is not complete + if (i > referenceData.getRowCount() - 1) + { + break; + } + + // test computed Pass + const auto &referenceRow = referenceData[i]; + + const Instant referencePassStartInstant = + Instant::DateTime(DateTime::Parse(referenceRow[1].accessString()), Scale::UTC); + const Instant referencePassEndInstant = + Instant::DateTime(DateTime::Parse(referenceRow[2].accessString()), Scale::UTC); + + const Instant referencePassAscendingNodeInstant = + Instant::DateTime(DateTime::Parse(referenceRow[3].accessString()), Scale::UTC); + const Instant referencePassDescendingNodeInstant = + Instant::DateTime(DateTime::Parse(referenceRow[5].accessString()), Scale::UTC); + + const Instant referencePassNorthPointInstant = + Instant::DateTime(DateTime::Parse(referenceRow[4].accessString()), Scale::UTC); + const Instant referencePassSouthPointInstant = + Instant::DateTime(DateTime::Parse(referenceRow[6].accessString()), Scale::UTC); + + EXPECT_TRUE(pass.isDefined()); + + EXPECT_EQ(pass.getRevolutionNumber(), referenceRow[0].accessInteger()); + EXPECT_LT( + std::fabs((referencePassStartInstant - pass.accessInstantAtAscendingNode()).inSeconds()), 1e-6 + ); + EXPECT_LT(std::fabs((referencePassEndInstant - pass.accessInstantAtPassBreak()).inSeconds()), 1e-6); + EXPECT_LT( + std::fabs((referencePassAscendingNodeInstant - pass.accessInstantAtAscendingNode()).inSeconds()), + 1e-6 + ); + EXPECT_LT( + std::fabs((referencePassDescendingNodeInstant - pass.accessInstantAtDescendingNode()).inSeconds()), + 1e-6 + ); + EXPECT_LT( + std::fabs((referencePassNorthPointInstant - pass.accessInstantAtNorthPoint()).inSeconds()), 3.0 + ); + EXPECT_LT( + std::fabs((referencePassSouthPointInstant - pass.accessInstantAtSouthPoint()).inSeconds()), 3.0 + ); + + ++i; + } + } + } + + // SGP4 model + { + const Instant startInstant = Instant::DateTime(DateTime::Parse("2020-07-18 20:38:11.102.455"), Scale::UTC); + const Instant endInstant = Instant::DateTime(DateTime::Parse("2020-07-18 22:26:11.102.455"), Scale::UTC); + + const TLE tle = + TLE("1 43890U 18111Q 20195.55622117 +.00000241 +00000-0 +28512-4 0 9991", + "2 43890 097.6899 099.9703 0003551 181.1072 179.0140 14.92932318084236"); + const SGP4 sgp4 = SGP4(tle); + + const Array passes = Orbit::ComputePassesWithModel(sgp4, startInstant, endInstant, 8423); + + for (const auto &pass : passes) + { + EXPECT_TRUE(pass.isDefined()); + } + + EXPECT_EQ(passes.getSize(), 2); + } + + // Propagated model + { + const Environment environment = {Instant::J2000(), {std::make_shared(Earth::EGM96(10, 10))}}; + const Propagator propagator = Propagator::FromEnvironment(NumericalSolver::Default(), environment); + + const Array states = { + { + Instant::DateTime(DateTime::Parse("2023-01-30T18:30:00.184", DateTime::Format::ISO8601), Scale::UTC), + Position::Meters({-5981016.371280898340, 2114677.607544674072, 2683244.080375305377}, Frame::GCRF()), + Velocity::MetersPerSecond({3142.190382984703, 50.122115740761, 6934.877387038766}, Frame::GCRF()), + }, + }; + + const Propagated propagated = {propagator, states, 5}; + + const Instant startInstant = states.accessFirst().accessInstant() + Duration::Hours(1.0); + const Instant endInstant = states.accessFirst().accessInstant() + Duration::Hours(3.0); + + const Array passes = Orbit::ComputePassesWithModel(propagated, startInstant, endInstant, 1); + + for (const auto &pass : passes) + { + EXPECT_TRUE(pass.isDefined()); + } + + EXPECT_EQ(passes.getSize(), 2); + } +} + TEST(OpenSpaceToolkit_Astrodynamics_Trajectory_Orbit, GetPassWithRevolutionNumber) { { @@ -884,6 +1147,108 @@ TEST(OpenSpaceToolkit_Astrodynamics_Trajectory_Orbit, GetPassWithRevolutionNumbe } } +TEST(OpenSpaceToolkit_Astrodynamics_Trajectory_Orbit, GetPassesWithinInterval) +{ + // undefined + { + // orbit + { + const Interval interval = Interval::Closed(Instant::Now(), Instant::Now() + Duration::Days(1.0)); + EXPECT_THROW(Orbit::Undefined().getPassesWithinInterval(interval), ostk::core::error::runtime::Undefined); + } + + // interval + { + const Orbit orbit = Orbit::SunSynchronous( + Instant::DateTime(DateTime::Parse("2018-01-01 00:00:00"), Scale::UTC), + Length::Kilometers(500.0), + Time::Parse("12:00:00"), + Environment::Default().accessCelestialObjectWithName("Earth") + ); + EXPECT_THROW(orbit.getPassesWithinInterval(Interval::Undefined()), ostk::core::error::runtime::Undefined); + } + } + + // Kepler model + { + // Environment setup + + const Environment environment = Environment::Default(); + + // Orbit setup + + const Length semiMajorAxis = Length::Kilometers(7000.0); + const Real eccentricity = 0.0; + const Angle inclination = Angle::Degrees(45.0); + const Angle raan = Angle::Degrees(0.0); + const Angle aop = Angle::Degrees(0.0); + const Angle trueAnomaly = Angle::Degrees(0.0); + + const COE coe = {semiMajorAxis, eccentricity, inclination, raan, aop, trueAnomaly}; + + const Instant epoch = Instant::DateTime(DateTime(2018, 1, 1, 0, 0, 0), Scale::UTC); + const Derived gravitationalParameter = EarthGravitationalModel::EGM2008.gravitationalParameter_; + const Length equatorialRadius = EarthGravitationalModel::EGM2008.equatorialRadius_; + const Real J2 = EarthGravitationalModel::EGM2008.J2_; + const Real J4 = EarthGravitationalModel::EGM2008.J4_; + + const Kepler keplerianModel = { + coe, epoch, gravitationalParameter, equatorialRadius, J2, J4, Kepler::PerturbationType::None + }; + + const Orbit orbit = {keplerianModel, environment.accessCelestialObjectWithName("Earth")}; + + // Reference data setup + + const Table referenceData = Table::Load( + File::Path(Path::Parse("/app/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Kepler/Test_1/" + "Satellite Passes.csv")), + Table::Format::CSV, + true + ); + + // Pass test + + const Array passes = orbit.getPassesWithinInterval(Interval::Closed( + Instant::DateTime(DateTime::Parse("2018-01-01 00:00:00"), Scale::UTC), + Instant::DateTime(DateTime::Parse("2018-01-01 23:00:00"), Scale::UTC) + )); + + EXPECT_TRUE(passes.getSize() > 0); + EXPECT_EQ(passes.getSize() - 1, referenceData.getRowCount()); + + // Test regenerating with cached passes + EXPECT_NO_THROW(orbit.getPassesWithinInterval(Interval::Closed( + Instant::DateTime(DateTime::Parse("2018-01-01 00:00:00"), Scale::UTC), + Instant::DateTime(DateTime::Parse("2018-01-01 23:00:00"), Scale::UTC) + ))); + + for (const auto &referenceRow : referenceData) + { + const Integer referenceRevolutionNumber = referenceRow[0].accessInteger(); + const Instant referencePassStartInstant = + Instant::DateTime(DateTime::Parse(referenceRow[1].accessString()), Scale::UTC); + const Instant referencePassEndInstant = + Instant::DateTime(DateTime::Parse(referenceRow[2].accessString()), Scale::UTC); + + const Pass pass = orbit.getPassWithRevolutionNumber(referenceRevolutionNumber); + + EXPECT_TRUE(pass.isDefined()); + + EXPECT_EQ(Pass::Type::Complete, pass.getType()); + + EXPECT_GT( + Duration::Microseconds(1.0), + Duration::Between(referencePassStartInstant, pass.accessInstantAtAscendingNode()).getAbsolute() + ); + EXPECT_GT( + Duration::Microseconds(1.0), + Duration::Between(referencePassEndInstant, pass.accessInstantAtPassBreak()).getAbsolute() + ); + } + } +} + TEST(OpenSpaceToolkit_Astrodynamics_Trajectory_Orbit, GetOrbitalFrame) { using ostk::mathematics::geometry::d3::transformation::rotation::RotationMatrix; diff --git a/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Propagated.test.cpp b/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Propagated.test.cpp index 3ca2333b0..9fc364794 100644 --- a/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Propagated.test.cpp +++ b/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Propagated.test.cpp @@ -712,7 +712,7 @@ TEST_F(OpenSpaceToolkit_Astrodynamics_Trajectory_Orbit_Model_Propagated, SetCach TEST_F(OpenSpaceToolkit_Astrodynamics_Trajectory_Orbit_Model_Propagated, PropagationAtExactInstants) { { - Array states = { + const Array states = { { Instant::DateTime(DateTime::Parse("2023-01-30T18:30:00.184", DateTime::Format::ISO8601), Scale::UTC), Position::Meters({-5981016.371280898340, 2114677.607544674072, 2683244.080375305377}, gcrfSPtr_), @@ -730,9 +730,9 @@ TEST_F(OpenSpaceToolkit_Astrodynamics_Trajectory_Orbit_Model_Propagated, Propaga } }; - Propagated propagatedModel(propagator_, states); + const Propagated propagatedModel(propagator_, states); - State estimatedState = propagatedModel.calculateStateAt(states[1].getInstant()); + const State estimatedState = propagatedModel.calculateStateAt(states[1].getInstant()); ASSERT_EQ(estimatedState, states[1]); } diff --git a/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Pass.test.cpp b/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Pass.test.cpp index 91924da9a..de4580153 100644 --- a/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Pass.test.cpp +++ b/test/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Pass.test.cpp @@ -682,6 +682,35 @@ TEST_F(OpenSpaceToolkit_Astrodynamics_Trajectory_Orbit_Pass, GetEndInstant) } } +TEST_F(OpenSpaceToolkit_Astrodynamics_Trajectory_Orbit_Pass, GetInterval) +{ + { + EXPECT_THROW(Pass::Undefined().getInterval(), ostk::core::error::runtime::Undefined); + } + + { + EXPECT_TRUE(defaultPass_.getInterval().isDefined()); + } + + { + const Pass pass = { + defaultRevolutionNumber_, + defaultInstantAtAscendingNode_, + defaultInstantAtNorthPoint_, + defaultInstantAtDescendingNode_, + Instant::Undefined(), + defaultInstantAtPassBreak_, + }; + EXPECT_FALSE(pass.getInterval().isDefined()); + } + + { + EXPECT_EQ( + defaultPass_.getInterval(), Interval::Closed(defaultInstantAtAscendingNode_, defaultInstantAtPassBreak_) + ); + } +} + TEST_F(OpenSpaceToolkit_Astrodynamics_Trajectory_Orbit_Pass, AccessInstantAtAscendingNode) { {