anchor.diff:         &                                          6, 13, 14, 17, 18,  0,  0, 23, 24, 25, 26,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0/)
   
   contains
  +
  +  logical pure function equal_nulp(x, y, n)
  +    implicit none
  +    real(8),intent(in) :: x, y
  +    integer,intent(in) :: n
  +    if(x == y) then
  +       equal_nulp = .true.
  +    else
  +       equal_nulp = abs(x - y) <= n * spacing(max(x, y))
  +    end if
  +  end function equal_nulp
   
     real(dp) function cell_width(cell, idir)
       implicit none
  
positive.code: 
      ! Find whether the photon is on a radial wall
      if((p%r%x * p%v%x + p%r%y * p%v%y) >= 0._dp) then
         if(r2 == geo%wr2(p%icell%i1)) then  ! TODO: shouldn't use equal, not precise enough
            p%on_wall_id = combine_wall(p%on_wall_id, 1)
         else if(r2 == geo%wr2(p%icell%i1 + 1)) then
            p%on_wall_id = combine_wall(p%on_wall_id, 1)
            p%icell%i1 = p%icell%i1 + 1
         end if
      else
         if(r2 == geo%wr2(p%icell%i1)) then
            p%on_wall_id = combine_wall(p%on_wall_id, 2)
            p%icell%i1 = p%icell%i1 - 1
         else if(r2 == geo%wr2(p%icell%i1 + 1)) then
            p%on_wall_id = combine_wall(p%on_wall_id, 2)
         end if
      end if
  
      ! Find whether the photon is on a vertical wall
      if(p%v%z > 0._dp) then
         if(p%r%z == geo%w2(p%icell%i2)) then
            p%on_wall_id = combine_wall(p%on_wall_id, 3)
         else if(p%r%z == geo%w2(p%icell%i2 + 1)) then
            p%on_wall_id = combine_wall(p%on_wall_id, 3)
            p%icell%i2 = p%icell%i2 + 1
         end if
      else if(p%v%z < 0._dp) then
         if(p%r%z == geo%w2(p%icell%i2)) then
            p%on_wall_id = combine_wall(p%on_wall_id, 4)
            p%icell%i2 = p%icell%i2 - 1
         else if(p%r%z == geo%w2(p%icell%i2 + 1)) then
            p%on_wall_id = combine_wall(p%on_wall_id, 4)
         end if
      end if
  
      ! Find whether the photon is on an azimuthal wall
      phiv = atan2(p%v%y, p%v%x)
      if(phi == geo%w3(p%icell%i3)) then
         dphi = phiv - geo%w3(p%icell%i3)
         if(dphi < -pi) dphi = dphi + twopi
         if(dphi > 0._dp) then
  
negative.code:           p%icell%i3 = p%icell%i3 - 1
            if(p%icell%i3==0) p%icell%i3 = geo%n3
         end if
      else if(equal_nulp(phi, geo%w3(p%icell%i3 + 1), 10)) then
         dphi = phiv - geo%w3(p%icell%i3 + 1)
         if(dphi < -pi) dphi = dphi + twopi
         if(dphi > 0._dp) then
  