diff --git a/src/mesh.F90 b/src/mesh.F90 index 3e489903b4..37847d58ac 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -218,7 +218,9 @@ contains end subroutine count_fission_sites !=============================================================================== -! MESH_INTERSECT +! MESH_INTERSECT determines if a line between xyz0 and xyz1 intersects the outer +! boundary of the given mesh. This is important for determining whether a track +! will score to a mesh tally. !=============================================================================== function mesh_intersects(m, xyz0, xyz1) result(intersects) @@ -228,9 +230,11 @@ contains real(8), intent(in) :: xyz1(3) logical :: intersects - real(8) :: x0, y0, z0 ! track start point - real(8) :: x1, y1, z1 ! track end point - real(8) :: xi, yi, zi ! track intersection point with mesh + real(8) :: x0, y0, z0 ! track start point + real(8) :: x1, y1, z1 ! track end point + real(8) :: xi, yi, zi ! track intersection point with mesh + real(8) :: xm0, ym0, zm0 ! lower-left coordinates of mesh + real(8) :: xm1, ym1, zm1 ! upper-right coordinates of mesh ! Copy coordinates of starting point x0 = xyz0(1) @@ -242,62 +246,66 @@ contains y1 = xyz1(2) z1 = xyz1(3) + ! Copy coordinates of mesh lower_left + xm0 = m % lower_left(1) + ym0 = m % lower_left(2) + zm0 = m % lower_left(3) + + ! Copy coordinates of mesh upper_right + xm1 = m % upper_right(1) + ym1 = m % upper_right(2) + zm1 = m % upper_right(3) + ! Check if line intersects bottom surface -- calculate the intersection ! point (y,z) - xi = m % lower_left(1) - yi = y0 + (xi - x0) * (y1 - y0) / (x1 - x0) - zi = z0 + (xi - x0) * (z1 - z0) / (x1 - x0) - if (yi >= y0 .and. yi < y1 .and. zi >= z0 .and. zi < z1) then + yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0) + zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0) + if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then intersects = .true. return end if ! Check if line intersects left surface -- calculate the intersection point ! (x,z) - yi = m % lower_left(2) - xi = x0 + (yi - y0) * (x1 - x0) / (y1 - y0) - zi = z0 + (yi - y0) * (z1 - z0) / (y1 - y0) - if (xi >= x0 .and. xi < x1 .and. zi >= z0 .and. zi < z1) then + xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0) + zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0) + if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then intersects = .true. return end if ! Check if line intersects front surface -- calculate the intersection point ! (x,y) - zi = m % lower_left(3) - xi = x0 + (zi - z0) * (x1 - x0) / (z1 - z0) - yi = y0 + (zi - z0) * (y1 - y0) / (z1 - z0) - if (xi >= x0 .and. xi < x1 .and. yi >= y0 .and. yi < y1) then + xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0) + yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0) + if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then intersects = .true. return end if ! Check if line intersects top surface -- calculate the intersection ! point (y,z) - xi = m % upper_right(1) - yi = y0 + (xi - x0) * (y1 - y0) / (x1 - x0) - zi = z0 + (xi - x0) * (z1 - z0) / (x1 - x0) - if (yi >= y0 .and. yi < y1 .and. zi >= z0 .and. zi < z1) then + yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0) + zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0) + if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then intersects = .true. return end if ! Check if line intersects right surface -- calculate the intersection point ! (x,z) - yi = m % upper_right(2) - xi = x0 + (yi - y0) * (x1 - x0) / (y1 - y0) - zi = z0 + (yi - y0) * (z1 - z0) / (y1 - y0) - if (xi >= x0 .and. xi < x1 .and. zi >= z0 .and. zi < z1) then + xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0) + zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0) + if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then intersects = .true. return end if ! Check if line intersects back surface -- calculate the intersection point ! (x,y) - zi = m % upper_right(3) - xi = x0 + (zi - z0) * (x1 - x0) / (z1 - z0) - yi = y0 + (zi - z0) * (y1 - y0) / (z1 - z0) - if (xi >= x0 .and. xi < x1 .and. yi >= y0 .and. yi < y1) then + xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0) + yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0) + if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then intersects = .true. return end if