diff --git a/src/geometry.f90 b/src/geometry.f90
index ee81c9d9b1..bdcc6c7415 100644
--- a/src/geometry.f90
+++ b/src/geometry.f90
@@ -101,11 +101,15 @@ contains
     integer :: y                    ! y-index for lattice
     integer :: n                    ! number of cells to search
     integer :: index_cell           ! index in cells array
+    real(8) :: xyz(3)               ! temporary location
     logical :: use_search_cells     ! use cells provided as argument
     type(Cell),     pointer :: c    ! pointer to cell
     type(Lattice),  pointer :: lat  ! pointer to lattice
     type(Universe), pointer :: univ ! universe to search in
 
+    ! Remove coordinates for any lower levels
+    call deallocate_coord(p % coord % next)
+
     ! set size of list to search
     if (present(search_cells)) then
        use_search_cells = .true.
@@ -128,8 +132,6 @@ contains
        ! get pointer to cell
        c => cells(index_cell)
 
-       ! print *, 'searching cell ', c % id, cell_contains(c,p)
-       
        if (cell_contains(c, p)) then
           ! Set cell on this level
           p % coord % cell = index_cell
@@ -164,8 +166,9 @@ contains
              lat => lattices(c % fill)
 
              ! determine universe based on lattice position
-             x = ceiling((p % coord % xyz(1) - lat % x0)/lat % width_x)
-             y = ceiling((p % coord % xyz(2) - lat % y0)/lat % width_y)
+             xyz = p % coord % xyz + TINY_BIT * p % coord % uvw
+             x = ceiling((xyz(1) - lat % x0)/lat % width_x)
+             y = ceiling((xyz(2) - lat % y0)/lat % width_y)
 
              ! Create new level of coordinates
              p % in_lower_universe = .true.
@@ -455,7 +458,10 @@ contains
        dist = INFINITY
 
        ! left and right sides
-       if (u > 0) then
+       if (u == ZERO) then
+          d_left = INFINITY
+          d_right = INFINITY
+       elseif (u > 0) then
           d_left = INFINITY
           d_right = (x0 - x)/u
        else
@@ -464,7 +470,10 @@ contains
        end if
 
        ! top and bottom sides
-       if (v > 0) then
+       if (v == ZERO) then
+          d_bottom = INFINITY
+          d_top = INFINITY
+       elseif (v > 0) then
           d_bottom = INFINITY
           d_top = (y0 - y)/v
        else
@@ -558,6 +567,7 @@ contains
 
     ! inialize distance to infinity (huge)
     dist = INFINITY
+    lattice_crossed = .false.
 
     ! Get pointer to top-level coordinates
     coord => p % coord0
@@ -846,6 +856,7 @@ contains
           if (d < dist) then
              dist = d
              surface_crossed = -cl % surfaces(i)
+             lattice_crossed = .false.
              final_coord => coord
           end if
 
@@ -854,19 +865,22 @@ contains
        ! =======================================================================
        ! FIND MINIMUM DISTANCE TO LATTICE SURFACES
 
-       lattice_crossed = .false.
-
-       if (p % coord % lattice /= NONE) then
-          lat => lattices(p % coord % lattice)
+       if (coord % lattice /= NONE) then
+          lat => lattices(coord % lattice)
           if (lat % type == LATTICE_RECT) then
-             x = p % coord % xyz(1)
-             y = p % coord % xyz(2)
-             z = p % coord % xyz(3)
+             ! copy local coordinates
+             x = coord % xyz(1)
+             y = coord % xyz(2)
+             z = coord % xyz(3)
+
+             ! determine oncoming edge
              x0 = lat % width_x * 0.5_8
              y0 = lat % width_y * 0.5_8
 
              ! left and right sides
-             if (u > 0) then
+             if (u == ZERO) then
+                d = INFINITY
+             elseif (u > 0) then
                 d = (x0 - x)/u
              else
                 d = -(x + x0)/u
@@ -888,7 +902,9 @@ contains
              end if
 
              ! top and bottom sides
-             if (v > 0) then
+             if (v == ZERO) then
+                d = INFINITY
+             elseif (v > 0) then
                 d = (y0 - y)/v
              else
                 d = -(y + y0)/v
