From d335606220db98ce68e11b59684d8425654c6104 Mon Sep 17 00:00:00 2001
From: Paul Romano <romano7@gmail.com>
Date: Sun, 4 Dec 2011 18:28:36 -0500
Subject: [PATCH] Fixed geometry plotting when bounding box is bigger than
 geometry.

---
 src/geometry.f90 |   5 ++-
 src/plot.f90     | 108 ++++++++++++++++++++++++-----------------------
 2 files changed, 58 insertions(+), 55 deletions(-)

diff --git a/src/geometry.f90 b/src/geometry.f90
index bdcc6c7415..da2a3febab 100644
--- a/src/geometry.f90
+++ b/src/geometry.f90
@@ -405,7 +405,7 @@ contains
     call find_cell(p, found)
 
     ! Couldn't find next cell anywhere!
-    if (.not. found) then
+    if ((.not. found) .and. (.not. plotting)) then
        message = "After particle crossed surface " // trim(int_to_str(p%surface)) &
             // ", it could not be located in any cell and it did not leak."
        call fatal_error()
@@ -568,6 +568,7 @@ contains
     ! inialize distance to infinity (huge)
     dist = INFINITY
     lattice_crossed = .false.
+    nullify(final_coord)
 
     ! Get pointer to top-level coordinates
     coord => p % coord0
@@ -928,7 +929,7 @@ contains
     end do LEVEL_LOOP
 
     ! Move particle to appropriate coordinate level
-    p % coord => final_coord
+    if (associated(final_coord)) p % coord => final_coord
 
   end subroutine distance_to_boundary
 
diff --git a/src/plot.f90 b/src/plot.f90
index 785b051353..8eb0a10823 100644
--- a/src/plot.f90
+++ b/src/plot.f90
@@ -6,7 +6,8 @@ module plot
                              cross_lattice, cell_contains
   use geometry_header, only: Universe, BASE_UNIVERSE
   use global
-  use particle_header, only: Particle, initialize_particle, LocalCoord
+  use particle_header, only: Particle, initialize_particle, LocalCoord,      &
+                             deallocate_coord
 
   implicit none
 
@@ -76,58 +77,59 @@ contains
        ! =======================================================================
        ! MOVE PARTICLE FORWARD TO NEXT CELL
 
-!!$       if (.not. found_cell) then
-!!$          univ => universes(BASE_UNIVERSE)
-!!$          do i = 1, univ % n_cells
-!!$             p % xyz = coord
-!!$             p % xyz_local = coord
-!!$             p % cell = univ % cells(i)
-!!$
-!!$             distance = INFINITY
-!!$             ! call distance_to_boundary(p, d, surf, in_lattice)
-!!$             if (d < distance) then
-!!$                ! Move particle forward to next surface
-!!$                p % xyz = p % xyz + d * p % uvw
-!!$
-!!$                ! Check to make sure particle is actually going into this cell
-!!$                ! by moving it slightly forward and seeing if the cell contains
-!!$                ! that coordinate
-!!$
-!!$                p % xyz = p % xyz + 1e-4 * p % uvw
-!!$                p % xyz_local = p % xyz
-!!$
-!!$                c => cells(p % cell)
-!!$                if (.not. cell_contains(c, p)) cycle
-!!$
-!!$                ! Reset coordinate to surface crossing
-!!$                p % xyz = p % xyz - 1e-4 * p % uvw
-!!$                p % xyz_local = p % xyz
-!!$
-!!$                ! Set new distance and retain pointer to this cell
-!!$                distance = d
-!!$                last_cell = p % cell
-!!$             end if
-!!$          end do
-!!$
-!!$          ! No cell was found on this horizontal ray
-!!$          if (distance == INFINITY) then
-!!$             p % xyz(1) = last_x_coord
-!!$             p % cell = 0
-!!$             write(UNIT_PLOT) p % xyz, p % cell
-!!$
-!!$             ! Move to next horizontal ray
-!!$             xyz(2) = xyz(2) - pixel
-!!$             cycle
-!!$          end if
-!!$
-!!$          ! Write coordinate where next cell begins
-!!$          write(UNIT=UNIT_PLOT) p % xyz, 0
-!!$
-!!$          ! Process surface crossing for next cell
-!!$          p % cell = 0
-!!$          p % surface = -surf
-!!$          call cross_surface(p, last_cell)
-!!$       end if
+       if (.not. found_cell) then
+          ! Clear any coordinates beyond first level
+          call deallocate_coord(p % coord0 % next)
+          p % coord => p % coord0
+
+          univ => universes(BASE_UNIVERSE)
+          do i = 1, univ % n_cells
+             p % coord0 % xyz = xyz
+             p % coord0 % cell = univ % cells(i)
+
+             distance = INFINITY
+             call distance_to_boundary(p, d, surface_crossed, lattice_crossed)
+             if (d < distance) then
+                ! Move particle forward to next surface
+                ! Advance particle
+                p % coord0 % xyz = p % coord0 % xyz + d * p % coord0 % uvw
+
+                ! Check to make sure particle is actually going into this cell
+                ! by moving it slightly forward and seeing if the cell contains
+                ! that coordinate
+
+                p % coord0 % xyz = p % coord0 % xyz + 1e-4 * p % coord0 % uvw
+
+                c => cells(p % coord0 % cell)
+                if (.not. cell_contains(c, p)) cycle
+
+                ! Reset coordinate to surface crossing
+                p % coord0 % xyz = p % coord0 % xyz - 1e-4 * p % coord0 % uvw
+
+                ! Set new distance and retain pointer to this cell
+                distance = d
+                last_cell = p % coord0 % cell
+             end if
+          end do
+
+          ! No cell was found on this horizontal ray
+          if (distance == INFINITY) then
+             p % coord0 % xyz(1) = last_x_coord
+             write(UNIT_PLOT) p % coord0 % xyz, 0
+
+             ! Move to next horizontal ray
+             xyz(2) = xyz(2) - pixel
+             cycle
+          end if
+
+          ! Write coordinate where next cell begins
+          write(UNIT=UNIT_PLOT) p % coord0 % xyz, 0
+
+          ! Process surface crossing for next cell
+          p % coord0 % cell = NONE
+          p % surface = -surface_crossed
+          call cross_surface(p, last_cell)
+       end if
 
        ! =======================================================================
        ! MOVE PARTICLE ACROSS HORIZONTAL TRACK