Line data    Source code

       1             : // SPDX-FileCopyrightText: 2024 PairInteraction Developers
       2             : // SPDX-License-Identifier: LGPL-3.0-or-later
       3             : 
       4             : #pragma once
       5             : 
       6             : #include "pairinteraction/utils/eigen_assertion.hpp"
       7             : #include "pairinteraction/utils/eigen_compat.hpp"
       8             : 
       9             : #include <Eigen/Dense>
      10             : #include <array>
      11             : #include <type_traits>
      12             : 
      13             : namespace pairinteraction::euler {
      14             : /**
      15             :  * @function get_rotation_matrix
      16             :  *
      17             :  * @brief Build a matrix that rotates the coordinate system to the new z-axis and y-axis.
      18             :  *
      19             :  * @tparam Real The type of a real number.
      20             :  *
      21             :  * @param to_z_axis The new z-axis.
      22             :  *
      23             :  * @param to_y_axis The new y-axis.
      24             :  *
      25             :  * @return The passive rotation matrix.
      26             :  */
      27             : 
      28             : template <typename Real>
      29           5 : inline Eigen::Matrix3<Real> get_rotation_matrix(std::array<Real, 3> to_z_axis,
      30             :                                                 std::array<Real, 3> to_y_axis) {
      31             :     static_assert(std::is_floating_point_v<Real>);
      32             : 
      33           5 :     auto to_z_axis_mapped = Eigen::Map<Eigen::Vector3<Real>>(to_z_axis.data()).normalized();
      34           5 :     auto to_y_axis_mapped = Eigen::Map<Eigen::Vector3<Real>>(to_y_axis.data()).normalized();
      35             : 
      36           5 :     auto scale = to_z_axis_mapped.norm() * to_y_axis_mapped.norm();
      37           5 :     auto numerical_precision = 100 * scale * std::numeric_limits<Real>::epsilon();
      38             : 
      39           5 :     if (std::abs(to_z_axis_mapped.dot(to_y_axis_mapped)) > numerical_precision) {
      40           1 :         throw std::runtime_error("The z-axis and the y-axis are not orhogonal.");
      41             :     }
      42             : 
      43           4 :     Eigen::Matrix3<Real> rotator;
      44           4 :     rotator << to_y_axis_mapped.cross(to_z_axis_mapped), to_y_axis_mapped, to_z_axis_mapped;
      45             : 
      46           8 :     return rotator;
      47             : }
      48             : 
      49             : /**
      50             :  * @function get_euler_angles
      51             :  *
      52             :  * @brief Extract the Euler angles alpha, beta, gamma
      53             :  *
      54             :  * @tparam Real The type of a real number.
      55             :  *
      56             :  * @param to_z_axis The new z-axis.
      57             :  *
      58             :  * @param to_y_axis The new y-axis.
      59             :  *
      60             :  * @return The Euler angles.
      61             :  */
      62             : 
      63             : template <typename Real>
      64           1 : inline std::array<Real, 3> get_euler_angles(std::array<Real, 3> to_z_axis,
      65             :                                             std::array<Real, 3> to_y_axis) {
      66             :     static_assert(std::is_floating_point_v<Real>);
      67             : 
      68           1 :     auto rotator = get_rotation_matrix(to_z_axis, to_y_axis);
      69           1 :     std::array<Real, 3> euler_zyz{};
      70           1 :     Eigen::Map<Eigen::Vector3<Real>>(euler_zyz.data()) = rotator.eulerAngles(2, 1, 2);
      71           2 :     return euler_zyz;
      72             : }
      73             : 
      74             : } // namespace pairinteraction::euler