Rotation and translation on filled cells. Closes #46 on github.

src/geometry.F90

@@ -142,6 +142,18 @@ contains
              p % coord => p % coord % next
              p % coord % universe = c % fill
 
+             ! Apply translation
+             if (allocated(c % translation)) then
+                p % coord % xyz = p % coord % xyz - c % translation
+             end if
+
+             ! Apply rotation
+             if (allocated(c % rotation)) then
+                p % coord % xyz = matmul(c % rotation, p % coord % xyz)
+                p % coord % uvw = matmul(c % rotation, p % coord % uvw)
+                p % coord % rotated = .true.
+             end if
+
              call find_cell(found)
              if (.not. found) exit
 

src/geometry_header.F90

@@ -64,6 +64,10 @@ module geometry_header
           & surfaces(:)    ! List of surfaces bounding cell -- note that
                            ! parentheses, union, etc operators will be listed
                            ! here too
+
+     ! Rotation matrix and translation vector
+     real(8), allocatable :: rotation(:,:)
+     real(8), allocatable :: translation(:)
   end type Cell
 
   ! array index of universe 0

src/input_xml.F90

@@ -493,6 +493,7 @@ contains
     integer :: universe_num
     integer :: n_cells_in_univ
     integer :: coeffs_reqd
+    real(8) :: phi, theta, psi
     logical :: file_exists
     character(MAX_LINE_LEN) :: filename
     character(MAX_WORD_LEN) :: word
@@ -561,7 +562,7 @@ contains
        end select
 
        ! Check to make sure that either material or fill was specified
-       if (c % material == 0 .and. c % fill == 0) then
+       if (c % material == NONE .and. c % fill == NONE) then
           message = "Neither material nor fill was specified for cell " // & 
                trim(to_str(c % id))
           call fatal_error()
@@ -569,7 +570,7 @@ contains
 
        ! Check to make sure that both material and fill haven't been
        ! specified simultaneously
-       if (c % material /= 0 .and. c % fill /= 0) then
+       if (c % material /= NONE .and. c % fill /= NONE) then
           message = "Cannot specify material and fill simultaneously"
           call fatal_error()
        end if
@@ -587,6 +588,64 @@ contains
        allocate(c % surfaces(n))
        c % surfaces = cell_(i) % surfaces
 
+       ! Rotation matrix
+       if (associated(cell_(i) % rotation)) then
+          ! Rotations can only be applied to cells that are being filled with
+          ! another universe
+          if (c % fill == NONE) then
+             message = "Cannot apply a rotation to cell " // trim(to_str(&
+                  c % id)) // " because it is not filled with another universe"
+             call fatal_error()
+          end if
+
+          ! Read number of rotation parameters
+          n = size(cell_(i) % rotation)
+          if (n /= 3) then
+             message = "Incorrect number of rotation parameters on cell " // &
+                  to_str(c % id)
+             call fatal_error()
+          end if
+
+          ! Copy rotation angles in x,y,z directions
+          phi   = -cell_(i) % rotation(1) * PI/180.0
+          theta = -cell_(i) % rotation(2) * PI/180.0
+          psi   = -cell_(i) % rotation(3) * PI/180.0
+          
+          ! Calculate rotation matrix based on angles given
+          allocate(c % rotation(3,3))
+          c % rotation = 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
+
+       ! Translation vector
+       if (associated(cell_(i) % translation)) then
+          ! Translations can only be applied to cells that are being filled with
+          ! another universe
+          if (c % fill == NONE) then
+             message = "Cannot apply a translation to cell " // trim(to_str(&
+                  c % id)) // " because it is not filled with another universe"
+             call fatal_error()
+          end if
+
+          ! Read number of translation parameters
+          n = size(cell_(i) % translation)
+          if (n /= 3) then
+             message = "Incorrect number of translation parameters on cell " &
+                  // to_str(c % id)
+             call fatal_error()
+          end if
+
+          ! Copy translation vector
+          allocate(c % translation(3))
+          c % translation = cell_(i) % translation
+       end if
+
        ! Add cell to dictionary
        call dict_add_key(cell_dict, c % id, i)
 

src/particle_header.F90

@@ -22,6 +22,9 @@ module particle_header
      real(8) :: xyz(3)
      real(8) :: uvw(3)
 
+     ! Is this level rotated?
+     logical :: rotated = .false.
+
      ! Pointer to next (more local) set of coordinates
      type(LocalCoord), pointer :: next => null()
   end type LocalCoord

src/physics.F90

@@ -149,8 +149,17 @@ contains
           ! Set all uvws to base level -- right now, after a collision, only the
           ! base level uvws are changed
           coord => p % coord0
-          do while(associated(coord))
-             coord % uvw = p % coord0 % uvw
+          do while(associated(coord % next))
+             if (coord % next % rotated) then
+                ! If next level is rotated, apply rotation matrix
+                coord % next % uvw = matmul(cells(coord % cell) % &
+                     rotation, coord % uvw)
+             else
+                ! Otherwise, copy this level's direction
+                coord % next % uvw = coord % uvw
+             end if
+
+             ! Advance coordinate level
              coord => coord % next
           end do
        end if

src/xml-fortran/templates/geometry_t.xml

@@ -11,6 +11,8 @@
     <component name="material" type="word" length="12" default="''" />
     <component name="fill" type="integer" default="0" />
     <component name="surfaces" type="integer-array" />
+    <component name="rotation" type="double-array" />
+    <component name="translation" type="double-array" />
   </typedef>
 
   <typedef name="surface_xml">