SUBROUTINE ROTSUM (R11,R12,R13,R21,R22,R23,R31,R32,R33)
C  THIS ROUTINE WILL KEEP A SUMMARY OF ALL ROTATIONS PERFORMED
C   IT IS CALLED FROM ALL THREE ROTATION ROUTINES
      IMPLICIT REAL (A-H,O-Z)
      CHARACTER*80 DUMMY
      LOGICAL DEBUG, DEBUGL, DEBUGN, DEBUGO, DEBUGP, DEBUGI
      COMMON /DEBCOM/ DEBUG, DEBUGL, DEBUGN, DEBUGO, DEBUGP, DEBUGI
      COMMON /ALLROT/ ROTPRD(3,3)
C
      IF (DEBUG) THEN
         WRITE (DUMMY, 3000) ((ROTPRD(I,J),J=1,3),I=1,3)
3000     FORMAT ( 1X,'ROTSUM:',9F7.3)
         CALL DEBUGR( DUMMY( 1: 75) )
      ENDIF
C
      DO 10 I=1,3
         X=ROTPRD(1,I)
         Y=ROTPRD(2,I)
         Z=ROTPRD(3,I)
         ROTPRD(1,I)= X*R11 + Y*R12 + Z*R13
         ROTPRD(2,I)= X*R21 + Y*R22 + Z*R23
         ROTPRD(3,I)= X*R31 + Y*R32 + Z*R33
 10   CONTINUE
C
      IF (DEBUG) THEN
         WRITE (DUMMY, 3000) ((ROTPRD(I,J),J=1,3),I=1,3)
         CALL DEBUGR( DUMMY( 1: 75) )
      ENDIF

      CALL PREULR( THETA, PHI, PSI )
      IF (DEBUG) THEN
        WRITE ( DUMMY,
     .   '('' EULERIAN EQUIVALENT: THETA='',F10.5,'' PHI='',F10.5,
     .   '' PSI='',F10.5)') THETA, PHI, PSI
        CALL DEBUGR( DUMMY( 1: 79) )
      ENDIF

      RETURN
      END

      SUBROUTINE PREULR( THETA, PHI, PSI )

* THIS ROUTINE TAKES THE CUMULATIVE ROTATION MATRIX: ROTPRD
*  AND CONVERTS IT TO THE EQUIVALENT EULERIAN ROTATION ANGLES

      IMPLICIT REAL (A-H,O-Z)
      CHARACTER*80 DUMMY
      LOGICAL DEBUG, DEBUGL, DEBUGN, DEBUGO, DEBUGP, DEBUGI
      COMMON /DEBCOM/ DEBUG, DEBUGL, DEBUGN, DEBUGO, DEBUGP, DEBUGI
      COMMON /ALLROT/ ROTPRD(3,3)
C
      IF (DEBUG) THEN
         WRITE (DUMMY, 3000) ((ROTPRD(I,J),J=1,3),I=1,3)
3000     FORMAT ( 1X,'PREULR:',9F7.3)
         CALL DEBUGR( DUMMY( 1: 75) )
      ENDIF
C


* LET JAMES (JIMMY) J.P. STEWART
*   PRINT THE EULERIAN REPRESENTATION OF THE MATIRX

      CCA = ROTPRD( 3, 3)
      IF ( ABS( CCA) .GT. 0.999999) THEN
         SSA = SQRT(( 1.D0 - ABS( CCA)) * 2.0D0)
      ELSE
         SSA = SQRT( 1.0D0 - CCA * CCA)
      ENDIF

      IF ( SSA .LT. 1.0D-5 ) THEN
* MUST USE 2 X 2  SUB-MATRIX
         CCB = ROTPRD( 1, 1) * ROTPRD( 3, 3)
         SSB = ROTPRD( 1, 2) * ROTPRD( 3, 3)
         CCC = 1.0D0
         SSC = 0.0D0
      ELSE
         CCB = ROTPRD( 3, 1) / SSA
         SSB = ROTPRD( 3, 2) / SSA
         CCC = -ROTPRD( 1, 3) / SSA
         SSC =  ROTPRD( 2, 3) / SSA

         IF ( ABS(ROTPRD(1,2)) .GT. 1.0D-12 ) THEN
           IF((SSB*CCA*CCC+CCB*SSC)/ROTPRD(1,2).GT.0.0D0) GOTO 20
         ELSEIF( ABS(ROTPRD( 1,1)) .GT. 1.0D-12) THEN
           IF((CCB*CCA*CCC-SSB*SSC)/ROTPRD(1,1).GT.0.0D0) GOTO 20
         ENDIF
* SIGN OF SQRT(SSA) IS INCORRECT
         SSA = -SSA
         SSB = -SSB
         CCB = -CCB
         CCC = -CCC
      ENDIF

 20   CONTINUE

      CCA =MAX( MIN( CCA, 1.0D0), -1.0D0)
      CCB =MAX( MIN( CCB, 1.0D0), -1.0D0)
      CCC =MAX( MIN( CCC, 1.0D0), -1.0D0)

      THETA = SIGN( ACOS( CCA), SSA) / 0.01745329252D0
      PHI   = SIGN( ACOS( CCB), SSB) / 0.01745329252D0
      PSI   = SIGN( ACOS( CCC), SSC) / 0.01745329252D0

      IF ( DEBUG ) THEN
         WRITE ( DUMMY, '(''EULERIAN EQUIVALENT: THETA='',
     .  F10.4, '' PHI='', F10.4,'' PSI='', F10.4)') THETA, PHI, PSI
         CALL DEBUGR( DUMMY( 1: 79) )
      ENDIF
C
      RETURN
      END

      SUBROUTINE ROTCLR
* CLEARS ROTATION MATRIX
      IMPLICIT REAL (A-H,O-Z)
      COMMON /ALLROT/ ROTPRD(3,3)

      ROTPRD(1,1)=1.0D0
      ROTPRD(1,2)=0.0D0
      ROTPRD(1,3)=0.0D0

      ROTPRD(2,1)=0.0D0
      ROTPRD(2,2)=1.0D0
      ROTPRD(2,3)=0.0D0

      ROTPRD(3,1)=0.0D0
      ROTPRD(3,2)=0.0D0
      ROTPRD(3,3)=1.0D0

      RETURN
      END