      SUBROUTINE THREED
C  THIS SUBROUTINE TRIES TO FIND THE ANGLE BETWEEN TWO INTERSECTING
C     ZONES (FOUND IN THE PRECEDING PART OF THE PROGRAM)
      DIMENSION SOM(10),ISOM(10),BE(10),IBE(10)
      INTEGER H
      LOGICAL SELF,ENL,TEST
      COMMON /LINES / IDOC,INP,IOUT,MAX,OBS(40),CEL(10)
      COMMON /PARMS / TOL2,TOL3,WAVEL,FMINEQ,FMINIM,SELF,ENL,TEST,LIST,
     *                NSYST(3),OBSMT(40),OBSPT(40),TWTHET(40),TOLG
      COMMON /TEMRES/ LINCO,NQ2,NINDEX,INDEX(40),ROSTER(150,12)
      COMMON /COUNTM/ ITL(-2500:2500)
      F3=6.0/TOL3
      RF3=TOL3/6.0
      IF (TEST) WRITE(IOUT,1)  F3
      NEWNR=0
      DO 1000 J=4,LINCO,4
      LASTNR=NEWNR
      A=ROSTER(J,1)
      B=ROSTER(J,2)
      C=ROSTER(J,3)
      E=ROSTER(J,5)
      F=ROSTER(J,6)
      NEWNR=NINT( ROSTER(J,8) )
      NC=0
      W=SQRT( (4.*A*C-E*E)*(4.*A*B-F*F) )
      BIGD=(W-E*F)/(2.05*A)
      LIMH=NINT(F3*BIGD)
      IF ( LIMH .GT. 2500) LIMH=2500
      IF (TEST) WRITE(IOUT,5) NEWNR,A,B,C,E,F,LIMH
      DO 10 I=-2500,2500
   10 ITL(I)=0
      DO 20 I=1,10
         ISOM(I)=0
         SOM(I)=0.0
         IBE(I)=0
         BE(I)=0.0
   20 CONTINUE
      MIN=4+MAX/10
      NR=1
      IF (ABS(E)+ABS(F) .GT. 1.0) GO TO  200
C  TRICLINIC LATTICES ARE TREATED BETWEEN LABELS 200 AND 300,
C  THE 'MONOCLINIC' LATTICES CAN BE TREATED IN A SIMPLE  WAY:
      DO 100 L=1,2
         DO 90 K=1,2
            X=B*FLOAT(K*K)+C*FLOAT(L*L)
            Z=F3/FLOAT(K*L)
            DO 80 H=0,2
               Y=X+A*FLOAT(H*H)
               DO 70 M=1,MAX
                  IT=NINT( Z*ABS(OBS(M)-Y) )
                  IF (IT.LT.LIMH) ITL(IT)=ITL(IT)+1
   70          CONTINUE
   80       CONTINUE
   90    CONTINUE
  100 CONTINUE
      KG=LIMH-5
      KST=2
      DO 110 I=0,2
  110 ITL(0)=ITL(0)+ITL(I)
      ITL(0)=2*ITL(0)
      IF (ITL(0) .LT. MIN) GO TO 300
      ISOM(NR)=ITL(0)
      SOM(NR)=0.0
      GO TO 300
C  TRICLINIC LATTICES ARE TREATED HERE.
  200 KG=INT(F3*(W+E*F)/(2.05*A) )
      IF (KG .GT. 2500) KG=2500
      DO 290 L=1,2
         RL=L
         DO 280 K=-2,+2
            IF (K .EQ. 0) GO TO 280
            RK=K
            X=RK*RK*B+RL*RL*C
            V=RK*F+RL*E
            Z=F3/(RK*RL)
            DO 270 H=-2,+2
               Y=X+A*FLOAT(H*H)+V*FLOAT(H)
               DO 260 M=1,MAX
                  IT=NINT( Z*(OBS(M)-Y) )
                  IF (ABS(IT) .LE. KG)  ITL(IT)=ITL(IT)+1
  260          CONTINUE
  270       CONTINUE
  280    CONTINUE
  290 CONTINUE
      KST=-KG
      KG=KG-5
  300 CONTINUE
C  NOW FIND THE MOST FREQUENT VALUES OF  * D *
      IF (TEST) WRITE(IOUT,2)
      DO 400 H=KST,KG
         K=0
         M=H+5
         DO 310 L=H,M
  310    K=K+ITL(L)
         IF (K .LE. MIN ) GO TO 400
         IF (TEST) WRITE(IOUT,3) H,K
         DO 350 N=1,NR
            IF (K .LE. ISOM(N) ) GO TO  350
            DO 340 NC=N,NR
               I=NR+N-NC
               ISOM(I+1)=ISOM(I)
  340          SOM(I+1)  =SOM(I)
            ISOM(N)=K
            NR=NR+1
            IF (NR .GT. 8) NR=8
            V=0.
            DO 345 L=H,M
  345       V=V+L*ITL(L)
            SOM(N)=RF3*( V/FLOAT(K) )
            GO TO 380
  350    CONTINUE
  380    IF (ISOM(4)  .GT. MIN) MIN=ISOM(4)
  400 CONTINUE
      IF (TEST) WRITE(IOUT,4) (SOM(I),ISOM(I),I=1,NR)
