From b0fe887aebbe5c79f834f6bf2c61f70ef3c7c649 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 1 Feb 2013 22:06:50 -0500 Subject: [PATCH] Fix distance to surface for cones, again. Closes #143. --- src/geometry.F90 | 118 +++++++++++++----------- tests/test_reflective_cone/geometry.xml | 1 - 2 files changed, 64 insertions(+), 55 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 506a3a9d85..cc3ca0e5b1 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -715,6 +715,10 @@ contains ! get pointer to surface surf => surfaces(index_surf) + ! TODO: Can probably combines a lot of the cases to reduce repetition + ! since the algorithm is the same for (x-plane, y-plane, z-plane), + ! (x-cylinder, y-cylinder, z-cylinder), etc. + select case (surf % type) case (SURF_PX) if (on_surface .or. u == ZERO) then @@ -964,39 +968,41 @@ contains c = y*y + z*z - r*x*x quad = k*k - a*c - if (quad < ZERO .or. a == ZERO) then + if (quad < ZERO) then ! no intersection with cone d = INFINITY elseif (on_surface) then - ! particle is on the cone, thus one distance is - ! positive/negative and the other is zero. The sign of k - ! determines if we are facing in or out + ! particle is on the cone, thus one distance is positive/negative + ! and the other is zero. The sign of k determines which distance is + ! zero and which is not. if (k >= ZERO) then - d = INFINITY + d = (-k - sqrt(quad))/a else d = (-k + sqrt(quad))/a end if - elseif (c < ZERO) then - ! particle is inside the cone, thus one distance must be - ! negative and one must be positive. The positive distance will - ! be the one with negative sign on sqrt(quad) - - d = (-k + sqrt(quad))/a - else - ! particle is outside the cone, thus both distances are either - ! positive or negative. If positive, the smaller distance is the - ! one with positive sign on sqrt(quad) - - d = (-k - sqrt(quad))/a - if (d <= ZERO) d = INFINITY + ! calculate both solutions to the quadratic + quad = sqrt(quad) + d = (-k - quad)/a + b = (-k + quad)/a + ! determine the smallest positive solution + if (d < ZERO) then + if (b > ZERO) then + d = b + end if + else + if (b > ZERO) d = min(d, b) + end if end if + ! If the distance was negative, set boundary distance to infinity + if (d <= ZERO) d = INFINITY + case (SURF_CONE_Y) x0 = surf % coeffs(1) y0 = surf % coeffs(2) @@ -1011,39 +1017,41 @@ contains c = x*x + z*z - r*y*y quad = k*k - a*c - if (quad < ZERO .or. a == ZERO) then + if (quad < ZERO) then ! no intersection with cone d = INFINITY elseif (on_surface) then - ! particle is on the cone, thus one distance is - ! positive/negative and the other is zero. The sign of k - ! determines if we are facing in or out + ! particle is on the cone, thus one distance is positive/negative + ! and the other is zero. The sign of k determines which distance is + ! zero and which is not. if (k >= ZERO) then - d = INFINITY + d = (-k - sqrt(quad))/a else d = (-k + sqrt(quad))/a end if - elseif (c < ZERO) then - ! particle is inside the cone, thus one distance must be - ! negative and one must be positive. The positive distance will - ! be the one with negative sign on sqrt(quad) - - d = (-k + sqrt(quad))/a - else - ! particle is outside the cone, thus both distances are either - ! positive or negative. If positive, the smaller distance is the - ! one with positive sign on sqrt(quad) - - d = (-k - sqrt(quad))/a - if (d <= ZERO) d = INFINITY + ! calculate both solutions to the quadratic + quad = sqrt(quad) + d = (-k - quad)/a + b = (-k + quad)/a + ! determine the smallest positive solution + if (d < ZERO) then + if (b > ZERO) then + d = b + end if + else + if (b > ZERO) d = min(d, b) + end if end if + ! If the distance was negative, set boundary distance to infinity + if (d <= ZERO) d = INFINITY + case (SURF_CONE_Z) x0 = surf % coeffs(1) y0 = surf % coeffs(2) @@ -1058,39 +1066,41 @@ contains c = x*x + y*y - r*z*z quad = k*k - a*c - if (quad < ZERO .or. a == ZERO) then + if (quad < ZERO) then ! no intersection with cone d = INFINITY elseif (on_surface) then - ! particle is on the cone, thus one distance is - ! positive/negative and the other is zero. The sign of k - ! determines if we are facing in or out + ! particle is on the cone, thus one distance is positive/negative + ! and the other is zero. The sign of k determines which distance is + ! zero and which is not. if (k >= ZERO) then - d = INFINITY + d = (-k - sqrt(quad))/a else d = (-k + sqrt(quad))/a end if - elseif (c < ZERO) then - ! particle is inside the cone, thus one distance must be - ! negative and one must be positive. The positive distance will - ! be the one with negative sign on sqrt(quad) - - d = (-k + sqrt(quad))/a - else - ! particle is outside the cone, thus both distances are either - ! positive or negative. If positive, the smaller distance is the - ! one with positive sign on sqrt(quad) - - d = (-k - sqrt(quad))/a - if (d <= ZERO) d = INFINITY + ! calculate both solutions to the quadratic + quad = sqrt(quad) + d = (-k - quad)/a + b = (-k + quad)/a + ! determine the smallest positive solution + if (d < ZERO) then + if (b > ZERO) then + d = b + end if + else + if (b > ZERO) d = min(d, b) + end if end if + ! If the distance was negative, set boundary distance to infinity + if (d <= ZERO) d = INFINITY + case (SURF_GQ) message = "Surface distance not yet implement for general quadratic." call fatal_error() diff --git a/tests/test_reflective_cone/geometry.xml b/tests/test_reflective_cone/geometry.xml index 979f18f7a7..5453e16155 100644 --- a/tests/test_reflective_cone/geometry.xml +++ b/tests/test_reflective_cone/geometry.xml @@ -1,7 +1,6 @@ -