C+
C NAME:
C	rotateq
C CATEGORY:
C	Math: coordinate transformation
C CALLING SEQUENCE:
C	call rotateq_xyz(Q,X,Y,Z)
C PROCEDURE:
C	given a unit quaternion, it specifies a rotation about a unit vector.  This rotation 
C	can be put into matrix form and multiplied by cartesian coordinates to get rotated
C	cartesian coordinates.
C MODIFICATION HISTORY:
C	2003 March, Aaron(UCSD/CASS)
C
C-

	subroutine rotateq_xyz(q, xnew, ynew, znew)

	real*8 q(4), xold, yold, zold, xnew, ynew, znew
	real*8 r(3,3), q1sq, q2sq, q3sq, q4sq

	xold = xnew
	yold = ynew
	zold = znew

	!!set up rotation matrix r(3,3)
	q1sq=q(1)*q(1)
	q2sq=q(2)*q(2)
	q3sq=q(3)*q(3)
	q4sq=q(4)*q(4)
	r(1,1) = q1sq + q2sq - q3sq - q4sq
	r(1,2) = 2.D0*( q(2)*q(3) + q(1)*q(4) )
	r(1,3) = 2.D0*( q(2)*q(4) - q(1)*q(3) )
	r(2,1) = 2.D0*( q(2)*q(3) - q(1)*q(4) )
	r(2,2) = q1sq - q2sq + q3sq - q4sq
	r(2,3) = 2.D0*( q(3)*q(4) + q(1)*q(2) )
	r(3,1) = 2.D0*( q(2)*q(4) + q(1)*q(3) )
	r(3,2) = 2.D0*( q(3)*q(4) - q(1)*q(2) )
	r(3,3) = q1sq - q2sq - q3sq + q4sq

	!!rotate X,Y,Z old by matrix multiplication
	xnew = r(1,1)*xold + r(1,2)*yold + r(1,3)*zold
	ynew = r(2,1)*xold + r(2,2)*yold + r(2,3)*zold
	znew = r(3,1)*xold + r(3,2)*yold + r(3,3)*zold

	return
	end
