C TRACER - PROGRAM TO CHECK CELL
C===============================
C
C
C TRACER:
C-------
C
C Tracer by S.L. Lawton, A Lattice Transformation - Cell Reduction Program.
C
CThe reduced cell is defined to be the unit cell with the smallest volume
Cwhose axes are the three shortest non-coplanar translations in the lattice.
CTRACER converts the cell parameters which you put in as input to the reduced
Ccell. Such a reduction to a triclinic lattice with c < a < b and the angles
Calpha and beta as non-acute allows for the identification of the Bravais
Clattice of highest symmetry which it represents. This is accomplished by
Ccalculation of the scalar products
C
C  A.A, B.B, C.C, B.C, C.A, A.B and comparison to matrices within the
Cprogram.
C
CThis is facilitated by using the reduced cell and considering all possible
Ccombinations of axial lengths in the 14 Bravais lattices. There exists
Conly 43 independent cells whose axes correspond to the 3 shortest non-coplanar
Ctranslations in the lattice.
C
CThe one input parameter which might require clarification is the value of
CDEL to be used. Lawton reports that for lattice parameters with errors
Cas low as 1 part in 10,000 Angstrom that a DEL = 0.10 is sufficient, while
Cerrors of the order of one part in 1000 Angstrom might require DEL values
Cin the range of 0.50 to 1.00 (for Proteins use at least 1.0). The DEL value
Cimposes the upper limit on the equivalence of special relationships between
Cthe scalar products used in the program. It essentially defines the absolute
Cvalue of the largest allowable difference between any two of the six reduced
Ccell scalar products. In as much as a DEL value to any extreme is undesirable
Cit has been our practice to run TRACER using DEL values of 0.10, 0.50 and
C1.00 in the same run and then examine the output as a function of these
Cdifferent DEL values. Usually the higher the standard deviation in the cell
Cparameters the higher the DEL value that should be used. 
C
C
C INPUT
C------
C
C
C----line 1 RUn number (terminates if equal to -1)
C
C---- line 2 TITLE card
C
C---- line 3 CONTROL input DEL
C            CONTROL = 1 for direct space
C                    = 2 for reciprocal space cell dimensions
C
C            DEL = Absolute value of the largest permissible difference
C                  between any two of the six reduced cell scalar products.
C                  
C---- line 4 cell dimensions (free format)
C            a, b, c, alpha, beta, gamma (Angstrom /degrees)
C
C---- line 5 INSTRUCTION (format I3  !!!!!!!!!)
C
C       This card(s) function as the commands to the program in
C       transforming one cell to another. Each card contains an instruction.
C       As many cards as necessary may be used to complete the desired
C       successive transformations for any one compound.
C       Their order specifies to the program the order in which the
C       transformations are to be executed. The last instruction for a
C       particular sequence of transformations must have either 7 or 8!!!!
C
C       Both may not be present and neither may precede any with
C       1,,2,3,4,5,6.
C      
C       NQ = 7 triggers cell reduction.
C       NQ = 8 merely identifies the last cell in a sequence to which a
C              transformation is to be made with no cell reduction is performed.
C
CINSTRUCTIONS
C=============
C
C  NQ = 1   :Transform this cell to a new cell using a matrix supplied by
C            the user on the corresponding transformation card.
C
C  NQ = 2   :The initial cell is A-centred, transform to a primitive cell
C            using a matrix held in subroutine CENTER
C
C  NQ = 3   :..... B-Centred........
C
C  NQ = 4   :...... C-Centred .......
C
C  NQ = 5   :...... I-Centred ......
C
C  NQ = 6   :........ F-Centred .......
C
C  NQ = 7   :Get reduced cell and then identify and generate from it the
C            unit cell of highest symmetry and its corresponding lattice
C            parameters
C
C  NQ=8   : This is a non-executive instruction notifying the program
C           that no transformation is to be applied to this cell.
C           This may be used for;
C           (i) to calculate the inverse cell data of the input parameters,
C           or
C           (ii) to identify the name of the cell to which the previous cell
C           is transformed when the previous cell in a transformation
C           sequence involves any of NQ = 1,2,3,4,5,6)
C
C
C TRANSFORMATION cards
C
C each contains a 3 x 3 matrix one per card  for each NQ = 1
C
C in the order
C 1,1  1,2  1,3  2,1  2,2  2,3  3,1  3,2  3,3
C
C
C
C    create input file testTR.DAT
C  with
C      line 1 : run number
C      line 2 : title
C      line 3 : control card TYPE and DEL
C            TYPE=1 direct space input cell
C            TYPE=2 reciprocal space input cell
C            DEL= allowed error for transformation match
C            DEL = 0.10, 0.50, 1.00 ie
C                  use three consective runs in one file
C 
C      line 4 : cell parameters
C      line 5 : INSTRUCTION number (FORMAT I3 !!!!!!!)
C            =1 you supply matrix
C            =2 A to P
C            =3 B to P
C            =4 C to P
C            =5 I to P
C            =6 F to P
C            =7 find reduced cell (normal value)
C       if you use 1 then next line has a 3x3 matrix on it
C       
C      next line can be next run number and repeat input
C      last line : -1 in columns 1 and 2 to terminate job
C      all other lines FREE FORMAT
C 
C to use
C            RUN TRACER
C      input file name testTR.DAT
C      results in file testTR.OUT
C      Sample TESTtr.dat
C       ------------
C       1                            !run number
C       test tracer                  ! title
C       1  0.1                       ! real cell/ del
C       10.51 15.15 6.54 90 151.7 90 ! cell
C         4                          ! supplied cell is C-faced
C         7                          ! find reduced cell
C       -1                           ! terminate run
C       ----------
C
C
C   TRACER program
C     A2       CELL PARAMETER, A, IN NEW UNIT CELL                       
C     ABC      A*B*COS(GAMMA)                                            
C     ACC      A*C*COS(BETA)                                             
C     ALPHA2   CELL PARAMETER, ALPHA, IN NEW UNIT CELL                   
C     B2       CELL PARAMETER, B, IN NEW UNIT CELL                       
C     BCC      B*C*COS(ALPHA)                                            
C     BETA2    CELL PARAMETER, BETA, IN NEW UNIT CELL                    
C     C2       CELL PARAMETER, C, IN NEW UNIT CELL                       
C     DEL      LARGEST EST  ERROR BETWEEN REDUCED CELL SCALAR PRODUCTS   
C     DETEU    MODULUS OF MATRIX EU                                      
C     DETEUI   MODULUS OF MATRIX EUI                                     
C     DETFN    MODULUS OF MATRIX F                                       
C     DETFNI   MODULUS OF MATRIX FI                                      
C     EU(9)    MATRIX ELEM  FOR TRANSF  OF CONV  RED  CELL TO UNIT CELL  
C     EUI(9)   INVERSE MATRIX ELEMENTS OF EU                             
C     F(9)     MATRIX ELEM  FOR TRANSF  OF ORIG  CELL TO UNIT CELL       
C     FI(9)    INVERSE MATRIX ELEMENTS OF F                              
C     GAMMA2   CELL PARAMETER, GAMMA, IN NEW UNIT CELL                   
C     ISYM     CODE IDENTIFYING EQUIVALENCY AMONG SYMMETRICAL SCALARS    
C     JSUM     ISUM + 1                                                  
C     JU       CODE IDENTIFYING SET OF UNSYMMETRICAL SCALAR PRODUCTS     
C     JU12     CODE FOR AB IN TEST FOR ZERO ELEMENT                      
C     JU23     CODE FOR BC IN TEST FOR ZERO ELEMENT                      
C     JU31     CODE FOR CA IN TEST FOR ZERO ELEMENT                      
C     ITT      CODE SPECIFYING SECTIONS OF SUBROUTINE TEST               
C               
C
      PROGRAM TRACER
C     ==============
C
C
C
C     .. Scalars in Common ..
      REAL A,AL,ALPHA,ALS,ALSTAR,AS,B,BE,BES,BESTAR,BETA,BS,C,COSAL,
     +     COSALS,COSBE,COSBES,COSGA,COSGAS,CS,DEL,GA,GAMMA,GAS,GASTAR,
     +     SINAL,SINALS,SINBE,SINBES,SINGA,SINGAS,V,VS
      INTEGER ICYCLE,IOUTPT,ISUM,JSUM,K,L,L2,LUNIN,LUNOUT,MATRIX
      CHARACTER TITLE*70
C     ..
C     .. Arrays in Common ..
      REAL BET(15),DET(15),DIS(6,12),P(3,3,12),S(3,3,15)
      INTEGER IHKL(12),MM(3,3,4),NN(3,3,3)
C     ..
C     .. Local Scalars ..
      REAL AA,AB,BB,BC,CA,CC,SUM
      INTEGER I,ICONT,INPUT,IS,ITL,ITOTAL,ITT,J,J1,J2,JB2,K4,KX,LINE,MA
      CHARACTER FNAME*10
C     ..
C     .. Local Arrays ..
      REAL CELL(25,12)
      INTEGER NQ(12)
C     ..
C     .. External Subroutines ..
      EXTERNAL CENTER,DIMEN,LATCON,MATIN,MATOFC,RCELL,STORE,TEST
C     ..
C     .. Common blocks ..
      COMMON /INOUT/LUNIN,LUNOUT
      COMMON /TCELL/A,B,C,ALPHA,BETA,GAMMA,COSAL,COSBE,COSGA,SINAL,
     +       SINBE,SINGA,ALS,BES,GAS,SINALS,SINBES,SINGAS,COSALS,COSBES,
     +       COSGAS,AS,BS,CS,AL,BE,GA,ALSTAR,BESTAR,GASTAR,V,VS
      COMMON /TCHAR/TITLE
      COMMON /WORK/DIS,K,IHKL,DEL,ISUM,JSUM,P,S,DET,BET,L,L2,MATRIX,NN,
     +       MM,ICYCLE,IOUTPT
C     ..
C
C
      LUNIN = 5
      LUNOUT = 6
C
C---- open files
C
      WRITE (LUNOUT,FMT=6080)
      READ (LUNIN,FMT=6082) FNAME
      OPEN (UNIT=21,ACCESS='SEQUENTIAL',FILE=FNAME,FORM='FORMATTED',
     +     STATUS='OLD')
      FNAME = FNAME(1:6)//'.OUT'
      OPEN (UNIT=20,FILE=FNAME,ACCESS='SEQUENTIAL',FORM='FORMATTED',
     +     STATUS='NEW')
      WRITE (20,FMT=6084)
      WRITE (20,FMT=6086)
      WRITE (20,FMT=6088)
      WRITE (20,FMT=6090)
C
C---- input
C
      LUNIN = 21
      LUNOUT = 20
   10 CONTINUE
      READ (LUNIN,FMT=*) ICONT
      IF (ICONT.NE.-1) THEN
        READ (LUNIN,FMT=6000) TITLE
        READ (LUNIN,FMT=*) INPUT,DEL
        GO TO (20,30) INPUT
   20   READ (LUNIN,FMT=*) A,B,C,ALPHA,BETA,GAMMA
        IS = 1
        GO TO 40
   30   READ (LUNIN,FMT=*) AS,BS,CS,ALSTAR,BESTAR,GASTAR
        IS = 3
   40   ITOTAL = 0
        MATRIX = 0
   50   CONTINUE
        ITOTAL = ITOTAL + 1
        READ (LUNIN,FMT=6008) NQ(ITOTAL), (CELL(I,ITOTAL),I=1,25)
        IF (NQ(ITOTAL).EQ.1) MATRIX = MATRIX + 1
        I = NQ(ITOTAL)
        GO TO (50,50,50,50,50,50,70,60) I
        GO TO 70
   60   ISUM = ITOTAL
        ITL = ITOTAL
        GO TO 80
   70   ISUM = ITOTAL + 2
        ITL = ITOTAL + 1
   80   IF (MATRIX.GT.0) THEN
          DO 90 K = 1,ITOTAL
            IF (NQ(K).EQ.1) READ (LUNIN,FMT=*) ((P(I,J,K),J=1,3),I=1,3)
   90     CONTINUE
          K = 1
          K4 = 1
          CALL MATIN(K4)
        END IF
C
C---- Preliminary output
C
        IOUTPT = 1
        ICYCLE = 0
        K = 1
        IHKL(1) = 1
        CALL LATCON(IS)
        CALL STORE
        IF (NQ(1).NE.8) THEN
  100     CONTINUE
          WRITE (LUNOUT,FMT=6010) TITLE
          WRITE (LUNOUT,FMT=6012)
          DO 110 J = 1,ITOTAL
            WRITE (LUNOUT,FMT=6014) J, (CELL(I,J),I=1,25)
  110     CONTINUE
          IF (NQ(ITOTAL).EQ.7) THEN
            J1 = ITOTAL + 1
            J2 = ITOTAL + 2
            WRITE (LUNOUT,FMT=6016) J1
            WRITE (LUNOUT,FMT=6018) J2
          END IF
  120     CONTINUE
C
C---- Call desired transformation subroutine
C
          I = NQ(K)
          GO TO (140,130,130,130,130,130,150) I
  130     CALL CENTER(K,NQ,P)
          K4 = 1
          CALL MATIN(K4)
  140     CALL DIMEN
          IS = 1
          CALL LATCON(IS)
          CALL STORE
          IF ((K-1).LE.0) THEN
            GO TO 100
          ELSE IF ((ISUM-K).LE.0) THEN
            GO TO 160
          ELSE
            GO TO 120
          END IF
  150     CALL RCELL
  160     MA = 1
          CALL MATOFC(MA,K4,ITL)
