From dce591b037e8bd48dc2df907fad729595143d1a5 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 23 Sep 2015 12:56:13 +0700
Subject: [PATCH] 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                 | 22 ++++++++++------------
 tests/test_plot_basis/results_true.dat |  2 +-
 tests/test_plot_mask/results_true.dat  |  2 +-
 3 files changed, 12 insertions(+), 14 deletions(-)

diff --git a/src/surface_header.F90 b/src/surface_header.F90
index 402e028ea1..8941fed725 100644
--- a/src/surface_header.F90
+++ b/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)
diff --git a/tests/test_plot_basis/results_true.dat b/tests/test_plot_basis/results_true.dat
index ae588fd4f3..b1d5dc8534 100644
--- a/tests/test_plot_basis/results_true.dat
+++ b/tests/test_plot_basis/results_true.dat
@@ -1 +1 @@
-6a2400a95ea5baee432dd04a85252818d682004dc80b225b99c4bbc7cfd20d8523bd8e5df8a272df288cb2f300e5eee55f7ffa7b4d28ea64f92b78839ab0e86b
\ No newline at end of file
+368e0135c136d5c8a2dabb4c8085279dc7ac0bd81b2ec905bdf11ecb5fe99803868631cdff0b3ddec941323bcc661747d4c16edfd4f8d38582155bd6fd7e82e8
\ No newline at end of file
diff --git a/tests/test_plot_mask/results_true.dat b/tests/test_plot_mask/results_true.dat
index 1076324545..a1e323203d 100644
--- a/tests/test_plot_mask/results_true.dat
+++ b/tests/test_plot_mask/results_true.dat
@@ -1 +1 @@
-d888a734d6e69c5ca92457dc865c4defd8a1d8a1d4a01e0c8e7c528dde25d96df1283da37889dac0857329051958be2f2b491f2fd2bc22c8424a60ec6cda04bd
\ No newline at end of file
+a7cb65bf40c84c0540d45ff292c398f9ae51b3d9396e88b9b4e5cdf05e8730f409bddb53aec6d396058194c6293c5bd3ef39efd0b0f30f2423f696193c85176c
\ No newline at end of file