Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <array>
       7             : #include <boost/math/quaternion.hpp>
       8             : #include <cmath>
       9             : #include <limits>
      10             : #include <pup.h>
      11             : #include <string>
      12             : #include <vector>
      13             : 
      14             : #include "DataStructures/DataVector.hpp"
      15             : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
      16             : #include "Domain/FunctionsOfTime/FunctionOfTimeHelpers.hpp"
      17             : #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp"
      18             : #include "Utilities/Gsl.hpp"
      19             : #include "Utilities/Serialization/CharmPupable.hpp"
      20             : 
      21             : namespace domain::FunctionsOfTime {
      22             : /// \ingroup ComputationalDomainGroup
      23             : /// \brief A FunctionOfTime that stores quaternions for the rotation map
      24             : ///
      25             : /// \details This FunctionOfTime stores quaternions that will be used in the
      26             : /// time-dependent rotation map as well as the orbital angular velocity that
      27             : /// will be controlled by the rotation control sytem. To get the quaternion, an
      28             : /// ODE is solved of the form \f$ \dot{q} = \frac{1}{2} q \times \omega \f$
      29             : /// where \f$ \omega \f$ is the orbital angular velocity which is stored
      30             : /// internally as the derivative of an angle `PiecewisePolynomial`, and
      31             : /// \f$ \times \f$ here is quaternion multiplication.
      32             : ///
      33             : /// Different from a `PiecewisePolynomial`, only the quaternion
      34             : /// itself is stored, not any of the derivatives because the derivatives must be
      35             : /// calculated from the solved ODE at every function call. Because
      36             : /// derivatives of the quaternion are not stored, the template parameter
      37             : /// `MaxDeriv` refers to both the max derivative of the stored angle
      38             : /// PiecewisePolynomial and the max derivative returned by the
      39             : /// QuaternionFunctionOfTime. The `update` function is then just a wrapper
      40             : /// around the internal `PiecewisePolynomial::update` function with the addition
      41             : /// that it then updates the stored quaternions as well.
      42             : ///
      43             : /// The angle PiecewisePolynomial is accessible through the `angle_func`,
      44             : /// `angle_func_and_deriv`, and `angle_func_and_2_derivs` functions which
      45             : /// correspond to the function calls of a normal PiecewisePolynomial except
      46             : /// without the `angle_` prefix.
      47             : ///
      48             : /// It is encouraged to use `quat_func` and `angle_func` when you want the
      49             : /// specific values of the functions to avoid ambiguity in what you are
      50             : /// calling. However, the original three `func` functions inherited from the
      51             : /// FunctionOfTime base class are necessary because the maps use the generic
      52             : /// `func` functions, thus they return the quaternion and its derivatives (which
      53             : /// are needed for the map). This is all to keep the symmetry of naming
      54             : /// `angle_func` and `quat_func` so that function calls won't be ambiguous.
      55             : template <size_t MaxDeriv>
      56           1 : class QuaternionFunctionOfTime : public FunctionOfTime {
      57             :  public:
      58           0 :   QuaternionFunctionOfTime() = default;
      59           0 :   QuaternionFunctionOfTime(
      60             :       double t, std::array<DataVector, 1> initial_quat_func,
      61             :       std::array<DataVector, MaxDeriv + 1> initial_angle_func,
      62             :       double expiration_time);
      63             : 
      64           0 :   ~QuaternionFunctionOfTime() override = default;
      65           0 :   QuaternionFunctionOfTime(QuaternionFunctionOfTime&&) = default;
      66           0 :   QuaternionFunctionOfTime& operator=(QuaternionFunctionOfTime&&) = default;
      67           0 :   QuaternionFunctionOfTime(const QuaternionFunctionOfTime&) = default;
      68           0 :   QuaternionFunctionOfTime& operator=(const QuaternionFunctionOfTime&) =
      69             :       default;
      70             : 
      71             :   // LCOV_EXCL_START
      72           0 :   explicit QuaternionFunctionOfTime(CkMigrateMessage* /*unused*/) {}
      73             :   // LCOV_EXCL_STOP
      74             : 
      75           0 :   auto get_clone() const -> std::unique_ptr<FunctionOfTime> override;
      76             : 
      77             :   // clang-tidy: google-runtime-references
      78             :   // clang-tidy: cppcoreguidelines-owning-memory,-warnings-as-errors
      79           0 :   WRAPPED_PUPable_decl_template(QuaternionFunctionOfTime<MaxDeriv>);  // NOLINT
      80             : 
      81           1 :   void reset_expiration_time(const double next_expiration_time) override {
      82             :     angle_f_of_t_.reset_expiration_time(next_expiration_time);
      83             :   }
      84             : 
      85             :   /// Returns domain of validity for the function of time
      86           1 :   std::array<double, 2> time_bounds() const override {
      87             :     return angle_f_of_t_.time_bounds();
      88             :   }
      89             : 
      90             :   /// Updates the `MaxDeriv`th derivative of the angle piecewisepolynomial at
      91             :   /// the given time, then updates the stored quaternions.
      92             :   ///
      93             :   /// `updated_max_deriv` is a datavector of the `MaxDeriv`s for each component.
      94             :   /// `next_expiration_time` is the next expiration time.
      95           1 :   void update(double time_of_update, DataVector updated_max_deriv,
      96             :               double next_expiration_time) override;
      97             : 
      98             :   // NOLINTNEXTLINE(google-runtime-references)
      99           0 :   void pup(PUP::er& p) override;
     100             : 
     101             :   /// Returns the quaternion at an arbitrary time `t`.
     102           1 :   std::array<DataVector, 1> func(const double t) const override {
     103             :     return quat_func(t);
     104             :   }
     105             : 
     106             :   /// Returns the quaternion and its first derivative at an arbitrary time `t`.
     107           1 :   std::array<DataVector, 2> func_and_deriv(const double t) const override {
     108             :     return quat_func_and_deriv(t);
     109             :   }
     110             : 
     111             :   /// Returns the quaternion and the first two derivatives at an arbitrary
     112             :   /// time `t`.
     113           1 :   std::array<DataVector, 3> func_and_2_derivs(const double t) const override {
     114             :     return quat_func_and_2_derivs(t);
     115             :   }
     116             : 
     117             :   /// Returns the quaternion at an arbitrary time `t`.
     118           1 :   std::array<DataVector, 1> quat_func(double t) const;
     119             : 
     120             :   /// Returns the quaternion and its first derivative at an arbitrary time `t`.
     121           1 :   std::array<DataVector, 2> quat_func_and_deriv(double t) const;
     122             : 
     123             :   /// Returns the quaternion and the first two derivatives at an arbitrary
     124             :   /// time `t`.
     125           1 :   std::array<DataVector, 3> quat_func_and_2_derivs(double t) const;
     126             : 
     127             :   /// Returns stored angle at an arbitrary time `t`.
     128           1 :   std::array<DataVector, 1> angle_func(const double t) const {
     129             :     return angle_f_of_t_.func(t);
     130             :   }
     131             : 
     132             :   /// Returns stored angle and its first derivative (omega) at an arbitrary time
     133             :   /// `t`.
     134           1 :   std::array<DataVector, 2> angle_func_and_deriv(const double t) const {
     135             :     return angle_f_of_t_.func_and_deriv(t);
     136             :   }
     137             : 
     138             :   /// Returns stored angle and the first two derivatives at an arbitrary
     139             :   /// time `t`.
     140           1 :   std::array<DataVector, 3> angle_func_and_2_derivs(const double t) const {
     141             :     return angle_f_of_t_.func_and_2_derivs(t);
     142             :   }
     143             : 
     144             :  private:
     145             :   template <size_t LocalMaxDeriv>
     146           0 :   friend bool operator==(  // NOLINT(readability-redundant-declaration)
     147             :       const QuaternionFunctionOfTime<LocalMaxDeriv>& lhs,
     148             :       const QuaternionFunctionOfTime<LocalMaxDeriv>& rhs);
     149             : 
     150             :   template <size_t LocalMaxDeriv>
     151           0 :   friend std::ostream& operator<<(  // NOLINT(readability-redundant-declaration)
     152             :       std::ostream& os,
     153             :       const QuaternionFunctionOfTime<LocalMaxDeriv>& quaternion_f_of_t);
     154             : 
     155             :   std::vector<FunctionOfTimeHelpers::StoredInfo<1, false>>
     156           0 :       stored_quaternions_and_times_;
     157             : 
     158           0 :   domain::FunctionsOfTime::PiecewisePolynomial<MaxDeriv> angle_f_of_t_;
     159             : 
     160             :   /// Integrates the ODE \f$ \dot{q} = \frac{1}{2} q \times \omega \f$ from time
     161             :   /// `t0` to time `t`. On input, `quaternion_to_integrate` is the initial
     162             :   /// quaternion at time `t0` and on output, it stores the result at time `t`
     163           1 :   void solve_quaternion_ode(
     164             :       gsl::not_null<boost::math::quaternion<double>*> quaternion_to_integrate,
     165             :       double t0, double t) const;
     166             : 
     167             :   /// Updates the `std::vector<StoredInfo>` to have the same number of stored
     168             :   /// quaternions as the `angle_f_of_t_ptr` has stored angles. This is necessary
     169             :   /// to ensure we can solve the ODE at any time `t`
     170           1 :   void update_stored_info();
     171             : 
     172             :   /// Does common operations to all the `func` functions such as updating stored
     173             :   /// info, solving the ODE, and returning the normalized quaternion as a boost
     174             :   /// quaternion for easy calculations
     175           1 :   boost::math::quaternion<double> setup_func(double t) const;
     176             : };
     177             : 
     178             : template <size_t MaxDeriv>
     179           0 : bool operator!=(const QuaternionFunctionOfTime<MaxDeriv>& lhs,
     180             :                 const QuaternionFunctionOfTime<MaxDeriv>& rhs);
     181             : 
     182             : /// \cond
     183             : template <size_t MaxDeriv>
     184             : PUP::able::PUP_ID QuaternionFunctionOfTime<MaxDeriv>::my_PUP_ID = 0;  // NOLINT
     185             : /// \endcond
     186             : }  // namespace domain::FunctionsOfTime