C WE HAVE FOUND THE VALUES OF D AND STORED THEM IN  SOM IN DESCENDING
C ORDER OF FREQUENCY  (K).
C  WE NOW TRY TO ELIMINATE 'EQUAL' VALUES OF D
      IF (ISOM(1) .LT.MIN) GO TO 1000
C  IF NO ANSWER HAS BEEN FOUND WE ARE READY
      DO 450 H=1,NR-1
         DO 440 K=H+1,NR
            IF (ABS(SOM(H)  -SOM(K)  ).GT. FMINEQ ) GO TO 440
            IF (ISOM(K) .LT. 1) GO TO 440
            X=ISOM(H)
            Y=ISOM(K)
            SOM(H)=(X*SOM(H)+Y*SOM(K))/(X+Y)
            ISOM(K)=0
  440    CONTINUE
  450 CONTINUE
C THE BEST FOUR VALUES OF D ARE STORED,FIRST IN BE
      N=0
      MIN=ISOM(1)-5
      DO 480 H=1,NR
         IF (ISOM(H) .LE. MIN) GO TO 480
         N=N+1
         IBE(N)=ISOM(H)
         BE(N)  =SOM(H)
         IF (N.EQ.4) GO TO 500
  480 CONTINUE
  500 LN=0
      IF (NEWNR .EQ. LASTNR) LN=5
      L=10*NEWNR+LN
C  WE HAVE ADDED 5 TO THE NUMBERS OF THE 'EQUAL' ZONES.
      DO 550 H=1,N
         M=J+H-1
         DO 540 K=1,6
 540     ROSTER(M,K)=ROSTER(J,K)
         ROSTER(M,4)=BE(H)
         ROSTER(M,7)=IBE(H)
         ROSTER(M,8)=L+H
  550 CONTINUE
 1000 CONTINUE
      RETURN
    1 FORMAT(32H0TEST OUTPUT FROM  ** THREED **  ,'THE FACTOR F3=',F6.3)
    2 FORMAT(12H0POSITION NR )
    3 FORMAT(1H ,I8,I3)
    4 FORMAT(20H0ARRAYS SOM AND ISOM /  (F11.2,I4 ) )
    5 FORMAT(18H0ZONE COMBINATION ,I2,3F8.2,5X,3HXXX,2F8.2,
     1  5X,'LIMH=',I4  )
      END

      SUBROUTINE NORM3
C-----------------------------------------------------------------------
C  THIS SUBROUTINE 'NORMALIZES' THE LATTICES AND PUTS THEM INTO A
C  STANDARD SHAPE
C-----------------------------------------------------------------------
C
      DIMENSION  OLD(6),          IZON(3,6),QVAL(5,5,5)
      INTEGER H,H1
      LOGICAL SELF,ENL,TEST,CHANGE
      COMMON /LINES / IDOC,INP,IOUT,MAX,OBS(40),CEL(10)
      COMMON /PARMS / TOL2,TOL3,WAVEL,FMINEQ,FMINIM,SELF,ENL,TEST,LIST,
     *                NSYST(3),OBSMT(40),OBSPT(40),TWTHET(40),TOLG
      EQUIVALENCE (CEL(1),A),(CEL(2),B),(CEL(3),C),(CEL(4),D),
     *(CEL(5),E),(CEL(6),F),(CEL(8),PED)
C-----------------------------------------------------------------------
C
      DO 80 I=1,6
   80 OLD(I)=CEL(I)
C  INITIALIZING *OLD*
      GO TO 100
C
C-----------------------------------------------------------------------
C  THE NEXT PART IS THE OLD PROCEDURE 'MAZEN'
C  MAZEN FINDS THE  VALUES OF D,E AND F IN A LATTICE, WHEN THE COORDINAT
C  OF A,B,AND C HAVE BEEN GIVEN.
C
   10 I1=IZON(H,1)+IZON(K,1)-3
      I2=IZON(H,2)+IZON(K,2)-3
      I3=IZON(H,3)+IZON(K,3)-3
      X=QVAL(I1,I2,I3)
      I1=IZON(H,1)-IZON(K,1)+3
      I2=IZON(H,2)-IZON(K,2)+3
      I3=IZON(H,3)-IZON(K,3)+3
      Y=QVAL(I1,I2,I3)
      Z=(X-Y)/2.0
      GO TO (601,602,603),ISW
