diff --git a/src/common/geometry.hh b/src/common/geometry.hh
index 72c58ff..cb0e281 100644
--- a/src/common/geometry.hh
+++ b/src/common/geometry.hh
@@ -1,507 +1,507 @@
 /**
  * @file   geometry.hh
  *
  * @author Till Junge <till.junge@altermail.ch>
  *
  * @date   18 Apr 2018
  *
  * @brief  Geometric calculation helpers
  *
  * Copyright © 2018 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include "common/common.hh"
 #include "common/tensor_algebra.hh"
 #include "common/eigen_tools.hh"
 
 #include <Eigen/Dense>
 #include <Eigen/Geometry>
 
 #include <array>
 
 namespace muSpectre {
 
   /**
    * The rotation matrices depend on the order in which we rotate
    * around different axes. See [[
    * https://en.wikipedia.org/wiki/Euler_angles#Rotation_matrix ]] to
    * find the matrices
    */
   enum class RotationOrder {Z,
       XZXEuler, XYXEuler, YXYEuler, YZYEuler, ZYZEuler, ZXZEuler,
       XZYTaitBryan, XYZTaitBryan, YXZTaitBryan, YZXTaitBryan, ZYXTaitBryan,
       ZXYTaitBryan};
 
   namespace internal {
 
     template <Dim_t Dim>
     struct DefaultOrder
     {
       constexpr static RotationOrder value{RotationOrder::ZXYTaitBryan};
     };
 
     template <>
     struct DefaultOrder<twoD>
     {
       constexpr static RotationOrder value{RotationOrder::Z};
     };
 
   }  // internal
 
   template <Dim_t Dim>
   class RotatorBase{
   public:
     using RotMat_t = Eigen::Matrix<Real, Dim, Dim>;
     //! Default constructor
     RotatorBase() = delete;
 
     RotatorBase(const RotMat_t rotation_matrix_input):
       rot_mat{rotation_matrix_input}
     {}
     //! Copy constructor
     RotatorBase(const RotatorBase &other) = default;
 
     //! Move constructor
     RotatorBase(RotatorBase &&other) = default;
 
     //! Destructor
     virtual ~RotatorBase() = default;
 
     //! Copy assignment operator
     RotatorBase& operator=(const RotatorBase &other) = default;
 
     //! Move assignment operator
     RotatorBase& operator=(RotatorBase &&other) = default;
 
     /**
      * Applies the rotation into the frame define my the rotation
      * matrix
      *
      * @param input is a first-, second-, or fourth-rank tensor
      * (column vector, square matrix, or T4Matrix, or a Eigen::Map of
      * either of these, or an expression that evaluates into any of
      * these)
      */
     template<class In_t>
     inline decltype(auto) rotate(In_t && input);
 
     /**
      * Applies the rotation back out from the frame define my the
      * rotation matrix
      *
      * @param input is a first-, second-, or fourth-rank tensor
      * (column vector, square matrix, or T4Matrix, or a Eigen::Map of
      * either of these, or an expression that evaluates into any of
      * these)
      */
     template<class In_t>
     inline decltype(auto) rotate_back(In_t && input);
 
     const RotMat_t & get_rot_mat() const {return this->rot_mat;}
 
   protected:
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
     RotMat_t rot_mat;
 
   };
 
 
 
 
   template <Dim_t Dim, RotationOrder Order=internal::DefaultOrder<Dim>::value>
   class RotatorAngle: public RotatorBase<Dim>
   {
     static_assert(((Dim == twoD) and (Order == RotationOrder::Z)) or
                   ((Dim == threeD) and (Order != RotationOrder::Z)),
                   "In 2d, only order 'Z' makes sense. In 3d, it doesn't");
 
   public:
     using Parent = RotatorBase<Dim>;
     using Angles_t = Eigen::Matrix<Real, (Dim == twoD) ? 1 :3, 1>;
     using RotMat_t = Eigen::Matrix<Real, Dim, Dim>;
 
 
     //! Default constructor
     RotatorAngle() = delete;
 
     RotatorAngle(const Eigen::Ref<const Angles_t> & angles_inp)
       :Parent(this->compute_rotation_matrix_angle(angles_inp)),
        angles{angles_inp}
     {}
 
     //! Copy constructor
     RotatorAngle(const RotatorAngle &other) = default;
 
     //! Move constructor
     RotatorAngle(RotatorAngle &&other) = default;
 
     //! Destructor
     virtual ~RotatorAngle() = default;
 
     //! Copy assignment operator
     RotatorAngle& operator=(const RotatorAngle &other) = default;
 
     //! Move assignment operator
     RotatorAngle& operator=(RotatorAngle &&other) = default;
 
   protected:
     inline RotMat_t compute_rotation_matrix_angle(Angles_t angles);
     inline RotMat_t compute_this_rotation_matrix_angle(){
       return compute_rotation_matrix(this->angles);
     };
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
     Angles_t angles;
 
   private:
   };
 
     /* ---------------------------------------------------------------------- */
   namespace internal {
 
     template <RotationOrder Order, Dim_t Dim>
     struct RotationMatrixComputerAngle
     {};
 
     template <RotationOrder Order>
     struct RotationMatrixComputerAngle<Order, twoD> {
       constexpr static Dim_t Dim{twoD};
       using RotMat_t = typename RotatorAngle<Dim, Order>::RotMat_t;
       using Angles_t = typename RotatorAngle<Dim, Order>::Angles_t;
 
       inline static decltype(auto) compute(const Eigen::Ref<Angles_t> & angles) {
         static_assert(Order == RotationOrder::Z,
                       "Two-d rotations can only be around the z axis");
         return RotMat_t(Eigen::Rotation2Dd(angles(0)));
       }
     };
 
     template <RotationOrder Order>
     struct RotationMatrixComputerAngle<Order, threeD> {
       constexpr static Dim_t Dim{threeD};
       using RotMat_t = typename RotatorAngle<Dim, Order>::RotMat_t;
       using Angles_t = typename RotatorAngle<Dim, Order>::Angles_t;
 
       inline static decltype(auto) compute(const Eigen::Ref<Angles_t> & angles) {
         static_assert(Order != RotationOrder::Z,
                       "three-d rotations cannot only be around the z axis");
 
         switch (Order) {
         case RotationOrder::ZXZEuler: {
           return RotMat_t((Eigen::AngleAxisd(angles(0), Eigen::Vector3d::UnitZ()) *
                            Eigen::AngleAxisd(angles(1), Eigen::Vector3d::UnitX()) *
                            Eigen::AngleAxisd(angles(2), Eigen::Vector3d::UnitZ())));
           break;
         }
         case RotationOrder::ZXYTaitBryan: {
           return RotMat_t((Eigen::AngleAxisd(angles(0), Eigen::Vector3d::UnitZ()) *
                            Eigen::AngleAxisd(angles(1), Eigen::Vector3d::UnitX()) *
                            Eigen::AngleAxisd(angles(2), Eigen::Vector3d::UnitY())));
         }
         default: {
           throw std::runtime_error("not yet implemented.");
         }
         }
       }
     };
 
   }  // internal
   /* ---------------------------------------------------------------------- */
   template <Dim_t Dim, RotationOrder Order>
   auto RotatorAngle<Dim, Order>::
   compute_rotation_matrix_angle(Angles_t angles) -> RotMat_t {
     return internal::RotationMatrixComputerAngle<Order, Dim>::compute(angles);
   }
 
   /* ---------------------------------------------------------------------- */
   /**
    * this class is used to make the vector a aligned to the vec b by means of a
    rotation system, the input for the constructor is the vector itself and the
    functions
    rotate and rotate back would be available as they exist in the parent class
    (RotatorBase) nad can be used in order to do the functionality of the class
   */
   template <Dim_t Dim>
   class RotatorTwoVec: public RotatorBase<Dim>
   {
   public:
     using Parent = RotatorBase<Dim>;
     using Vec_t = Eigen::Matrix<Real, (Dim == twoD) ? 2 :3, 1>;
     using RotMat_t = Eigen::Matrix<Real, Dim, Dim>;
 
 
     //! Default constructor
     RotatorTwoVec() = delete;
 
     RotatorTwoVec(Vec_t  vec_a_inp,
                   Vec_t  vec_b_inp)
       :Parent(this->compute_rotation_matrix_TwoVec(vec_a_inp,
                                                    vec_b_inp)),
        vec_ref{vec_a_inp},
        vec_des{vec_b_inp}
     {}
 
     //! Copy constructor
     RotatorTwoVec(const RotatorTwoVec &other) = default;
 
     //! Move constructor
     RotatorTwoVec(RotatorTwoVec &&other) = default;
 
     //! Destructor
     virtual ~RotatorTwoVec() = default;
 
     //! Copy assignment operator
     RotatorTwoVec& operator=(const RotatorTwoVec &other) = default;
 
     //! Move assignment operator
     RotatorTwoVec& operator=(RotatorTwoVec &&other) = default;
 
   protected:
     inline RotMat_t compute_rotation_matrix_TwoVec(Vec_t vec_ref, Vec_t vec_des);
     inline RotMat_t compute_this_rotation_matrix_TwoVec(){
       return compute_rotation_matrix(this->vec_ref, this->vec_des);
     }
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
     Vec_t vec_ref, vec_des;
   private:
   };
   /* ---------------------------------------------------------------------- */
   namespace internal {
     template <Dim_t Dim>
     struct RotationMatrixComputerTwoVec
     {};
     template<>
     struct RotationMatrixComputerTwoVec<twoD> {
       constexpr static Dim_t Dim{twoD};
       using RotMat_t = typename RotatorTwoVec<Dim>::RotMat_t;
       using Vec_t    = typename RotatorTwoVec<Dim>::Vec_t;
 
       inline static RotMat_t compute(Vec_t vec_ref,
                                      Vec_t vec_des) {
         Real v_ref_norm = sqrt(vec_ref(0)*vec_ref(0) + vec_ref(1)*vec_ref(1));
         Real v_des_norm = sqrt(vec_des(0)*vec_des(0) + vec_des(1)*vec_des(1));
         RotMat_t ret_mat ;
         ret_mat(0,0) = ret_mat(1,1) = (((vec_ref(0)/v_ref_norm)*
                                         (vec_des(0)/v_des_norm))+
                                        ((vec_des(1)/v_des_norm)*
                                         (vec_ref(1)/v_ref_norm)));
 
         ret_mat(1,0) =   (((vec_ref(0)/v_ref_norm)*
                            (vec_des(1)/v_des_norm))-
                           ((vec_des(0)/v_des_norm)*
                            (vec_ref(1)/v_ref_norm)));
         ret_mat(0,1 ) = -ret_mat(1,0);
         return ret_mat;
       }
     };
 
     template<>
     struct RotationMatrixComputerTwoVec<threeD> {
       constexpr static Dim_t Dim{threeD};
       using RotMat_t = typename RotatorTwoVec<Dim>::RotMat_t;
       using Vec_t    = typename RotatorTwoVec<Dim>::Vec_t;
       inline static RotMat_t compute(Vec_t vec_ref,
                                      Vec_t vec_des) {
         return Eigen::Quaternion<double>::FromTwoVectors(vec_ref, vec_des).
           normalized().toRotationMatrix();
       }
-
     };
+
   } //internal
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t Dim>
   auto RotatorTwoVec<Dim>::
   compute_rotation_matrix_TwoVec(Vec_t vec_ref,
                                  Vec_t vec_des) -> RotMat_t {
     return internal::RotationMatrixComputerTwoVec<Dim>::compute(vec_ref, vec_des);
   }
 
   /* ---------------------------------------------------------------------- */
   /**
    * this class is used to make a vector aligned to x-axis of the coordinate
    system, the input for the constructor is the vector itself and the functions
    rotate and rotate back would be available as they exist in the parent class
    (RotatorBase) nad can be used in order to do the functionality of the class
    */
   template <Dim_t Dim>
   class RotatorNormal: public RotatorBase<Dim>
   {
   public:
     using Parent = RotatorBase<Dim>;
     using Vec_t = Eigen::Matrix<Real, Dim, 1>;
     using RotMat_t = Eigen::Matrix<Real, Dim, Dim>;
 
     //! Default constructor
     RotatorNormal() = delete;
 
     RotatorNormal(Vec_t vec)
       :Parent(this->compute_rotation_matrix_normal(vec)),
        vec{vec}
     {}
 
     //! Copy constructor
     RotatorNormal(const RotatorNormal &other) = default;
 
     //! Move constructor
     RotatorNormal(RotatorNormal &&other) = default;
 
     //! Destructor
     virtual ~RotatorNormal() = default;
 
     //! Copy assignment operator
     RotatorNormal& operator=(const RotatorNormal &other) = default;
 
     //! Move assignment operator
     RotatorNormal& operator=(RotatorNormal &&other) = default;
 
   protected:
     inline RotMat_t compute_rotation_matrix_normal(Vec_t vec);
     inline RotMat_t compute_this_rotation_matrix_normal(){
       return compute_rotation_matrix_normal(this->vec);
     }
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
     Vec_t vec;
 
   private:
   };
 
 
     /* ---------------------------------------------------------------------- */
   namespace internal {
     template <Dim_t Dim>
     struct RotationMatrixComputerNormal
     {};
     template<>
     struct RotationMatrixComputerNormal<twoD> {
       constexpr static Dim_t Dim{twoD};
       using RotMat_t = typename RotatorTwoVec<Dim>::RotMat_t;
       using Vec_t    = typename RotatorTwoVec<Dim>::Vec_t;
 
       inline static RotMat_t compute(Vec_t vec) {
         Real vec_norm = sqrt(vec(0)*vec(0) + vec(1)*vec(1));
         Vec_t x; x<<1.0,
                    0.0;
 
         RotMat_t ret_mat;
         ret_mat(0,0) = ret_mat(1,1) = ((vec(0)/vec_norm)* x(0));
         ret_mat(1,0) = -(-(vec(1)/vec_norm)* x(0));
         ret_mat(0,1) = - ret_mat(1,0);
         return ret_mat;
       }
     };
 
     template<>
     struct RotationMatrixComputerNormal<threeD> {
       constexpr static Dim_t Dim{threeD};
       using RotMat_t = typename RotatorTwoVec<Dim>::RotMat_t;
       using Vec_t    = typename RotatorTwoVec<Dim>::Vec_t;
       inline static RotMat_t compute(Vec_t vec) {
         Real eps = 0.1 ;
         Vec_t vec1 = vec/vec.norm();
         Vec_t x(Vec_t::UnitX());
         Vec_t y(Vec_t::UnitY());
         Vec_t n_x = vec1.cross(x);
         Vec_t vec2 =  ((n_x.norm()>eps)   * n_x +
                       (1-(n_x.norm()>eps))*(vec1.cross(y)));
         Vec_t vec3 = vec1.cross(vec2);
         RotMat_t ret_mat;
         ret_mat <<vec1(0), vec2(0)/vec2.norm(), vec3(0)/vec3.norm(),
                   vec1(1), vec2(1)/vec2.norm(), vec3(1)/vec3.norm(),
                   vec1(2), vec2(2)/vec2.norm(), vec3(2)/vec3.norm();
         return ret_mat;
 
       }
 
     };
   } //internal
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t Dim>
   auto RotatorNormal<Dim>::
   compute_rotation_matrix_normal(Vec_t vec) -> RotMat_t{
     return internal::RotationMatrixComputerNormal<Dim>::compute(vec);
   }
 
   /* ---------------------------------------------------------------------- */
 
   namespace internal {
 
     template<Dim_t Rank>
     struct RotationHelper {
     };
 
     /* ---------------------------------------------------------------------- */
     template <>
     struct RotationHelper<firstOrder> {
       template <class In_t, class Rot_t>
       inline static decltype(auto) rotate(In_t && input,
                                    Rot_t && R) {
         return R * input;
       }
     };
 
     /* ---------------------------------------------------------------------- */
     template <>
     struct RotationHelper<secondOrder> {
       template <class In_t, class Rot_t>
       inline static decltype(auto) rotate(In_t && input,
                                    Rot_t && R) {
         return R * input * R.transpose();
       }
     };
 
     /* ---------------------------------------------------------------------- */
     template <>
     struct RotationHelper<fourthOrder> {
       template <class In_t, class Rot_t>
       inline static decltype(auto) rotate(In_t && input,
                                    Rot_t && R) {
         constexpr Dim_t Dim{EigenCheck::tensor_dim<Rot_t>::value};
         auto && rotator_forward{Matrices::outer_under(R.transpose(),
                                                       R.transpose())};
         auto && rotator_back = Matrices::outer_under(R, R);
 
         // unclear behaviour. When I return this value as an
         // expression, clange segfaults or returns an uninitialised
         // tensor
         return T4Mat<Real, Dim>(rotator_back * input * rotator_forward);
       }
     };
   }
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t Dim>
   template <class In_t>
   auto RotatorBase<Dim>::rotate(In_t && input) -> decltype(auto) {
     constexpr Dim_t tensor_rank{EigenCheck::tensor_rank<In_t, Dim>::value};
 
     return internal::RotationHelper<tensor_rank>::
       rotate(std::forward<In_t>(input), this->rot_mat);
   }
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t Dim>
   template <class In_t>
   auto RotatorBase<Dim>::rotate_back(In_t && input) -> decltype(auto) {
     constexpr Dim_t tensor_rank{EigenCheck::tensor_rank<In_t, Dim>::value};
 
     return internal::RotationHelper<tensor_rank>::
       rotate(std::forward<In_t>(input), this->rot_mat.transpose());
   }
 
 }  // muSpectre
