diff --git a/include/OpenSpaceToolkit/Astrodynamics/Access/Generator.hpp b/include/OpenSpaceToolkit/Astrodynamics/Access/Generator.hpp index 67dbeba50..a5e9685e9 100644 --- a/include/OpenSpaceToolkit/Astrodynamics/Access/Generator.hpp +++ b/include/OpenSpaceToolkit/Astrodynamics/Access/Generator.hpp @@ -9,6 +9,7 @@ #include #include +#include #include #include @@ -34,6 +35,7 @@ using ostk::core::type::Real; using ostk::core::type::Shared; using ostk::mathematics::object::Interval; +using ostk::mathematics::object::Matrix3d; using ostk::physics::coordinate::Position; using ostk::physics::coordinate::spherical::AER; @@ -43,6 +45,7 @@ using ostk::physics::time::Duration; using ostk::physics::time::Instant; using ostk::physics::unit::Angle; using ostk::physics::unit::Length; +using EarthGravitationalModel = ostk::physics::environment::gravitational::Earth; using ostk::astrodynamics::Access; using ostk::astrodynamics::Trajectory; @@ -51,6 +54,117 @@ using ostk::astrodynamics::trajectory::State; #define DEFAULT_STEP Duration::Minutes(1.0) #define DEFAULT_TOLERANCE Duration::Microseconds(1.0) +/// @brief Represents the configuration for a ground target, including azimuth, elevation, and range intervals, as well +/// as position and LLA (Latitude, Longitude, Altitude). +/// +/// @details This class provides methods to retrieve the trajectory, position, LLA, and intervals for azimuth, +/// elevation, and range. It also includes a method to get the SEZ (South-East-Zenith) rotation matrix. +struct GroundTargetConfiguration +{ + /// @brief Constructor + /// @param anAzimuthInterval An azimuth interval [deg] + /// @param anElevationInterval An elevation interval [deg] + /// @param aRangeInterval A range interval [m] + /// @param aPosition A position + GroundTargetConfiguration( + const Interval& anAzimuthInterval, + const Interval& anElevationInterval, + const Interval& aRangeInterval, + const Position& aPosition + ); + + /// @brief Constructor + /// @param anAzimuthInterval An azimuth interval [deg] + /// @param anElevationInterval An elevation interval [deg] + /// @param aRangeInterval A range interval [m] + /// @param aLLA An LLA + GroundTargetConfiguration( + const Interval& anAzimuthInterval, + const Interval& anElevationInterval, + const Interval& aRangeInterval, + const LLA& aLLA + ); + + /// @brief Get the trajectory + /// + /// @code{.cpp} + /// GroundStationConfiguration groundStationConfiguration = { ... } ; + /// Trajectory trajectory = groundStationConfiguration.getTrajectory(); + /// @endcode + /// + /// @return The trajectory + Trajectory getTrajectory() const; + + /// @brief Get the position + /// + /// @code{.cpp} + /// GroundStationConfiguration groundStationConfiguration = { ... } ; + /// Position position = groundStationConfiguration.getPosition(); + /// @endcode + /// + /// @return The position + Position getPosition() const; + + /// @brief Get the latitude, longitude, and altitude (LLA) + /// + /// @code{.cpp} + /// GroundStationConfiguration groundStationConfiguration = { ... } ; + /// LLA lla = groundStationConfiguration.getLLA(); + /// @endcode + /// + /// @return The latitude, longitude, and altitude (LLA) + LLA getLLA() const; + + /// @brief Get the azimuth interval + /// + /// @code{.cpp} + /// GroundStationConfiguration groundStationConfiguration = { ... } ; + /// Interval groundStationConfiguration = generator.getAzimuthInterval(); + /// @endcode + /// + /// @return The azimuth interval + Interval getAzimuthInterval() const; + + /// @brief Get the elevation interval + /// + /// @code{.cpp} + /// GroundStationConfiguration groundStationConfiguration = { ... } ; + /// Interval groundStationConfiguration = generator.getElevationInterval(); + /// @endcode + /// + /// @return The elevation interval + Interval getElevationInterval() const; + + /// @brief Get the range interval + /// + /// @code{.cpp} + /// GroundStationConfiguration groundStationConfiguration = { ... } ; + /// Interval groundStationConfiguration = generator.getRangeInterval(); + /// @endcode + /// + /// @return The range interval + Interval getRangeInterval() const; + + /// @brief Get the rotation matrix (Matrix3d) from ECEF (Earth-Centered-Earth-Fixed) to SEZ (South-East-Zenith) + /// frame + /// + /// @code{.cpp} + /// GroundStationConfiguration groundStationConfiguration = { ... } ; + /// Matrix3d sezRotation = groundStationConfiguration.getR_SEZ_ECEF(); + /// @endcode + /// + /// @return The SEZ rotation matrix + Matrix3d getR_SEZ_ECEF() const; + + private: + Interval azimuthInterval_; + Interval elevationInterval_; + Interval rangeInterval_; + + Position position_; + LLA lla_; +}; + class Generator { public: @@ -89,6 +203,12 @@ class Generator const physics::time::Interval& anInterval, const Trajectory& aFromTrajectory, const Trajectory& aToTrajectory ) const; + Array> computeAccessesWithGroundTargets( + const physics::time::Interval& anInterval, + const Array& someGroundTargetConfigurations, + const Trajectory& aToTrajectory + ) const; + void setStep(const Duration& aStep); void setTolerance(const Duration& aTolerance); @@ -137,6 +257,13 @@ class Generator std::function accessFilter_; std::function stateFilter_; + Array generateAccessesFromIntervals( + const Array& someIntervals, + const physics::time::Interval& anInterval, + const Trajectory& aFromTrajectory, + const Trajectory& aToTrajectory + ) const; + static Access GenerateAccess( const physics::time::Interval& anAccessInterval, const physics::time::Interval& aGlobalInterval, diff --git a/src/OpenSpaceToolkit/Astrodynamics/Access/Generator.cpp b/src/OpenSpaceToolkit/Astrodynamics/Access/Generator.cpp index ac72c30b9..106472ff3 100644 --- a/src/OpenSpaceToolkit/Astrodynamics/Access/Generator.cpp +++ b/src/OpenSpaceToolkit/Astrodynamics/Access/Generator.cpp @@ -1,5 +1,6 @@ /// Apache License 2.0 +#include #include #include @@ -9,17 +10,23 @@ #include #include +#include #include using ostk::mathematics::geometry::d3::object::Point; using ostk::mathematics::geometry::d3::object::Segment; +using ostk::mathematics::object::MatrixXd; +using ostk::mathematics::object::Vector3d; +using ostk::mathematics::object::VectorXd; -using ostk::astrodynamics::solver::TemporalConditionSolver; using ostk::physics::coordinate::Frame; using ostk::physics::coordinate::spherical::LLA; using ostk::physics::environment::Object; using ostk::physics::environment::object::celestial::Earth; +using ostk::astrodynamics::RootSolver; +using ostk::astrodynamics::solver::TemporalConditionSolver; + namespace ostk { namespace astrodynamics @@ -27,6 +34,92 @@ namespace astrodynamics namespace access { +GroundTargetConfiguration::GroundTargetConfiguration( + const Interval& anAzimuthInterval, + const Interval& anElevationInterval, + const Interval& aRangeInterval, + const Position& aPosition +) + : azimuthInterval_(anAzimuthInterval), + elevationInterval_(anElevationInterval), + rangeInterval_(aRangeInterval), + position_(aPosition), + lla_(LLA::Cartesian( + aPosition.getCoordinates(), + EarthGravitationalModel::EGM2008.equatorialRadius_, + EarthGravitationalModel::EGM2008.flattening_ + )) +{ + if (aPosition.isDefined() && aPosition.accessFrame() != Frame::ITRF()) + { + throw ostk::core::error::RuntimeError("The position frame must be ITRF."); + } +} + +GroundTargetConfiguration::GroundTargetConfiguration( + const Interval& anAzimuthInterval, + const Interval& anElevationInterval, + const Interval& aRangeInterval, + const LLA& aLLA +) + : azimuthInterval_(anAzimuthInterval), + elevationInterval_(anElevationInterval), + rangeInterval_(aRangeInterval), + position_(Position::Meters( + aLLA.toCartesian( + EarthGravitationalModel::EGM2008.equatorialRadius_, EarthGravitationalModel::EGM2008.flattening_ + ), + Frame::ITRF() + )), + lla_(aLLA) +{ +} + +Trajectory GroundTargetConfiguration::getTrajectory() const +{ + return Trajectory::Position(position_); +} + +Position GroundTargetConfiguration::getPosition() const +{ + return position_; +} + +LLA GroundTargetConfiguration::getLLA() const +{ + return lla_; +} + +Interval GroundTargetConfiguration::getAzimuthInterval() const +{ + return azimuthInterval_; +} + +Interval GroundTargetConfiguration::getElevationInterval() const +{ + return elevationInterval_; +} + +Interval GroundTargetConfiguration::getRangeInterval() const +{ + return rangeInterval_; +} + +Matrix3d GroundTargetConfiguration::getR_SEZ_ECEF() const +{ + // TBM: Move this to OSTk physics as a utility function + const double sinLat = std::sin(lla_.getLatitude().inRadians()); + const double cosLat = std::cos(lla_.getLatitude().inRadians()); + const double sinLon = std::sin(lla_.getLongitude().inRadians()); + const double cosLon = std::cos(lla_.getLongitude().inRadians()); + + Matrix3d SEZRotation; + SEZRotation << sinLat * cosLon, sinLat * sinLon, -cosLat, -sinLon, cosLon, 0.0, cosLat * cosLon, cosLat * sinLon, + sinLat; + + return SEZRotation; +} + Generator::Generator(const Environment& anEnvironment, const Duration& aStep, const Duration& aTolerance) : environment_(anEnvironment), step_(aStep), @@ -152,23 +245,235 @@ Array Generator::computeAccesses( const Shared earthSPtr = this->environment_.accessCelestialObjectWithName("Earth"); - return accessIntervals - .map( - [&anInterval, &aFromTrajectory, &aToTrajectory, &earthSPtr, this]( - const physics::time::Interval& anAccessInterval - ) -> Access + return generateAccessesFromIntervals(accessIntervals, anInterval, aFromTrajectory, aToTrajectory); +} + +Array> Generator::computeAccessesWithGroundTargets( + const physics::time::Interval& anInterval, + const Array& someGroundTargetConfigurations, + const Trajectory& aToTrajectory +) const +{ + // create a stacked matrix of SEZ rotations for all ground targets + const Index targetCount = someGroundTargetConfigurations.getSize(); + Eigen::Matrix SEZRotations(3, 3 * targetCount); + + for (Index i = 0; i < targetCount; ++i) + { + SEZRotations.block(0, 3 * i, 3, 3) = someGroundTargetConfigurations[i].getR_SEZ_ECEF(); + } + + // create a stacked matrix of ITRF positions for all ground targets + MatrixXd fromPositionCoordinates_ITRF = MatrixXd::Zero(3, targetCount); + + for (Index i = 0; i < targetCount; ++i) + { + fromPositionCoordinates_ITRF.col(i) = someGroundTargetConfigurations[i].getPosition().getCoordinates(); + } + + // create a stacked matrix of azimuth, elevation, and range bounds for all ground targets + MatrixXd aerLowerBounds = MatrixXd::Zero(targetCount, 3); + MatrixXd aerUpperBounds = MatrixXd::Zero(targetCount, 3); + + for (Index i = 0; i < targetCount; ++i) + { + aerLowerBounds(i, 0) = + someGroundTargetConfigurations[i].getAzimuthInterval().accessLowerBound() * Real::Pi() / 180.0; + aerLowerBounds(i, 1) = + someGroundTargetConfigurations[i].getElevationInterval().accessLowerBound() * Real::Pi() / 180.0; + aerLowerBounds(i, 2) = + someGroundTargetConfigurations[i].getRangeInterval().accessLowerBound() * Real::Pi() / 180.0; + + aerUpperBounds(i, 0) = + someGroundTargetConfigurations[i].getAzimuthInterval().accessUpperBound() * Real::Pi() / 180.0; + aerUpperBounds(i, 1) = + someGroundTargetConfigurations[i].getElevationInterval().accessUpperBound() * Real::Pi() / 180.0; + aerUpperBounds(i, 2) = + someGroundTargetConfigurations[i].getRangeInterval().accessUpperBound() * Real::Pi() / 180.0; + } + + // initialize solver and condition + const RootSolver rootSolver = RootSolver(100, this->tolerance_.inSeconds()); + + const auto createAccessCondition = + [&fromPositionCoordinates_ITRF, &SEZRotations, &aToTrajectory, &aerLowerBounds, &aerUpperBounds]( + const Index& targetIdx + ) + { + const Vector3d fromPositionCoordinate_ITRF = fromPositionCoordinates_ITRF.col(targetIdx); + + const Matrix3d SEZRotation = SEZRotations.block(0, 3 * targetIdx, 3, 3); + + const double azimuthLowerBound = aerLowerBounds(targetIdx, 0); + const double azimuthUpperBound = aerUpperBounds(targetIdx, 0); + const double elevationLowerBound = aerLowerBounds(targetIdx, 1); + const double elevationUpperBound = aerUpperBounds(targetIdx, 1); + const double rangeLowerBound = aerLowerBounds(targetIdx, 2); + const double rangeUpperBound = aerUpperBounds(targetIdx, 2); + + return [fromPositionCoordinate_ITRF, + SEZRotation, + azimuthLowerBound, + azimuthUpperBound, + elevationLowerBound, + elevationUpperBound, + rangeLowerBound, + rangeUpperBound, + &aToTrajectory](const Instant& instant) -> bool + { + const auto toPositionCoordinates_ITRF = + aToTrajectory.getStateAt(instant).inFrame(Frame::ITRF()).getPosition().getCoordinates(); + + const VectorXd dx = toPositionCoordinates_ITRF - fromPositionCoordinate_ITRF; + + const MatrixXd dx_SEZ = SEZRotation * dx; + + const double range = dx_SEZ.norm(); + const double elevation_rad = std::asin(dx_SEZ(2) / range); + double azimuth_rad = std::atan2(dx_SEZ(1), dx_SEZ(0)); + + azimuth_rad = azimuth_rad < 0.0 ? azimuth_rad + 2.0 * M_PI : azimuth_rad; + + return azimuth_rad > azimuthLowerBound && azimuth_rad < azimuthUpperBound && + elevation_rad > elevationLowerBound && elevation_rad < elevationUpperBound && + range > rangeLowerBound && range < rangeUpperBound; + }; + }; + + // initialize containers + const Array instants = anInterval.generateGrid(this->step_); + + Array> conditionArray(targetCount, nullptr); + for (Index i = 0; i < targetCount; ++i) + { + conditionArray[i] = createAccessCondition(i); + } + + Array> accessIntervals = + Array>(targetCount, Array::Empty()); + Array previousAccessStates = Array(targetCount, false); + + for (Index index = 0; index < instants.getSize(); ++index) + { + const Instant& instant = instants[index]; + + // calculate target to satellite vector in ITRF + const Vector3d toPositionCoordinates_ITRF = + aToTrajectory.getStateAt(instant).inFrame(Frame::ITRF()).getPosition().getCoordinates(); + + const MatrixXd dx = + toPositionCoordinates_ITRF.replicate(1, fromPositionCoordinates_ITRF.cols()) - fromPositionCoordinates_ITRF; + + // calculate target to satellite vector in SEZ + MatrixXd dx_SEZ = MatrixXd::Zero(3, dx.cols()); + for (Index i = 0; i < (Index)dx.cols(); ++i) + { + dx_SEZ.col(i) = SEZRotations.block<3, 3>(0, 3 * i) * dx.col(i); + } + + // calculate azimuth, elevation, and range + const VectorXd range = dx_SEZ.colwise().norm(); + const VectorXd elevation_rad = (dx_SEZ.row(2).transpose().array() / range.array()).asin(); + VectorXd azimuth_rad = dx_SEZ.row(0).array().binaryExpr( + dx_SEZ.row(1).array(), + [](double x, double y) { - return Generator::GenerateAccess( - anAccessInterval, anInterval, aFromTrajectory, aToTrajectory, earthSPtr, this->tolerance_ - ); + return std::atan2(y, x); } - ) - .getWhere( - [this](const Access& anAccess) -> bool + ); + + azimuth_rad = azimuth_rad.unaryExpr( + [](double x) { - return this->accessFilter_ ? this->accessFilter_(anAccess) : true; + return x < 0.0 ? x + 2.0 * M_PI : x; } ); + + // check if satellite is in access + const auto inAccess = + (azimuth_rad.array() > aerLowerBounds.col(0).array() && + azimuth_rad.array() < aerUpperBounds.col(0).array() && + elevation_rad.array() > aerLowerBounds.col(1).array() && + elevation_rad.array() < aerUpperBounds.col(1).array() && range.array() > aerLowerBounds.col(2).array() && + range.array() < aerUpperBounds.col(2).array()) + .eval(); + + // First instant handling + if (index == 0) + { + for (Index i = 0; i < (Index)inAccess.size(); ++i) + { + previousAccessStates[i] = inAccess[i]; + if (inAccess[i]) + { + accessIntervals[i].add(physics::time::Interval::Closed(instant, instant)); + } + } + continue; + } + + const Instant& previousInstant = instants[index - 1]; + + // Check for state changes and find exact crossings + + for (Index i = 0; i < (Index)inAccess.size(); ++i) + { + const bool currentAccess = inAccess[i]; + + if (currentAccess != previousAccessStates[i]) + { + // Find exact crossing time + const auto condition = conditionArray[i]; + + const auto result = rootSolver.solve( + [&previousInstant, &condition](double aDurationInSeconds) -> double + { + return condition(previousInstant + Duration::Seconds(aDurationInSeconds)) ? +1.0 : -1.0; + }, + 0.0, + Duration::Between(previousInstant, instant).inSeconds() + ); + + const Instant crossingTime = previousInstant + Duration::Seconds(result.root); + + if (currentAccess) // Rising edge + { + accessIntervals[i].add(physics::time::Interval::Closed(crossingTime, instant)); + } + else // Falling edge + { + if (!accessIntervals[i].isEmpty()) + { + const physics::time::Interval& lastInterval = accessIntervals[i].accessLast(); + accessIntervals[i].accessLast() = + physics::time::Interval::Closed(lastInterval.getStart(), crossingTime); + } + } + } + else if (currentAccess) // Continuing access + { + if (!accessIntervals[i].isEmpty()) + { + const physics::time::Interval& lastInterval = accessIntervals[i].accessLast(); + accessIntervals[i].accessLast() = physics::time::Interval::Closed(lastInterval.getStart(), instant); + } + } + + previousAccessStates[i] = currentAccess; + } + } + + const Shared earthSPtr = this->environment_.accessCelestialObjectWithName("Earth"); + + Array> accesses = Array>(targetCount, Array::Empty()); + for (Index i = 0; i < accessIntervals.getSize(); ++i) + { + const Trajectory fromTrajectory = someGroundTargetConfigurations[i].getTrajectory(); + accesses[i] = + this->generateAccessesFromIntervals(accessIntervals[i], anInterval, fromTrajectory, aToTrajectory); + } + + return accesses; } void Generator::setStep(const Duration& aStep) @@ -224,16 +529,6 @@ Generator Generator::AerRanges( const Interval elevationRange_deg = anElevationRange; const Interval rangeRange_m = aRangeRange; - // const Interval azimuthRange_deg = anAzimuthRange; // anAzimuthRange.isDefined() ? - // Interval(anAzimuthRange.accessLowerBound().inDegrees(0.0, +360.0), - // anAzimuthRange.accessUpperBound().inDegrees(0.0, +360.0), anAzimuthRange.getType()) : Interval::Undefined() - //; const Interval elevationRange_deg = anElevationRange; // anElevationRange.isDefined() ? - // Interval(anElevationRange.accessLowerBound().inDegrees(-180.0, +180.0), - // anElevationRange.accessUpperBound().inDegrees(-180.0, +180.0), anElevationRange.getType()) : - // Interval::Undefined(); const Interval rangeRange_m = aRangeRange; // aRangeRange.isDefined() ? - // Interval(aRangeRange.accessLowerBound().inMeters(), aRangeRange.accessUpperBound().inMeters(), - // aRangeRange.getType()) : Interval::Undefined(); - const std::function aerFilter = [azimuthRange_deg, elevationRange_deg, rangeRange_m](const AER& anAER) -> bool { @@ -311,6 +606,34 @@ Generator Generator::AerMask( return {anEnvironment, aerFilter}; } +Array Generator::generateAccessesFromIntervals( + const Array& someIntervals, + const physics::time::Interval& anInterval, + const Trajectory& aFromTrajectory, + const Trajectory& aToTrajectory +) const +{ + const Shared earthSPtr = this->environment_.accessCelestialObjectWithName("Earth"); + + return someIntervals + .map( + [&anInterval, &aFromTrajectory, &aToTrajectory, &earthSPtr, this]( + const physics::time::Interval& anAccessInterval + ) -> Access + { + return Generator::GenerateAccess( + anAccessInterval, anInterval, aFromTrajectory, aToTrajectory, earthSPtr, this->tolerance_ + ); + } + ) + .getWhere( + [this](const Access& anAccess) -> bool + { + return this->accessFilter_ ? this->accessFilter_(anAccess) : true; + } + ); +} + Access Generator::GenerateAccess( const physics::time::Interval& anAccessInterval, const physics::time::Interval& aGlobalInterval, @@ -474,6 +797,7 @@ bool GeneratorContext::isAccessActive(const Instant& anInstant) const auto [fromPosition, toPosition] = GeneratorContext::GetPositionsFromStates(fromState, toState); // Line of sight + // TBI: Remove this check as it is redundant static const Shared commonFrameSPtr = Frame::GCRF(); @@ -542,12 +866,10 @@ AER GeneratorContext::CalculateAer( ) { // [TBM] This logic is Earth-specific - const Point referencePoint_ITRF = - Point::Vector(aFromPosition.inFrame(Frame::ITRF(), anInstant).accessCoordinates()); + const Vector3d referenceCoordinates_ITRF = aFromPosition.inFrame(Frame::ITRF(), anInstant).accessCoordinates(); - const LLA referencePoint_LLA = LLA::Cartesian( - referencePoint_ITRF.asVector(), anEarthSPtr->getEquatorialRadius(), anEarthSPtr->getFlattening() - ); + const LLA referencePoint_LLA = + LLA::Cartesian(referenceCoordinates_ITRF, anEarthSPtr->getEquatorialRadius(), anEarthSPtr->getFlattening()); const Shared nedFrameSPtr = anEarthSPtr->getFrameAt(referencePoint_LLA, Earth::FrameType::NED); diff --git a/test/OpenSpaceToolkit/Astrodynamics/Access/Generator.test.cpp b/test/OpenSpaceToolkit/Astrodynamics/Access/Generator.test.cpp index 134761821..06d8e4807 100644 --- a/test/OpenSpaceToolkit/Astrodynamics/Access/Generator.test.cpp +++ b/test/OpenSpaceToolkit/Astrodynamics/Access/Generator.test.cpp @@ -7,6 +7,8 @@ #include #include +#include + #include #include #include @@ -17,6 +19,8 @@ #include #include #include +#include +#include #include @@ -25,9 +29,13 @@ using ostk::core::container::Table; using ostk::core::container::Tuple; using ostk::core::filesystem::File; using ostk::core::filesystem::Path; +using ostk::core::type::Index; using ostk::core::type::Real; using ostk::core::type::String; +using ostk::mathematics::object::Matrix3d; +using ostk::mathematics::object::Vector3d; + using ostk::physics::coordinate::Frame; using ostk::physics::coordinate::Position; using ostk::physics::coordinate::spherical::AER; @@ -46,6 +54,7 @@ using ostk::physics::unit::Length; using ostk::astrodynamics::Access; using ostk::astrodynamics::access::Generator; +using ostk::astrodynamics::access::GroundTargetConfiguration; using ostk::astrodynamics::Trajectory; using ostk::astrodynamics::trajectory::Orbit; using ostk::astrodynamics::trajectory::orbit::model::Kepler; @@ -646,6 +655,226 @@ TEST(OpenSpaceToolkit_Astrodynamics_Access_Generator, ComputeAccesses_4) } } +TEST(OpenSpaceToolkit_Astrodynamics_Access_Generator, GroundTargetConfiguration) +{ + const ostk::mathematics::object::Interval azimuthInterval = + ostk::mathematics::object::Interval::Closed(0.0, 360.0); + const ostk::mathematics::object::Interval elevationInterval = + ostk::mathematics::object::Interval::Closed(0.0, 90.0); + const ostk::mathematics::object::Interval rangeInterval = + ostk::mathematics::object::Interval::Closed(0.0, 1.0e10); + const LLA lla = LLA::Vector({0.0, 0.0, 0.0}); + + // Constructor + { + { + EXPECT_NO_THROW(GroundTargetConfiguration(azimuthInterval, elevationInterval, rangeInterval, lla)); + } + + { + { + const Position position = Position::Meters({0.0, 0.0, 0.0}, Frame::ITRF()); + EXPECT_NO_THROW(GroundTargetConfiguration(azimuthInterval, elevationInterval, rangeInterval, position)); + } + + { + const Position position = Position::Meters({0.0, 0.0, 0.0}, Frame::GCRF()); + EXPECT_THROW( + GroundTargetConfiguration(azimuthInterval, elevationInterval, rangeInterval, position), + ostk::core::error::RuntimeError + ); + } + } + } + + // Getters + { + const GroundTargetConfiguration groundTargetConfiguration = + GroundTargetConfiguration(azimuthInterval, elevationInterval, rangeInterval, lla); + EXPECT_EQ(groundTargetConfiguration.getAzimuthInterval(), azimuthInterval); + EXPECT_EQ(groundTargetConfiguration.getElevationInterval(), elevationInterval); + EXPECT_EQ(groundTargetConfiguration.getRangeInterval(), rangeInterval); + EXPECT_EQ(groundTargetConfiguration.getLLA(), lla); + EXPECT_EQ( + groundTargetConfiguration.getPosition(), + Position::Meters( + lla.toCartesian(Earth::EGM2008.equatorialRadius_, Earth::EGM2008.flattening_), Frame::ITRF() + ) + ); + EXPECT_NO_THROW(groundTargetConfiguration.getTrajectory()); + + Matrix3d r_SEZ_ECEF; + r_SEZ_ECEF.row(0) = Vector3d {0.0, 0.0, 1.0}; + r_SEZ_ECEF.row(1) = Vector3d {0.0, 1.0, 0.0}; + r_SEZ_ECEF.row(2) = Vector3d {-1.0, 0.0, 0.0}; + + EXPECT_EQ(groundTargetConfiguration.getR_SEZ_ECEF(), r_SEZ_ECEF); + } +} + +TEST(OpenSpaceToolkit_Astrodynamics_Access_Generator, ComputeAccessesWithGroundTargets) +{ + { + const Environment environment = Environment::Default(); + + const TLE tle = { + "1 60504U 24149AN 24293.10070306 .00000000 00000-0 58313-3 0 08", + "2 60504 97.4383 7.6998 0003154 274.9510 182.9597 15.19652001 9607", + }; + const SGP4 sgp4 = SGP4(tle); + const Orbit aToTrajectory = Orbit(sgp4, environment.accessCelestialObjectWithName("Earth")); + + const Instant startInstant = Instant::Parse("2024-10-19 02:25:00.744.384", Scale::UTC); + const Instant endInstant = startInstant + Duration::Days(1.0); + const Interval interval = Interval::Closed(startInstant, endInstant); + + const Array LLAs = { + LLA(Angle::Degrees(53.406), Angle::Degrees(-6.225), Length::Meters(50.5)), + LLA(Angle::Degrees(13.51), Angle::Degrees(144.82), Length::Meters(46)), + LLA(Angle::Degrees(42.77), Angle::Degrees(141.62), Length::Meters(100)), + LLA(Angle::Degrees(47.2393), Angle::Degrees(-119.88515), Length::Meters(392.5)), + LLA(Angle::Degrees(78.22702), Angle::Degrees(15.38624), Length::Meters(493)), + LLA(Angle::Degrees(60.674), Angle::Degrees(17.142), Length::Meters(100)), + LLA(Angle::Degrees(37.7716), Angle::Degrees(-122.4135), Length::Meters(100)), + LLA(Angle::Degrees(69.09), Angle::Degrees(17.6986), Length::Meters(12)), + LLA(Angle::Degrees(78.23164), Angle::Degrees(15.37725), Length::Meters(483)), + LLA(Angle::Degrees(-72.0021), Angle::Degrees(2.5251), Length::Meters(1401)), + LLA(Angle::Degrees(-25.89), Angle::Degrees(27.71), Length::Meters(1562.66)), + LLA(Angle::Degrees(-46.53), Angle::Degrees(168.38), Length::Meters(13.65)), + LLA(Angle::Degrees(71.275), Angle::Degrees(-156.806), Length::Meters(24)), + LLA(Angle::Degrees(-52.9351), Angle::Degrees(-70.8713), Length::Meters(23)) + }; + + const ostk::mathematics::object::Interval azimuthInterval = + ostk::mathematics::object::Interval::Closed(0.0, 360.0); + const ostk::mathematics::object::Interval elevationInterval = + ostk::mathematics::object::Interval::Closed(0.0, 90.0); + const ostk::mathematics::object::Interval rangeInterval = + ostk::mathematics::object::Interval::Closed(0.0, 1.0e10); + + Array groundTargets = LLAs.map( + [&azimuthInterval, &elevationInterval, &rangeInterval](const LLA& lla) -> GroundTargetConfiguration + { + return GroundTargetConfiguration(azimuthInterval, elevationInterval, rangeInterval, lla); + } + ); + + const Generator generator = + Generator::AerRanges(azimuthInterval, elevationInterval, rangeInterval, environment); + + const Array> accessesPerTarget = + generator.computeAccessesWithGroundTargets(interval, groundTargets, aToTrajectory); + + ASSERT_EQ(accessesPerTarget.getSize(), groundTargets.getSize()); + + for (Index i = 0; i < accessesPerTarget.getSize(); ++i) + { + const Array accesses = accessesPerTarget.at(i); + const GroundTargetConfiguration groundTarget = groundTargets.at(i); + + const Array expectedAccesses = + generator.computeAccesses(interval, groundTarget.getTrajectory(), aToTrajectory); + + ASSERT_EQ(accesses.getSize(), expectedAccesses.getSize()); + + for (Index j = 0; j < accesses.getSize(); ++j) + { + const Access& access = accesses.at(j); + const Access& expectedAccess = expectedAccesses.at(j); + + EXPECT_TRUE(access.getAcquisitionOfSignal().isNear( + expectedAccess.getAcquisitionOfSignal(), Duration::Microseconds(1.0) + )) << access.getAcquisitionOfSignal().toString() + << " ~ " << expectedAccess.getAcquisitionOfSignal().toString(); + EXPECT_TRUE(access.getTimeOfClosestApproach().isNear( + expectedAccess.getTimeOfClosestApproach(), Duration::Microseconds(1.0) + )) << access.getTimeOfClosestApproach().toString() + << " ~ " << expectedAccess.getTimeOfClosestApproach().toString(); + EXPECT_TRUE( + access.getLossOfSignal().isNear(expectedAccess.getLossOfSignal(), Duration::Microseconds(1.0)) + ) << access.getLossOfSignal().toString() + << " ~ " << expectedAccess.getLossOfSignal().toString(); + EXPECT_TRUE(access.getDuration().isNear(expectedAccess.getDuration(), Duration::Microseconds(1.0))) + << access.getDuration().toString() << " ~ " << expectedAccess.getDuration().toString(); + } + } + } +} + +TEST(OpenSpaceToolkit_Astrodynamics_Access_Generator, SpeedTest) +{ + { + using ostk::astrodynamics::trajectory::orbit::model::Tabulated; + using ostk::mathematics::curvefitting::Interpolator; + const Environment environment = Environment::Default(); + + const TLE tle = { + "1 60504U 24149AN 24293.10070306 .00000000 00000-0 58313-3 0 08", + "2 60504 97.4383 7.6998 0003154 274.9510 182.9597 15.19652001 9607", + }; + const SGP4 sgp4 = SGP4(tle); + const Orbit aToTrajectory = Orbit(sgp4, environment.accessCelestialObjectWithName("Earth")); + + const Instant startInstant = Instant::Parse("2024-10-19 02:25:00.744.384", Scale::UTC); + const Instant endInstant = startInstant + Duration::Days(14.0); + const Interval interval = Interval::Closed(startInstant, endInstant); + + const Array instants = interval.generateGrid(Duration::Seconds(20.0)); + const Array states = aToTrajectory.getStatesAt(instants); + const Orbit tabulatedOrbit = Orbit( + Tabulated(states, 1, Interpolator::Type::CubicSpline), environment.accessCelestialObjectWithName("Earth") + ); + + const Array LLAs = { + LLA(Angle::Degrees(53.406), Angle::Degrees(-6.225), Length::Meters(50.5)), + LLA(Angle::Degrees(13.51), Angle::Degrees(144.82), Length::Meters(46)), + LLA(Angle::Degrees(42.77), Angle::Degrees(141.62), Length::Meters(100)), + LLA(Angle::Degrees(47.2393), Angle::Degrees(-119.88515), Length::Meters(392.5)), + LLA(Angle::Degrees(78.22702), Angle::Degrees(15.38624), Length::Meters(493)), + LLA(Angle::Degrees(60.674), Angle::Degrees(17.142), Length::Meters(100)), + LLA(Angle::Degrees(37.7716), Angle::Degrees(-122.4135), Length::Meters(100)), + LLA(Angle::Degrees(69.09), Angle::Degrees(17.6986), Length::Meters(12)), + LLA(Angle::Degrees(78.23164), Angle::Degrees(15.37725), Length::Meters(483)), + LLA(Angle::Degrees(-72.0021), Angle::Degrees(2.5251), Length::Meters(1401)), + LLA(Angle::Degrees(-25.89), Angle::Degrees(27.71), Length::Meters(1562.66)), + LLA(Angle::Degrees(-46.53), Angle::Degrees(168.38), Length::Meters(13.65)), + LLA(Angle::Degrees(71.275), Angle::Degrees(-156.806), Length::Meters(24)), + LLA(Angle::Degrees(-52.9351), Angle::Degrees(-70.8713), Length::Meters(23)), + }; + + const ostk::mathematics::object::Interval azimuthInterval = + ostk::mathematics::object::Interval::Closed(0.0, 360.0); + const ostk::mathematics::object::Interval elevationInterval = + ostk::mathematics::object::Interval::Closed(0.0, 90.0); + const ostk::mathematics::object::Interval rangeInterval = + ostk::mathematics::object::Interval::Closed(0.0, 1.0e10); + + Array groundTargets = LLAs.map( + [&azimuthInterval, &elevationInterval, &rangeInterval](const LLA& lla) -> GroundTargetConfiguration + { + return GroundTargetConfiguration(azimuthInterval, elevationInterval, rangeInterval, lla); + } + ); + + std::cout << "Speed test: ComputeAccessesWithGroundTargets" << std::endl; + // start timer + std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now(); + + const Generator generator = + Generator::AerRanges(azimuthInterval, elevationInterval, rangeInterval, environment); + + const Array> accessesPerTarget = + generator.computeAccessesWithGroundTargets(interval, groundTargets, tabulatedOrbit); + + // end timer + std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now(); + std::cout << "Time elapsed: " << std::chrono::duration_cast(end - start).count() + << " ms" << std::endl; + + ASSERT_EQ(accessesPerTarget.getSize(), groundTargets.getSize()); + } +} + TEST(OpenSpaceToolkit_Astrodynamics_Access_Generator, SetStep) { { @@ -835,7 +1064,8 @@ TEST(OpenSpaceToolkit_Astrodynamics_Access_Generator, AerRanges) << String::Format( "{} ~ {}", reference_acquisitionOfSignal.toString(), access.getAcquisitionOfSignal().toString() ); - // EXPECT_TRUE(access.getTimeOfClosestApproach().isNear(reference_timeOfClosestApproach, toleranceDuration)) + // EXPECT_TRUE(access.getTimeOfClosestApproach().isNear(reference_timeOfClosestApproach, + // toleranceDuration)) // << String::Format("{} ~ {}", reference_timeOfClosestApproach.toString(), // access.getTimeOfClosestApproach().toString()); EXPECT_TRUE(access.getLossOfSignal().isNear(reference_lossOfSignal, toleranceDuration)) @@ -940,7 +1170,8 @@ TEST(OpenSpaceToolkit_Astrodynamics_Access_Generator, AerMask) << String::Format( "{} ~ {}", reference_acquisitionOfSignal.toString(), access.getAcquisitionOfSignal().toString() ); - // EXPECT_TRUE(access.getTimeOfClosestApproach().isNear(reference_timeOfClosestApproach, toleranceDuration)) + // EXPECT_TRUE(access.getTimeOfClosestApproach().isNear(reference_timeOfClosestApproach, + // toleranceDuration)) // << String::Format("{} ~ {}", reference_timeOfClosestApproach.toString(), // access.getTimeOfClosestApproach().toString()); EXPECT_TRUE(access.getLossOfSignal().isNear(reference_lossOfSignal, toleranceDuration))