      PROGRAM REDUCE
C
C     BUERGER CELL REDUCTION
C
C     L. V. AZAROFF AND M. J. BUERGER, THE POWDER METHOD IN X-RAY
C     CRYSTALLOGRAPHY, MC GRAW-HILL, NEW YORK, 1958, PP. 125-159.
C
      DIMENSION A(3),ALPHA(3),U(3),G(3,3),TT(0:2,-2:+2,-2:+2)
      DIMENSION T(3,3),V(3),B(3),BETA(3)
      DIMENSION DATA(75),INDEX(75)
      PI=3.141593
C
C     START WITH AN ARBITRARY PRIMITIVE CELL.
C                             ---------
C
      IO1=10
      IO5=5
      IO6=6
      OPEN(UNIT=IO1,FILE='cell.dat',STATUS='UNKNOWN')
      WRITE (IO6,999)
 999  FORMAT (/1X,'CELL PARAMETERS FROM "cell.dat" FILE?  Y OR N.')
      READ (IO5,899) ANS
 899  FORMAT (A)
 100  IF (ANS.NE.'Y') THEN
        WRITE (IO6,896)
        READ (IO5,*) A,ALPHA
        REWIND IO1
        WRITE (IO1,895) A,ALPHA
      ELSE
        READ (IO1,*) A,ALPHA
      END IF
 896  FORMAT (/1X,'ENTER A, B, C, ALPHA, BETA, GAMMA.')
 895  FORMAT (3F10.5,3F10.4)
      CA=COS(ALPHA(1)*PI/180)
      CB=COS(ALPHA(2)*PI/180)
      CG=COS(ALPHA(3)*PI/180)
      VOL=A(1)*A(2)*A(3)*(1-CA**2-CB**2-CG**2+2*CA*CB*CG)
      WRITE (IO6,898) A,ALPHA,VOL
 898  FORMAT (/1X,'         A         B         C     ALPHA      BETA ',
     &'    GAMMA         V'/1X,3F10.4,3F10.3,F10.2)
      WRITE (IO6,897)
 897  FORMAT (1X,'CORRECT?  Y OR N.')
      READ (IO5,899) ANS
      IF (ANS.NE.'Y') GO TO 100
C
C     STORE SCALAR PRODUCTS OF LATTICE TRANSLATION VECTORS, I.E., THE
C     METRIC TENSOR G(I,J).
C
      G(1,1)=A(1)**2
      G(1,2)=A(1)*A(2)*COS(ALPHA(3)*PI/180)
      G(1,3)=A(1)*A(3)*COS(ALPHA(2)*PI/180)
      G(2,1)=G(1,2)
      G(2,2)=A(2)**2
      G(2,3)=A(2)*A(3)*COS(ALPHA(1)*PI/180)
      G(3,1)=G(1,3)
      G(3,2)=G(2,3)
      G(3,3)=A(3)**2
C
C     STORE IN MATRIX FORM TABLES OF T*T VALUES, WHERE
C     T = U*A + V*B + W*C, FOR U = 0, 1, 2 AND V, W = -2, -1, 0, 1, 2.
C
      DO 20 IU= 0,+2
      DO 20 IV=-2,+2
      DO 20 IW=-2,+2
 20   TT(IU,IV,IW)=0
      DO 21 IU= 0,+2
      U(1)=IU
      DO 21 IV=-2,+2
      U(2)=IV
      DO 21 IW=-2,+2
      U(3)=IW
C
C     FUNCTION VTMV CALCULATES THE BILINEAR FORM,
C     SCALAR = (VECTOR TRANSPOSE)*(MATRIX)*(VECTOR).
C
 21   TT(IU,IV,IW)=VTMV(U,G,U)