@@ -911,9 +927,6 @@ contains
 
     end do LEVEL_LOOP
 
-    ! Remove coordinates for any lower levels
-    call deallocate_coord(final_coord % next)
-
     ! Move particle to appropriate coordinate level
     p % coord => final_coord
 
diff --git a/src/particle_header.f90 b/src/particle_header.f90
index 64fbc2e6d5..fcd6068039 100644
--- a/src/particle_header.f90
+++ b/src/particle_header.f90
@@ -39,9 +39,9 @@ module particle_header
      integer    :: type          ! Particle type (n, p, e, etc)
 
      ! Particle coordinates
-     logical :: in_lower_universe        ! is particle in lower universe?
-     type(LocalCoord), pointer :: coord0 ! coordinates on universe 0
-     type(LocalCoord), pointer :: coord  ! coordinates on lowest universe
+     logical :: in_lower_universe                  ! is particle in lower universe?
+     type(LocalCoord), pointer :: coord0 => null() ! coordinates on universe 0
+     type(LocalCoord), pointer :: coord  => null() ! coordinates on lowest universe
 
      ! Other physical data
      real(8)    :: wgt           ! particle weight
diff --git a/src/physics.f90 b/src/physics.f90
index 2b911329d5..6c5d4ec902 100644
--- a/src/physics.f90
+++ b/src/physics.f90
@@ -10,7 +10,7 @@ module physics
   use geometry_header,      only: Universe, BASE_UNIVERSE
   use global
   use interpolation,        only: interpolate_tab1
-  use output,               only: write_message, print_particle
+  use output,               only: write_message
   use particle_header,      only: Particle, LocalCoord
   use random_lcg,           only: prn
   use search,               only: binary_search
@@ -104,7 +104,7 @@ contains
           call collision(p)
 
           ! Save coordinates at collision for tallying purposes
-          p % last_xyz = p % coord % xyz
+          p % last_xyz = p % coord0 % xyz
 
           ! Set all uvws to base level -- right now, after a collision, only the
           ! base level uvws are changed
diff --git a/src/plot.f90 b/src/plot.f90
index 79a08cd70d..785b051353 100644
--- a/src/plot.f90
+++ b/src/plot.f90
@@ -13,7 +13,9 @@ module plot
 contains
 
 !===============================================================================
-! RUN_PLOT
+! RUN_PLOT generates a binary stream file containing a list of surface/lattice
+! crossings and what cell was traveled through. A Python script can then be used
+! to generate a plot based on the recorded crossings and cells
 !===============================================================================
 
   subroutine run_plot()
@@ -131,6 +133,8 @@ contains
        ! MOVE PARTICLE ACROSS HORIZONTAL TRACK
 
        do while (p % alive)
+          ! save particle's current cell
+          last_cell = p % coord % cell
 
           ! Calculate distance to next boundary
           call distance_to_boundary(p, distance, surface_crossed, lattice_crossed)
@@ -152,14 +156,13 @@ contains
              if (distance == INFINITY) p % coord % cell = 0
 
              ! Write ending coordinates to file
-             write(UNIT=UNIT_PLOT) p % coord0 % xyz, p % coord % cell
+             write(UNIT=UNIT_PLOT) p % coord0 % xyz, last_cell
              cycle
           end if
 
           ! Write boundary crossing coordinates to file
-          write(UNIT=UNIT_PLOT) p % coord0 % xyz, p % coord % cell
+          write(UNIT=UNIT_PLOT) p % coord0 % xyz, last_cell
 
-          last_cell = p % coord % cell
           p % coord % cell = 0
           if (lattice_crossed) then
              p % surface = NONE
@@ -170,7 +173,7 @@ contains
 
              ! Since boundary conditions are disabled in plotting mode, we need
              ! to manually add the last segment
-             if (surfaces(surface_crossed) % bc == BC_VACUUM) then
+             if (surfaces(abs(surface_crossed)) % bc == BC_VACUUM) then
                 p % coord0 % xyz(1) = last_x_coord
                 write(UNIT=UNIT_PLOT) p % coord0 % xyz, 0
                 exit