feat: exponential drag acceleration (#99)

* feat: add constant drag accel just to test

* feat: add drag-defined flag

* test: trying to modify current tests to ignore atmospheric model

* test: get the tests to pass by setting the drag coefficient to 0 for existing tests

* feat: add RungeKutta4 solver from boost so we can compare to OREkit

* test: add validation tests against OREKit for exponential atmospheric drag

* test: add perf tests for new RK4 integrator wrapper

* test: add high-drag exponential drag test comparison with OREKit

* fix: remove accidental git rebase markers

* fix: set drag coefficient for non-drag performance tests to 0

* test: explicitly set Exponential atmosphere instead of relying on it as the default in tests

* chore: remove some debug stuff

* fix: forgot an import

* feat: python bindings for rk4 numerical solver

* Apply suggestions from code review

style: use full atmosphericDensity naming instead  of rho

Co-authored-by: vishwa2710 <vishwa2710@gmail.com>

* style: formatting, weird header changes. May have to revert

* test: re-add test that was removed during rebase

* chore: bump to OSTk Physics requirement to 0.8.3

* style: re-order numerical solvers to put RK4 first

* style: formatting

* test: set drag coefficient to 0 in existing python bindings tests

* refactor: modify atmosphere query in drag calculation to use isDefined() function

* refactor: suggestions from PR

* style: ordering for rk4

* fix: Celestial models cannot be changed after instantiation apparently

* fix: get back to a working state

* fix: add atmosphere model explicitly to other drag tests

* style: formatting

* fix: had the frame rotation backwards

* style: formatting/comments

* fix: bindings typo

* style: comment clarity

* refactor: let there be const

* build: use physics 1.0.0 version

* build: temporarily adapt version pattern for find package in cmake

* chore: propagate changes from new physics version in astro code and tests

---------

Co-authored-by: vishwa2710 <vishwa2710@gmail.com>
Co-authored-by: Remy Derollez <remy@loftorbital.com>

CMakeLists.txt

@@ -277,7 +277,7 @@ ENDIF ()
 
 ### Open Space Toolkit ▸ Physics
 
-FIND_PACKAGE ("OpenSpaceToolkitPhysics" "0.8.2" REQUIRED)
+FIND_PACKAGE ("OpenSpaceToolkitPhysics" "1.0" REQUIRED)
 
 IF (NOT OpenSpaceToolkitPhysics_FOUND)
 

bindings/python/requirements.txt

@@ -1,3 +1,3 @@
 # Apache License 2.0
 
-open-space-toolkit-physics~=0.8.1
+open-space-toolkit-physics~=1.0.0

bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/NumericalSolver.cpp

@@ -106,6 +106,7 @@ inline void OpenSpaceToolkitAstrodynamicsPy_NumericalSolver(pybind11::module& aM
 
         enum_<NumericalSolver::StepperType>(numericalSolver, "StepperType")
 
+            .value("RungeKutta4", NumericalSolver::StepperType::RungeKutta4)
             .value("RungeKuttaCashKarp54", NumericalSolver::StepperType::RungeKuttaCashKarp54)
             .value("RungeKuttaFehlberg78", NumericalSolver::StepperType::RungeKuttaFehlberg78)
 

bindings/python/test/test_numerical_solver.py

@@ -80,6 +80,12 @@ def test_get_types(
         assert numerical_solver.get_absolute_tolerance() == absolute_tolerance
 
     def test_get_string_from_types(self):
+        assert (
+            NumericalSolver.string_from_stepper_type(
+                NumericalSolver.StepperType.RungeKutta4
+            )
+            == "RungeKutta4"
+        )
         assert (
             NumericalSolver.string_from_stepper_type(
                 NumericalSolver.StepperType.RungeKuttaCashKarp54

bindings/python/test/trajectory/orbit/models/test_propagated.py

@@ -52,7 +52,7 @@ def propagated_default_inputs():
     )
     inertia_tensor = np.ndarray(shape=(3, 3))
     surface_area = 0.8
-    drag_coefficient = 2.2
+    drag_coefficient = 0.0
 
     satellitesystem = SatelliteSystem(
         mass, satellite_geometry, inertia_tensor, surface_area, drag_coefficient

bindings/python/test/trajectory/test_propagator.py

@@ -52,7 +52,7 @@ def propagator_default_inputs():
     )
     inertia_tensor = np.identity(3)
     surface_area = 0.8
-    drag_coefficient = 2.2
+    drag_coefficient = 0.0
 
     satellitesystem = SatelliteSystem(
         mass, satellite_geometry, inertia_tensor, surface_area, drag_coefficient
@@ -77,7 +77,7 @@ def propagator(propagator_default_inputs) -> Propagator:
     return Propagator(*propagator_default_inputs[:2])
 
 
-class TestPropagated:
+class TestPropagator:
     def test_constructors(self, propagator: Propagator):
         assert propagator is not None
         assert isinstance(propagator, Propagator)

docker/development/Dockerfile

@@ -127,8 +127,8 @@ RUN mkdir -p /tmp/open-space-toolkit-mathematics \
 
 ## Open Space Toolkit ▸ Physics
 
-ARG OSTK_PHYSICS_MAJOR="0"
-ARG OSTK_PHYSICS_MINOR="8"
+ARG OSTK_PHYSICS_MAJOR="1"
+ARG OSTK_PHYSICS_MINOR="0"
 
 RUN mkdir -p /tmp/open-space-toolkit-physics \
  && cd /tmp/open-space-toolkit-physics \

include/OpenSpaceToolkit/Astrodynamics/Flight/System/Dynamics/SatelliteDynamics.hpp

@@ -11,6 +11,7 @@
 
 #include <OpenSpaceToolkit/Mathematics/Objects/Vector.hpp>
 
+#include <OpenSpaceToolkit/Physics/Coordinate/Frame.hpp>
 #include <OpenSpaceToolkit/Physics/Data/Vector.hpp>
 #include <OpenSpaceToolkit/Physics/Environment.hpp>
 #include <OpenSpaceToolkit/Physics/Environment/Objects/CelestialBodies/Earth.hpp>
@@ -163,15 +164,7 @@ class SatelliteDynamics : public Dynamics
     SatelliteSystem satelliteSystem_;
     Instant instant_;
 
-    // Only force model currently used that incorporates solely Earth's gravity
     void DynamicalEquations(const Dynamics::StateVector& x, Dynamics::StateVector& dxdt, const double t);
-
-    // // Atmospheric perturbations only
-    // void                    Exponential_Dynamics                        (   const   SatelliteDynamics::StateVector&
-    // x,
-    //                                                                                 SatelliteDynamics::StateVector&
-    //                                                                                 dxdt,
-    //                                                                         const   double ) const ;
 };
 
 }  // namespace dynamics

include/OpenSpaceToolkit/Astrodynamics/NumericalSolver.hpp

@@ -33,6 +33,7 @@ class NumericalSolver
    public:
     enum class StepperType
     {
+        RungeKutta4,
         RungeKuttaCashKarp54,
         RungeKuttaFehlberg78
     };

src/OpenSpaceToolkit/Astrodynamics/Flight/System/Dynamics/SatelliteDynamics.cpp

@@ -3,6 +3,8 @@
 #include <OpenSpaceToolkit/Core/Error.hpp>
 #include <OpenSpaceToolkit/Core/Utilities.hpp>
 
+#include <OpenSpaceToolkit/Physics/Environment/Objects/Celestial.hpp>
+
 #include <OpenSpaceToolkit/Astrodynamics/Flight/System/Dynamics/SatelliteDynamics.hpp>
 
 namespace ostk
@@ -19,6 +21,8 @@ namespace dynamics
 using ostk::physics::units::Derived;
 using ostk::physics::units::Length;
 using ostk::physics::units::Time;
+using ostk::physics::coord::Frame;
+using ostk::physics::env::obj::Celestial;
 
 static const Derived::Unit GravitationalParameterSIUnit =
     Derived::Unit::GravitationalParameter(Length::Unit::Meter, Time::Unit::Second);
@@ -134,88 +138,73 @@ void SatelliteDynamics::DynamicalEquations(const Dynamics::StateVector& x, Dynam
     const Instant currentInstant = instant_ + Duration::Seconds(t);
     environment_.setInstant(currentInstant);
 
+    // GRAVITY
     // Initialize gravitational acceleration vector
     Vector3d totalGravitationalAcceleration_SI = {0.0, 0.0, 0.0};
 
     // Access all objects in the environment and loop through them
     for (const auto& objectName : environment_.getObjectNames())
     {
+        const Shared<const Celestial> object = environment_.accessCelestialObjectWithName(objectName);
+
         if (objectName != "Earth")
         {
             // Obtain 3rd body effect on center of Earth (origin in GCRF) aka 3rd body correction
             const Vector gravitationalAcceleration3rdBodyCorrection =
-                environment_.accessCelestialObjectWithName(objectName)
-                    ->getGravitationalFieldAt(Position::Meters({0.0, 0.0, 0.0}, gcrfSPtr_));
+                object->getGravitationalFieldAt(Position::Meters({0.0, 0.0, 0.0}, gcrfSPtr_), currentInstant);
 
             // Subtract 3rd body correct from total gravitational acceleration
             totalGravitationalAcceleration_SI -=
                 gravitationalAcceleration3rdBodyCorrection.inFrame(gcrfSPtr_, currentInstant).getValue();
         }
 
         // Obtain gravitational acceleration from current object
-        const Vector gravitationalAcceleration =
-            environment_.accessCelestialObjectWithName(objectName)->getGravitationalFieldAt(currentPosition);
+        const Vector gravitationalAcceleration = object->getGravitationalFieldAt(currentPosition, currentInstant);
 
         // Add object's gravity to total gravitational acceleration
         totalGravitationalAcceleration_SI += gravitationalAcceleration.inFrame(gcrfSPtr_, currentInstant).getValue();
     }
 
-    // Integrate position and velocity states
+    // DRAG
+    Vector3d totalDragAcceleration = {0.0, 0.0, 0.0};
+
+    const Real mass = satelliteSystem_.getMass().inKilograms();
+    const Real dragCoefficient = satelliteSystem_.getDragCoefficient();
+    const Real surfaceArea = satelliteSystem_.getCrossSectionalSurfaceArea();
+
+    for (const auto& objectName : environment_.getObjectNames())
+    {
+        const Shared<const Celestial> object = environment_.accessCelestialObjectWithName(objectName);
+
+        // TBI: currently only defined for Earth
+        if (object->atmosphericModelIsDefined())
+        {
+            const Real atmosphericDensity = object->getAtmosphericDensityAt(currentPosition, currentInstant).getValue();
+
+            const Vector3d earthAngularVelocity =
+                Frame::GCRF()->getTransformTo(Frame::ITRF(), currentInstant).getAngularVelocity();  // rad/s
+
+            const Vector3d relativeVelocity =
+                Vector3d(x[3], x[4], x[5]) - earthAngularVelocity.cross(Vector3d(x[0], x[1], x[2]));
+
+            // Calculate drag acceleration
+            const Vector3d dragAcceleration = -(0.5 / mass) * dragCoefficient * surfaceArea * atmosphericDensity *
+                                              relativeVelocity.norm() * relativeVelocity;
+
+            totalDragAcceleration += dragAcceleration;
+        }
+    }
+
+    // Propagate velocity
     dxdt[0] = x[3];
     dxdt[1] = x[4];
     dxdt[2] = x[5];
-    dxdt[3] = totalGravitationalAcceleration_SI[0];
-    dxdt[4] = totalGravitationalAcceleration_SI[1];
-    dxdt[5] = totalGravitationalAcceleration_SI[2];
-}
 
-// void                            SatelliteDynamics::Exponential_Dynamics
-//                                                                             (   const SatelliteDynamics::StateVector&
-//                                                                             x,
-//                                                                                         SatelliteDynamics::StateVector&
-//                                                                                         dxdt,
-//                                                                                 const   double ) const
-// {
-//     // Find position magnitude
-//     const double positionMagnitude = std::sqrt(std::pow(x[0],2) + std::pow(x[1],2) + std::pow(x[2],2));
-
-//     // Integrate velocity states
-//     const double positionMagnitudeCubed = std::pow(positionMagnitude,3) ;
-
-//     // Access constants
-//     const double mu_SI =
-//     ((environment_.accessCelestialObjectWithName("Earth"))->getGravitationalParameter()).in(GravitationalParameterSIUnit)
-//     ; const double Re = ((environment_.accessCelestialObjectWithName("Earth"))->getEquatorialRadius()).inMeters() ;
-//     const double dragCoefficient = satelliteSystem_.getDragCoefficient() ;
-//     const double surfaceArea = satelliteSystem_.getCrossSectionalSurfaceArea() ;
-//     const double mass = satelliteSystem_.getMass().inKilograms() ;
-//     const Vector3d earthAngularVelocity = { 0, 0, 7.2921159e-5 } ; // rad/s
-//     const double rho0 = 6.967e-13 ; // kg/m^3
-//     const double h0 = 500000 ;
-//     const double H = 63822 ;
-
-//     // Calculate relative velocity to wind
-//     const Vector3d relativeVelocity = Vector3d({ x[3], x[4], x[5] }) -
-//     Vector3d({earthAngularVelocity.cross(Vector3d({ x[0], x[1], x[2] }))}) ; const double relativeVelocityMagnitude =
-//     std::sqrt(std::pow(relativeVelocity[0],2) + std::pow(relativeVelocity[1],2) + std::pow(relativeVelocity[2],2)) ;
-
-//     // Calculate density parameters
-//     const double hEllipse = positionMagnitude - Re ;
-//     const double rho = rho0 * std::exp( - (hEllipse - h0) / H ) ;
-
-//     // Calculate drag acceleration
-//     const Vector3d dragAcceleration = -0.5 * dragCoefficient * (surfaceArea / mass) * rho * relativeVelocityMagnitude
-//     * relativeVelocity ;
-
-//     // Integrate position and velocity states
-//     dxdt[0] = x[3] ;
-//     dxdt[1] = x[4] ;
-//     dxdt[2] = x[5] ;
-//     dxdt[3] = -(mu_SI / positionMagnitudeCubed) * x[0] + dragAcceleration[0] ;
-//     dxdt[4] = -(mu_SI / positionMagnitudeCubed) * x[1] + dragAcceleration[1] ;
-//     dxdt[5] = -(mu_SI / positionMagnitudeCubed) * x[2] + dragAcceleration[2] ;
-
-// }
+    // Propagate acceleration
+    dxdt[3] = totalGravitationalAcceleration_SI[0] + totalDragAcceleration[0];
+    dxdt[4] = totalGravitationalAcceleration_SI[1] + totalDragAcceleration[1];
+    dxdt[5] = totalGravitationalAcceleration_SI[2] + totalDragAcceleration[2];
+}
 
 }  // namespace dynamics
 }  // namespace system

