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_sky_to_ccd(Q,PHI,RLAT,X,Y,Z)
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	X,Y,Z	real	new cartesian coordinates(by passing transformation to spherical)
C OUTPUTS:
C	X,Y,Z	real	new cartesian coordinates(by passing transformation to spherical)
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_sky_to_ccd(q, phi, rlat, xnew, ynew, znew)

	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.0) znew=1.D0
	if(znew.lt.-1.0) znew=-1.D0

	return
	end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c	subroutine rotateq_negpi_to_pi(q, phi, rlat)
c
c	real*8 q(4), phi, rlat, xold, yold, zold, xnew, ynew, znew
c	real*8 r(3,3), q1sq, q2sq, q3sq, q4sq
c
c	!!find cartesian coordinates from ra/dec
c	xold = dcosd(rlat)*dcosd(phi)
c	yold = dcosd(rlat)*dsind(phi)
c	zold = dsind(rlat)
c
c	!!set up rotation matrix r(3,3)
c	q1sq=q(1)*q(1)
c	q2sq=q(2)*q(2)
c	q3sq=q(3)*q(3)
c	q4sq=q(4)*q(4)
c	r(1,1) = q1sq + q2sq - q3sq - q4sq
c	r(1,2) = 2.D0*( q(2)*q(3) + q(1)*q(4) )
c	r(1,3) = 2.D0*( q(2)*q(4) - q(1)*q(3) )
c	r(2,1) = 2.D0*( q(2)*q(3) - q(1)*q(4) )
c	r(2,2) = q1sq - q2sq + q3sq - q4sq
c	r(2,3) = 2.D0*( q(3)*q(4) + q(1)*q(2) )
c	r(3,1) = 2.D0*( q(2)*q(4) + q(1)*q(3) )
c	r(3,2) = 2.D0*( q(3)*q(4) - q(1)*q(2) )
c	r(3,3) = q1sq - q2sq - q3sq + q4sq
c
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)
c
c	!!rotate X,Y,Z old by matrix multiplication
c	xnew = r(1,1)*xold + r(1,2)*yold + r(1,3)*zold
c	ynew = r(2,1)*xold + r(2,2)*yold + r(2,3)*zold
c	znew = r(3,1)*xold + r(3,2)*yold + r(3,3)*zold
c
c	!!at poles |znew| becomes greater than 1 (by less than 10**-6)
c	if(znew.gt.1.0) znew=1.D0
c	if(znew.lt.-1.0) znew=-1.D0
c
c	phi = datan2d(ynew,xnew)
c	rlat = dasind(znew)
c
c	return
c	end