diff --git a/tests/test_geometry.cc b/tests/test_geometry.cc
new file mode 100644
index 0000000..51444b1
--- /dev/null
+++ b/tests/test_geometry.cc
@@ -0,0 +1,340 @@
+/**
+ * file   test_geometry.cc
+ *
+ * @author Till Junge <till.junge@epfl.ch>
+ *
+ * @date   19 Apr 2018
+ *
+ * @brief  Tests for tensor rotations
+ *
+ * @section LICENSE
+ *
+ * Copyright © 2018 Till Junge
+ *
+ * µSpectre is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License as
+ * published by the Free Software Foundation, either version 3, or (at
+ * your option) any later version.
+ *
+ * µSpectre is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with GNU Emacs; see the file COPYING. If not, write to the
+ * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+ * Boston, MA 02111-1307, USA.
+ */
+
+#include "common/geometry.hh"
+#include "tests.hh"
+#include "test_goodies.hh"
+#include "common/T4_map_proxy.hh"
+
+#include <Eigen/Dense>
+#include <boost/mpl/list.hpp>
+
+#include <cmath>
+
+namespace muSpectre {
+
+  enum class IsAligned{
+    yes,
+    no
+  };
+
+  BOOST_AUTO_TEST_SUITE(geometry);
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t Dim_>
+  struct RotationFixture {
+    static constexpr Dim_t Dim{Dim_};
+    using Vec_t = Eigen::Matrix<Real, Dim, 1>;
+    using Mat_t = Eigen::Matrix<Real, Dim, Dim>;
+    using Ten_t = T4Mat<Real, Dim>;
+    static constexpr Dim_t get_Dim() {return Dim_;}
+    using Rot_t = RotatorBase<Dim>;
+    RotationFixture(Rot_t rot_mat_inp):
+      rot_mat{rot_mat_inp},
+      rotator(rot_mat)
+    {}
+    Vec_t v{Vec_t::Random()};
+    Mat_t m{Mat_t::Random()};
+    Ten_t t{Ten_t::Random()};
+    Mat_t rot_mat;
+    Rot_t rotator;
+    testGoodies::RandRange<Real> rr{};
+
+  };
+
+
+  template <Dim_t Dim_, RotationOrder Rot>
+  struct RotationAngleFixture {
+    static constexpr Dim_t Dim{Dim_};
+    using Parent = RotationFixture<Dim>;
+    using Vec_t = Eigen::Matrix<Real, Dim, 1>;
+    using Mat_t = Eigen::Matrix<Real, Dim, Dim>;
+    using Ten_t = T4Mat<Real, Dim>;
+    using Angles_t = Eigen::Matrix<Real, (Dim == threeD ? 3: 1), 1>;
+    using RotAng_t = RotatorAngle<Dim, Rot>;
+    static constexpr RotationOrder EulerOrder{Rot};
+    static constexpr Dim_t get_Dim() {return Dim_;}
+    RotationAngleFixture():
+      rotator{euler}
+      {}
+    Vec_t v{Vec_t::Random()};
+    Mat_t m{Mat_t::Random()};
+    Ten_t t{Ten_t::Random()};
+    Angles_t euler{2*pi*Angles_t::Random()};
+    RotatorAngle<Dim, Rot> rotator;
+  };
+
+  template <Dim_t Dim_, IsAligned is_aligned = IsAligned::no>
+  struct RotationTwoVecFixture {
+    static constexpr Dim_t Dim{Dim_};
+    using Parent = RotationFixture<Dim>;
+    using Vec_t = Eigen::Matrix<Real, Dim, 1>;
+    using Mat_t = Eigen::Matrix<Real, Dim, Dim>;
+    using Ten_t = T4Mat<Real, Dim>;
+    static constexpr Dim_t get_Dim() {return Dim_;}
+    RotationTwoVecFixture():
+      vec_ref{this->ref_vec_maker()},
+      vec_des{this->des_vec_maker()},
+      rotator(vec_ref, vec_des)
+    {}
+
+    Vec_t ref_vec_maker(){
+      Vec_t ret_vec{Vec_t::Random()};
+      return ret_vec/ret_vec.norm();
+    }
+    Vec_t des_vec_maker(){
+      if (is_aligned == IsAligned::yes) {
+        return -this->vec_ref;
+      } else {
+        Vec_t ret_vec{Vec_t::Random()};
+        return ret_vec/ret_vec.norm();
+      }
+    }
+    Vec_t v{Vec_t::Random()};
+    Mat_t m{Mat_t::Random()};
+    Ten_t t{Ten_t::Random()};
+    Vec_t vec_ref{Vec_t::Random()};
+    Vec_t vec_des{Vec_t::Random()};
+    RotatorTwoVec<Dim> rotator;
+  };
+
+
+  template <Dim_t Dim_, IsAligned is_aligned = IsAligned::no>
+  struct RotationNormalFixture {
+    static constexpr Dim_t Dim{Dim_};
+    using Parent = RotationFixture<Dim>;
+    using Vec_t = Eigen::Matrix<Real, Dim, 1>;
+    using Mat_t = Eigen::Matrix<Real, Dim, Dim>;
+    using Ten_t = T4Mat<Real, Dim>;
+    static constexpr Dim_t get_Dim() {return Dim_;}
+    RotationNormalFixture():
+      vec_norm{this->vec_maker()},
+      rotator(vec_norm)
+    {}
+
+    Vec_t vec_maker(){
+      if (is_aligned == IsAligned::yes) {
+        return -Vec_t::UnitX();
+      } else {
+        Vec_t ret_vec{Vec_t::Random()};
+        return ret_vec/ret_vec.norm();
+      }
+  }
+    Vec_t v{Vec_t::Random()};
+    Mat_t m{Mat_t::Random()};
+    Ten_t t{Ten_t::Random()};
+    Vec_t vec_norm;
+    RotatorNormal<Dim> rotator;
+  };
+  /* ---------------------------------------------------------------------- */
+  using fix_list = boost::mpl::list
+    <RotationAngleFixture<  twoD, RotationOrder::Z>,
+     RotationAngleFixture<threeD, RotationOrder::ZXYTaitBryan>,
+     RotationNormalFixture<  twoD>,
+     RotationNormalFixture<threeD>,
+     RotationTwoVecFixture<  twoD>,
+     RotationTwoVecFixture<threeD>>;
+
+  /* ---------------------------------------------------------------------- */
+  BOOST_FIXTURE_TEST_CASE_TEMPLATE(rotation_test, Fix, fix_list, Fix) {
+
+    using Vec_t = typename Fix::Vec_t;
+    using Mat_t = typename Fix::Mat_t;
+    using Ten_t = typename Fix::Ten_t;
+
+    constexpr const Dim_t Dim{Fix::get_Dim()};
+
+    Vec_t &v{Fix::v};
+    Mat_t &m{Fix::m};
+    Ten_t &t{Fix::t};
+    const Mat_t &R{Fix::rotator.get_rot_mat()};
+
+
+    Vec_t v_ref{R * v};
+    Mat_t m_ref{R * m * R.transpose()};
+    Ten_t t_ref{Ten_t::Zero()};
+    for (int i = 0; i < Dim; ++i) {
+      for (int a = 0; a < Dim; ++a) {
+        for (int l = 0; l < Dim; ++l) {
+          for (int b = 0; b < Dim; ++b) {
+            for (int m = 0; m < Dim; ++m) {
+              for (int n = 0; n < Dim; ++n) {
+                for (int o = 0; o < Dim; ++o) {
+                  for (int p = 0; p < Dim; ++p) {
+                    get(t_ref, a, b, o, p) +=
+                      R(a, i) * R(b, l) * get(t, i, l, m, n) * R(o, m) * R(p, n);
+                  }
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+
+    Vec_t v_rotator(Fix::rotator.rotate(v));
+    Mat_t m_rotator(Fix::rotator.rotate(m));
+    Ten_t t_rotator(Fix::rotator.rotate(t));
+
+    auto v_error{(v_rotator-v_ref).norm()/v_ref.norm()};
+    BOOST_CHECK_LT(v_error, tol);
+
+    auto m_error{(m_rotator-m_ref).norm()/m_ref.norm()};
+    BOOST_CHECK_LT(m_error, tol);
+
+    auto t_error{(t_rotator-t_ref).norm()/t_ref.norm()};
+    BOOST_CHECK_LT(t_error, tol);
+    if (t_error >= tol) {
+      std::cout << "t4_reference:" << std::endl
+                << t_ref << std::endl;
+      std::cout << "t4_rotator:" << std::endl
+                << t_rotator << std::endl;
+    }
+
+    Vec_t v_back{Fix::rotator.rotate_back(v_rotator)};
+    Mat_t m_back{Fix::rotator.rotate_back(m_rotator)};
+    Ten_t t_back{Fix::rotator.rotate_back(t_rotator)};
+
+    v_error = (v_back-v).norm()/v.norm();
+    BOOST_CHECK_LT(v_error, tol);
+
+    m_error = (m_back-m).norm()/m.norm();
+    BOOST_CHECK_LT(m_error, tol);
+
+    t_error = (t_back-t).norm()/t.norm();
+    BOOST_CHECK_LT(t_error, tol);
+
+  }
+
+  /* ---------------------------------------------------------------------- */
+  using threeD_list = boost::mpl::list
+    <RotationAngleFixture<threeD, RotationOrder::ZXYTaitBryan>>;
+
+  /* ---------------------------------------------------------------------- */
+  BOOST_FIXTURE_TEST_CASE_TEMPLATE(rotation_matrix_test,
+                                   Fix, threeD_list, Fix) {
+    using Mat_t = typename Fix::Mat_t;
+    auto c{Eigen::cos(Fix::euler.array())};
+    Real c_1{c[0]}, c_2{c[1]}, c_3{c[2]};
+    auto s{Eigen::sin(Fix::euler.array())};
+    Real s_1{s[0]}, s_2{s[1]}, s_3{s[2]};
+    Mat_t rot_ref;
+
+
+    switch (Fix::EulerOrder) {
+    case RotationOrder::ZXYTaitBryan: {
+      rot_ref <<  c_1 * c_3 - s_1 * s_2 * s_3, - c_2 * s_1, c_1 * s_3 + c_3 * s_1 * s_2 ,
+        c_3 * s_1 + c_1 * s_2 * s_3 , c_1 * c_2, s_1 * s_3 - c_1 * c_3 * s_2 ,
+        - c_2 * s_3, s_2, c_2 * c_3;
+
+      break;
+    }
+    default: {
+      BOOST_CHECK(false);
+      break;
+    }
+    }
+    auto err{(rot_ref - Fix::rotator.get_rot_mat()).norm()};
+    BOOST_CHECK_LT(err, tol);
+    if (err >= tol) {
+      std::cout << "Reference:" << std::endl << rot_ref << std::endl;
+      std::cout << "Rotator:" << std::endl << Fix::rotator.get_rot_mat()
+                << std::endl;
+    }
+  }
+
+
+  /* ---------------------------------------------------------------------- */
+  using twovec_list = boost::mpl::list<
+    RotationTwoVecFixture<threeD>,
+     RotationTwoVecFixture<  twoD>,
+    RotationTwoVecFixture<threeD, IsAligned::yes>,
+    RotationTwoVecFixture<  twoD, IsAligned::yes>>;
+
+  /* ---------------------------------------------------------------------- */
+  BOOST_FIXTURE_TEST_CASE_TEMPLATE(rotation_twovec_test,
+                                   Fix, twovec_list, Fix) {
+    using Vec_t = typename Fix::Vec_t;
+    Vec_t vec_ref{Fix::vec_ref};
+    Vec_t vec_des{Fix::vec_des};
+    Vec_t vec_res{Fix::rotator.rotate(vec_ref)};
+    Vec_t vec_back{Fix::rotator.rotate_back(vec_res)};
+
+    auto err_f{(vec_res - vec_des).norm()};
+    BOOST_CHECK_LT(err_f, tol);
+    if (err_f >= tol) {
+      std::cout << "Destination:" << std::endl << vec_des << std::endl;
+      std::cout << "Rotated:" << std::endl << vec_res
+                << std::endl;
+    }
+    auto err_b{(vec_back - vec_ref).norm()};
+    BOOST_CHECK_LT(err_b, tol);
+    if (err_b >= tol) {
+      std::cout << "Refrence:" << std::endl << vec_ref << std::endl;
+      std::cout << "Rotated Back:" << std::endl << vec_back
+                << std::endl;
+    }
+  }
+
+  /* ---------------------------------------------------------------------- */
+  using normal_list = boost::mpl::list<
+    RotationNormalFixture<threeD>,
+    RotationNormalFixture<  twoD>,
+    RotationNormalFixture<threeD, IsAligned::yes>,
+    RotationNormalFixture<  twoD, IsAligned::yes>>;
+
+  /* ---------------------------------------------------------------------- */
+  BOOST_FIXTURE_TEST_CASE_TEMPLATE(rotation_normal_test,
+                                   Fix, normal_list, Fix) {
+    using Vec_t = typename Fix::Vec_t;
+    Vec_t vec_ref{Fix::vec_norm};
+    Vec_t vec_des{Vec_t::UnitX()};
+    Vec_t vec_res{Fix::rotator.rotate_back(vec_ref)};
+    Vec_t vec_back{Fix::rotator.rotate(vec_res)};
+
+    auto err_f{(vec_res - vec_des).norm()};
+    BOOST_CHECK_LT(err_f, tol);
+    if (err_f >= tol) {
+      std::cout << "Destination:" << std::endl << vec_des << std::endl;
+      std::cout << "Rotated:" << std::endl << vec_res
+                << std::endl;
+    }
+    auto err_b{(vec_back - vec_ref).norm()};
+    BOOST_CHECK_LT(err_b, tol);
+    if (err_b >= tol) {
+      std::cout << "Refrence:" << std::endl << vec_ref << std::endl;
+      std::cout << "Rotated Back:" << std::endl << vec_back
+                << std::endl;
+    }
+  }
+
+
+  BOOST_AUTO_TEST_SUITE_END();
+
+}  // muSpectre