SUBROUTINE ROTREF(NUM,NAXIS,ANGLE,ROT)                                0270
C***********************************************************************    0280
C*                                                                     *    0290
C*    COMPUTE PRODUCT MATRIX OF SEQUENCE OF ROTATIONS AND REFLECTIONS  *    0300
C*                                                                     *    0310
C*    INPUTS                                                           *    0320
C*        NUM = NUMBER OF ROTATIONS AND REFLECTIONS IN SEQUENCE        *    0330
C*        NAXIS = SEQUENCE OF ROTATION AND REFLECTION AXES             *    0340
C*                FOR ROTATIONS USE    1, 2, OR 3                      *    0350
C*                FOR REFLECTIONS USE   -1, -2, OR -3                  *    0360
C*        ANGLE = SEQUENCE OF ROTATION ANGLES IN RADIANS               *    0370
C*                FOR REFLECTIONS THIS ANGLE IS IGNORED (ASSUMED ZERO) *    0380
C*                                                                     *    0390
C*     OUTPUT                                                          *    0400
C*        ROT = 3X3 PRODUCT MATRIX                                     *    0410
C*                                                                     *    0420
C***********************************************************************    0430
      IMPLICIT REAL*8(A-H,O-Z)                                              0440
      DIMENSION ROT(3,3),R1(3,3),R2(3),ANGLE(NUM),NAXIS(NUM)                0450
      EPS = 1D-15                                                           0460
C                                                                           0470
C     ...INITIALIZE (SET 'ROT' = IDENTITY MATRIX)                           0480
C                                                                           0490
      DO 10 I=1,3                                                           0500
      DO 10 J=1,3                                                           0510
      ROT(I,J)=0.                                                           0520
      IF(I.EQ.J)  ROT(I,J)=1.                                               0530
 10   CONTINUE                                                              0540
C                                                                           0550
C      ...PROCESS SEQUENCE OF ROTATIONS AND REFLECTIONS ONE AT A TIME       0560
C                                                                           0570
      DO 100 N=1,NUM                                                        0580
C                                                                           0590
C     ...DEFINE 3 AXES FOR CURRENT ROTATION OR REFLECTION                   0600
C                                                                           0610
      N1=IABS(NAXIS(N))                                                     0620
      N2=N1+1                                                               0630
      IF(N2.GT.3)  N2=N2-3                                                  0640
      N3=N2+1                                                               0650
      IF(N3.GT.3)  N3=N3-3                                                  0660
C                                                                           0670
C     ...DEFINE DIAGONAL ELEMENTS                                           0680
C                                                                           0690
      R1(N1,N1) = 1.                                                        0700
      IF(NAXIS(N).LT.0.)  R1(N1,N1) = -1.                                   0710
      R1(N2,N2) = DCOS(ANGLE(N))                                            0720
      IF(NAXIS(N).LT.0.)  R1(N2,N2) = 1.                                    0730
      R1(N3,N3) = R1(N2,N2)                                                 0740
C                                                                           0750
C     ...DEFINE NONZERO OFF-DIAGONAL ELEMENTS                               0760
C                                                                           0770
      R1(N2,N3) = DSIN(ANGLE(N))                                            0780
      IF(NAXIS(N).LT.0.)  R1(N2,N3) = 0.                                    0790
      R1(N3,N2) = -R1(N2,N3)                                                0800
C                                                                           0810
C     ...DEFINE ZERO OFF-DIAGONAL ELEMENTS                                  0820
C                                                                           0830
      R1(N1,N2) = 0.                                                        0840
      R1(N1,N3) = 0.                                                        0850
      R1(N2,N1) = 0.                                                        0860
      R1(N3,N1) = 0.                                                        0870
C                                                                           0880
C     ...FORM PRODUCT (SET 'ROT' = 'R1' * 'ROT')                            0890
C                                                                           0900
      DO 100 J=1,3                                                          0910
      DO 30 I=1,3                                                           0920
      R2(I) = 0.                                                            0930
      DO 30 K=1,3                                                           0940
 30   R2(I) = R2(I) + R1(I,K)*ROT(K,J)                                      0950
      DO 100 I=1,3                                                          0960
      ROT(I,J) = R2(I)                                                      0970
      IF(DABS(ROT(I,J)).LT.EPS)  ROT(I,J) = 0.                              0980
 100  CONTINUE                                                              0990
      RETURN                                                                1000
      END                                                                   1010