anchor.diff:               p%on_wall_id%w2 = +1
             end if
          end if
  +
       end if
   
       ! Find whether the photon is on an azimuthal wall
  -    phiv = atan2(p%v%y, p%v%x)
  -    if(equal_nulp(phi, geo%w3(p%icell%i3), eps)) then
  -       dphi = phiv - geo%w3(p%icell%i3)
  +
  +    if(equal_nulp(phi, geo%w3(p%icell%i3), eps)) then  ! photon is on inner wall
  +
  +       ! Find the angle between the direction of motion and the wall.
  +       phi_v = atan2(p%v%y, p%v%x)
  +       dphi = phi_v - geo%w3(p%icell%i3)
          if(dphi < -pi) dphi = dphi + twopi
  -       if(dphi > 0._dp) then
  +
  +       if(dphi > 0._dp) then  ! photon is moving towards upper cells
             p%on_wall_id%w3 = -1
  -       else
  +       else  ! photon is moving towards lower cells
             p%on_wall_id%w3 = +1
             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), eps)) then
  -       dphi = phiv - geo%w3(p%icell%i3 + 1)
  +
  +    else if(equal_nulp(phi, geo%w3(p%icell%i3 + 1), eps)) then  ! photon is on outer wall
  +
  +       ! Find the angle between the direction of motion and the wall.
  +       phi_v = atan2(p%v%y, p%v%x)
  +       dphi = phi_v - geo%w3(p%icell%i3 + 1)
          if(dphi < -pi) dphi = dphi + twopi
  -       if(dphi > 0._dp) then
  +
  +       if(dphi > 0._dp) then  ! photon is moving towards upper cells
             p%on_wall_id%w3 = -1
             p%icell%i3 = p%icell%i3 + 1
             if(p%icell%i3==geo%n3+1) p%icell%i3 = 1
  -       else
  +       else  ! photon is moving towards lower cells
             p%on_wall_id%w3 = +1
          end if
       end if
  
positive.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) < 1.e-3_dp
         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) < 1.e-3_dp
         else
            in_correct_cell = in_correct_cell .and. icell_actual%i3 == p%icell%i3
         end if
  
negative.code:        phi = atan2(p%r%y,p%r%x)
         if(phi < 0._dp) phi = phi + twopi
      end if
  
      ! Find the cells that the photon falls in
      i1 = locate(geo%wr2,r_sq)
      i2 = locate(geo%w2,theta)
      i3 = locate(geo%w3,phi)
  
      ! We now check whether the photon falls outside the grid in any of the
      ! dimensions.
  
      if(i1<1.or.i1>geo%n1) then
         call warn("find_cell","photon not in cell (in r direction)")
  