Improve check for coincidence in sense().

The previous solution was to move the particle forward artificially and then
call sense() recursively. A better way is to look at the projection of the
particle's direction on the surface normal from which the sign determines which
side the particle will be on when it moves forward. This enables us to make
sense() a pure non-recursive function.

src/surface_header.F90

@@ -1,6 +1,6 @@
 module surface_header
 
-  use constants, only: ONE, TWO, ZERO, INFINITY, FP_COINCIDENT, TINY_BIT
+  use constants, only: ONE, TWO, ZERO, INFINITY, FP_COINCIDENT
 
   implicit none
 
@@ -178,12 +178,13 @@ contains
 !===============================================================================
 ! SENSE determines whether a point is on the 'positive' or 'negative' side of a
 ! surface. This routine is crucial for determining what cell a particular point
-! is in.
+! is in. The positive side is indicated by a returned value of .true. and the
+! negative side is indicated by a returned value of .false.
 !===============================================================================
 
-  recursive function sense(this, xyz, uvw) result(s)
+  pure function sense(this, xyz, uvw) result(s)
     class(Surface), intent(in) :: this  ! surface
-    real(8), intent(inout) :: xyz(3)
+    real(8), intent(in) :: xyz(3)
     real(8), intent(in) :: uvw(3)
     logical :: s  ! sense of particle
 
@@ -195,16 +196,13 @@ contains
 
     ! Check which side of surface the point is on
     if (abs(f) < FP_COINCIDENT) then
-      ! Particle may be coincident with this surface. Artifically move the
-      ! particle forward a tiny bit.
-      xyz(:) = xyz + TINY_BIT * uvw
-      s = this%sense(xyz, uvw)
-    elseif (f > 0) then
-      s = .true.
+      ! Particle may be coincident with this surface. To determine the sense, we
+      ! look at the direction of the particle relative to the surface normal (by
+      ! default in the positive direction) via their dot product.
+      s = (dot_product(uvw, this%normal(xyz)) > ZERO)
     else
-      s = .false.
+      s = (f > ZERO)
     end if
-
   end function sense
 
   subroutine reflect(this, xyz, uvw)

tests/test_plot_basis/results_true.dat

@@ -1 +1 @@
-6a2400a95ea5baee432dd04a85252818d682004dc80b225b99c4bbc7cfd20d8523bd8e5df8a272df288cb2f300e5eee55f7ffa7b4d28ea64f92b78839ab0e86b
+368e0135c136d5c8a2dabb4c8085279dc7ac0bd81b2ec905bdf11ecb5fe99803868631cdff0b3ddec941323bcc661747d4c16edfd4f8d38582155bd6fd7e82e8

tests/test_plot_mask/results_true.dat

@@ -1 +1 @@
-d888a734d6e69c5ca92457dc865c4defd8a1d8a1d4a01e0c8e7c528dde25d96df1283da37889dac0857329051958be2f2b491f2fd2bc22c8424a60ec6cda04bd
+a7cb65bf40c84c0540d45ff292c398f9ae51b3d9396e88b9b4e5cdf05e8730f409bddb53aec6d396058194c6293c5bd3ef39efd0b0f30f2423f696193c85176c