C
C-----------------------------------------------------------------------
C  THE FOLLOWING IS THE PROCEDURE 'CHANGE'
C  IF ANY CONSTANT HAS BEEN  CHANGED,THE PROGRAM REPEATS UNTIL NO
C  MORE CHANGES ARE DETECTED.
C-----------------------------------------------------------------------
   20 U=0.0
      DO 30 I=1,6
      U=U+ ABS(CEL(I)   -OLD(I))
      OLD(I)=CEL(I)
   30 CONTINUE
      CHANGE=.FALSE.
      IF (U.GT.0.1 ) CHANGE=.TRUE.
      GO TO (310,610),ISW
C-----------------------------------------------------------------------
C
  100 IF (TEST) WRITE(IOUT,9930) A,B,C,D,E,F,PED
      QMAX=1.0E+6
C--CHECK ON  CERTAIN RELATIONS BETWEEN A,B,C,D,E,F  WHICH MEAN THAT
C--THE LATTICE HAS HIGHER SYMMETRY.
  200 DO 300 II=1,3
      X=A
      A=B
      B=C
      C=X
      X=D
      D=E
      E=F
      F=X
      IF (B .LT. FMINIM) GO TO 1000
      IF (D.GT. 0.0)  GO TO 205
      D=-D
      E=-E
  205 IF (E.GT. 0.0)  GO TO 210
      E=-E
      F=-F
  210 IF (ABS(A-C) .GT. FMINEQ) GO TO 230
      IF (E .LT. FMINEQ) GO TO 230
C  IF A=C .AND. E>0 THEN WE CAN SIMPLIFY
      A=(A+C+E)/4.
      C=A-E/2.
      F=(D+F)/2.
      D=ABS(D-F)
      E=0
  230 IF (ABS(E-A).GT.FMINEQ ) GO TO 240
      A=A/4.
      C=C-A
      E=0.0
      F=F/2.
      D=ABS(D-F)
  240 IF (ABS(E-C).GT.FMINEQ) GO TO 250
      C=C/4.0
      A=A-C
      E=0.0
      D=D/2.
      F=ABS(D-F)
  250 IF (ABS(D) .LE. FMINEQ) D=0.0
      IF (ABS(E) .LE. FMINEQ) E=0.0
      IF (ABS(F) .LE. FMINEQ) F=0.0
  300 CONTINUE
      IF (TEST) WRITE(IOUT,9931) A,B,C,D,E,F,PED
C  REDUCED    ......  REDUCED2 IS ABOUT AT CARD 200
      ISW=1
      GO TO 20
C--TEST ON CHANGES.
  310 IF (CHANGE) GO TO 200
C--CALCULATE  THE Q-VALUES NEAR THE ORIGIN
      DO 400 L1=3,4
      L=L1-3
      DO 390 K1=2,4
      K=K1-3
      DO 380 H1=2, 4
      H=H1-3
      IF (L  .EQ.0 .AND. K .LT.1) GO TO 330
      X=A*FLOAT(H*H)+B*FLOAT(K*K)+C*FLOAT(L*L)+D*FLOAT(K*L)+E*FLOAT(H*L)
     1  +F*FLOAT(H*K)
      IF (X .LT.(-0.1))GO TO 1000
  320 QVAL (H1,K1,L1)=X
      GO TO 380
  330 QVAL(H1,K1,L1)=QMAX
  380 CONTINUE
  390 CONTINUE
  400 CONTINUE
      QVAL(4,3,3) =A
C--NOW FIND THE SMALLEST Q-VALUES.
      N1=0
  410 N1=N1+1
      U=QMAX-10.0
      DO 500 L=3,4
      DO 490 K=2,4
      DO 480 H=2,4
      IF (QVAL(H,K,L).GE.U) GO TO 480
      U=QVAL(H,K,L)
      IZON(N1,3)=L
      IZON(N1,1)=H
      IZON(N1,4)=H
      IZON(N1,2)=K
      IZON(N1,5)=K
  480 CONTINUE
  490 CONTINUE
  500 CONTINUE
      I1=IZON(N1,1)
      I2=IZON(N1,2)
      I3=IZON(N1,3)
      QVAL(I1,I2,I3)=QMAX
      IF (N1.LT.3) GO TO 410
C--CHECK WHETHER THEY ARE IN ONE PLANE.
      II=0
      DO  510 I=1,3
      K=IZON(1,I+1)-3
      L=IZON(2,I+2)-3
      M=IZON(2,I+1)-3
      N=IZON(1,I+2)-3
      J=K*L-M*N
      II=II+J*(IZON(3,I)-3)
  510 CONTINUE
      IF (II .NE. 0) GO TO 520
      N1=2
      GO TO  410
