! routine to rotate vector about a selected axis and by s specified angle
!Input: p_in(3) - real
!       axis - character*1, 'x','y', or 'z'
!       angle - real, rotation angle in radians
!
!Output: p_out(3) - real
!
subroutine rotate(p_in, axis, angle, p_out)

 use precision_def
 implicit none
 real(rp) p_in(3), angle, p_out(3)

 real(rp) matrix(3,3),s,c

 character*1 axis

! construct matrix

  s=sin(angle)
  c=cos(angle)

 matrix(1:3,1:3) = 0.
 if(axis == 'x')then
   matrix(1,1)=1.
   matrix(2,2) = c
   matrix(2,3) = s
   matrix(3,2)= -s
   matrix(3,3) = c
 endif 
 if(axis == 'y')then
   matrix(1,1)= c
   matrix(1,3) = s
   matrix(2,2) = 1
   matrix(3,1)= -s
   matrix(3,3) = c
 endif 
 if(axis == 'x')then
   matrix(1,1)= c
   matrix(1,2) = s
   matrix(2,1) = -s
   matrix(2,2)= c
   matrix(3,3) = 1.
 endif 

 p_out = matmul(matrix,p_in)

return
end