diff --git a/src/geometry.F90 b/src/geometry.F90 index 5e6daac25e..07ce177020 100644 --- a/src/geometry.F90 +++ b/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 diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index fbcf0353a9..58f7d02878 100644 --- a/src/geometry_header.F90 +++ b/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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index c056d2cc90..8eb26a57b1 100644 --- a/src/input_xml.F90 +++ b/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) diff --git a/src/particle_header.F90 b/src/particle_header.F90 index 7ce921bca1..8527419792 100644 --- a/src/particle_header.F90 +++ b/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 diff --git a/src/physics.F90 b/src/physics.F90 index 73eb40cce1..18ae23496d 100644 --- a/src/physics.F90 +++ b/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 diff --git a/src/xml-fortran/templates/geometry_t.xml b/src/xml-fortran/templates/geometry_t.xml index fcae61b709..ae8105a6bd 100644 --- a/src/xml-fortran/templates/geometry_t.xml +++ b/src/xml-fortran/templates/geometry_t.xml @@ -11,6 +11,8 @@ + +