anchor.diff:  
       type(photon), intent(inout) :: p
   
  -    real(dp) :: r_sq,w_sq,theta,phi, dphi,phiv, thetav
  -    logical :: radial
  +    real(dp) :: r_sq,w_sq
  +    real(dp) :: theta, theta_v
  +    real(dp) :: phi, phi_v, dphi
       integer,parameter :: eps = 3
   
  -    ! Initialize values
  +    ! Initialize values. We start off by assuming the photon is not on any
  +    ! walls.
       p%on_wall = .false.
       p%on_wall_id = no_wall
   
  +    ! Pre-compute useful quantities
       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
  +    ! Find whether the photon is on a radial wall.
  +
  +    if((p%r .dot. p%v) >= 0._dp) then  ! photon is moving outwards
  +       if(equal_nulp(r_sq, geo%wr2(p%icell%i1), eps)) then ! photon is on inner wall
             p%on_wall_id%w1 = -1
  -       else if(equal_nulp(r_sq, geo%wr2(p%icell%i1 + 1), eps)) then
  +       else if(equal_nulp(r_sq, geo%wr2(p%icell%i1 + 1), eps)) then ! photon is on outer wall
             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
  +    else  ! photon is moving inwards
  +       if(equal_nulp(r_sq, geo%wr2(p%icell%i1), eps)) then ! photon is on inner wall
             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
  +       else if(equal_nulp(r_sq, geo%wr2(p%icell%i1 + 1), eps)) then ! photon is on outer wall
             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
  +    ! Find whether the photon is on a theta wall.
  +
  +    ! We first have to check whether we are at the origin, because in that
  +    ! case we cannot rely on the position of the photon to check which wall it
  +    ! is on.
  +    if(r_sq == 0._dp) then  ! photon is at origin
  +
  +       ! We need to ensure that the photon is not moving straight up or
  +       ! straight down. If it is, then we do not do anything since we do not
  +       ! treat the theta=0 and theta=pi walls as real walls. If the magnitude
  +       ! of the vertical velocity is not 1, then the photon is not moving
  +       ! vertically (since v is normalized).
          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
  +
  +          ! Find direction of motion of the photon.
  +          theta_v = atan2(sqrt(p%v%x*p%v%x + p%v%y*p%v%y), p%v%z)
  +
  +          ! Compare to the position of the lower and upper walls, and set
  +          ! the wall ID accordingly. We don't need to worry about changing the
  +          ! cell of the photon, since it is moving along the wall, so it
  +          ! doesn't matter whether it is set to one cell or the other.
  +          if(equal_nulp(theta_v, geo%w2(p%icell%i2), eps)) then
                p%on_wall_id%w2 = -1
  -          else if(equal_nulp(thetav, geo%wtant(p%icell%i2 + 1), eps)) then
  +          else if(equal_nulp(theta_v, geo%w2(p%icell%i2 + 1), eps)) then
                p%on_wall_id%w2 = +1
             end if
  -       end if
  +
  +       end if
  +
  +       ! If the photon is not at the origin, then we check whether it is on the
  +       ! lower wall, but we only do that if the lower wall is not theta=0.
       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%icell%i2 == geo%midplane) then  ! lower wall is grid midplane
  +
  +          ! The lower wall is the midplane of the grid, so we can just use
  +          ! v_x to determine the direction of motion relative to the wall.
             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
  +
  +       else  ! lower wall is not grid midplane
  +
  +          ! We need to do some vector arithmetic and dot products to
  +          ! determine whether the photon is moving towards the lower or upper
  +          ! cells.
             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
  +
  +       end if
  +
  +       ! If the photon is not at the origin, or on the lower wall, then we
  +       ! check whether it is on the upper wall, but we only do that if the lower
  +       ! wall is not theta=pi.
       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%icell%i2 + 1 == geo%midplane) then  ! upper wall is grid midplane
  +
  +          ! The lower wall is the midplane of the grid, so we can just use
  +          ! v_x to determine the direction of motion relative to the wall.
             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
  +
  +       else  ! upper wall is not grid midplane
  +
  +          ! We need to do some vector arithmetic and dot products to
  +          ! determine whether the photon is moving towards the lower or upper
  +          ! cells.
             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
  
positive.code:              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(dphi < -pi) dphi = dphi + twopi
         if(dphi > 0._dp) then
            p%on_wall_id%w3 = -1
         else
            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)
         if(dphi < -pi) dphi = dphi + twopi
         if(dphi > 0._dp) then
            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
            p%on_wall_id%w3 = +1
         end if
      end if
  
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
  