src/OpenSpaceToolkit/Astrodynamics/NumericalSolver.cpp

@@ -14,6 +14,7 @@ namespace astro
 
 using namespace boost::numeric::odeint;
 
+typedef runge_kutta4<NumericalSolver::StateVector> stepper_type_4;
 typedef runge_kutta_cash_karp54<NumericalSolver::StateVector> error_stepper_type_54;
 typedef runge_kutta_fehlberg78<NumericalSolver::StateVector> error_stepper_type_78;
 
@@ -195,6 +196,23 @@ Array<NumericalSolver::StateVector> NumericalSolver::integrateStatesAtSortedInst
 
     switch (stepperType_)
     {
+        case NumericalSolver::StepperType::RungeKutta4:
+        {
+            integrate_times(
+                stepper_type_4(),
+                aSystemOfEquations,
+                aStateVector,
+                anIntegrationDurationInSecsArray.begin(),
+                anIntegrationDurationInSecsArray.end(),
+                adjustedTimeStep,
+                [&](const NumericalSolver::StateVector& x, double t) -> void
+                {
+                    this->observeNumericalIntegration(x, t);
+                }
+            );
+            break;
+        }
+
         case NumericalSolver::StepperType::RungeKuttaCashKarp54:
         {
             integrate_times(
@@ -261,6 +279,36 @@ NumericalSolver::StateVector NumericalSolver::integrateStateForDuration(
 
     switch (stepperType_)
     {
+        case NumericalSolver::StepperType::RungeKutta4:
+        {
+            // Integrate_adaptive uses constant step size under the hood
+            // for a stepper without error control like RK4.
+            // Therefore, just use integrate_const for simplicity.
+            switch (logType_)
+            {
+                case NumericalSolver::LogType::NoLog:
+                case NumericalSolver::LogType::LogAdaptive:
+                case NumericalSolver::LogType::LogConstant:
+                {
+                    integrate_const(
+                        stepper_type_4(),
+                        aSystemOfEquations,
+                        aStateVector,
+                        (0.0),
+                        integrationDurationInSecs,
+                        adjustedTimeStep,
+                        [&](const NumericalSolver::StateVector& x, double t) -> void
+                        {
+                            this->observeNumericalIntegration(x, t);
+                        }
+                    );
+                    return aStateVector;
+                }
+                default:
+                    throw ostk::core::error::runtime::Wrong("Log type");
+            }
+        }
+
         case NumericalSolver::StepperType::RungeKuttaCashKarp54:
         {
             switch (logType_)
@@ -388,6 +436,9 @@ String NumericalSolver::StringFromStepperType(const NumericalSolver::StepperType
 {
     switch (aStepperType)
     {
+        case NumericalSolver::StepperType::RungeKutta4:
+            return "RungeKutta4";
+
         case NumericalSolver::StepperType::RungeKuttaCashKarp54:
             return "RungeKuttaCashKarp54";
 

src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit.cpp

@@ -874,7 +874,7 @@ Orbit Orbit::SunSynchronous(
         // Sun direction in GCRF
 
         const Vector3d sunDirection_GCRF = environment.accessCelestialObjectWithName("Sun")
-                                               ->getPositionIn(Frame::GCRF())
+                                               ->getPositionIn(Frame::GCRF(), anEpoch)
                                                .getCoordinates()
                                                .normalized();