C
C---- Final output
C
          IOUTPT = 3
          IF (NQ(ITOTAL).NE.8) THEN
            LINE = 0
            ITT = 1
            CALL TEST(ITT,LINE)
          END IF
          WRITE (LUNOUT,FMT=6020)
          K = 1
          IF ((ISUM-5).GE.0) THEN
            LINE = 3
          ELSE
            LINE = 2
          END IF
  170     CONTINUE
          I = NQ(K)
          GO TO (180,180,180,180,180,180,190) I
  180     KX = K + 1
          WRITE (LUNOUT,FMT=6022) K,KX,KX,K
          WRITE (LUNOUT,FMT=6024) P(1,1,K),P(1,2,K),P(1,3,K),S(1,1,K),
     +      S(1,2,K),S(1,3,K)
          WRITE (LUNOUT,FMT=6024) P(2,1,K),P(2,2,K),P(2,3,K),S(2,1,K),
     +      S(2,2,K),S(2,3,K)
          WRITE (LUNOUT,FMT=6024) P(3,1,K),P(3,2,K),P(3,3,K),S(3,1,K),
     +      S(3,2,K),S(3,3,K)
          WRITE (LUNOUT,FMT=6074) DET(K),BET(K)
          GO TO 200
  190     L = 1
          KX = K + 1
          WRITE (LUNOUT,FMT=6022) K,KX,KX,K
          WRITE (LUNOUT,FMT=6070) NN(1,1,L),NN(1,2,L),NN(1,3,L),
     +      MM(1,1,L),MM(1,2,L),MM(1,3,L)
          WRITE (LUNOUT,FMT=6070) NN(2,1,L),NN(2,2,L),NN(2,3,L),
     +      MM(2,1,L),MM(2,2,L),MM(2,3,L)
          WRITE (LUNOUT,FMT=6070) NN(3,1,L),NN(3,2,L),NN(3,3,L),
     +      MM(3,1,L),MM(3,2,L),MM(3,3,L)
          WRITE (LUNOUT,FMT=6074) DET(K),BET(K)
          LINE = LINE + 1
          IF ((LINE-5).GE.0) THEN
            LINE = 1
            WRITE (LUNOUT,FMT=6010) TITLE
          END IF
          L = 2
          KX = K + 2
          WRITE (LUNOUT,FMT=6022) K,KX,KX,K
          WRITE (LUNOUT,FMT=6070) NN(1,1,L),NN(1,2,L),NN(1,3,L),
     +      MM(1,1,L),MM(1,2,L),MM(1,3,L)
          WRITE (LUNOUT,FMT=6070) NN(2,1,L),NN(2,2,L),NN(2,3,L),
     +      MM(2,1,L),MM(2,2,L),MM(2,3,L)
          WRITE (LUNOUT,FMT=6070) NN(3,1,L),NN(3,2,L),NN(3,3,L),
     +      MM(3,1,L),MM(3,2,L),MM(3,3,L)
          WRITE (LUNOUT,FMT=6074) DET(K+1),BET(K+1)
  200     LINE = LINE + 1
          IF ((LINE-5).GE.0) THEN
            LINE = 1
            WRITE (LUNOUT,FMT=6010) TITLE
          END IF
          K = K + 1
          IF ((ITL-K).GT.0) GO TO 170
          IF ((ITL-3).GE.0) THEN
            IF (NQ(ITOTAL).EQ.8) GO TO 210
          ELSE IF (NQ(1).NE.7) THEN
            GO TO 250
          END IF
          ITT = 2
          CALL TEST(ITT,LINE)
          IF (NQ(1).EQ.7) GO TO 220
  210     WRITE (LUNOUT,FMT=6026) ITL,ITL
          WRITE (LUNOUT,FMT=6024) P(1,1,L2),P(1,2,L2),P(1,3,L2),
     +      S(1,1,L2),S(1,2,L2),S(1,3,L2)
          WRITE (LUNOUT,FMT=6024) P(2,1,L2),P(2,2,L2),P(2,3,L2),
     +      S(2,1,L2),S(2,2,L2),S(2,3,L2)
          WRITE (LUNOUT,FMT=6024) P(3,1,L2),P(3,2,L2),P(3,3,L2),
     +      S(3,1,L2),S(3,2,L2),S(3,3,L2)
          WRITE (LUNOUT,FMT=6074) DET(L2),BET(L2)
          IF (ICYCLE.LE.0) GO TO 250
  220     K4 = 4
          JB2 = ITL - 1
          DO 240 I = 1,3
            DO 230 J = 1,3
              P(I,J,JB2) = NN(I,J,2)
  230       CONTINUE
  240     CONTINUE
          MA = 2
          IF (NQ(1).EQ.7) L2 = 2
          IF (NQ(1).NE.7) THEN
            CALL MATOFC(MA,K4,ITL)
            LINE = LINE + 1
            IF ((LINE-5).GE.0) THEN
              LINE = 1
              WRITE (LUNOUT,FMT=6010) TITLE
            END IF
            WRITE (LUNOUT,FMT=6026) ISUM,ISUM
            WRITE (LUNOUT,FMT=6024) P(1,1,L2),P(1,2,L2),P(1,3,L2),
     +        S(1,1,L2),S(1,2,L2),S(1,3,L2)
            WRITE (LUNOUT,FMT=6024) P(2,1,L2),P(2,2,L2),P(2,3,L2),
     +        S(2,1,L2),S(2,2,L2),S(2,3,L2)
            WRITE (LUNOUT,FMT=6024) P(3,1,L2),P(3,2,L2),P(3,3,L2),
     +        S(3,1,L2),S(3,2,L2),S(3,3,L2)
            WRITE (LUNOUT,FMT=6074) DET(L2),BET(L2)
          END IF
          ITT = 3
          CALL TEST(ITT,LINE)
        END IF
C
  250   J = 0
  260   CONTINUE
        J = J + 1
        WRITE (LUNOUT,FMT=6010) TITLE
        IF ((IHKL(J)-2).EQ.0) THEN
          WRITE (LUNOUT,FMT=6076) J
        ELSE IF ((IHKL(J)-2).GT.0) THEN
          WRITE (LUNOUT,FMT=6072) J
        ELSE
          WRITE (LUNOUT,FMT=6028) J, (CELL(I,J),I=1,25)
        END IF
        A = DIS(1,J)
        B = DIS(2,J)
        C = DIS(3,J)
        ALPHA = DIS(4,J)
        BETA = DIS(5,J)
        GAMMA = DIS(6,J)
        IS = 1
        CALL LATCON(IS)
  270   CONTINUE
        WRITE (LUNOUT,FMT=6030) A,ALPHA
        WRITE (LUNOUT,FMT=6032) B,BETA
        WRITE (LUNOUT,FMT=6034) C,GAMMA
        WRITE (LUNOUT,FMT=6036) AS,ALSTAR
        WRITE (LUNOUT,FMT=6038) BS,BESTAR
        WRITE (LUNOUT,FMT=6040) CS,GASTAR
        WRITE (LUNOUT,FMT=6042) V
        WRITE (LUNOUT,FMT=6044) VS
C
        WRITE (LUNOUT,FMT=6046)
        WRITE (LUNOUT,FMT=6048)
        WRITE (LUNOUT,FMT=6050) ALPHA,SINAL,COSAL
        WRITE (LUNOUT,FMT=6052) BETA,SINBE,COSBE
        WRITE (LUNOUT,FMT=6054) GAMMA,SINGA,COSGA
        WRITE (LUNOUT,FMT=6056) ALSTAR,SINALS,COSALS
        WRITE (LUNOUT,FMT=6058) BESTAR,SINBES,COSBES
        WRITE (LUNOUT,FMT=6060) GASTAR,SINGAS,COSGAS
C
        IF (NQ(1).EQ.8) THEN
          GO TO 10
        ELSE
          I = IHKL(J)
          GO TO (290,280,280,10) I
C
  280     WRITE (LUNOUT,FMT=6062)
          AA = A*A
          BB = B*B
          CC = C*C
          AB = A*B*COSGA
          BC = B*C*COSAL
          CA = C*A*COSBE
          SUM = AB + BC + CA
          WRITE (LUNOUT,FMT=6064) AA,BB,CC
          WRITE (LUNOUT,FMT=6066) BC,CA,AB
          WRITE (LUNOUT,FMT=6068) SUM
C
  290     IF ((ISUM-J).GT.0) THEN
            GO TO 260
          ELSE IF (NQ(ITOTAL).EQ.8) THEN
            GO TO 10
          ELSE
            ITT = 4
            J = JSUM
            IHKL(J) = 4
            CALL TEST(ITT,LINE)
            IS = 1
            CALL LATCON(IS)
            GO TO 270
          END IF
        END IF
      END IF
      STOP
