Add Position::rotate method

include/openmc/position.h

@@ -55,13 +55,16 @@ struct Position {
   //! Dot product of two vectors
   //! \param[in] other Vector to take dot product with
   //! \result Resulting dot product
-  inline double dot(Position other) {
+  inline double dot(Position other) const {
     return x*other.x + y*other.y + z*other.z;
   }
-  inline double norm() {
+  inline double norm() const {
     return std::sqrt(x*x + y*y + z*z);
   }
 
+  //! Rotate the position based on a rotation matrix
+  Position rotate(const std::vector<double>& rotation) const;
+
   // Data members
   double x = 0.;
   double y = 0.;

src/particle.cpp

@@ -337,9 +337,7 @@ Particle::transport()
           // If next level is rotated, apply rotation matrix
           const auto& m {model::cells[coord_[j].cell]->rotation_};
           const auto& u {coord_[j].u};
-          coord_[j + 1].u.x = m[3]*u.x + m[4]*u.y + m[5]*u.z;
-          coord_[j + 1].u.y = m[6]*u.x + m[7]*u.y + m[8]*u.z;
-          coord_[j + 1].u.z = m[9]*u.x + m[10]*u.y + m[11]*u.z;
+          coord_[j + 1].u = u.rotate(m);
         } else {
           // Otherwise, copy this level's direction
           coord_[j+1].u = coord_[j].u;

src/position.cpp

@@ -84,4 +84,14 @@ Position::operator-() const
   return {-x, -y, -z};
 }
 
+Position
+Position::rotate(const std::vector<double>& rotation) const
+{
+  return {
+    x*rotation[3] + y*rotation[4] + z*rotation[5],
+    x*rotation[6] + y*rotation[7] + z*rotation[8],
+    x*rotation[9] + y*rotation[10] + z*rotation[11]
+  };
+}
+
 } // namespace openmc