C+
C NAME:
C	rotatexyz
C PURPOSE:
C	Calculate cartesian coordinates in rotated coordinate system
C CATEGORY:
C	Math: coordinate transformation
C CALLING SEQUENCE:
C	call Rotate(ALFA,BETA,GAMMA,X,Y,Z)
C CALLS:
C	DACOSD,DATAN2D,DCOSD,DSIND,DSQRT,DMOD
C INPUTS:
C	ALFA	real*8		phase angle of the pole Z(new) in old coordinate system
C	BETA	real*8		polar angle of the pole Z(new) in old coordinate system
C	GAMMA	real*8		phase angle of X(new) in coordinate system (2)
C	X	real*8		X-axis coordinate old coordinate system
C	Y	real*8		Y-axis coordinate old coordinate system
C	Z	real*8		Z-axis coordinate old coordinate system
C OUTPUTS:
C	X	real*8		X-axis coordinate new coordinate system
C	Y	real*8		Y-axis coordinate new coordinate system
C	Z	real*8		Z-axis coordinate new coordinate system
C PROCEDURE:
C	New cartesian coordinates are calculated in the new coordinate system 
C	rotated with respect to the original one over the Euler angles ALFA, BETA, GAMMA.
C	The rotation from old to new system is realized by
C	(1) rotation around Z(old) over ALFA  ---> X(1),Y(1),Z(1)=Z(old)
C	(2) rotation around Y(1)   over BETA  ---> X(2),Y(2)=Y(1),Z(2)=Z(new)
C	(3) rotation around Z(new) over GAMMA ---> X(new),Y(new),Z(new).
C	All angles are in degrees.
C MODIFICATION HISTORY:
C	2003 March, Aaron (UCSD/CASS)
C-
	subroutine rotatexyz(alfa, beta, gamma, x, y, z)

	real*8 sina, cosa, sinb, cosb, sing, cosg
	real*8 alfa, beta, gamma, x, y, z, xnew, ynew, znew

	sina = dsind(alfa)
	cosa = dcosd(beta)
	sinb = dsind(beta)
	cosb = dcosd(beta)
	sing = dsind(gamma)
	cosg = dcosd(gamma)

	xnew = x*(cosg*cosa - cosb*sina*sing) + y*(cosg*sina + cosb*cosa*sing) + z*sing*sinb
	ynew = x*(-sing*cosa - cosb*sina*cosg) + y*(-sing*sina + cosb*cosa*cosg) + z*cosg*sinb
	znew = x*sinb*sina - y*sinb*cosa + z*cosb
	x = xnew
	y = ynew
	z = znew

	return
	end