c
c=====================================================================
c
c	ROTATION.F
c
c	A collection of utility subroutines for doing coordinate rotations
c
c	rotate -- rotation about an arbitrary axis
c	rotate_x -- rotation about the X axis
c	rotate_y -- rotation about the Y axis
c	rotate_z -- rotation about the Z axis
c
c=====================================================================
c=====================================================================
c
	subroutine rotate(x1, y1, z1, ax, ay, az, angle, x2, y2, z2)
c
	implicit none
c
	real x1, y1, z1		! input vector
	real ax, ay, az		! rotation axis
	real angle		! rotation angle in degrees
	real x2, y2, z2		! output vector
c
	real theta
	real ax2, ay2, az2
	real xtemp, ytemp, ztemp
	real thetax, thetay
c
	double precision axd, ayd, azd
c
c=====================================================================
c
c----- Copy input rotation angle into local variable
c
	theta = angle
c
c----- Get rotation angle in range (-180, 180)
c
	do while (theta .gt. 180.0) 
	  theta = theta - 360.0
	enddo
c
	do while (theta .lt. -180.0)
	  theta = theta + 360.0
	enddo
c
c----- Check for null case
c
	if (theta .eq. 0.0) then
	  x2 = x1
	  y2 = y1
	  z2 = z1
c
c		COULD ALSO CHECK HERE FOR INPUT VECTOR
c		PARALLEL TO ROTATION AXIS
c
	  return
	endif
c
c=====================================================================
c
c	Here's the solution in this code ...
c
c	-- Rotate the axis vector into the positive Z axis
c	    -- first rotate about X axis into XZ plane
c	    -- then rotate about Y axis into YZ plane
c	-- Rotate by angle aound the Z axis (now the axis vector)
c	-- Reverse the rotations to restore the axis vector
c	    -- reverse rotation about Y axis out of the YZ plane
c	    -- then reverse rotation about X axis out of the XZ plane
c
c=====================================================================
c
c----- Rotate the axis vector into the XZ plane
c
	if (az .eq. 0.0) then
c
	  if (ay .eq. 0.0) then
	    thetax = 0.0
	  else
	    thetax = 90.0*sign(1.0,ay)
	  endif
c
	else
c
	  ayd = dble(ay)
	  azd = dble(az)
	  thetax = sngl(atan(ayd/abs(azd)))*57.29578
c
ccc	  thetax = atan(ay/abs(az))
c
	  if (az .lt. 0.0) thetax = sign(1.0,ay)*180.0 - thetax
c
	endif
c
	if (thetax .eq. 0.0) then
	  x2 = x1
	  y2 = y1
	  z2 = z1
	  ax2 = ax
	  ay2 = ay
	  az2 = az
	else
	  call rotate_x(x1, y1, z1, thetax, x2, y2, z2)
	  call rotate_x(ax, ay, az, thetax, ax2, ay2, az2)
	endif
c
c=====================================================================
c
c----- Rotate the axis vector into the YZ plane --> hence the Z axis
c	after the above rotation
c
	if (az2 .eq. 0.0) then
c
	  if (ax2 .eq. 0.0) then
	    thetay = 0.0
	  else
	    thetay = -90.0*sign(1.0,ax2)
	  endif
c
	else
c
	  axd = dble(ax2)
	  azd = dble(az2)
	  thetay = -1.0*sngl(atan(axd/abs(azd)))*57.29578
c
ccc	  thetay = -1.0*atan(ax2/abs(az2))
c
	  if (az2 .lt. 0.0) thetay = -1.0*sign(1.0,ay2)*180.0 - thetay
c
	endif
c
	if (thetay .eq. 0.0) then
	  xtemp = x2
	  ytemp = y2
	  ztemp = z2
	else
	  call rotate_y(x2, y2, z2, thetay, xtemp, ytemp, ztemp)
	endif
c
c=====================================================================
c
c----- Rotate around the Z axis - now the rotation axis
c
	call rotate_z(xtemp, ytemp, ztemp, theta, x2, y2, z2)
