C+
C NAME:
C	rotate
C PURPOSE:
C	Calculate polar angles in rotated coordinate system
C CATEGORY:
C	Math: coordinate transformation
C CALLING SEQUENCE:
	subroutine rotate(ALFA,BETA,GAMMA,PHI,RLAT)
C INPUTS:
C	ALFA	real	phase angle of the pole Z(new) in old coordinate system
C	BETA	real	polar angle of the pole Z(new) in old coordinate system
C	GAMMA	real	phase angle of X(new) in coordinate system (2)
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 CALLS:
C	dsind, dcosd, dacosd, datan2d
C PROCEDURE:
C	Given phase and polar angles of a point in some coordinate system
C	[axes: X(old),Y(old),Z(old)], phase and polar angles are calculated in
C	the new coordinate system [axes: X(new),Y(new),Z(new)] rotated with
C	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	1989, Paul Hick (MPAE,UCSD/CASS; pphick@ucsd.edu)
C-
	    real	ALFA
	    real	BETA
	    real	GAMMA
	    real	PHI
	    real	RLAT

	double precision	SINB
	double precision	COSB
	double precision	SINL
	double precision	COSL
	double precision	SINP
	double precision	COSP
	double precision	COST

	! Needed for g77 when using goniod.f
	! What happens if these are available as intrinsic functions??

	double precision	dsind
	double precision	dcosd
	double precision	dacosd
	double precision	datan2d

	if (abs(RLAT) .gt. 90.0) call Say('rotate','E','invalid','latitude')

	PHI  = PHI-ALFA			! Rotation (1)

	SINB = dsind(dble(BETA))	! Rotation (2)
	COSB = dcosd(dble(BETA))
	SINL = dsind(dble(RLAT))
	COSL = dcosd(dble(RLAT))
	SINP = dsind(dble(PHI ))
	COSP = dcosd(dble(PHI ))

	COST = SINL*COSB+COSL*SINB*COSP	! cos(polar angle)

	if (dabs(COST) .ne. 1.0d0) then
	    SINP = COSL*SINP
	    COSP = COSL*COSB*COSP-SINL*SINB

	else				!Point on new Z-axis; PHI not defined

	    !-------
	    ! The following statement is inserted for the benefit of the
	    ! synoptic map drawing program.
	    ! If BETA=0 then SINB=0, COSB=+/-1, i.e. the z-axis stays the same
	    ! or is flipped by 180 degrees. The input phase angle for points on
	    ! the z-axis will remain the same or is changed to 180-(phase angle)
	    ! (just as for points not on the z-axis).

	    if (SINB .eq. 0.0d0) COSP = COSB*COSP

	end if

	PHI  = datan2d(SINP,COSP)
	PHI  = PHI-GAMMA		! Rotation (3)
	PHI  = amod(amod(PHI,360.0)+360.0,360.0)!0<=PHI<360
	RLAT = 90.0d0-dacosd(COST)	! acosd returns between 0 and 180 deg

	return
	end