anchor.diff:  
       ! 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
  +       if(equal_nulp(r2, geo%wr2(p%icell%i1), 10)) 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
  +       else if(equal_nulp(r2, geo%wr2(p%icell%i1 + 1), 10)) 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
  +       if(equal_nulp(r2, geo%wr2(p%icell%i1), 10)) 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
  +       else if(equal_nulp(r2, geo%wr2(p%icell%i1 + 1), 10)) 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
  +       if(equal_nulp(p%r%z, geo%w2(p%icell%i2), 10)) then
             p%on_wall_id = combine_wall(p%on_wall_id, 3)
  -       else if(p%r%z == geo%w2(p%icell%i2 + 1)) then
  +       else if(equal_nulp(p%r%z, geo%w2(p%icell%i2 + 1), 10)) 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
  +       if(equal_nulp(p%r%z, geo%w2(p%icell%i2), 10)) 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
  +       else if(equal_nulp(p%r%z, geo%w2(p%icell%i2 + 1), 10)) 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
  +    if(equal_nulp(phi, geo%w3(p%icell%i3), 10)) then
          dphi = phiv - geo%w3(p%icell%i3)
          if(dphi < -pi) dphi = dphi + twopi
          if(dphi > 0._dp) then
  
positive.code:        &                                          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
  
    real(dp) function cell_width(cell, idir)
      implicit none
  
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
  