From 417d36b9c7d00f074182510cc7edf02a1b554e8b Mon Sep 17 00:00:00 2001
From: Paul Romano <romano7@gmail.com>
Date: Mon, 17 Oct 2011 13:10:41 -0400
Subject: [PATCH] Fixed cross_surface to handle lattices when searching
 neighbor lists. Geometry plotter can now handle lattices.

---
 src/geometry.f90 | 57 ++++++++++++++++++++++++++++++++++++++++++++++--
 src/plot.f90     | 12 ++++++++--
 2 files changed, 65 insertions(+), 4 deletions(-)

diff --git a/src/geometry.f90 b/src/geometry.f90
index d48419d58b..6b2254cf8f 100644
--- a/src/geometry.f90
+++ b/src/geometry.f90
@@ -176,6 +176,8 @@ contains
 
     integer                 :: i          ! index of neighbors
     integer                 :: index_cell ! index in cells array
+    integer                 :: i_x        ! x index in lattice
+    integer                 :: i_y        ! y index in lattice
     real(8)                 :: x          ! x-x0 for sphere
     real(8)                 :: y          ! y-y0 for sphere
     real(8)                 :: z          ! z-z0 for sphere
@@ -190,8 +192,9 @@ contains
     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
-    type(Cell),     pointer :: c
+    type(Surface),  pointer :: surf => null()
+    type(Cell),     pointer :: c => null()
+    type(Lattice),  pointer :: lat => null()
     type(Universe), pointer :: lower_univ => null()
 
     surf => surfaces(abs(p % surface))
@@ -332,6 +335,31 @@ contains
                    msg = "Could not locate particle in universe: "
                    call fatal_error(msg)
                 end if
+             elseif (c % type == CELL_LATTICE) then
+                ! Set current lattice
+                lat => lattices(c % fill)
+                p % lattice = c % fill
+
+                ! determine universe based on lattice position
+                i_x = ceiling((p%xyz(1) - lat%x0)/lat%width_x)
+                i_y = ceiling((p%xyz(2) - lat%y0)/lat%width_y)
+                lower_univ => universes(lat % element(i_x,i_y))
+
+                ! adjust local position of particle
+                p%xyz_local(1) = p%xyz(1) - (lat%x0 + (i_x-0.5)*lat%width_x)
+                p%xyz_local(2) = p%xyz(2) - (lat%y0 + (i_y-0.5)*lat%width_y)
+                p%xyz_local(3) = p%xyz(3)
+             
+                ! set particle lattice indices
+                p % index_x = i_x
+                p % index_y = i_y
+
+                call find_cell(lower_univ, p, found)
+                if (.not. found) then
+                   msg = "Could not locate particle in lattice: " & 
+                        & // int_to_str(lat % uid)
+                   call fatal_error(msg)
+                end if
              else
                 ! set current pointers
                 p % cell = index_cell
@@ -354,6 +382,31 @@ contains
                    msg = "Could not locate particle in universe: "
                    call fatal_error(msg)
                 end if
+             elseif (c % type == CELL_LATTICE) then
+                ! Set current lattice
+                lat => lattices(c % fill)
+                p % lattice = c % fill
+
+                ! determine universe based on lattice position
+                i_x = ceiling((p%xyz(1) - lat%x0)/lat%width_x)
+                i_y = ceiling((p%xyz(2) - lat%y0)/lat%width_y)
+                lower_univ => universes(lat % element(i_x,i_y))
+
+                ! adjust local position of particle
+                p%xyz_local(1) = p%xyz(1) - (lat%x0 + (i_x-0.5)*lat%width_x)
+                p%xyz_local(2) = p%xyz(2) - (lat%y0 + (i_y-0.5)*lat%width_y)
+                p%xyz_local(3) = p%xyz(3)
+             
+                ! set particle lattice indices
+                p % index_x = i_x
+                p % index_y = i_y
+
+                call find_cell(lower_univ, p, found)
+                if (.not. found) then
+                   msg = "Could not locate particle in lattice: " & 
+                        & // int_to_str(lat % uid)
+                   call fatal_error(msg)
+                end if
              else
                 ! set current pointers
                 p % cell = index_cell
diff --git a/src/plot.f90 b/src/plot.f90
index 953e7cc326..d247f7536c 100644
--- a/src/plot.f90
+++ b/src/plot.f90
@@ -142,7 +142,7 @@ contains
 
           ! If next boundary crossing is out of range of the plot, only include
           ! the visible portion and move to next horizontal ray
-          if (p % xyz(1) > last_x_coord) then
+          if (p % xyz(1) >= last_x_coord) then
              p % alive = .false.
              p % xyz(1) = last_x_coord
 
@@ -165,8 +165,16 @@ contains
           else
              p % surface = surf
              call cross_surface(p, last_cell)
+
+             ! Since boundary conditions are disabled in plotting mode, we need
+             ! to manually add the last segment
+             if (surfaces(surf) % bc == BC_VACUUM) then
+                p % xyz(1) = last_x_coord
+                write(UNIT=UNIT_PLOT) p % xyz, 0
+                exit
+             end if
           end if
-       
+
        end do
 
        ! Move y-coordinate to next position