C
C---- Format statements
C
 6000 FORMAT (A)
 6008 FORMAT (I3,25A3)
 6010 FORMAT ('1',//5X,A,/)
 6012 FORMAT ('0',/4X,'Lattices used and generated by TRACER',/)
 6014 FORMAT (10X,'CELL (',I1,')',3X,25A3)
 6016 FORMAT (10X,'CELL (',I1,')   Reduced Cell')
 6018 FORMAT (10X,'CELL (',I1,')   Reduced Cell (Conventional Orientat',
     +       'ion)')
 6020 FORMAT ('0',//4X,'Transformation Matrices',/)
 6022 FORMAT ('0',/13X,'CELL (',I1,' To CELL (',I1,')',12X,'CELL (',I1,
     +       ' To CELL (',I1,')',//)
 6024 FORMAT (12X,3F7.2,12X,3F7.2,/)
 6026 FORMAT ('0',/13X,'CELL (1) To CELL (',I1,')',12X,'CELL (',I1,') ',
     +       'To CELL (1)',//)
 6028 FORMAT ('0',/10X,'Lattice Parameters of CELL (',I1,') --',25A3)
 6030 FORMAT ('0',/16X,'A  = ',F9.4,15X,'ALPHA  = ',F10.4)
 6032 FORMAT (16X,'B  = ',F9.4,15X,'BETA   = ',F10.4)
 6034 FORMAT (16X,'C  = ',F9.4,15X,'GAMMA  = ',F10.4)
 6036 FORMAT ('0',15X,'A* = ',F9.6,15X,'ALPHA* = ',F10.4)
 6038 FORMAT (16X,'B* = ',F9.6,15X,'BETA * = ',F10.4)
 6040 FORMAT (16X,'C* = ',F9.6,15X,'GAMMA* = ',F10.4)
 6042 FORMAT ('0',15X,'V  = ',F15.3)
 6044 FORMAT (16X,'V* = ',F15.6,//)
 6046 FORMAT ('0',/10X,'Natural trigonometric values of the ANGLES')
 6048 FORMAT ('0',/29X,'ANGLE',8X,'Sin',9X,'Cos ')
 6050 FORMAT ('0',15X,'ALPHA ',F13.4,2F12.5)
 6052 FORMAT (16X,'BETA  ',F13.4,2F12.5)
 6054 FORMAT (16X,'GAMMA ',F13.4,2F12.5)
 6056 FORMAT ('0',/16X,'ALPHA*',F13.4,2F12.5)
 6058 FORMAT (16X,'BETA *',F13.4,2F12.5)
 6060 FORMAT (16X,'GAMMA*',F13.4,2F12.5)
 6062 FORMAT ('0',//10X,'Scalars corresponding to this REDUCED CELL')
 6064 FORMAT ('0',/17X,'R(1,1) = ',F7.2,6X,'R(2,2) = ',F7.2,6X,'R(3,3)',
     +       ' = ',F7.2)
 6066 FORMAT ('0',16X,'R(2,3) = ',F7.2,6X,'R(3,1) = ',F7.2,6X,'R(1,2) ',
     +       '= ',F7.2)
 6068 FORMAT ('0',/17X,'R(2,3) + R(3,1) + R(1,2)  =  ',F7.2)
 6070 FORMAT (17X,3I4,21X,3I4,/)
 6072 FORMAT ('0',/11X,'Lattice parameters of CELL (',I1,') -- Reduced',
     +       ' cell (conventional orientation)')
 6074 FORMAT (18X,'MOD = ',F8.4,20X,'MOD = ',F8.4,/)
 6076 FORMAT ('0',/11X,'Lattice parameters of CELL (',I1,') -- Reduced',
     +       ' CELL')
 6080 FORMAT (1X,'ENTER FILENAME')
 6082 FORMAT (A10)
 6084 FORMAT (5X,'TRACER RUN ON',2A5,'AT',A5)
 6086 FORMAT (5X,'ACKNOWLEDGEMENTS')
 6088 FORMAT (5X,'PROGRAM TRACER FROM S.L. LAWTON')
 6090 FORMAT (5X,'J. APPLIED CRYSTALLOGRAPHY 6,309-346,(1973)')
C
C
      END
C
C
C
      SUBROUTINE RCELL
C     ================
C
C
C---- Transform given cell to reduced cell   called when instruction
C     nq = 7
C
C
C---- part 1   transformation
C
C
C
C     .. Scalars in Common ..
      REAL A,AL,ALPHA,ALS,ALSTAR,AS,B,BE,BES,BESTAR,BETA,BS,C,COSAL,
     +     COSALS,COSBE,COSBES,COSGA,COSGAS,CS,DEL,GA,GAMMA,GAS,GASTAR,
     +     SINAL,SINALS,SINBE,SINBES,SINGA,SINGAS,V,VS
      INTEGER ICYCLE,IOUTPT,ISUM,JSUM,K,L,L2,LUNIN,LUNOUT,MATRIX
      CHARACTER TITLE*70
C     ..
C     .. Arrays in Common ..
      REAL BET(15),DET(15),DIS(6,12),P(3,3,12),S(3,3,15)
      INTEGER IHKL(12),MM(3,3,4),NN(3,3,3)
C     ..
C     .. Local Scalars ..
      REAL AB,ABSAB,ABSBC,ABSCA,BC,CA,PN,QMAX,QQMAX,SAVE,X,XX,XY,XY2,Y,
     +     YY,YZ,Z,ZX,ZX2,ZZ
      INTEGER I,IAXIS,IS,ISIGS,J,JA,JB,JC,JF,JJ,JL,K2,K4,KK,LL,LM,LOC,M,
     +        MAXT,NCHECK,NSTOP,NTURN
C     ..
C     .. Local Arrays ..
      REAL ANGLE(3),CELL(25,12),D(3)
      INTEGER N(9),NA(9),NB(9),NQ(12),NTEMP(9)
C     ..
C     .. External Subroutines ..
      EXTERNAL LATCON,MATIN,STORE
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,ACOS,MAX,MIN,SIN,SQRT
C     ..
C     .. Common blocks ..
      COMMON /INOUT/LUNIN,LUNOUT
      COMMON /TCELL/A,B,C,ALPHA,BETA,GAMMA,COSAL,COSBE,COSGA,SINAL,
     +       SINBE,SINGA,ALS,BES,GAS,SINALS,SINBES,SINGAS,COSALS,COSBES,
     +       COSGAS,AS,BS,CS,AL,BE,GA,ALSTAR,BESTAR,GASTAR,V,VS
      COMMON /TCHAR/TITLE
      COMMON /WORK/DIS,K,IHKL,DEL,ISUM,JSUM,P,S,DET,BET,L,L2,MATRIX,NN,
     +       MM,ICYCLE,IOUTPT
C     ..
C
C
      K = K + 1
      IHKL(K) = 2
      ICYCLE = 1
      AB = A*B*COSGA
      BC = B*C*COSAL
      CA = C*A*COSBE
      DO 10 J = 1,9
        N(J) = 0
        NA(J) = 0
   10 CONTINUE
      N(1) = 1
      N(5) = 1
      N(9) = 1
      NA(1) = 1
      NA(5) = 1
      NA(9) = 1
   20 NSTOP = 1
      K2 = 0
      ABSAB = ABS(AB)
      ABSBC = ABS(BC)
      ABSCA = ABS(CA)
      QQMAX = MAX(ABSAB,ABSBC,ABSCA)
      IF ((QQMAX-ABSAB).EQ.0.0) GO TO 50
      IF ((QQMAX-ABSBC).EQ.0.0) THEN
        GO TO 60
      ELSE
        IF ((QQMAX-ABSBC).GT.0.0) GO TO 70
      END IF
      LOC = 7
      GO TO 1530
   50 QMAX = MAX(A,B)
      IF ((QMAX-A).LT.0.0) THEN
        GO TO 350
      ELSE IF ((QMAX-A).EQ.0.0) THEN
        GO TO 80
      ELSE
        GO TO 90
      END IF
   60 QMAX = MAX(B,C)
      IF ((QMAX-B).LT.0.0) THEN
        GO TO 350
      ELSE IF ((QMAX-B).EQ.0.0) THEN
        GO TO 100
      ELSE
        GO TO 110
      END IF
   70 QMAX = MAX(C,A)
      IF ((QMAX-C).LT.0.0) THEN
        GO TO 350
      ELSE IF ((QMAX-C).EQ.0.0) THEN
        GO TO 120
      ELSE
        GO TO 130
      END IF
   80 IAXIS = 2
      GO TO 140
   90 IAXIS = 1
      GO TO 210
  100 IAXIS = 3
      GO TO 210
  110 IAXIS = 2
      GO TO 280
  120 IAXIS = 1
      GO TO 280
  130 IAXIS = 3
      GO TO 140
C
  140 MAXT = 1
      GO TO 380
  150 IAXIS = 1
      QMAX = MAX(B,C)
      IF ((QMAX-B).LT.0.0) THEN
        GO TO 360
      ELSE IF ((QMAX-B).EQ.0.0) THEN
        GO TO 160
      ELSE
        GO TO 170
      END IF
  160 K2 = 1
      GO TO 450
  170 K2 = 2
      GO TO 500
  180 IF ((K2-1).LT.0) THEN
        GO TO 370
      ELSE IF ((K2-1).EQ.0) THEN
        GO TO 190
      ELSE
        GO TO 200
      END IF
  190 IAXIS = 2
      GO TO 500
  200 IAXIS = 3
  210 MAXT = 4
      GO TO 450
  220 IAXIS = 2
      QMAX = MAX(A,C)
      IF ((QMAX-A).LT.0.0) THEN
        GO TO 360
      ELSE IF ((QMAX-A).EQ.0.0) THEN
        GO TO 230
      ELSE
        GO TO 240
      END IF
  230 K2 = 1
      GO TO 380
  240 K2 = 2
      GO TO 500
  250 IF ((K2-1).LT.0) THEN
        GO TO 370
      ELSE IF ((K2-1).EQ.0) THEN
        GO TO 260
      ELSE
        GO TO 270
      END IF
  260 IAXIS = 1
      GO TO 500
  270 IAXIS = 3
      GO TO 380
C
  280 MAXT = 7
      GO TO 500
  290 IAXIS = 3
      QMAX = MAX(A,B)
      IF ((QMAX-A).LT.0.0) THEN
        GO TO 360
      ELSE IF ((QMAX-A).EQ.0.0) THEN
        GO TO 300
      ELSE
        GO TO 310
      END IF
  300 K2 = 1
      GO TO 380
  310 K2 = 2
      GO TO 450
  320 IF ((K2-1).LT.0) THEN
        GO TO 370
      ELSE IF ((K2-1).EQ.0) THEN
        GO TO 330
      ELSE
        GO TO 340
      END IF
  330 IAXIS = 1
      GO TO 450
  340 IAXIS = 2
      GO TO 380
  350 LOC = 1
      GO TO 1530
  360 LOC = 2
      GO TO 1530
  370 LOC = 3
      GO TO 1530
C
C---- reduction of a    (qmax = a)
C
  380 X = A
      Y = B
      Z = C
      XY = AB
      YZ = BC
      ZX = CA
C
      NCHECK = 1
      GO TO (550,390,400) IAXIS
  390 I = 2
      GO TO 570
  400 I = 3
      GO TO 720
C
  410 A = X
      B = Y
      C = Z
      AB = XY
      BC = YZ
      CA = ZX
C
      NCHECK = NCHECK - 1
      IF (NCHECK.LT.0) GO TO 430
      IF ((I-2).LT.0) THEN
        GO TO 560
      ELSE IF ((I-2).EQ.0) THEN
        GO TO 400
      ELSE
        GO TO 390
      END IF
  430 MAXT = MAXT + 1
      NSTOP = NSTOP + 1
      IF ((NSTOP-3).GT.0) GO TO 920
      GO TO (140,150,180,210,220,250,280,290,320) MAXT
C
C---- reduction of b    (qmax = b)
C
  450 X = B
      Y = C
      Z = A
      XY = BC
      YZ = CA
      ZX = AB
C
      NCHECK = 1
      GO TO (470,550,460) IAXIS
  460 I = 6
      GO TO 570
  470 I = 4
      GO TO 720
C
  480 A = Z
      B = X
      C = Y
      AB = ZX
      BC = XY
      CA = YZ
C
      NCHECK = NCHECK - 1
      IF (NCHECK.LT.0) GO TO 430
      IF ((I-4).LT.0) THEN
        GO TO 560
      ELSE IF ((I-4).EQ.0) THEN
        GO TO 460
      ELSE
        GO TO 470
      END IF
C
C---- Reduction of c    (qmax = c)
C
  500 X = C
      Y = A
      Z = B
      XY = CA
      YZ = AB
      ZX = BC
C
      NCHECK = 1
      GO TO (510,520,550) IAXIS
  510 I = 7
      GO TO 570
  520 I = 8
      GO TO 720
C
  530 A = Y
      B = Z
      C = X
      AB = YZ
      BC = ZX
      CA = XY
C
      NCHECK = NCHECK - 1
      IF (NCHECK.LT.0) GO TO 430
      IF ((I-7).LT.0) THEN
        GO TO 560
      ELSE IF ((I-7).EQ.0) THEN
        GO TO 520
      ELSE
        GO TO 510
      END IF
C
  550 LOC = 4
      GO TO 1530
  560 LOC = 5
      GO TO 1530
C
C---- reduction of x   absf(x y) greater than (y**2)/2
C
  570 XX = X*X
      YY = Y*Y
      PN = 0.0
      IF ((ABS(XY)-0.5*YY).LE.0.0) GO TO 910
  580 IF (XY.LT.0.0) THEN
        GO TO 590
      ELSE IF (XY.EQ.0.0) THEN
        GO TO 910
      ELSE
        GO TO 670
      END IF
  590 XY2 = PN*YY + XY
      IF (XY2.GE.0.0) GO TO 610
      SAVE = XY2
      PN = PN + 1.0
      GO TO 580
  610 IF ((ABS(XY2)-ABS(SAVE)).LT.0.0) GO TO 630
      PN = PN - 1.0
  630 NA(I) = PN
      X = SQRT(2.0*PN*XY+XX+PN*PN*YY)
      ZX = PN*YZ + ZX
  640 IF ((ABS(XY2)-ABS(SAVE)).GE.0.0) GO TO 660
      XY = XY2
      GO TO 870
  660 XY = SAVE
      GO TO 870
C
  670 XY2 = XY - PN*YY
      IF (XY2.LE.0.0) GO TO 690
      SAVE = XY2
      PN = PN + 1.0
      GO TO 670
  690 IF ((ABS(XY2)-ABS(SAVE)).LE.0.0) GO TO 710
      PN = PN - 1.0
  710 NA(I) = -PN
      X = SQRT(XX-2.0*PN*XY+PN*PN*YY)
      ZX = ZX - PN*YZ
      GO TO 640
C
C---- Reduction of x   absf(z x) greater than (z**2)/2
C
  720 XX = X*X
      ZZ = Z*Z
      PN = 0.0
      IF ((ABS(ZX)-0.5*ZZ).LE.0.0) GO TO 910
  730 IF (ZX.LT.0.0) THEN
        GO TO 740
      ELSE IF (ZX.EQ.0.0) THEN
        GO TO 910
      ELSE
        GO TO 820
      END IF
  740 ZX2 = PN*ZZ + ZX
      IF (ZX2.GE.0.0) GO TO 760
      SAVE = ZX2
      PN = PN + 1.0
      GO TO 730
  760 IF ((ABS(ZX2)-ABS(SAVE)).LT.0.0) GO TO 780
      PN = PN - 1.0
  780 NA(I) = PN
      X = SQRT(2.0*PN*ZX+XX+PN*PN*ZZ)
      XY = PN*YZ + XY
  790 IF ((ABS(ZX2)-ABS(SAVE)).GE.0.0) GO TO 810
      ZX = ZX2
      GO TO 870
  810 ZX = SAVE
      GO TO 870
C
  820 ZX2 = ZX - PN*ZZ
      IF (ZX2.LE.0.0) GO TO 840
      SAVE = ZX2
      PN = PN + 1.0
      GO TO 820
  840 IF ((ABS(ZX2)-ABS(SAVE)).LE.0.0) GO TO 860
      PN = PN - 1.0
  860 NA(I) = -PN
      X = SQRT(XX-2.0*PN*ZX+PN*PN*ZZ)
      XY = XY - PN*YZ
      GO TO 790
C
  870 NB(1) = NA(1)*N(1) + NA(2)*N(4) + NA(3)*N(7)
      NB(2) = NA(1)*N(2) + NA(2)*N(5) + NA(3)*N(8)
      NB(3) = NA(1)*N(3) + NA(2)*N(6) + NA(3)*N(9)
      NB(4) = NA(4)*N(1) + NA(5)*N(4) + NA(6)*N(7)
      NB(5) = NA(4)*N(2) + NA(5)*N(5) + NA(6)*N(8)
      NB(6) = NA(4)*N(3) + NA(5)*N(6) + NA(6)*N(9)
      NB(7) = NA(7)*N(1) + NA(8)*N(4) + NA(9)*N(7)
      NB(8) = NA(7)*N(2) + NA(8)*N(5) + NA(9)*N(8)
      NB(9) = NA(7)*N(3) + NA(8)*N(6) + NA(9)*N(9)
      DO 900 J = 1,9
        N(J) = NB(J)
        NA(J) = 0
  900 CONTINUE
      NA(1) = 1
      NA(5) = 1
      NA(9) = 1
  910 GO TO (920,410,410,480,920,480,530,530,920) I
C
C---- output of part 1     (reduced cell)
C
  920 COSAL = BC/ (B*C)
      COSBE = CA/ (C*A)
      COSGA = AB/ (A*B)
      AL = ACOS(COSAL)
      BE = ACOS(COSBE)
      GA = ACOS(COSGA)
      SINAL = SIN(AL)
      SINBE = SIN(BE)
      SINGA = SIN(GA)
      ALPHA = 57.295780*AL
      BETA = 57.295780*BE
      GAMMA = 57.295780*GA
      LL = 1
      GO TO 1390
C
C---- A test is now made to determine if the cell just obtained
C     is the reduced cell   if one or more of the axes can be further
C     reduced, the program returns to stmnt 382 for another cycle
C
  930 IF ((ABS(AB)-0.5*A*A).GT.0.0) GO TO 990
      IF ((ABS(AB)-0.5*B*B).GT.0.0) GO TO 990
      IF ((ABS(BC)-0.5*B*B).GT.0.0) GO TO 990
      IF ((ABS(BC)-0.5*C*C).GT.0.0) GO TO 990
      IF ((ABS(CA)-0.5*C*C).GT.0.0) GO TO 990
      IF ((ABS(CA)-0.5*A*A).LE.0.0) GO TO 1000
  990 ICYCLE = ICYCLE + 1
      GO TO 20
C
C
C---- Part 2    transformation of reduced cell to conventional
C     reduced cell   the convention is c less than a less than b,
C     with alpha, beta obtuse
C
C
 1000 JA = 0
      JB = 0
      JC = 0
      K = K + 1
      IHKL(K) = 3
      NTURN = 1
      KK = 3
      IF ((MIN(A,B,C)-A).LT.0.0) THEN
        GO TO 1020
      ELSE IF ((MIN(A,B,C)-A).EQ.0.0) THEN
        GO TO 1010
      ELSE
        GO TO 1180
      END IF
 1010 JJ = 7
      M = -6
      GO TO 1190
 1020 IF ((MIN(A,B,C)-B).LT.0.0) THEN
        GO TO 1030
      ELSE IF ((MIN(A,B,C)-B).EQ.0.0) THEN
        GO TO 1040
      ELSE
        GO TO 1180
      END IF
 1030 JC = 1
      JJ = 7
      M = 0
      GO TO 1230
 1040 JJ = 7
      M = -3
      GO TO 1220
 1050 NTURN = 2
      KK = 2
      IF ((MAX(A,B,C)-A).LT.0.0) THEN
        GO TO 1180
      ELSE IF ((MAX(A,B,C)-A).EQ.0.0) THEN
        GO TO 1060
      ELSE
        GO TO 1070
      END IF
 1060 IF (MIN(A,B,C).EQ.A) GO TO 1070
      JJ = 4
      M = -3
      GO TO 1190
 1070 IF ((MAX(A,B,C)-B).LT.0.0) THEN
        GO TO 1180
      ELSE IF ((MAX(A,B,C)-B).EQ.0.0) THEN
        GO TO 1080
      ELSE
        GO TO 1090
      END IF
 1080 JB = 1
      JJ = 4
      M = 0
      GO TO 1220
 1090 JJ = 4
      M = 3
      GO TO 1230
 1100 NTURN = 3
      KK = 1
      IF ((D(2)-A).EQ.0.0) GO TO 1120
      IF ((D(2)-B).LT.0.0) THEN
        GO TO 1180
      ELSE IF ((D(2)-B).EQ.0.0) THEN
        GO TO 1150
      ELSE
        GO TO 1170
      END IF
 1120 IF ((MIN(B,C)-B).LT.0.0) THEN
        GO TO 1140
      ELSE IF ((MIN(B,C)-B).EQ.0.0) THEN
        GO TO 1130
      ELSE
        GO TO 1180
      END IF
 1130 JJ = 1
      M = 6
      GO TO 1230
 1140 JJ = 1
      M = 3
      GO TO 1220
 1150 IF ((MIN(A,C)-A).LT.0.0) THEN
        GO TO 1160
      ELSE IF ((MIN(A,C)-A).EQ.0.0) THEN
        GO TO 1130
      ELSE
        GO TO 1180
      END IF
 1160 JA = 1
      JJ = 1
      M = 0
      GO TO 1190
 1170 IF ((MIN(A,B)-A).LT.0.0) THEN
        GO TO 1160
      ELSE IF ((MIN(A,B)-A).EQ.0.0) THEN
        GO TO 1140
      ELSE
        GO TO 1180
      END IF
C
 1180 LOC = 6
      GO TO 1530
C
 1190 D(KK) = A
      ANGLE(KK) = ALPHA
 1200 JL = JJ + 2
      DO 1210 L = JJ,JL
        LM = L + M
        NTEMP(L) = N(LM)
 1210 CONTINUE
      GO TO (1050,1100,1240) NTURN
 1220 D(KK) = B
      ANGLE(KK) = BETA
      GO TO 1200
 1230 D(KK) = C
      ANGLE(KK) = GAMMA
      GO TO 1200
C
 1240 DO 1250 I = 1,9
        N(I) = NTEMP(I)
 1250 CONTINUE
      A = D(1)
      B = D(2)
      C = D(3)
      ALPHA = ANGLE(1)
      BETA = ANGLE(2)
      GAMMA = ANGLE(3)
C
C---- The following section checks to make sure alpha and beta are
C     each equal to or greater than 90 degrees   if one or both are
C     less than 90 degrees, the proper adjustments are made
C
      LL = 2
      ISIGS = 0
      IF ((ALPHA-90.0).GE.0.0) GO TO 1280
      ISIGS = ISIGS + 1
 1280 IF ((BETA-90.0).GE.0.0) GO TO 1300
      ISIGS = ISIGS + 1
 1300 IF ((ISIGS-1).LT.0) THEN
        GO TO 1390
      ELSE IF ((ISIGS-1).EQ.0) THEN
        GO TO 1310
      ELSE
        GO TO 1370
      END IF
C
C---- isigs = 1
C     either alpha or beta is less than 90 degrees   redefine that
C     angle as its supplement and set gamma = 180 0 - gamma
C     (two angles must always be transformed simultaneously )
C
 1310 GAMMA = 180.0 - GAMMA
      DO 1320 I = 7,9
        N(I) = -N(I)
 1320 CONTINUE
      IF ((ALPHA-90.0).GE.0.0) GO TO 1350
      ALPHA = 180.0 - ALPHA
      DO 1340 I = 1,3
        N(I) = -N(I)
 1340 CONTINUE
      GO TO 1390
 1350 BETA = 180.0 - BETA
      DO 1360 I = 4,6
        N(I) = -N(I)
 1360 CONTINUE
      GO TO 1390
C
C---- isigs = 2
C     redefine alpha and beta as its own supplement and reverse the
C     directions of a and b
C
 1370 ALPHA = 180.0 - ALPHA
      BETA = 180.0 - BETA
      DO 1380 I = 1,6
        N(I) = -N(I)
 1380 CONTINUE
C
 1390 GO TO (1450,1400) LL
C
C---- A check is made to determine if the system has gone from
C     right-handed to left-handed or vice versa   if jf = 1 the system
C     has reversed   if jf = 0 or 3, it has not reversed
C
 1400 JF = JA + JB + JC
      IF ((JF-1).NE.0) GO TO 1450
      DO 1440 I = 1,9
        N(I) = -N(I)
 1440 CONTINUE
C
 1450 K = K - 1
      L = 0
      DO 1470 I = 1,3
        DO 1460 J = 1,3
          L = L + 1
          NN(I,J,LL) = N(L)
          P(I,J,K) = N(L)
 1460   CONTINUE
 1470 CONTINUE
      K4 = 2
      CALL MATIN(K4)
      L = 0
      DO 1490 I = 1,3
        DO 1480 J = 1,3
          L = L + 1
          MM(I,J,LL) = S(I,J,K)
 1480   CONTINUE
 1490 CONTINUE
      K = K + 1
      IF ((IHKL(K)-2).LE.0) GO TO 1510
      IS = 1
      GO TO 1520
 1510 IS = 2
 1520 CALL LATCON(IS)
      CALL STORE
      I = IHKL(K)
      IF (I.EQ.2) GO TO 930
      GO TO 1540
C
 1530 WRITE (LUNOUT,FMT=6000) LOC
 1540 RETURN
C
C---- Format statements
C
 6000 FORMAT ('0',//'An ERROR has occurred in subroutine RCELL, LOC = ',
     +       I2)
C
C
      END
C
C
C
      SUBROUTINE LATCON(IS)
C     =====================
C
C---- Calculate cell parameter data in direct and reciprocal space
C
C
C     .. Scalar Arguments ..
      INTEGER IS
C     ..
C     .. Scalars in Common ..
      REAL A,AL,ALPHA,ALS,ALSTAR,AS,B,BE,BES,BESTAR,BETA,BS,C,COSAL,
     +     COSALS,COSBE,COSBES,COSGA,COSGAS,CS,DEL,GA,GAMMA,GAS,GASTAR,
     +     SINAL,SINALS,SINBE,SINBES,SINGA,SINGAS,V,VS
      INTEGER ICYCLE,IOUTPT,ISUM,JSUM,K,L,L2,LUNIN,LUNOUT,MATRIX
      CHARACTER TITLE*70
C     ..
C     .. Arrays in Common ..
      REAL BET(15),DET(15),DIS(6,12),P(3,3,12),S(3,3,15)
      INTEGER IHKL(12),MM(3,3,4),NN(3,3,3)
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ACOS,COS,SIN
C     ..
C     .. Common blocks ..
      COMMON /INOUT/LUNIN,LUNOUT
      COMMON /TCELL/A,B,C,ALPHA,BETA,GAMMA,COSAL,COSBE,COSGA,SINAL,
     +       SINBE,SINGA,ALS,BES,GAS,SINALS,SINBES,SINGAS,COSALS,COSBES,
     +       COSGAS,AS,BS,CS,AL,BE,GA,ALSTAR,BESTAR,GASTAR,V,VS
      COMMON /TCHAR/TITLE
      COMMON /WORK/DIS,K,IHKL,DEL,ISUM,JSUM,P,S,DET,BET,L,L2,MATRIX,NN,
     +       MM,ICYCLE,IOUTPT
C     ..
C
C
      GO TO (10,30,20) IS
   10 SINAL = SIN(ALPHA*1.745329E-02)
      SINBE = SIN(BETA*1.745329E-02)
      SINGA = SIN(GAMMA*1.745329E-02)
      COSAL = COS(ALPHA*1.745329E-02)
      COSBE = COS(BETA*1.745329E-02)
      COSGA = COS(GAMMA*1.745329E-02)
      GO TO 30
C
   20 SINALS = SIN(ALSTAR*1.745329E-02)
      SINBES = SIN(BESTAR*1.745329E-02)
      SINGAS = SIN(GASTAR*1.745329E-02)
      COSALS = COS(ALSTAR*1.745329E-02)
      COSBES = COS(BESTAR*1.745329E-02)
      COSGAS = COS(GASTAR*1.745329E-02)
      AL = ACOS((COSBES*COSGAS-COSALS)/ (SINBES*SINGAS))
      BE = ACOS((COSALS*COSGAS-COSBES)/ (SINALS*SINGAS))
      GA = ACOS((COSALS*COSBES-COSGAS)/ (SINALS*SINBES))
      SINAL = SIN(AL)
      SINBE = SIN(BE)
      SINGA = SIN(GA)
      COSAL = COS(AL)
      COSBE = COS(BE)
      COSGA = COS(GA)
      A = 1.0/ (AS*SINBE*SINGAS)
      B = 1.0/ (BS*SINAL*SINGAS)
      C = 1.0/ (CS*SINAL*SINBES)
      ALPHA = AL*57.295780
      BETA = BE*57.295780
      GAMMA = GA*57.295780
      GO TO 40
   30 ALS = ACOS((COSBE*COSGA-COSAL)/ (SINBE*SINGA))
      BES = ACOS((COSAL*COSGA-COSBE)/ (SINAL*SINGA))
      GAS = ACOS((COSAL*COSBE-COSGA)/ (SINAL*SINBE))
      SINALS = SIN(ALS)
      SINBES = SIN(BES)
      SINGAS = SIN(GAS)
      COSALS = COS(ALS)
      COSBES = COS(BES)
      COSGAS = COS(GAS)
      AS = 1.0/ (A*SINBES*SINGA)
      BS = 1.0/ (B*SINALS*SINGA)
      CS = 1.0/ (C*SINALS*SINBE)
      ALSTAR = ALS*57.295780
      BESTAR = BES*57.295780
      GASTAR = GAS*57.295780
C
   40 V = A*B*C*SINAL*SINBE*SINGAS
      VS = 1.0/V
C
C
      END
C
C
C
      SUBROUTINE STORE
C     ================
C
C---- Store cell parameters into variable dis(i,k)
C
C
C     .. Scalars in Common ..
      REAL A,AL,ALPHA,ALS,ALSTAR,AS,B,BE,BES,BESTAR,BETA,BS,C,COSAL,
     +     COSALS,COSBE,COSBES,COSGA,COSGAS,CS,DEL,GA,GAMMA,GAS,GASTAR,
     +     SINAL,SINALS,SINBE,SINBES,SINGA,SINGAS,V,VS
      INTEGER ICYCLE,IOUTPT,ISUM,JSUM,K,L,L2,LUNIN,LUNOUT,MATRIX
      CHARACTER TITLE*70
C     ..
C     .. Arrays in Common ..
      REAL BET(15),DET(15),DIS(6,12),P(3,3,12),S(3,3,15)
      INTEGER IHKL(12),MM(3,3,4),NN(3,3,3)
C     ..
C     .. Common blocks ..
      COMMON /INOUT/LUNIN,LUNOUT
      COMMON /TCELL/A,B,C,ALPHA,BETA,GAMMA,COSAL,COSBE,COSGA,SINAL,
     +       SINBE,SINGA,ALS,BES,GAS,SINALS,SINBES,SINGAS,COSALS,COSBES,
     +       COSGAS,AS,BS,CS,AL,BE,GA,ALSTAR,BESTAR,GASTAR,V,VS
      COMMON /TCHAR/TITLE
      COMMON /WORK/DIS,K,IHKL,DEL,ISUM,JSUM,P,S,DET,BET,L,L2,MATRIX,NN,
     +       MM,ICYCLE,IOUTPT
C     ..            
C
C
      DIS(1,K) = A
      DIS(2,K) = B
      DIS(3,K) = C
      DIS(4,K) = ALPHA
      DIS(5,K) = BETA
      DIS(6,K) = GAMMA
C
C
      END
C
C
C
      SUBROUTINE DIMEN
C     ================
C
C---- Transform given cell (defined by vectors in direct space) to a
C     new cell using a vectorial transformation matrix p(i,j,k)
C     called if nq le 6
C
C
C---- calculation of the new a, b, c
C
C     .. Scalars in Common ..
      REAL A,AL,ALPHA,ALS,ALSTAR,AS,B,BE,BES,BESTAR,BETA,BS,C,COSAL,
     +     COSALS,COSBE,COSBES,COSGA,COSGAS,CS,DEL,GA,GAMMA,GAS,GASTAR,
     +     SINAL,SINALS,SINBE,SINBES,SINGA,SINGAS,V,VS
      INTEGER ICYCLE,IOUTPT,ISUM,JSUM,K,L,L2,LUNIN,LUNOUT,MATRIX
      CHARACTER TITLE*70
C     ..
C     .. Arrays in Common ..
      REAL BET(15),DET(15),DIS(6,12),P(3,3,12),S(3,3,15)
      INTEGER IHKL(12),MM(3,3,4),NN(3,3,3)
C     ..
C     .. Local Scalars ..
      INTEGER I,J
C     ..
C     .. Local Arrays ..
      REAL ANGLE(3),AXIS(3),CSE(3)
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ACOS,SQRT
C     ..            
C     .. Common blocks ..
      COMMON /INOUT/LUNIN,LUNOUT
      COMMON /TCELL/A,B,C,ALPHA,BETA,GAMMA,COSAL,COSBE,COSGA,SINAL,
     +       SINBE,SINGA,ALS,BES,GAS,SINALS,SINBES,SINGAS,COSALS,COSBES,
     +       COSGAS,AS,BS,CS,AL,BE,GA,ALSTAR,BESTAR,GASTAR,V,VS
      COMMON /TCHAR/TITLE
      COMMON /WORK/DIS,K,IHKL,DEL,ISUM,JSUM,P,S,DET,BET,L,L2,MATRIX,NN,
     +       MM,ICYCLE,IOUTPT
C     ..
C
C
      DO 10 I = 1,3
        AXIS(I) = SQRT((P(I,1,K)*A)**2+ (P(I,2,K)*B)**2+
     +            (P(I,3,K)*C)**2+ (P(I,1,K)*P(I,2,K)*A*B*COSGA+P(I,1,
     +            K)*P(I,3,K)*A*C*COSBE+P(I,2,K)*P(I,3,K)*B*C*COSAL)*
     +            2.0)
   10 CONTINUE
C
C---- calculation of the new alpha, beta, gamma
C
      J = 2
      L = 3
      DO 40 I = 1,3
        CSE(I) = (P(J,1,K)*P(L,1,K)*A*A+P(J,2,K)*P(L,2,K)*B*B+
     +           P(J,3,K)*P(L,3,K)*C*C+ (P(J,1,K)*P(L,2,K)+P(J,2,K)*P(L,
     +           1,K))*A*B*COSGA+ (P(J,1,K)*P(L,3,K)+P(J,3,K)*P(L,1,K))*
     +           A*C*COSBE+ (P(J,2,K)*P(L,3,K)+P(J,3,K)*P(L,2,K))*B*C*
     +           COSAL)/ (AXIS(J)*AXIS(L))
        ANGLE(I) = ACOS(CSE(I))*57.295780
        GO TO (20,30,40) I
   20   J = 1
        GO TO 40
   30   L = 2
   40 CONTINUE
      A = AXIS(1)
      B = AXIS(2)
      C = AXIS(3)
      ALPHA = ANGLE(1)
      BETA = ANGLE(2)
      GAMMA = ANGLE(3)
C
      K = K + 1
      IHKL(K) = 1
C
C
      END
C
C
C
      SUBROUTINE MATIN(K4)
C     ====================
C
C---- Generate matrix elements of inverse transformation matrix
C
C
C     .. Scalar Arguments ..
      INTEGER K4
C     ..
C     .. Scalars in Common ..
      REAL A,AL,ALPHA,ALS,ALSTAR,AS,B,BE,BES,BESTAR,BETA,BS,C,COSAL,
     +     COSALS,COSBE,COSBES,COSGA,COSGAS,CS,DEL,GA,GAMMA,GAS,GASTAR,
     +     SINAL,SINALS,SINBE,SINBES,SINGA,SINGAS,V,VS
      INTEGER ICYCLE,IOUTPT,ISUM,JSUM,K,L,L2,LUNIN,LUNOUT,MATRIX
      CHARACTER TITLE*70
C     ..
C     .. Arrays in Common ..
      REAL BET(15),DET(15),DIS(6,12),P(3,3,12),S(3,3,15)
      INTEGER IHKL(12),MM(3,3,4),NN(3,3,3)
C     ..
C     .. Local Scalars ..
      INTEGER I,J
C     ..
C     .. Common blocks ..
      COMMON /INOUT/LUNIN,LUNOUT
      COMMON /TCELL/A,B,C,ALPHA,BETA,GAMMA,COSAL,COSBE,COSGA,SINAL,
     +       SINBE,SINGA,ALS,BES,GAS,SINALS,SINBES,SINGAS,COSALS,COSBES,
     +       COSGAS,AS,BS,CS,AL,BE,GA,ALSTAR,BESTAR,GASTAR,V,VS
      COMMON /TCHAR/TITLE
      COMMON /WORK/DIS,K,IHKL,DEL,ISUM,JSUM,P,S,DET,BET,L,L2,MATRIX,NN,
     +       MM,ICYCLE,IOUTPT
C     ..
C
C
      GO TO (10,10,20,20) K4
   10 L = K
      GO TO 30
   20 L = L2
   30 CONTINUE
      DET(L) = (P(2,2,L)*P(3,3,L)-P(3,2,L)*P(2,3,L))*P(1,1,L) -
     +         (P(2,1,L)*P(3,3,L)-P(3,1,L)*P(2,3,L))*P(1,2,L) +
     +         (P(2,1,L)*P(3,2,L)-P(3,1,L)*P(2,2,L))*P(1,3,L)
      S(1,1,L) = (P(2,2,L)*P(3,3,L)-P(3,2,L)*P(2,3,L))/DET(L)
      S(1,2,L) = - (P(1,2,L)*P(3,3,L)-P(3,2,L)*P(1,3,L))/DET(L)
      S(1,3,L) = (P(1,2,L)*P(2,3,L)-P(2,2,L)*P(1,3,L))/DET(L)
      S(2,1,L) = - (P(2,1,L)*P(3,3,L)-P(3,1,L)*P(2,3,L))/DET(L)
      S(2,2,L) = (P(1,1,L)*P(3,3,L)-P(3,1,L)*P(1,3,L))/DET(L)
      S(2,3,L) = - (P(1,1,L)*P(2,3,L)-P(2,1,L)*P(1,3,L))/DET(L)
      S(3,1,L) = (P(2,1,L)*P(3,2,L)-P(3,1,L)*P(2,2,L))/DET(L)
      S(3,2,L) = - (P(1,1,L)*P(3,2,L)-P(3,1,L)*P(1,2,L))/DET(L)
      S(3,3,L) = (P(1,1,L)*P(2,2,L)-P(2,1,L)*P(1,2,L))/DET(L)
C
C
      BET(L) = 1.0/DET(L)
      IF ((K4-1).EQ.0) THEN
        IF ((L-MATRIX).LT.0) THEN
          L = L + 1
          GO TO 30
        END IF
      END IF
      RETURN
C
C
      END
C
C
C
      SUBROUTINE MATOFC(MA,K4,ITL)
C     ===========================
C
C     generate matrix for vectorial transformation of lattice axes of
C     original cell to final cell when both cells are separated by
C     intermediate cells   called only if 3 or more cells are involved
C     in a sequence of transformations   if a series of instructions
C     ends with nq = 8, matofc called only once   if the series ends
C     with nq = 7, matofc called twice, once for each of the two
C     orientations of the reduced cell   matrix obtained by matrix mult
C
C
C     .. Scalar Arguments ..
      INTEGER ITL,K4,MA
C     ..
C     .. Scalars in Common ..
      REAL A,AL,ALPHA,ALS,ALSTAR,AS,B,BE,BES,BESTAR,BETA,BS,C,COSAL,
     +     COSALS,COSBE,COSBES,COSGA,COSGAS,CS,DEL,GA,GAMMA,GAS,GASTAR,
     +     SINAL,SINALS,SINBE,SINBES,SINGA,SINGAS,V,VS
      INTEGER ICYCLE,IOUTPT,ISUM,JSUM,K,L,L2,LUNIN,LUNOUT,MATRIX
      CHARACTER TITLE*70
C     ..
C     .. Arrays in Common ..
      REAL BET(15),DET(15),DIS(6,12),P(3,3,12),S(3,3,15)
      INTEGER IHKL(12),MM(3,3,4),NN(3,3,3)
C     ..
C     .. Local Scalars ..
      INTEGER I,J,K8,LL,M
C     ..
C     .. Local Arrays ..
      REAL Q(3,3,12),W(3,3,12)
C     ..
C     .. External Subroutines ..
      EXTERNAL MATIN
C     ..            
C     .. Common blocks ..
      COMMON /INOUT/LUNIN,LUNOUT
      COMMON /TCELL/A,B,C,ALPHA,BETA,GAMMA,COSAL,COSBE,COSGA,SINAL,
     +       SINBE,SINGA,ALS,BES,GAS,SINALS,SINBES,SINGAS,COSALS,COSBES,
     +       COSGAS,AS,BS,CS,AL,BE,GA,ALSTAR,BESTAR,GASTAR,V,VS
      COMMON /TCHAR/TITLE
      COMMON /WORK/DIS,K,IHKL,DEL,ISUM,JSUM,P,S,DET,BET,L,L2,MATRIX,NN,
     +       MM,ICYCLE,IOUTPT
C     ..
C
C
      GO TO (10,20) MA
   10 IF ((ITL-3).LT.0) THEN
        GO TO 130
      ELSE
        K4 = 3
        K = K + 1
        IHKL(K) = 4
      END IF
   20 DO 50 L = 1,ISUM
        DO 40 J = 1,3
          DO 30 I = 1,3
            Q(I,J,L) = P(I,J,L)
   30     CONTINUE
   40   CONTINUE
   50 CONTINUE
      IF ((K4-4).GE.0) THEN
        K = K8
      ELSE
        K8 = K
      END IF
      IF (ICYCLE.LE.0) K = K + 1
      L = K - 3
      M = K - 4
      LL = L
   60 CONTINUE
      DO 80 I = 1,3
        DO 70 J = 1,3
          W(I,J,L) = Q(I,1,L)*Q(1,J,M) + Q(I,2,L)*Q(2,J,M) +
     +               Q(I,3,L)*Q(3,J,M)
   70   CONTINUE
   80 CONTINUE
      IF ((M-1).GT.0) THEN
        DO 100 I = 1,3
          DO 90 J = 1,3
            Q(I,J,L) = W(I,J,L)
   90     CONTINUE
  100   CONTINUE
        M = M - 1
        GO TO 60
      END IF
      IF (K4.EQ.3) THEN
        L2 = ISUM
      ELSE
        L2 = L2 + 1
      END IF
      DO 120 I = 1,3
        DO 110 J = 1,3
          P(I,J,L2) = W(I,J,L)
  110   CONTINUE
  120 CONTINUE
      CALL MATIN(K4)
  130 RETURN
C
C
      END
C
C
C
      SUBROUTINE CENTER(K,NQ,P)
C     =========================
C
C---- Define matrix elements for transformation of an a-, b-, c-, i- and
C     f-centered cell to a primitive cell
C
C
C     .. Scalar Arguments ..
      INTEGER K
C     ..
C     .. Array Arguments ..
      REAL P(3,3,12)
      INTEGER NQ(12)
C     ..
C     .. Local Scalars ..
      INTEGER I,J,L
C     ..
C     .. Local Arrays ..
      REAL E(9)
C     ..
      J = NQ(K)
      DO 140 I = 1,9
        GO TO (10,20,30,40,50,60,70,80,90) I
C
   10   GO TO (170,100,120,110,100,110) J
   20   GO TO (170,120,120,110,120,110) J
   30   GO TO (170,120,130,120,120,120) J
   40   GO TO (170,120,120,130,120,120) J
   50   GO TO (170,110,100,120,100,110) J
   60   GO TO (170,110,120,120,120,110) J
   70   GO TO (170,120,110,120,110,110) J
   80   GO TO (170,130,120,120,110,120) J
   90   GO TO (170,120,110,100,110,110) J
        GO TO 130
  100   E(I) = 1.0
        GO TO 140
  110   E(I) = 0.5
        GO TO 140
C    2 E(I) = -0.5                                                   
C      GO TO 100                                                     
  120   E(I) = 0.0
        GO TO 140
C
  130   E(I) = -1.0
  140 CONTINUE
C
      L = 0
      DO 160 I = 1,3
        DO 150 J = 1,3
          L = L + 1
          P(I,J,K) = E(L)
  150   CONTINUE
  160 CONTINUE
  170 RETURN
C
      END
C
C
C
      SUBROUTINE TEST(ITT,LINE)
C     ========================
C
C   !!!!! this subroutine cannot be toolpacked (too many cases)
C
C--- Analyze scalar products of reduced cell (in conventional orienta-
C     tion c lt a lt b, alpha, beta obtuse) for unit cell of highest    
C     symmetry   matrix elements eu(i) taken from table 3 of - the
C     reduced cell and its crystallographic applications, u s a e c ,
C     ames lab report is-1141, april 1965, by lawton and jacobson
C
C---- Glossary of symbols   remaining symbols defined in glossary
C     preceding main program
C
C     A2       CELL PARAMETER, A, IN NEW UNIT CELL
C     ABC      A*B*COS(GAMMA)
C     ACC      A*C*COS(BETA)
C     ALPHA2   CELL PARAMETER, ALPHA, IN NEW UNIT CELL
C     B2       CELL PARAMETER, B, IN NEW UNIT CELL
C     BCC      B*C*COS(ALPHA)
C     BETA2    CELL PARAMETER, BETA, IN NEW UNIT CELL
C     C2       CELL PARAMETER, C, IN NEW UNIT CELL
C     DEL      LARGEST EST  ERROR BETWEEN REDUCED CELL SCALAR PRODUCTS
C     DETEU    MODULUS OF MATRIX EU
C     DETEUI   MODULUS OF MATRIX EUI
C     DETFN    MODULUS OF MATRIX F
C     DETFNI   MODULUS OF MATRIX FI
C     EU(9)    MATRIX ELEM  FOR TRANSF  OF CONV  RED  CELL TO UNIT CELL
C     EUI(9)   INVERSE MATRIX ELEMENTS OF EU
C     F(9)     MATRIX ELEM  FOR TRANSF  OF ORIG  CELL TO UNIT CELL
C     FI(9)    INVERSE MATRIX ELEMENTS OF F
C     GAMMA2   CELL PARAMETER, GAMMA, IN NEW UNIT CELL
C     ISYM     CODE IDENTIFYING EQUIVALENCY AMONG SYMMETRICAL SCALARS
C     JSUM     ISUM + 1
C     JU       CODE IDENTIFYING SET OF UNSYMMETRICAL SCALAR PRODUCTS
C     JU12     CODE FOR AB IN TEST FOR ZERO ELEMENT
C     JU23     CODE FOR BC IN TEST FOR ZERO ELEMENT
C     JU31     CODE FOR CA IN TEST FOR ZERO ELEMENT
C     ITT      CODE SPECIFYING SECTIONS OF SUBROUTINE TEST
C               
      CHARACTER TITLE*70
      REAL DIS(6,12),DET(15),BET(15),P(3,3,12),S(3,3,15)
      INTEGER NN(3,3,3),MM(3,3,4),IHKL(12)
      REAL EU(9),EUI(9),F(9),FI(9)                                  
      COMMON /INOUT/ LUNIN,LUNOUT           
      COMMON /TCELL/ A,B,C,ALPHA,BETA,GAMMA,COSAL,COSBE,COSGA,
     +               SINAL,SINBE,SINGA,ALS,BES,GAS,SINALS,SINBES,
     +               SINGAS,COSALS,COSBES,COSGAS,AS,BS,CS,AL,BE,
     +               GA,ALSTAR,BESTAR,GASTAR,V,VS         
      COMMON /TCHAR/ TITLE
      COMMON /WORK/  DIS,K,IHKL,DEL,ISUM,JSUM,P,S,DET,BET,L,L2,
     +               MATRIX,NN,MM,ICYCLE,IOUTPT
C 
C
C
      IF (ITT.EQ.1) THEN
C     
C---- test r(1,1), r(2,2), r(3,3)                                   
C
      AA = A * A
      BB = B * B
      CC = C * C
      BC = B * C * COSAL
      CA = C * A * COSBE
      AB = A * B * COSGA
      JU23 = 0
      JU31 = 0
      JU12 = 0
      JU = 0
C
      IF ((ABS(AA-BB) - DEL).GT.0.0)GOTO722
      IF ((ABS((AA+BB)/2. - CC) - DEL).GT.0.0) THEN
       ISYM = 2
       ELSE
       ISYM = 1
        ENDIF
      GO TO 100
C
C
  722 IF ((ABS(AA-CC) - DEL).LE.0.0)THEN
      ISYM = 3
      GO TO 100
       ENDIF
       IF ((ABS(BB-CC) - DEL).LE.0.0) WRITE (LUNOUT, 7001 )
  727 ISYM = 4
C     
C---- test r(2,3), r(3,1), r(1,2)                                   
C     1  test for zero elements                                     
C
  100 IF ((ABS(ALPHA-90.0) - DEL).LE.0.0)GOTO102
      IF ((ABS(BC) - DEL).GT.0.0)GOTO103
  102 JU23 = 1
  103 IF ((ABS(BETA-90.0) - DEL).LE.0.0)GOTO105
      IF ((ABS(CA) - DEL).GT.0.0)GOTO106
  105 JU31 = 1
  106 IF ((ABS(GAMMA-90.0) - DEL).LE.0.0)GOTO108
      IF ((ABS(AB) - DEL).GT.0.0)GOTO120
  108 JU12 = 1
C
  120 IF ((JU23+JU31+JU12).EQ.3) GO TO 130
      IF ((JU23+JU31).EQ.2) GO TO 140
      IF ((JU23+JU12).EQ.2) GO TO 150
      IF ((JU31+JU12).EQ.2) GO TO 160
      IF (JU23.EQ.1) GO TO 170
      IF (JU31.EQ.1) GO TO 180
      IF (JU12.EQ.1) GO TO 240
      GO TO 200
C     
C---- r(2,3) = r(3,1) = r(1,2) = 0.0                             
C     cubic(41), tetragonal(31), tetragonal(30) and              
C     orthorhombic(19)                                           
C
  130 IF (ISYM.LE.4) GO TO 11
C
C---- r(2,3) = r(3,1) = 0.0                                      
C     hexagonal(40a,b), orthorhombic(24a,b), orthorhombic(22a,b) 
C     and monoclinic(4a,b)                                       
C
  140 IF ((ABS(AB - AA/2.) - DEL).GT.0.0)GOTO142
      IF (ISYM.EQ.1) ISYM = 2
      IF (ISYM.EQ.1 .OR. ISYM.EQ.3) THEN
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
       ELSE
      WRITE (LUNOUT, 7004  )
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
        ENDIF
        ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.4) THEN
            GOTO 12
           END IF
  142 IF ((ABS(AB + AA/2.) - DEL).GT.0.0)GOTO144
      IF (ISYM.EQ.1) ISYM = 2
      IF (ISYM.EQ.1 .OR. ISYM.EQ.3) THEN
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
       ELSE
      WRITE (LUNOUT, 7004  )
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
        ENDIF
        ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.4) THEN
            GOTO 13
           END IF
  144 IF (AB.LT.0.0)THEN
         GOTO 146
       ELSE IF (AB.EQ.0.0) THEN
         GOTO 130
        ELSE
         GOTO 145
          ENDIF
  145 IF (ISYM.EQ.1) ISYM = 4
      IF (ISYM.EQ.1 .OR. ISYM.EQ.3) THEN
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
       ELSE
      WRITE (LUNOUT, 7004  )
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
        ENDIF
        ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.4) THEN
            GOTO 14
           END IF
  146 IF (ISYM.EQ.1) ISYM = 4
      IF (ISYM.EQ.1 .OR. ISYM.EQ.3) THEN
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
       ELSE
      WRITE (LUNOUT, 7004  )
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
        ENDIF
        ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.4) THEN
            GOTO 15
           END IF
C     
C---- r(2,3) = r(1,2) = 0.0                                      
C     hexagonal(39), orthorhombic(23), orthorhombic(20) and      
C     monoclinic(5)                                              
C
  150 IF ((ABS(CA + AA/2.) - DEL).GT.0.0)GOTO152
      IF (ISYM.EQ.1) ISYM = 3
      IF (ISYM.EQ.3) THEN
           GO TO 12
        ELSE 
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
       ELSE
      WRITE (LUNOUT, 7004  )
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
        ENDIF
           END IF
  152 IF (ISYM.EQ.3) GO TO 13
      IF (ISYM.NE.4) THEN
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
       ELSE
      WRITE (LUNOUT, 7004  )
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
        ENDIF
        ENDIF
      IF ((ABS(CA + CC/2.) - DEL).LE.0.0) THEN
        GOTO 16
        ELSE
        GOTO17
        ENDIF
C
C---- r(3,1) = r(1,2) = 0.0                                      
C     orthorhombic(21), monoclinic(3) and hexagonal(special      
C     case)                                                      
C
  160 IF ((ABS(BC + CC/2.) - DEL).GT.0.0)GOTO162
      IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
       ELSE
      WRITE (LUNOUT, 7004  )
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
        ENDIF
        ELSE IF (ISYM.EQ.1) THEN
            GOTO 23
         ELSE
           GOTO20
           END IF
162      IF (ISYM.EQ.4) THEN
           GO TO 21
        ELSE 
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
       ELSE
      WRITE (LUNOUT, 7004  )
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
        ENDIF
           END IF
C
C---- r(2,3) = 0.0     monoclinic(9a,b)                          
C
  170 IF (ISYM.NE.4) THEN
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
        ENDIF
        ENDIF
      IF ((ABS(CA+CC/2.) - DEL).GT.0.0)THEN
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
       ELSE
      WRITE (LUNOUT, 7004  )
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
        ENDIF
        ENDIF
      IF (AB.LT.0.0)THEN 
            GOTO19
          ELSE IF (AB.EQ.0.0) THEN
          GOTO150
          ELSE
           GOTO 18
           ENDIF
C     
C---- r(3,1) = 0.0                                               
C     tetragonal(32a,b), orthorhombic(27a,b), monoclinic(11a,b)  
C     and monoclinic(10a,b)                                      
C
  180 IF (ISYM.NE.3) GO TO 183
      IF ((ABS(BC+AA/2.) - DEL).GT.0.0)THEN
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
       ELSE
      WRITE (LUNOUT, 7004  )
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
        ENDIF
       ENDIF
      IF ((ABS(AB-AA/2.) - DEL).LE.0.0)GOTO14
      IF ((ABS(AB+AA/2.) - DEL).LE.0.0)THEN
          GOTO15
         ELSE
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
       ELSE
      WRITE (LUNOUT, 7004  )
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
        ENDIF
        ENDIF
  183 IF (ISYM.NE.4) THEN
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
        ENDIF
        ENDIF
      IF ((ABS(BC+CC/2.) - DEL).GT.0.0)GOTO187
      IF ((ABS(AB-AA/2.) - DEL).LE.0.0)GOTO22
      IF ((ABS(AB+AA/2.) - DEL).LE.0.0)GOTO23
      IF (AB.LT.0.0)THEN
          GOTO25
        ELSE IF (AB.EQ.0.0) THEN
          GOTO 160
          ELSE
         GOTO24
         ENDIF
  187 IF ((ABS(AB-AA/2.) - DEL).LE.0.0)GOTO26
      IF ((ABS(AB+AA/2.) - DEL).LE.0.0) THEN
        GOTO 27
         ELSE
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
       ELSE
      WRITE (LUNOUT, 7004  )
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
        ENDIF
         ENDIF
C     
C---- 2  test for cubic(42), tetragonal(34a,b,c) and                
C        orthorhombic(25a,b,c)                                      
C
  200 IF (ISYM.NE.1) GO TO 230
      IF ((ABS(BC+CA+AB + AA) - DEL).GT.0.0)GOTO230
      IF ((ABS(ABS(BC) - AA/3.) - DEL).GT.0.0)GOTO210
      IF ((ABS(ABS(CA) - AA/3.) - DEL).GT.0.0)GOTO210
      IF ((ABS(ABS(AB) - AA/3.) - DEL).LE.0.0)THEN
         GOTO 16
        ELSE
       GOTO 210
       ENDIF
C
  210 IF ((ABS(ABS(CA) - ABS(AB)) - DEL).GT.0.0)GOTO212
      IF ((ABS(ABS(CA) - ABS((-AA-BC)/2.)) - DEL).LE.0.0)THEN
       GOTO17
      ELSE
      GOTO230
       ENDIF
  212 IF ((ABS(ABS(BC) - ABS(AB)) - DEL).GT.0.0)GOTO214
      IF ((ABS(ABS(BC) - ABS(-AA-CA)/2.) - DEL).LE.0.0)THEN
        GOTO18
        ELSE
        GOTO230
        ENDIF
  214 IF ((ABS(ABS(BC) - ABS(CA)) - DEL).GT.0.0)GOTO220
      IF ((ABS(ABS(BC) - ABS((-AA-AB)/2.)) - DEL).LE.0.0)THEN
          GOTO 19
         ELSE
       GOTO230
       ENDIF
C
  220 IF ((ABS(ABS(BC) - ABS(-AA-CA-AB)) - DEL).LE.0.0)GOTO20
      IF ((ABS(ABS(CA) - ABS(-AA-BC-AB)) - DEL).LE.0.0)GOTO21
      IF ((ABS(ABS(AB) - ABS(-AA-BC-CA)) - DEL).LE.0.0)GOTO 22
C     
C---- 3  test for r(2,3) = r(3,1)                                   
C
  230 IF ((ABS(ABS(BC) - ABS(CA)) - DEL).GT.0.0)GOTO270
      IF ((ABS(BC+AA/2.) - DEL).LE.0.0)GOTO241
      IF ((ABS(BC+CC/2.) - DEL).LE.0.0)GOTO250
      IF ((ABS(BC+CC/3.) - DEL).LE.0.0)THEN
         GOTO256
          ELSE
          GOTO260
          ENDIF
C     
C---- r(2,3) = r(3,1) = -r(1,1)/2                                
C     cubic(43a,b) and rhombohedral(35)                          
C
240      IF (ISYM.EQ.1 ) THEN
           GO TO 12
        ELSE 
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
       ELSE
      WRITE (LUNOUT, 7004  )
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
        ENDIF
           END IF
  241 IF ((ABS(AB-AA/2.) - DEL).GT.0.0)GOTO260
      IF (ISYM.EQ.2 .OR. ISYM.EQ.4) THEN
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
       ELSE
      WRITE (LUNOUT, 7004  )
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
        ENDIF
        ELSE IF (ISYM.EQ.1) THEN
            GOTO 13
         ELSE
        GOTO 18
           END IF
C     
C---- r(2,3) = r(3,1) = -r(3,3)/2                                
C     tetragonal(33), orthorhombic(29), orthorhombic(26a,b)      
C     and monoclinic(8)                                          
C
  250 IF ((ABS(AB-CC/4.) - DEL).GT.0.0)GOTO252
      IF (ISYM.EQ.1 .OR. ISYM.EQ.3) THEN
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
       ELSE
      WRITE (LUNOUT, 7004  )
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
        ENDIF
        ELSE IF (ISYM.EQ.2 ) THEN
            GOTO 16
         ELSE
        GOTO 28
           END IF
  252 IF (AB.GE.0.0)GOTO254
      IF (ISYM.EQ.2) GO TO 18
      GO TO 260
254      IF (ISYM.EQ.1 .OR. ISYM.EQ.3) THEN
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
       ELSE
      WRITE (LUNOUT, 7004  )
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
        ENDIF
        ELSE IF (ISYM.EQ.2 ) THEN
            GOTO 17
           ELSE
           GOTO 29
           END IF
C     
C---- r(2,3) = r(3,1) = -r(3,3)/3                                
C     rhombohedral(38)                                           
C
  256 IF (ISYM.NE.2) GO TO 260
      IF ((ABS(ABS(AB) - ABS((-AA+CC/3.)/2.)) - DEL).LE.0.0)
     +   GOTO 19
C     
C---- r(2,3) = r(3,1) not necessarily = r(1,2)                   
C     rhombohedral(36), rhombohedral(37), monoclinic(14)         
C     and monoclinic(15)                                         
C
  260 IF ((ABS(ABS(BC) - ABS(AB)) - DEL).GT.0.0)GOTO262
      IF (ISYM.NE.1) GO TO 262
      IF (AB.LT.0.0)THEN
         GOTO 15
         ELSE IF (AB.EQ.0.0) THEN
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
       ELSE
      WRITE (LUNOUT, 7004  )
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
        ENDIF
          ELSE
          GOTO 14
         ENDIF
  262 IF (ISYM.NE.2) GO TO 270
      IF (AB.LT.0.0)THEN
           GOTO21
           ELSE IF (AB.EQ.0.0) THEN
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
       ELSE
      WRITE (LUNOUT, 7004  )
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
        ENDIF
           ELSE
            GOTO20
           ENDIF
C     
C----  test for abs(r(2,3)) = abs(r(1,2))                         
C      orthorhombic(28), monoclinic(12) and monoclinic(13)        
C
  270 IF ((ABS(ABS(BC) - ABS(AB)) - DEL).GT.0.0)GOTO280
      IF (ISYM.NE.3) GO TO 280
      IF ((ABS(CA - (-AA - 2.*BC)) - DEL).LE.0.0)GOTO19
      IF (AB.LT.0.0)THEN
         GOTO17
        ELSE IF (AB.EQ.0.0) THEN
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
       ELSE
      WRITE (LUNOUT, 7004  )
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
        ENDIF
         ELSE
         GOTO 16
         ENDIF
C     
C---- 5  test for remaining cells                                   
C        monoclinic(7)                                              
C
  280 IF (ISYM.NE.3) GO TO 300
      IF ((ABS(ABS(CA) - ABS(-AA-BC-AB)) - DEL).LE.0.0)GOTO20
C
C---- monoclinic(16)                                             
C
300      IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
      IF (ISYM.EQ.1) THEN
      ISYM = 3
      GO TO 100
      ELSE IF (ISYM.EQ.2 .OR. ISYM.EQ.3) THEN
         GOTO 727
       ELSE
      WRITE (LUNOUT, 7004  )
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
        ENDIF
        ELSE IF (ISYM.EQ.1) THEN
            GOTO 340
          ELSE
          GOTO301
           END IF
  301 IF ((ABS(BC + CC/2.) - DEL).GT.0.0)GOTO310
      IF ((ABS(AB + CA/2.) - DEL).LE.0.0)GOTO30
C
C---- monoclinic(17)                                             
C
  310 IF ((ABS(BC - CA/2.) - DEL).GT.0.0)GOTO320
      IF ((ABS(AB - AA/2.) - DEL).LE.0.0)GOTO31
C
C---- monoclinic(18)                                             
C
  320 IF ((ABS(AB + BC/2.) - DEL).GT.0.0)GOTO330
      IF ((ABS(CA + CC/2.) - DEL).LE.0.0)GOTO32
C
C---- monoclinic(6)                                              
C
  330 IF ((ABS(ABS(BC) - ABS((-CC - CA)/2.)) - DEL).GT.0.0)GOTO340
      IF ((ABS(ABS(AB) - ABS((-AA - CA)/2.)) - DEL).LE.0.0)GOTO33
C
C---- triclinic(1) and triclinic(2)                              
C
  340 ISYM = 4
      IF (AB.LT.0.0) THEN
         GOTO 35
         ELSE
         GOTO34
        ENDIF
C
   11 JU = 1
      GO TO 65
   12 JU = 2
      GO TO 65
   13 JU = 3
      GO TO 65
   14 JU = 4
      GO TO 65
   15 JU = 5
      GO TO 65
   16 JU = 6
      GO TO 65
   17 JU = 7
      GO TO 65
   18 JU = 8
      GO TO 65
   19 JU = 9
      GO TO 65
   20 JU = 10
      GO TO 65
   21 JU = 11
      GO TO 65
   22 JU = 12
      GO TO 65
   23 JU = 13
      GO TO 65
   24 JU = 14
      GO TO 65
   25 JU = 15
      GO TO 65
   26 JU = 16
      GO TO 65
   27 JU = 17
      GO TO 65
   28 JU = 18
      GO TO 65
   29 JU = 19
      GO TO 65
   30 JU = 20
      GO TO 65
   31 JU = 21
      GO TO 65
   32 JU = 22
      GO TO 65
   33 JU = 23
      GO TO 65
   34 JU = 24
      GO TO 65
   35 JU = 25
C
C---- tests completed   now execute assignment of matrix elements   
C     for transformation of reduced cell (in conventional form) to
C     unit cell
C
C
   65 DO 7 I = 1,9
      IF (ISYM.EQ.1) THEN
      GO TO (1001,1002,1003,1004,1005,1006,1007,1008,1009), I
      ELSE IF (ISYM.EQ.2) THEN
      GO TO (1011,1012,1013,1014,1015,1016,1017,1018,1019), I
      ELSE IF (ISYM.EQ.3) THEN
      GO TO (1021,1022,1023,1024,1025,1026,1027,1028,1029), I
      ELSE IF (ISYM.EQ.4) THEN
      GO TO (1031,1032,1033,1034,1035,1036,1037,1038,1039), I
      ENDIF
C
C---- isym = 1                                                      
C
C---- cell order = 41,43a,43b,36,37,42,34a,34b,34c,25a,25b,25c,
C                  hex(spec case)  
C
 1001 GO TO (1,1,1,4,1,1,1,1,4,1,4,1,4), JU
 1002 GO TO (4,5,1,1,4,1,4,1,1,1,1,4,1), JU
 1003 GO TO (4,4,1,4,4,4,1,4,1,4,1,1,4), JU
 1004 GO TO (4,1,1,1,4,4,1,4,1,4,1,1,4), JU
 1005 GO TO (1,1,5,4,1,1,1,1,4,1,4,1,4), JU
 1006 GO TO (4,4,5,4,4,1,4,1,1,1,1,4,1), JU
 1007 GO TO (4,1,5,4,4,1,4,1,1,1,1,4,1), JU
 1008 GO TO (4,1,1,4,4,4,1,4,1,4,1,1,4), JU
 1009 GO TO (1,2,5,5,1,1,1,1,4,1,4,1,4), JU
C
C---- isym = 2                                                      
C      
C---- cell order = 31,40a,40b,24a,24b,33,26a,26b,38,14,15 
C 
 1011 GO TO (1,1,1,5,1,1,1,1,1,1,1), JU
 1012 GO TO (4,4,4,1,1,5,1,5,1,1,1), JU
 1013 GO TO (4,4,4,4,4,4,1,4,1,4,4), JU
 1014 GO TO (4,5,4,1,5,1,5,1,5,5,5), JU                     
 1015 GO TO (1,1,1,1,1,1,1,1,4,1,1), JU
 1016 GO TO (4,4,4,4,4,1,4,1,4,4,4), JU
 1017 GO TO (4,4,4,4,4,4,4,4,4,4,4), JU
 1018 GO TO (4,4,4,4,4,4,4,4,5,4,4), JU
 1019 GO TO (1,1,1,1,1,1,1,1,4,1,1), JU
C
C---- isym = 3                                                      
C
C---- cell order = 30,39,23,32a,32b,12,13,35,28,7    
C
 1021 GO TO (4,4,1,5,4,5,1,4,5,1), JU
 1022 GO TO (4,4,4,4,4,4,4,1,4,4), JU
 1023 GO TO (1,1,1,4,1,1,1,1,1,5), JU
 1024 GO TO (1,1,1,4,1,1,1,4,1,1), JU
 1025 GO TO (4,4,4,4,4,4,4,1,2,4), JU
 1026 GO TO (4,4,5,1,4,1,5,4,1,1), JU
 1027 GO TO (4,4,4,5,1,4,4,5,5,5), JU
 1028 GO TO (1,1,1,2,2,1,1,1,4,5), JU
 1029 GO TO (4,4,4,1,1,4,4,4,5,4), JU
C      
C---- isym = 4                                                      
C
C---- cell order = 19,22a,22b,4a,4b,20,5,9a,9b,21,3,27a,27b,11a,11b,
C                  10a,10b,29,8,16,17,18,6,1,2
C
 1031 GO TO (1,1,1,5,1,4,4,2,2,4,4,5,1,4,4,5,1,2,6,4,1,4,1,1,1), JU
 1032 GO TO (4,4,4,4,4,4,4,4,4,4,4,4,4,2,2,2,2,4,4,4,4,4,4,4,4), JU
 1033 GO TO (4,4,4,4,4,1,1,1,1,1,1,4,4,1,1,4,4,1,5,1,4,1,1,4,4), JU
 1034 GO TO (4,5,5,4,4,6,4,4,4,4,1,5,1,4,4,5,1,4,4,4,5,2,1,4,4), JU
 1035 GO TO (1,2,6,4,4,4,1,4,4,6,4,2,2,4,4,4,4,2,4,2,2,4,2,1,1), JU
 1036 GO TO (4,4,4,1,5,5,4,1,1,5,4,1,1,1,1,4,4,1,1,1,4,1,1,4,4), JU
 1037 GO TO (4,4,4,4,4,4,1,4,4,1,4,4,4,5,1,4,4,4,1,1,4,4,4,4,4), JU
 1038 GO TO (4,4,4,1,1,1,4,5,1,4,1,4,4,4,4,4,4,4,5,4,4,1,4,4,4), JU
 1039 GO TO (1,1,1,4,4,4,4,4,4,4,4,1,1,4,4,1,1,1,4,4,1,4,5,1,1), JU
C
    1 EU(I) = 1
      GO TO 7
    2 EU(I) = 2
      GO TO 7
C   3 EU(I) = 3                        
C     GO TO 7                          
    4 EU(I) = 0
      GO TO 7
    5 EU(I) = -1
      GO TO 7
    6 EU(I) = -2
    7 CONTINUE
C 
C---- generate inverse transformation matrix of eu                  
C
      DETEU = EU(1)*(EU(5)*EU(9)-EU(8)*EU(6)) - EU(2)*(EU(4)*EU(9)-
     1 EU(7)*EU(6)) + EU(3)*(EU(4)*EU(8)-EU(7)*EU(5))
C
      IF (DETEU.EQ.0.0)THEN
         GOTO 96
         ELSE IF(DETEU.GT.0.0) THEN
          GOTO 54
          ENDIF    
      DETEU = -DETEU
      DO 53 I = 1,9    
      EU(I) = -EU(I)
   53 CONTINUE
C
   54 EUI(1) =  (EU(5)*EU(9) - EU(8)*EU(6)) / DETEU
      EUI(2) = -(EU(2)*EU(9) - EU(8)*EU(3)) / DETEU
      EUI(3) =  (EU(2)*EU(6) - EU(5)*EU(3)) / DETEU
      EUI(4) = -(EU(4)*EU(9) - EU(7)*EU(6)) / DETEU
      EUI(5) =  (EU(1)*EU(9) - EU(7)*EU(3)) / DETEU
      EUI(6) = -(EU(1)*EU(6) - EU(4)*EU(3)) / DETEU
      EUI(7) =  (EU(5)*EU(8) - EU(7)*EU(5)) / DETEU
      EUI(8) = -(EU(1)*EU(8) - EU(7)*EU(2)) / DETEU
      EUI(9) =  (EU(1)*EU(5) - EU(4)*EU(2)) / DETEU
      DETEUI = 1./DETEU
C     
C---- calculate the new cell parameters                             
C
      ABC = A*B*COSGA
      ACC = A*C*COSBE
      BCC = B*C*COSAL
      A2 = SQRT((EU(1)*A)**2 + (EU(2)*B)**2 + (EU(3)*C)**2
     1 + 2.*(EU(2)*EU(3)*BCC + EU(1)*EU(3)*ACC + EU(1)*EU(2)*ABC))
      B2 = SQRT((EU(4)*A)**2 + (EU(5)*B)**2 + (EU(6)*C)**2
     1 + 2.*(EU(5)*EU(6)*BCC + EU(4)*EU(6)*ACC + EU(4)*EU(5)*ABC))
      C2 = SQRT((EU(7)*A)**2 + (EU(8)*B)**2 + (EU(9)*C)**2
     1 + 2.*(EU(8)*EU(9)*BCC + EU(7)*EU(9)*ACC + EU(7)*EU(8)*ABC))
      ALPHA2 = 57.295780 * ACOS((EU(4)*EU(7)*A*A + EU(5)*EU(8)*B*B
     1 + EU(6)*EU(9)*C*C + (EU(4)*EU(8)+EU(5)*EU(7))*ABC + (EU(4)*EU(9)
     2 +EU(6)*EU(7))*ACC + (EU(5)*EU(9)+EU(6)*EU(8))*BCC) / (B2*C2))
      BETA2  = 57.295780 * ACOS((EU(1)*EU(7)*A*A + EU(2)*EU(8)*B*B
     1 + EU(3)*EU(9)*C*C + (EU(1)*EU(8)+EU(2)*EU(7))*ABC + (EU(1)*EU(9)
     2 +EU(3)*EU(7))*ACC + (EU(2)*EU(9)+EU(3)*EU(8))*BCC) / (A2*C2))
      GAMMA2 = 57.295780 * ACOS((EU(1)*EU(4)*A*A + EU(2)*EU(5)*B*B
     1 + EU(3)*EU(6)*C*C + (EU(1)*EU(5)+EU(2)*EU(4))*ABC + (EU(1)*EU(6)
     2 +EU(3)*EU(4))*ACC + (EU(2)*EU(6)+EU(3)*EU(5))*BCC) / (A2*B2))
C     
C---- output results                                                
C
      JSUM = ISUM + 1
      IF (ISYM.EQ.1) THEN
      GO TO (712,714,714,710,710,713,709,709,709,706,706,706,711), JU
      ELSE IF (ISYM.EQ.2) THEN
      GO TO (708,711,711,705,705,709,706,706,710,703,703), JU
      ELSE IF (ISYM.EQ.3)THEN
      GO TO (708,711,705,709,709,703,703,710,707,703), JU
      ELSE IF (ISYM.EQ.4) THEN
      GO TO (704,705,705,702,702,705,702,703,703,705,702,706,706,703,
     1   703,703,703,707,703,703,703,703,703,701,701), JU
       ENDIF
  701 WRITE (LUNOUT, 5001) JSUM
      GO TO 89
  702 WRITE (LUNOUT, 5002) JSUM
      GO TO 89
  703 WRITE (LUNOUT, 5003) JSUM
      GO TO 89
  704 WRITE (LUNOUT, 5004) JSUM
      GO TO 89
  705 WRITE (LUNOUT, 5005) JSUM
      GO TO 89
  706 WRITE (LUNOUT, 5006) JSUM
      GO TO 89
  707 WRITE (LUNOUT, 5007) JSUM
      GO TO 89
  708 WRITE (LUNOUT, 5008) JSUM
      GO TO 89
  709 WRITE (LUNOUT, 5009) JSUM
      GO TO 89
  710 WRITE (LUNOUT, 5010) JSUM
      GO TO 89
  711 WRITE (LUNOUT, 5011) JSUM
      GO TO 89
  712 WRITE (LUNOUT, 5012) JSUM
      GO TO 89
  713 WRITE (LUNOUT, 5013) JSUM
      GO TO 89
  714 WRITE (LUNOUT, 5014) JSUM
      GO TO 89
C
C---- ITT = 2
C
       ELSE IF (ITT.EQ.2) THEN
      WRITE (LUNOUT, 6019) ISUM,JSUM,JSUM,ISUM
      WRITE (LUNOUT, 6020) 
     +  EU(1),EU(2),EU(3),EUI(1),EUI(2),EUI(3),EU(4),EU(5),    
     1 EU(6),EUI(4),EUI(5),EUI(6),EU(7),EU(8),EU(9),EUI(7),EUI(8),EUI(9)
      WRITE (LUNOUT, 6046) DETEU,DETEUI
      LINE = LINE + 1
      IF ((LINE - 5).LT.0)GOTO89
      LINE = 1
      WRITE (LUNOUT, 4998) TITLE
      GO TO 89
C     
C---- calculate matrix elements, f, for transformation of original  
C     cell to final unit cell
C    
C
C---- ITT = 3
C
      ELSE IF (ITT.EQ.3) THEN
      F(1) = EU(1)*P(1,1,L2) + EU(2)*P(2,1,L2) + EU(3)*P(3,1,L2)
      F(2) = EU(1)*P(1,2,L2) + EU(2)*P(2,2,L2) + EU(3)*P(3,2,L2)
      F(3) = EU(1)*P(1,3,L2) + EU(2)*P(2,3,L2) + EU(3)*P(3,3,L2)
      F(4) = EU(4)*P(1,1,L2) + EU(5)*P(2,1,L2) + EU(6)*P(3,1,L2)
      F(5) = EU(4)*P(1,2,L2) + EU(5)*P(2,2,L2) + EU(6)*P(3,2,L2)
      F(6) = EU(4)*P(1,3,L2) + EU(5)*P(2,3,L2) + EU(6)*P(3,3,L2)
      F(7) = EU(7)*P(1,1,L2) + EU(8)*P(2,1,L2) + EU(9)*P(3,1,L2)
      F(8) = EU(7)*P(1,2,L2) + EU(8)*P(2,2,L2) + EU(9)*P(3,2,L2)
      F(9) = EU(7)*P(1,3,L2) + EU(8)*P(2,3,L2) + EU(9)*P(3,3,L2)
C
C---- generate inverse transformation matrix of f                   
C
      DETFN = F(1)*(F(5)*F(9)-F(8)*F(6)) - F(2)*(F(4)*F(9)-F(7)*F(6)) +
     1 F(3)*(F(4)*F(8)-F(7)*F(5))
      IF (DETFN.EQ.0.0) GO TO 97
      FI(1) =  (F(5)*F(9) - F(8)*F(6)) / DETFN
      FI(2) = -(F(2)*F(9) - F(8)*F(3)) / DETFN
      FI(3) =  (F(2)*F(6) - F(5)*F(3)) / DETFN
      FI(4) = -(F(4)*F(9) - F(7)*F(6)) / DETFN
      FI(5) =  (F(1)*F(9) - F(7)*F(3)) / DETFN
      FI(6) = -(F(1)*F(6) - F(4)*F(3)) / DETFN
      FI(7) =  (F(4)*F(8) - F(7)*F(5)) / DETFN
      FI(8) = -(F(1)*F(8) - F(7)*F(2)) / DETFN
      FI(9) =  (F(1)*F(5) - F(4)*F(2)) / DETFN
      DETFNI = 1./DETFN
C     
C
       LINE = LINE + 1
      IF ((LINE - 5).GE.0)THEN
      LINE = 1
      WRITE (LUNOUT, 4998) TITLE
      END IF
      WRITE (LUNOUT, 6021) JSUM,JSUM
      WRITE (LUNOUT, 6020) 
     +F(1),F(2),F(3),FI(1),FI(2),FI(3),F(4),F(5),F(6),FI(4), 
     1 FI(5),FI(6),F(7),F(8),F(9),FI(7),FI(8),FI(9)
      WRITE (LUNOUT, 6046) DETFN,DETFNI
      GO TO 89 
C      
C---- ITT = 4
C
      ELSE IF (ITT.EQ.4) THEN
        GOTO 94
        END IF
C
C---- end of first goto ITT
C
94    CONTINUE
      A = A2
      B = B2
      C = C2
      ALPHA = ALPHA2
      BETA  = BETA2
      GAMMA = GAMMA2
      WRITE (LUNOUT, 4998) TITLE
      IF (ISYM.EQ.1) THEN
      GO TO (812,814,814,810,810,813,809,809,809,806,806,806,811), JU
      ELSE IF (ISYM.EQ.2) THEN
      GO TO (808,811,811,805,805,809,806,806,810,803,803), JU 
      ELSE IF (ISYM.EQ.3) THEN
      GO TO (808,811,805,809,809,803,803,810,807,803), JU 
      ELSE IF (ISYM.EQ.4) THEN
      GO TO (804,805,805,802,802,805,802,803,803,805,802,806,806,803,
     1   803,803,803,807,803,803,803,803,803,801,801), JU
        END IF
C
  801 WRITE (LUNOUT, 6001) JSUM
      GO TO 910
  802 WRITE (LUNOUT, 6002) JSUM
      GO TO 901
  803 WRITE (LUNOUT, 6003) JSUM
      GO TO 901
  804 WRITE (LUNOUT, 6004) JSUM
      GO TO 907
  805 WRITE (LUNOUT, 6005) JSUM
      GO TO 907
  806 WRITE (LUNOUT, 6006) JSUM
      GO TO 907
  807 WRITE (LUNOUT, 6007) JSUM
      GO TO 907
  808 WRITE (LUNOUT, 6008) JSUM
      GO TO 902
  809 WRITE (LUNOUT, 6009) JSUM
      GO TO 902
  810 WRITE (LUNOUT, 6010) JSUM
      GO TO 903
  811 WRITE (LUNOUT, 6011) JSUM
      GO TO 904
  812 WRITE (LUNOUT, 6012) JSUM
      GO TO 905
  813 WRITE (LUNOUT, 6013) JSUM
      GO TO 905
  814 WRITE (LUNOUT, 6014) JSUM
      GO TO 905
C     
C---- reset cell edges and cell angles as prescribed by symmetry    
C     go to 901 (monoclinic), 902 (tetragonal), 903 (rhombohedral),
C     904 (hexagonal), 905 (cubic), 907 (orthorhombic)                  
C
  901 ALPHA = 90.0
      GO TO 909
  902 A = (A + B) / 2.
      GO TO 906
  903 A = (A + B + C) / 3.
      B = A
      C = A
      ALPHA = (ALPHA + BETA + GAMMA) / 3.
      GO TO 908
  904 A = (A + B) / 2.
      B = A
      ALPHA = 90.0
      BETA = 90.0
      GAMMA = 120.0
      GO TO 910
  905 A = (A + B + C) / 3.
      C = A
  906 B = A
  907 ALPHA = 90.0
  908 BETA  = ALPHA
  909 GAMMA = ALPHA
C
  910 WRITE (LUNOUT, 6030) DEL
 6030 FORMAT (' ',9X,' THIS CELL IS THE SUGGESTED UNIT CELL OF HIGHEST',
     1 ' SYMMETRY BASED ON DEL = ',F4.2)
      WRITE (LUNOUT, 6031) ISYM,JU
 6031 FORMAT (' ',9X,'ISYM = ',I2,4X,' JU = ',I3)
      IF ((ABS(AA - BB) - DEL).GT.0.0)GOTO922
      IF ((ABS((AA + BB)/2. - CC) - DEL).GT.0.0)GOTO924
      IF (JU.EQ.13) GO TO 931
      IF (ISYM.NE.1) GO TO 931
      GO TO 89
  924 IF (ISYM.NE.2) GO TO 932
      GO TO 89
  922 IF ((ABS(AA - CC) - DEL).GT.0.0)GOTO89
      IF (ISYM.NE.3) GO TO 933
      GO TO 89
C     
C---- special program messages                                      
C
  931 WRITE (LUNOUT, 6040) ISUM
      GO TO 89
  932 WRITE (LUNOUT, 6041) ISUM
      GO TO 89
  933 WRITE (LUNOUT, 6042) ISUM
      GO TO 89
C
C
C
   96 JSUM = ISUM + 1
      WRITE (LUNOUT, 7002) ISUM,JSUM                                
      GO TO 89
   97 I = 1
      WRITE (LUNOUT, 7002) I,JSUM
      LINE = LINE + 1
      IF ((LINE - 5).GE.0)THEN
      LINE = 1
      WRITE (LUNOUT, 4998) TITLE
      ENDIF
      WRITE (LUNOUT, 6021) JSUM,JSUM
      WRITE (LUNOUT, 6020) 
     +F(1),F(2),F(3),FI(1),FI(2),FI(3),F(4),F(5),F(6),FI(4), 
     1 FI(5),FI(6),F(7),F(8),F(9),FI(7),FI(8),FI(9)
      WRITE (LUNOUT, 6046) DETFN,DETFNI
C
   89 RETURN
C
C---- Format statements
C
 6001 FORMAT ('0',/,10X,' LATTICE PARAMETERS OF CELL (',I1,
     1 ') -- TRICLINIC (PRIMITIVE)')
 6002 FORMAT ('0',/,10X,' LATTICE PARAMETERS OF CELL (',I1,
     1 ') -- MONOCLINIC (PRIMITIVE)')
 6003 FORMAT ('0',/,10X,' LATTICE PARAMETERS OF CELL (',I1,
     1 ') -- MONOCLINIC (C-CENTERED)')                
 6004 FORMAT ('0',/,10X,' LATTICE PARAMETERS OF CELL (',I1,
     1 ') -- ORTHORHOMBIC (PRIMITIVE)')
 6005 FORMAT ('0',/,10X,' LATTICE PARAMETERS OF CELL (',I1,
     1 ') -- ORTHORHOMBIC (C-CENTERED)')
 6006 FORMAT ('0',/,10X,' LATTICE PARAMETERS OF CELL (',I1,
     1 ') -- ORTHORHOMBIC (BODY-CENTERED)')
 6007 FORMAT ('0',/,10X,' LATTICE PARAMETERS OF CELL (',I1,
     1 ') -- ORTHORHOMBIC (FACE-CENTERED)')
 6008 FORMAT ('0',/,10X,' LATTICE PARAMETERS OF CELL (',I1,
     1 ') -- TETRAGONAL (PRIMITIVE)')
 6009 FORMAT ('0',/,10X,' LATTICE PARAMETERS OF CELL (',I1,
     1 ') -- TETRAGONAL (BODY-CENTERED)')
 6010 FORMAT ('0',/,10X,' LATTICE PARAMETERS OF CELL (',I1,
     1 ') -- RHOMBOHEDRAL (PRIMITIVE)')
 6011 FORMAT ('0',/,10X,' LATTICE PARAMETERS OF CELL (',I1,
     1 ') -- HEXAGONAL (PRIMITIVE)')                                    
 6012 FORMAT ('0',/,10X,' LATTICE PARAMETERS OF CELL (',I1,
     1 ') -- CUBIC (PRIMITIVE)')
 6013 FORMAT ('0',/,10X,' LATTICE PARAMETERS OF CELL (',I1,
     1 ') -- CUBIC (BODY-CENTERED)')
 6014 FORMAT ('0',/,10X,' LATTICE PARAMETERS OF CELL (',I1,
     1 ') -- CUBIC (FACE-CENTERED)')
 5001 FORMAT (' ',9X,' CELL (',I1,')   TRICLINIC (PRIMITIVE)')       
 5002 FORMAT (' ',9X,' CELL (',I1,')   MONOCLINIC (PRIMITIVE)')      
 5003 FORMAT (' ',9X,' CELL (',I1,')   MONOCLINIC (C-CENTERED)')     
 5004 FORMAT (' ',9X,' CELL (',I1,')   ORTHORHOMBIC (PRIMITIVE)')    
 5005 FORMAT (' ',9X,' CELL (',I1,')   ORTHORHOMBIC (C-CENTERED)')   
 5006 FORMAT (' ',9X,' CELL (',I1,')   ORTHORHOMBIC (BODY-CENTERED)') 
 5007 FORMAT (' ',9X,' CELL (',I1,')   ORTHORHOMBIC (FACE-CENTERED)')
 5008 FORMAT (' ',9X,' CELL (',I1,')   TETRAGONAL (PRIMITIVE)')      
 5009 FORMAT (' ',9X,' CELL (',I1,')   TETRAGONAL (BODY-CENTERED)')  
 5010 FORMAT (' ',9X,' CELL (',I1,')   RHOMBOHEDRAL (PRIMITIVE)')    
 5011 FORMAT (' ',9X,' CELL (',I1,')   HEXAGONAL (PRIMITIVE)')       
 5012 FORMAT (' ',9X,' CELL (',I1,')   CUBIC (PRIMITIVE)')           
 5013 FORMAT (' ',9X,' CELL (',I1,')   CUBIC (BODY-CENTERED)')       
 5014 FORMAT (' ',9X,' CELL (',I1,')   CUBIC (FACE-CENTERED)')       
 7002 FORMAT ('0',9X,' A VALUE OF MOD = 0 00 HAS BEEN CALCULATED',
     +  ' IN SUBROUTINE TEST ',/,' FOR THE MATRIX OF CELL (',I1,
     +  ') TO CELL (',I1,')',/)  
 7004 FORMAT ('0',20X,' THIS REDUCED CELL CONTAINS ONE OR MORE',
     +  ' INTERRELATED OR SPECIAL SCALAR', /, 21X,
     +  ' PRODUCTS INCONSISTENT WITH THE UNIQUE REDUCED CELLS OF',/,
     +  ' THE TYPE ISYM = 4',  /, 21X,' IT WILL THEREFORE BE ',
     +  ' ARBITRARILY ASSIGNED AS TRICLINIC',/ )                      
 7001 FORMAT ('0',20X,' SUBROUTINE TEST HAS FOUND R(2,2) = R(3,3).NE.R',
     1'(1,1)',/,'   THIS RELATIONSHIP IS NOT ALLOWED', /,
     + 21X,' UPON ENTRY INTO THIS SUBROUTINE AND INDICATES',/,
     + 21X,' THE PRESENCE OF AN ERROR IN THE PROGRAM',/,
     +  21X,' ISYM HAS BEEN SET TO 4 BEFORE CONTINUING',/)      
 6041 FORMAT ('0',9X,' WARNING   THE REDUCED CELL, CELL (',I1,'), OF',
     1' THIS UNIT CELL IS A SPECIAL CASE HAVING PARAMETERS A = B' )     
 6042 FORMAT ('0',9X,' WARNING   THE REDUCED CELL, CELL (',I1,'), OF',
     1' THIS UNIT CELL IS A SPECIAL CASE HAVING PARAMETERS A = C' )     
 6040 FORMAT ('0',9X,' WARNING   THE REDUCED CELL, CELL (',I1,'), OF',
     1' THIS UNIT CELL IS A SPECIAL CASE HAVING PARAMETERS A = B = C ') 
 6019 FORMAT ('0',/,13X, 'CELL (',I1,') TO CELL (',I1,')',
     1             12X, 'CELL (',I1,') TO CELL (',I1,')', // )
 6020 FORMAT (' ', 11X, 3F7.2, 12X, 3F7.2,/ )
 6021 FORMAT ('0',/,13X, ' CELL (1) TO CELL (',I1,')',
     1             12X,  ' CELL (',I1,') TO CELL (1)', // )
 6046 FORMAT (' ',16X,' MOD = ',F8.4,19X,' MOD = ',F8.4 ,/ )
 4998 FORMAT ('1',//5X,A, / )                                          
C
C
      END
