From 47c19353e4f704038c34cf3562d90947c19e6fdb Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Thu, 14 Feb 2019 15:30:00 -0600
Subject: [PATCH] Fix rotation matrix in transport()

---
 src/geometry_header.F90 |  5 -----
 src/input_xml.F90       | 35 -----------------------------------
 src/particle.cpp        |  8 +++++---
 3 files changed, 5 insertions(+), 43 deletions(-)

diff --git a/src/geometry_header.F90 b/src/geometry_header.F90
index 563a60687..afb9e2c0d 100644
--- a/src/geometry_header.F90
+++ b/src/geometry_header.F90
@@ -82,11 +82,6 @@ module geometry_header
 
   type Cell
     type(C_PTR) :: ptr
-
-    ! Rotation matrix and translation vector
-    real(8), allocatable :: rotation(:)
-    real(8), allocatable :: rotation_matrix(:,:)
-
   contains
 
     procedure :: id => cell_id
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index 574d2949b..60c6485ac 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -217,41 +217,6 @@ contains
              // to_str(c % id()))
       end if
 
-      ! Rotation matrix
-      if (check_for_node(node_cell, "rotation")) then
-        ! Rotations can only be applied to cells that are being filled with
-        ! another universe
-        if (c % fill() == C_NONE) then
-          call fatal_error("Cannot apply a rotation to cell " // trim(to_str(&
-               &c % id())) // " because it is not filled with another universe")
-        end if
-
-        ! Read number of rotation parameters
-        n = node_word_count(node_cell, "rotation")
-        if (n /= 3) then
-          call fatal_error("Incorrect number of rotation parameters on cell " &
-               // to_str(c % id()))
-        end if
-
-        ! Copy rotation angles in x,y,z directions
-        allocate(c % rotation(3))
-        call get_node_array(node_cell, "rotation", c % rotation)
-        phi   = -c % rotation(1) * PI/180.0_8
-        theta = -c % rotation(2) * PI/180.0_8
-        psi   = -c % rotation(3) * PI/180.0_8
-
-        ! Calculate rotation matrix based on angles given
-        allocate(c % rotation_matrix(3,3))
-        c % rotation_matrix = reshape((/ &
-             cos(theta)*cos(psi), cos(theta)*sin(psi), -sin(theta), &
-             -cos(phi)*sin(psi) + sin(phi)*sin(theta)*cos(psi), &
-             cos(phi)*cos(psi) + sin(phi)*sin(theta)*sin(psi), &
-             sin(phi)*cos(theta), &
-             sin(phi)*sin(psi) + cos(phi)*sin(theta)*cos(psi), &
-             -sin(phi)*cos(psi) + cos(phi)*sin(theta)*sin(psi), &
-             cos(phi)*cos(theta) /), (/ 3,3 /))
-      end if
-
       ! Add cell to dictionary
       call cell_dict % set(c % id(), i)
 
diff --git a/src/particle.cpp b/src/particle.cpp
index 26c4e8a7f..73ec98ce6 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -357,9 +357,11 @@ Particle::transport()
       for (int j = 0; j < n_coord - 1; ++j) {
         if (coord[j + 1].rotated) {
           // If next level is rotated, apply rotation matrix
-          // FIXME
-          // coord[j + 1].uvw = matmul(cells(coord[j].cell + 1) % &
-          //       rotation_matrix, coord[j].uvw)
+          const auto& m {model::cells[coord[j].cell]->rotation_};
+          Direction u {coord[j].uvw};
+          coord[j + 1].uvw[0] = m[3]*u.x + m[4]*u.y + m[5]*u.z;
+          coord[j + 1].uvw[1] = m[6]*u.x + m[7]*u.y + m[8]*u.z;
+          coord[j + 1].uvw[2] = m[9]*u.x + m[10]*u.y + m[11]*u.z;
         } else {
           // Otherwise, copy this level's direction
           std::copy(coord[j].uvw, coord[j].uvw + 3, coord[j + 1].uvw);