anchor.diff:      cy = l(i)*uz(i)-ux(i)*n(i)
       cz = ux(i)*m(i)-l(i)*uy(i)
       !Rotate about cross vector an angle t2-t1
  -    call rotateaxis(l(i),m(i),n(i),t2-t1,cx,cy,cz)
  +    dt = t2-t1
  +    call rotateaxis(l(i),m(i),n(i),dt,cx,cy,cz)
       !Normalize
       alpha = sqrt(l(i)**2 + m(i)**2 + n(i)**2)
       l(i) = l(i)/alpha
  
positive.code:   real*8, intent(in) :: n1,n2
    real*8, intent(inout) :: l(num),m(num),n(num),ux(num),uy(num),uz(num)
    integer :: i
    real*8 :: dot,t1,t2,alpha,cx,cy,cz
  
    !Loop through rays
    !$omp parallel do private(dot,t1,t2,cx,cy,cz,alpha)
  
negative.code:   real*8, intent(inout) :: tx,ty,tz,rx,ry,rz
    real*8, intent(inout) :: x(num),y(num),z(num),l(num),m(num),n(num),ux(num),uy(num),uz(num)
    integer :: i
    real*8 :: tmp
  
    !Loop through rays
    !$omp parallel do
    do i=1,num
      !Perform z rotation
      tmp = -rz
      call rotatevector(x(i),y(i),z(i),tmp,3)
      call rotatevector(l(i),m(i),n(i),tmp,3)
      call rotatevector(ux(i),uy(i),uz(i),tmp,3)
      !Perform y rotation
      tmp = -ry
      call rotatevector(x(i),y(i),z(i),tmp,2)
      call rotatevector(l(i),m(i),n(i),tmp,2)
      call rotatevector(ux(i),uy(i),uz(i),tmp,2)
      !Perform x rotation
      tmp = -rx
      call rotatevector(x(i),y(i),z(i),tmp,1)
      call rotatevector(l(i),m(i),n(i),tmp,1)
      call rotatevector(ux(i),uy(i),uz(i),tmp,1)
      !Perform translation
      x(i) = x(i) - tx
      y(i) = y(i) - ty
  