C+
C NAME:
C	rotateq
C PURPOSE:
C	Calculate polar angles in rotated coordinate system
C CATEGORY:
C	Math: coordinate transformation
C CALLING SEQUENCE:
C	call rotateq(Q,PHI,RLAT)
C CALLS:
C	DACOSD,DATAN2D,DCOSD,DSIND,DSQRT,DMOD
C INPUTS:
C	Q(4)	real	quaternion by which rotation is to be performed
C	PHI	real	phase angle in old coordinate system
C	RLAT	real	latitude in old coordinate system (-90<=RLAT<=90)
C OUTPUTS:
C	PHI	real	phase angle in new coordinate system (0<=PHI<360)
C	RLAT	real	latitude in new coordinate system
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	All angles are in degrees.
C MODIFICATION HISTORY:
C	2003 March, Aaron(UCSD/CASS)
C
C-

	subroutine rotateq(q, phi, rlat)

	real*8 q(4), phi, rlat, xold, yold, zold, xnew, ynew, znew
	real*8 r(3,3), q1sq, q2sq, q3sq, q4sq
	
	!!find cartesian coordinates from ra/dec
	xold = dcosd(rlat)*dcosd(phi)
	yold = dcosd(rlat)*dsind(phi)
	zold = dsind(rlat)

	!!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
	
C	print*, 'DCM Matrix:'
C	print*, r(1,1), R(1,2), r(1,3)
C	print*, r(2,1), R(2,2), r(2,3)
C	print*, r(3,1), R(3,2), r(3,3)

	!!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

	!!at poles |znew| becomes greater than 1 (by less than 10**-6)
	if(znew.gt.1.D0) znew=1.D0
	if(znew.lt.-1.D0) znew=-1.D0

	phi = datan2d(ynew,xnew)
	phi = dmod(dmod(phi,360.D0)+360.D0,360.D0) !!0<phi<360
	rlat = dasind(znew)

	return
	end