anchor.diff:      ! 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
  +
  +    if(r_sq == 0._dp) then
  +       theta = atan2(sqrt(p%v%x*p%v%x+p%v%y*p%v%y),p%v%z)
  +    else
  +       theta = atan2(sqrt(p%r%x*p%r%x+p%r%y*p%r%y),p%r%z)
  +    end if
  +
  +    if(w_sq == 0._dp) then
  +       phi = atan2(p%v%y,p%v%x)
  +       if(phi < 0._dp) phi = phi + twopi
  +    else
  +       phi = atan2(p%r%y,p%r%x)
  +       if(phi < 0._dp) phi = phi + twopi
  +    end if
   
       ! Find whether the photon is on a radial wall.
   
  
positive.code:        ! make sure the photon is still considered to be on the wall
         ! afterwards. Therefore, we calculate dphi to find the angle between
         ! the direction of motion and the wall the photon is on.
         if(p%on_wall_id%w3 == -1) then
            dphi = atan2(p%v%y, p%v%x) - geo%w3(p%icell%i3)
            if(dphi > pi) dphi = dphi - twopi
            if(dphi < -pi) dphi = dphi + twopi
        end if
         if(p%on_wall_id%w3 == +1) then
            dphi = atan2(p%v%y, p%v%x) - geo%w3(p%icell%i3+1)
            if(dphi > pi) dphi = dphi - twopi
  
negative.code: 
         if(p%on_wall_id%w3 == +1 .and. abs(dphi) < geo%ew3(p%icell%i3+1)) then
  
            iext%w3 = +1
  
         else if(p%on_wall_id%w3 == -1 .and. abs(dphi) < geo%ew3(p%icell%i3)) then
  
            iext%w3 = -1
  
         ! Only check for intersections if the current position is not along
         ! the (x, y) = (0, 0) axis. If the photon is along this axis, then
  