anchor.diff:      type(photon),intent(in) :: p
       type(grid_cell) :: icell_actual
       real(dp) :: rad,theta,phi,frac,dphi
  +    real(dp),parameter :: threshold = 1.e-3_dp
   
       icell_actual = find_cell(p)
   
       if(p%on_wall) then
   
  +       in_correct_cell = .true.
  +
          rad = sqrt(p%r.dot.p%r)
  +
  +       ! If we are at the origin, then there isn't much point in checking,
  +       ! since this is difficult numerically, and the photon has not
  +       ! propagated anyway.
  +       if(rad == 0._dp) return
  +
          theta = atan2(sqrt(p%r%x*p%r%x+p%r%y*p%r%y),p%r%z)
          phi = atan2(p%r%y,p%r%x)
          if(phi < 0._dp) phi = phi + twopi
   
  -       in_correct_cell = .true.
  -
  -       if(rad == 0._dp) return
  -
          if(p%on_wall_id%w1 == -1) then
             if(geo%w1(p%icell%i1) .ne. rad) then
  -             frac = sqrt(rad / geo%w1(p%icell%i1)) - 1._dp
  -             in_correct_cell = in_correct_cell .and. abs(frac) < 1.e-3_dp
  +             frac = rad / geo%w1(p%icell%i1) - 1._dp
  +             in_correct_cell = in_correct_cell .and. abs(frac) < threshold
             end if
          else if(p%on_wall_id%w1 == +1) then
             if(geo%w1(p%icell%i1 + 1) .ne. rad) then
  -             frac = sqrt(rad / geo%w1(p%icell%i1+1)) - 1._dp
  -             in_correct_cell = in_correct_cell .and. abs(frac) < 1.e-3_dp
  +             frac = rad / geo%w1(p%icell%i1+1) - 1._dp
  +             in_correct_cell = in_correct_cell .and. abs(frac) < threshold
             end if
          else
             in_correct_cell = in_correct_cell .and. icell_actual%i1 == p%icell%i1
  
positive.code: 
      type(photon), intent(inout) :: p
  
      real(dp) :: r_sq,w_sq,theta,phi, dphi,phiv, thetav
      logical :: radial
      integer,parameter :: eps = 3
  
      ! Initialize values
      p%on_wall = .false.
      p%on_wall_id = no_wall
  
      r_sq = p%r.dot.p%r
      w_sq = p%r%x*p%r%x+p%r%y*p%r%y
      theta = atan2(sqrt(w_sq),p%r%z)
      phi = atan2(p%r%y,p%r%x)
      if(phi < 0._dp) phi = phi + twopi
  
      ! Find whether the photon is on a radial wall
      if((p%r .dot. p%v) >= 0._dp) then
         if(equal_nulp(r_sq, geo%wr2(p%icell%i1), eps)) then
            p%on_wall_id%w1 = -1
         else if(equal_nulp(r_sq, geo%wr2(p%icell%i1 + 1), eps)) then
            p%on_wall_id%w1 = -1
            p%icell%i1 = p%icell%i1 + 1
         end if
      else
         if(equal_nulp(r_sq, geo%wr2(p%icell%i1), eps)) then
            p%on_wall_id%w1 = +1
            p%icell%i1 = p%icell%i1 - 1
         else if(equal_nulp(r_sq, geo%wr2(p%icell%i1 + 1), eps)) then
            p%on_wall_id%w1 = +1
         end if
      end if
  
      ! Find whether the photon is on a theta wall
      if(r_sq == 0._dp) then
         if(abs(p%v%z) < 1._dp) then
            thetav = p%v%z / sqrt(p%v%x*p%v%x + p%v%y*p%v%y)
            if(equal_nulp(thetav, geo%wtant(p%icell%i2), eps)) then
               p%on_wall_id%w2 = -1
            else if(equal_nulp(thetav, geo%wtant(p%icell%i2 + 1), eps)) then
               p%on_wall_id%w2 = +1
            end if
         end if
      else if(p%icell%i2 > 1 .and. equal_nulp(theta, geo%w2(p%icell%i2), eps)) then
         if(p%icell%i2 == geo%midplane) then
            if(p%v%z > 0._dp) then
               p%on_wall_id%w2 = +1
               p%icell%i2 = p%icell%i2 - 1
            else
               p%on_wall_id%w2 = -1
            end if
         else
            if((sqrt(w_sq) * p%v%z * geo%wtant(p%icell%i2) -  (p%r%x * p%v%x + p%r%y * p%v%y) < 0._dp) .eqv. p%r%z > 0._dp) then
               p%on_wall_id%w2 = -1
            else
               p%on_wall_id%w2 = +1
               p%icell%i2 = p%icell%i2 - 1
            end if
         end if
      else if(p%icell%i2 + 1 < geo%n2 + 1 .and. equal_nulp(theta, geo%w2(p%icell%i2 + 1), eps)) then
         if(p%icell%i2 + 1 == geo%midplane) then
            if(p%v%z > 0._dp) then
               p%on_wall_id%w2 = +1
            else
               p%on_wall_id%w2 = -1
               p%icell%i2 = p%icell%i2 + 1
            end if
         else
            if((sqrt(w_sq) * p%v%z * geo%wtant(p%icell%i2 + 1)  -  (p%r%x * p%v%x + p%r%y * p%v%y) < 0._dp) .eqv. p%r%z > 0._dp) then
               p%on_wall_id%w2 = -1
               p%icell%i2 = p%icell%i2 + 1
  
negative.code:           if(dphi > pi) dphi = dphi - twopi
            if(dphi < -pi) dphi = dphi + twopi
            frac = dphi / (geo%w3(p%icell%i3+1) - geo%w3(p%icell%i3))
            in_correct_cell = in_correct_cell .and. abs(frac) < threshold
         else if(p%on_wall_id%w3 == +1) then
            dphi = phi - geo%w3(p%icell%i3+1)
            if(dphi > pi) dphi = dphi - twopi
            if(dphi < -pi) dphi = dphi + twopi
            frac = dphi / (geo%w3(p%icell%i3+1) - geo%w3(p%icell%i3))
            in_correct_cell = in_correct_cell .and. abs(frac) < threshold
         else
            in_correct_cell = in_correct_cell .and. icell_actual%i3 == p%icell%i3
         end if
  