c
c=====================================================================
c
c----- Reverse rotation around the Y axis
c
	if (thetay .eq. 0.0) then
	  xtemp = x2
	  ytemp = y2
	  ztemp = z2
	else
	  thetay = -1.0*thetay
	  call rotate_y(x2, y2, z2, thetay, xtemp, ytemp, ztemp)
	endif
c
c=====================================================================
c
c----- Reverse rotation around the X axis
c
	if (thetax .eq. 0.0) then
	  x2 = xtemp
	  y2 = ytemp
	  z2 = ztemp
	else
	  thetax = -1.0*thetax
	  call rotate_x(xtemp, ytemp, ztemp, thetax, x2, y2, z2)
	endif
c
c=====================================================================
c
	return
	end
c
c
c
c=====================================================================
c=====================================================================
c
	subroutine rotate_x(x1, y1, z1, angle, x2, y2, z2)
c
	implicit none
c
	real x1, y1, z1		! input vector
	real angle		! rotation angle in degrees
	real x2, y2, z2		! output vector
c
	double precision theta
c
c=====================================================================
c
c----- Copy input rotation angle into local variable
c
	theta = dble(angle)
c
c----- Get rotation angle in range (-180, 180)
c
	do while (theta .gt. 180.0) 
	  theta = theta - 360.0
	enddo
c
	do while (theta .lt. -180.0)
	  theta = theta + 360.0
	enddo
c
c----- Check for null case
c
	if (theta .eq. 0.0) then
	  x2 = x1
	  y2 = y1
	  z2 = z1
c
	  return
	endif
c
c=====================================================================
c
	theta = theta/57.29578
c
	x2 = x1
c
	y2 = y1*cos(theta) - z1*sin(theta)
c
	z2 = y1*sin(theta) + z1*cos(theta)
c
c=====================================================================
c
	return
	end
c
c
c
c=====================================================================
c=====================================================================
c
c
	subroutine rotate_y(x1, y1, z1, angle, x2, y2, z2)
c
	implicit none
c
	real x1, y1, z1		! input vector
	real angle		! rotation angle in degrees
	real x2, y2, z2		! output vector
c
	double precision theta
c
c=====================================================================
c
c----- Copy input rotation angle into local variable
c
	theta = dble(angle)
c
c----- Get rotation angle in range (-180, 180)
c
	do while (theta .gt. 180.0) 
	  theta = theta - 360.0
	enddo
c
	do while (theta .lt. -180.0)
	  theta = theta + 360.0
	enddo
c
c----- Check for null case
c
	if (theta .eq. 0.0) then
	  x2 = x1
	  y2 = y1
	  z2 = z1
c
	  return
	endif
c
c=====================================================================
c
	theta = theta/57.29578
c
	x2 = z1*sin(theta) + x1*cos(theta)
c
	y2 = y1
c
	z2 = z1*cos(theta) - x1*sin(theta)
c
c=====================================================================
c
	return
	end
c
c
c
c=====================================================================
c=====================================================================
c
c
	subroutine rotate_z(x1, y1, z1, angle, x2, y2, z2)
c
	implicit none
c
	real x1, y1, z1		! input vector
	real angle		! rotation angle in degrees
	real x2, y2, z2		! output vector
c
	double precision theta
c
c=====================================================================
c
c----- Copy input rotation angle into local variable
c
	theta = dble(angle)
c
c----- Get rotation angle in range (-180, 180)
c
	do while (theta .gt. 180.0) 
	  theta = theta - 360.0
	enddo
c
	do while (theta .lt. -180.0)
	  theta = theta + 360.0
	enddo
c
c----- Check for null case
c
	if (theta .eq. 0.0) then
	  x2 = x1
	  y2 = y1
	  z2 = z1
c
	  return
	endif
c
c=====================================================================
c
	theta = theta/57.29578
c
	x2 = x1*cos(theta) - y1*sin(theta)
c
	y2 = x1*sin(theta) + y1*cos(theta)
c
	z2 = z1
c
c=====================================================================
c
	return
	end