From 70e6eaf4583062221c90f587073bb30c36d99f20 Mon Sep 17 00:00:00 2001
From: Paul Romano <romano7@gmail.com>
Date: Thu, 8 Sep 2011 18:16:42 -0400
Subject: [PATCH] Initial code for reflective boundary conditions on planes --
 doesn't currently work.

---
 src/geometry.f90 | 59 ++++++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 57 insertions(+), 2 deletions(-)

diff --git a/src/geometry.f90 b/src/geometry.f90
index c4c6771704..4a1051ce8c 100644
--- a/src/geometry.f90
+++ b/src/geometry.f90
@@ -178,6 +178,13 @@ contains
 
     integer                 :: i          ! index of neighbors
     integer                 :: index_cell ! index in cells array
+    real(8)                 :: u          ! x-component of direction
+    real(8)                 :: v          ! y-component of direction
+    real(8)                 :: w          ! z-component of direction
+    real(8)                 :: n1         ! x-component of surface normal
+    real(8)                 :: n2         ! y-component of surface normal
+    real(8)                 :: n3         ! z-component of surface normal
+    real(8)                 :: norm       ! "norm" of surface normal
     logical                 :: found      ! particle found in universe?
     character(MAX_LINE_LEN) :: msg        ! output/error message?
     type(Surface),  pointer :: surf
@@ -190,14 +197,58 @@ contains
        call message(msg)
     end if
 
-    ! check for leakage
     if (surf%bc == BC_VACUUM) then
+       ! =======================================================================
+       ! PARTICLE LEAKS OUT OF PROBLEM
+
        p%alive = .false.
        if (verbosity >= 10) then
           msg = "    Leaked out of surface " // trim(int_to_str(surf%uid))
           call message(msg)
        end if
        return
+
+    elseif (surf%bc == BC_REFLECT) then
+       ! =======================================================================
+       ! PARTICLE REFLECTS FROM SURFACE
+
+       ! Copy particle's direction cosines
+       u = p % uvw(1)
+       v = p % uvw(2)
+       w = p % uvw(3)
+
+       select case (surf%type)
+       case (SURF_PX)
+          p % uvw = (/ -u, v, w /)
+       case (SURF_PY)
+          p % uvw = (/ u, -v, w /)
+       case (SURF_PZ)
+          p % uvw = (/ u, v, -w /)
+       case (SURF_PLANE)
+          ! Find surface coefficients and norm of vector normal to surface
+          n1 = surf % coeffs(1)
+          n2 = surf % coeffs(2)
+          n3 = surf % coeffs(3)
+          norm = n1*n1 + n2*n2 + n3*n3
+
+          ! Reflect direction according to normal
+          u = u - 2*u*n1*n1/norm
+          v = v - 2*v*n2*n2/norm
+          w = w - 2*w*n3*n3/norm
+
+          ! Set vector
+          p % uvw = (/ u, v, w /)
+       case default
+          msg = "Reflection not supported for surface " // &
+               trim(int_to_str(surf % uid))
+          call fatal_error(msg)
+       end select
+
+       if (verbosity >= 10) then
+          msg = "    Reflected from surface " // trim(int_to_str(surf%uid))
+          call message(msg)
+       end if
+       return
     end if
 
     if (p%surface > 0 .and. allocated(surf%neighbor_pos)) then
@@ -452,7 +503,11 @@ contains
           y = p % xyz(2)
           z = p % xyz(3)
        end if
-       ! check for coincident surface
+
+       ! check for coincident surface -- note that we can't skip the calculation
+       ! in general because a particle could be on one side of a cylinder and
+       ! still have a positive distance to the other
+
        index_surf = expression(i)
        if (index_surf == current_surf) then
           on_surface = .true.