C--CALCULATE THE RECIPROCAL LATTICE AGAIN., TO DETERMINE D,E, AND F.
  520 DO 580 L1=1,5
      L=L1-3
      DO 570 K1=1,5
      K=K1-3
      DO 560 H1=1,5
      H=H1-3
      Q=A*FLOAT(H*H)+B*FLOAT(K*K)+C*FLOAT(L*L)+D*FLOAT(K*L)+E*FLOAT(H*L)
     1  +F*FLOAT(H*K)
      IF (Q .LT.(-0.1)) GO TO 1000
C  IF ANY Q-VALUE IS NEGATIVE,THE LATTICE SHOULD BE DISCARDED
      QVAL(H1,K1,L1)=Q
  560 CONTINUE
  570 CONTINUE
  580 CONTINUE
      I11=IZON(1,1)
      I12=IZON(1,2)
      I13=IZON(1,3)
      I21=IZON(2,1)
      I22=IZON(2,2)
      I23=IZON(2,3)
      I31=IZON(3,1)
      I32=IZON(3,2)
      I33=IZON(3,3)
      B=QVAL(I11,I12,I13)
      A=QVAL(I21,I22,I23)
      C=QVAL(I31,I32,I33)
C  I KNOW, IKNOW,    BUT THIS PIDGIN-FORTRAN WORKS ON ALL MACHINES....
C--NOW FIND D,E AND F.
      ISW=1
      H=1
      K=3
      GO TO  10
  601 D=Z
      ISW=2
      K=2
      GO TO 10
  602 F=Z
      ISW=3
      H=3
      GO TO 10
  603 E=Z
      IF (TEST) WRITE(IOUT,9934) (CEL(I),I=1,6),PED
C  REDUCED2 ETC.
      ISW=2
      GO TO 20
C--TEST ON CHANGES.
  610 IF (CHANGE) GO TO 200
C  REARRANGE IN A CRYSTALLOGRAPHICAL ORDER. FIRST FIND ZERO TERMS.
      L=0
      J=1
      DO 620 I=4,6
      L=L+I
      IF (ABS(CEL(I)) .GT. FMINEQ) GO TO 620
      J=J+1
      K=I
      L=L-I
      CEL(I)=0.0
  620 CONTINUE
C  J COUNTS THE NUMBER OF ZERO ELEMENTS,STARTING AT 1
C  IF ONLY ONE ELEMENT IS ZERO, IT IS INDICATED BY K
C  IF ONLY ONE ELEMENT IS NON-ZERO, IT IS INDICATED BY L
C  IN THE OTHER CASES, ALL ELEMENTS CAN BE TREATED ALIKE.
C          TRICL    MONO ORTHO
      GO TO (750,700,650,800),J
C  THE MONOCLINIC CASE
  650 L=L-3
      IF (L .NE. 2) CALL WISSEL(L,2)
      IF (A .GT. C) CALL WISSEL(1,3)
      E=ABS(E)
      IF (E .LE.A) GO TO 800
      C=C+A-E
      E=A+A-E
      GO TO 200
C  THE SEMI-TRICLINIC CASE (ONE RECIPROCAL ANGLE 90 DEGREES)
  700 K=K-3
      IF (K .NE. 3) CALL WISSEL(K,3)
C  THE TRUE TRICLINIC CASE
  750 IF (D .GE. 0.) GO TO 760
      D=-D
      E=-E
  760 IF (E .GE. 0.) GO TO 770
      E=-E
      F=-F
  770 IF (ABS(A+B+F-D-E) .GT. FMINEQ) GO TO 800
C--THIS IS A MONOCLINIC LATTICE IN DISGUISE.
      E=ABS(A-B)*0.5
      A=(A+B+F)/4.
      B=C-A
      C=A-F/2.
      D=0.0
      F=0.0
      GO TO 200
C  THIS WAS THE TRANSFORMATION OF CERTAIN MONOCLINIC CENTERED LATTICES.
 800  CONTINUE
      IF (TEST) WRITE(IOUT,9935) (CEL(I),I=1,6),PED
C  REDUCED3  ETC.
      RETURN
 1000 A=0.0
      RETURN
 9930 FORMAT(11H0LATTICE IS, 6F8.1,F8.0)
 9931 FORMAT(11H REDUCED   , 6F8.1,F8.0)
 9934 FORMAT(11H REDUCED2  , 6F8.1,F8.0)
 9935 FORMAT(11H REDUCED3  , 6F8.1,F8.0)
      END

