Tracking

Using the awesome new particle restart capability I finally got around to figuring out why a tracking error occurs in my two-region lattice surrounded by reflective boundary conditions.  The two region lattice geometry is shown below.

![lattice](https://f.cloud.github.com/assets/1017998/363595/14d0e63c-a208-11e2-83be-8ec11dab4ffa.png)

A tracking error occurs when a neutron reflects off a surface that is close to a corner. An example of this is shown below (not to scale) where a neutron starts in a lattice cell and traverses to a reflective surface that is close to a corner.

![trackingerror](https://f.cloud.github.com/assets/1017998/363579/9fdc6298-a207-11e2-9a4b-229899ef6e27.png)

After the reflection math is performed, the following block of code is called from the cross_surface routine in geometry.F90:
```fortran
      if (associated(p % coord0 % next)) then
        call deallocate_coord(p % coord0 % next)
        call find_cell(found)
        if (.not. found) then
          call write_particle_restart()
          message = "Couldn't find particle after reflecting from surface."
          call fatal_error()
        end if
      end if
```
OpenMC will now try to locate the neutron again after the cell reflection by calling find_cell routine in geometry.F90.  In executing find_cell recursively, OpenMC will eventually determine that the neutron is in a lattice. OpenMC always moves the neutron a ```TINY_BIT``` in its direction so that it is not exactly on a surface and should have no problem computing the correct lattice indices ```(i_x, i_y, i_z)```.  This is operation is shown in the following code in find_cell:
```fortran
          ! now move lattice tiny bit along its unit vector
          xyz = p % coord % xyz + TINY_BIT * p % coord % uvw
          i_x = ceiling((xyz(1) - lat % lower_left(1))/lat % width(1))
          i_y = ceiling((xyz(2) - lat % lower_left(2))/lat % width(2))
          n_x = lat % dimension(1)
          n_y = lat % dimension(2)
          if (lat % n_dimension == 3) then
            i_z = ceiling((xyz(3) - lat % lower_left(3))/lat % width(3))
            n_z = lat % dimension(3)
          else
            i_z = 1
            n_z = 1
          end if
```
In the case that I am showing above, the neutron is so close to the corner that when it is moved that ```TINY_BIT```, it effectively is now outside the lattice and the index ```i_x``` will be greater than the max x index ```n_x```. When this occurs the code will not find the particle since it is outside of the lattice bounds.

To fix this tracking error I first compute the original indices of the neutron before moving it a ```TINY_BIT```.  This assumes that since OpenMC put that particle *exactly* on the surface, it will return indices that are inside the lattice. This is done in the following code:
```fortran
          ! determine lattice index based on position
          xyz = p % coord % xyz
          i_x0 = ceiling((xyz(1) - lat % lower_left(1))/lat % width(1))
          i_y0 = ceiling((xyz(2) - lat % lower_left(2))/lat % width(2))
          if (lat % n_dimension == 3) &
             i_z0 = ceiling((xyz(3) - lat % lower_left(3))/lat % width(3))
```

Then, I allow the particle to be moved the ```TINY_BIT``` along its unit vector and calculate the new indices. Then I go direction by direction and check and see if by moving this particle, we disturbed its lattice location by moving it out.  This is checked with the following code:
```fortran
          ! Check if neutron out of lattice bounds by moving TINY_BIT.
          ! This can happen if neutron is within tolerance of TINY_BIT
          ! to corner, this action will put neutron out of lattice.
          if (i_x0 >= 1 .and. i_x < 1) i_x = 1       ! -x direction
          if (i_y0 >= 1 .and. i_y < 1) i_y = 1       ! -y direction
          if (i_z0 >= 1 .and. i_z < 1) i_z = 1       ! -z direction
          if (i_x0 <= n_x .and. i_x > n_x) i_x = n_x ! +x direction
          if (i_y0 <= n_x .and. i_y > n_x) i_y = n_y ! +y direction
          if (i_z0 <= n_x .and. i_z > n_x) i_z = n_z ! +z direction
```

If any of these checks are satisfied, the lattice indices are either set to the max or the min lattice bounds.