C
C     PRINT THE T*T TABLES.
C
      ILP=60
      OPEN (UNIT=ILP,FILE='reduce.lp',STATUS='NEW')
      WRITE (ILP,997) (A(I),ALPHA(I),I=1,3),VOL
 997  FORMAT (1H1,'CELL REDUCTION'//1X,'A = ',F10.4,'    ALPHA = ',
     &F10.3/' B = ',F10.4,'     BETA = ',F10.3/' C = ',F10.4,'    GAMM',
     &'A = ',F10.3//' V = ',F10.2)
      WRITE (ILP,996)
 996  FORMAT (/1X,'MATRIX TABLES OF VALUES OF T*T, WHERE T = U*A + V*B +
     & W*C.'/1X,'TO FIND THE THREE SHORTEST NONCOPLANAR LATTICE TRANSLAT
     &ION VECTORS,'/' LOOK FOR THE U, V, W FOR THE THREE SMALLEST T*T VA
     &LUES.')
      DO 22 IU=0,2
      WRITE (ILP,995) IU
 995  FORMAT (/1X,'U = ',I2,8X,'V = -2',8X,'-1',9X,'0',9X,'1',9X,'2')
      IW=-2
      WRITE (ILP,994) IW,(TT(IU,IV,IW),IV=-2,+2)
 994  FORMAT (/1X,'W = ',I2,4X,5F10.3)
      DO 22 IW=-1,+2
      WRITE (ILP,993) IW,(TT(IU,IV,IW),IV=-2,+2)
 993  FORMAT (5X,I2,4X,5F10.3)
 22   CONTINUE
C
C     PRINT SORTED LIST OF THE T*T MATRIX ELEMENTS.
C
      WRITE (IO6,'(/1X,''PRINT SORTED LIST OF THE T*T MATRIX ELEMENTS'',
     &''?  Y OR N.'')')
      READ (IO5,'(A1)') ANS
      IF (ANS.EQ.'Y') THEN
        WRITE (ILP,970)
        N=0
        DO IU= 0,+2
        DO IV=-2,+2
        DO IW=-2,+2
          N=N+1
          DATA(N)=TT(IU,IV,IW)
        END DO
        END DO
        END DO
        CALL SORT (N,DATA,INDEX)
        DO I=1,N
          DO IU= 0,+2
          DO IV=-2,+2
          DO IW=-2,+2
            IF (DATA(INDEX(I)).EQ.TT(IU,IV,IW)) WRITE (ILP,971)
     &      IU,IV,IW,TT(IU,IV,IW)
          END DO
          END DO
          END DO
        END DO
      END IF
 970  FORMAT (/1X,' SORTED LIST OF T*T MATRIX ELEMENTS:'//1X,' IU IV I',
     &'W  T*T(IU,IV,IW)'/1X,' -- -- --  -------------')
 971  FORMAT (1X,3I3,F10.3)
C
C     TRANSFORM TO NEW UNIT CELL.
C
      WRITE (IO6,980)
 980  FORMAT (/1X,'ENTER THE TRANSFORMATION MATRIX?  Y OR N.')
      READ (IO5,981) ANS
 981  FORMAT (A1)
      IF (ANS.NE.'Y') STOP
      WRITE (IO6,992)
 992  FORMAT (/1X,'NEW A = T11*A + T12*B + T13*C'/' NEW B = T21*A + T22*
     &B + T23*C'/' NEW C = T31*A + T32*B + T33*C'//1X,'ENTER T11, T12, T
     &13.')
      READ (IO5,*) T(1,1),T(1,2),T(1,3)
      WRITE (IO6,991)
 991  FORMAT (' ENTER T21, T22, T23.')
      READ (IO5,*) T(2,1),T(2,2),T(2,3)
      WRITE (IO6,990)
 990  FORMAT (' ENTER T31, T32, T33.')
      READ (IO5,*) T(3,1),T(3,2),T(3,3)
      WRITE (ILP,989) ((T(I,J),J=1,3),I=1,3)
 989  FORMAT (/1X,'TRANSFORMATION FROM GIVEN TO TRANSFORMED LATTICE TRAN
     &SLATION VECTORS'//1X,'T11 T12 T13  ',3F5.1/1X,'T21 T22 T23  ',
     &3F5.1/' T31 T32 T33  ',3F5.1)
      DO 31 I=1,3
      DO 32 J=1,3
 32   U(J)=T(I,J)
 31   B(I)=SQRT(VTMV(U,G,U))
      DO 41 I=1,3
      J=I+1
      K=I+2
      IF (J.GT.3) J=J-3
      IF (K.GT.3) K=K-3
      DO 42 L=1,3
      U(L)=T(J,L)
 42   V(L)=T(K,L)
 41   BETA(I)=ACOS(VTMV(U,G,V)/SQRT(VTMV(U,G,U)*VTMV(V,G,V)))*180/PI
      CA=COS(BETA(1)*PI/180)
      CB=COS(BETA(2)*PI/180)
      CG=COS(BETA(3)*PI/180)
      VOL=B(1)*B(2)*B(3)*(1-CA**2-CB**2-CG**2+2*CA*CB*CG)
      WRITE (ILP,988) (B(I),BETA(I),I=1,3),VOL
 988  FORMAT (/1X,'TRANSFORMED CELL PARAMETERS'//1X,'A = ',F10.4,'    AL
     &PHA = ',F10.3/1x,'B = ',F10.4,'     BETA = ',F10.3/1x,'C = ',
     &F10.4,'    GAMMA = 'F10.3//1X,'V = ',F10.2)
      WRITE (IO1,895) B,BETA
      CLOSE (UNIT=IO1)
C
C     EVALUATE AND PRINT NIGGLI MATRIX OF SCALAR PRODUCTS.
C
      S11=B(1)**2
      S22=B(2)**2
      S33=B(3)**2
      S23=B(2)*B(3)*COS(BETA(1)*PI/180)
      S31=B(3)*B(1)*COS(BETA(2)*PI/180)
      S12=B(1)*B(2)*COS(BETA(3)*PI/180)
      WRITE (ILP,987) S11,S22,S33,S23,S31,S12
 987  FORMAT (/1X,'NIGGLI MATRIX FOR THE TRANSFORMED CELL'//1X,'S11 S2',
     &'2 S33  ',3F10.4/1X,'S23 S31 S12  ',3F10.4)
      STOP
      END
C-----------------------------------------------------------------------
      FUNCTION VTMV(U,G,V)
C
C     BILINEAR FORM = (VECTOR TRANSPOSE)*(MATRIX)*(VECTOR) = SCALAR
C     BILINEAR FORM = (ROW MATRIX)*(SQUARE MATRIX)*(COLUMN MATRIX)
C
      DIMENSION U(3),G(3,3),V(3)
      VTMV=0
      DO 1 I=1,3
      DO 1 J=1,3
    1 VTMV=VTMV+U(I)*G(I,J)*V(J)
      RETURN
      END
C------------------------------------------------------------------------
      SUBROUTINE SORT (N,DATA,INDEX)
C
C INDEXES THE ARRAY DATA(N) AND RETURNS THE ARRAY INDEX(N) SORTED SUCH
C THAT THE VALUES DATA(INDEX(I)) ARE IN ASCENDING ORDER FOR I = 1, 2,
C ..., N.  THE INPUT VARIABLES N AND DATA(N) ARE NOT CHANGED.
C
C EMPLOYS THE "HEAPSORT" ALGORITHM.
C
C FORTRAN CODE QUOTED FROM WILLIAM H. PRESS, BRIAN P. FLANNERY, SAUL A.
C TEUKOLSKY, AND WILLIAM T. VETTERLING (1986).  NUMERICAL RECIPIES:  THE
C                                               --------- --------   ---
C ART OF SCIENTIFIC COMPUTING, PP. 229-233.  CAMBRIDGE, ENGLAND:
C --- -- ---------- ---------
C CAMBRIDGE UNIVERSITY PRESS.
C
      REAL DATA,T
      DIMENSION DATA(N),INDEX(N)
      DO I=1,N
        INDEX(I)=I
      END DO
      IF (N.EQ.1) RETURN
      I=N/2+1
      J=N
 1    CONTINUE
      IF (I.GT.1) THEN
        I=I-1
        INEXT=INDEX(I)
        T=DATA(INEXT)
      ELSE
        INEXT=INDEX(J)
        T=DATA(INEXT)
        INDEX(J)=INDEX(1)
        J=J-1
        IF (J.EQ.1) THEN
          INDEX(1)=INEXT
          RETURN
        END IF
      END IF
      K=I
      L=I+I
      DO WHILE (L.LE.J)
        IF (L.LT.J) THEN
          IF (DATA(INDEX(L)).LT.DATA(INDEX(L+1))) L=L+1
        END IF
        IF (T.LT.DATA(INDEX(L))) THEN
          INDEX(K)=INDEX(L)
          K=L
          L=L+L
        ELSE
          L=J+1
        END IF
      END DO
      INDEX(K)=INEXT
      GO TO 1
      END
C-----------------------------------------------------------------------
