      SUBROUTINE ANITO2(ALL,NR,NSOLMX,WMOL,DOBS,PRNTMR,PRNTLN,NTST,
     1  NRRUN)
C-----------------------------------------------------------------------
C  THIS IS THE SECOND HALF OF THE MAIN PROGRAM.
C  THIS PART IS CALLED WHEN THE ZONES HAVE BEEN ESTABLISHED AND
C  COMBINATIONS HAVE BEEN MADE THAT COULD RESULT IN CORRECT LATTICES.
C-----------------------------------------------------------------------
C  .... CHANGED:  210585    ADDED SUMM OUTPUT+ /HISYM /
C  .... CHANGED:  111085    REPLACED OUTPUT FOR /HISYM/ BY *SYMPRT*
C  .... CHANGED:  271185    ADDED SRTSYM (SORTS HIGH-SYM LATTICES)
C  .... CHANGED:  080486    ADDED SOME 'NSUMM' OUTPUT
C  LAST CHANGED:  150986    REMOVED POSS. INFINITE LOOP 2000-2100
C-----------------------------------------------------------------------
      DIMENSION OLD(6)
      INTEGER  H
C-----------------------------------------------------------------------
C
C        ****************************
C        * ATTENTION   H = INTEGER  *
C        ****************************
C
C-----------------------------------------------------------------------
      LOGICAL SELF,ENL,TEST,ALL
      CHARACTER*1 BRAVTP(6)
      CHARACTER NAME*80,CARD*80,SNAME*70,INTSY*2
      COMMON /LINES / IDOC,INP,IOUT ,MAX,OBS(40),CEL(10)
      COMMON /CHARCM/ NAME,CARD,INTSY(50)
      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 /CONSTS/ NSUMM,RAD,TWORAD,RRAD,HLFLAM,WAVSQR,STORE(200,9)
      COMMON /THETAS/ ZERTOT,OBSTET(40),SHFTET(40),TETPT(40),TETMT(40)
      COMMON /HISYM / ISOL,SYMCON(12,6),NMATCH(45)
      DATA  BRAVTP(1)/'A'/,BRAVTP(2)/'B'/,BRAVTP(3)/'C'/
     *   ,  BRAVTP(4)/'I'/,BRAVTP(5)/'F'/,BRAVTP(6)/'P'/
c$$$$$$$$$
     	do 6457 i=1,nr
6457	print *, roster(i,1), roster(i,7)
      IF (ALL) THEN
         N4=NR
         GO TO 2000
      END IF
     	PRINT *,'$$$ START ANO2'
C  ---------------------------------------------------------------------
C  ALL IS  ABUSED HERE TO EXECUTE THE JUMP FROM AAITO1 TO 2000
C  WHEN COMPLETE LATTICES HAVE BEEN GIVEN AS INPUT.
C  ---------------------------------------------------------------------
      CALL NTEST(NTST,5)
     	PRINT *,'$$$ START THREED'
      CALL THREED
     	PRINT *,'$$$ END THREED'
C  ---------------------------------------------------------------------
C  THE PROCEDURE THREED TRIES TO FIND THE ANGLE BETWEEN TWO
C  INTERSECTING ZONES, DETERMINED BY THE LAST SECTION.(REAR),
C  THE POSSIBLE  LATTICES WERE ON  EVERY FOURTH POSITION IN ROSTER.
C  IF MORE THAN ONE VALUE OF D IS FOUND, THE ADDITIONAL POSSIBLE
C  LATTICES ARE STORED AFTER THE FIRST. A MAXIMUM OF FOUR DIFFERENT
C  VALUES OF D CAN THUS BE ACCOMODATED.  EMPTY POSITIONS ARE MARKED BY
C  A ZERO IN POSITION 7.
C-----------------------------------------------------------------------
      WRITE(IOUT ,9960) NAME
      N4=LINCO+4
      N=0
      DO 1910 I=1,N4
      IF (ROSTER(I,7).LT.1.)GO TO  1910
      N=N+1
      DO 1905 J=1,8
 1905 ROSTER(N,J)=ROSTER(I,J)
      WRITE(IOUT ,9937)  N,(ROSTER(N,J),J=1,8)
 1910 CONTINUE
      N4=N
C  ---------------------------------------------------------------------
C--THE LATTICES FOUND WILL BE PUT INTO A STANDARD SHAPE.(NORMALISED)
C  ---------------------------------------------------------------------
      CALL NTEST(NTST,6)
      IF (N4 .EQ. 0 ) THEN
         WRITE(NSUMM,9924)
         WRITE(IOUT ,9925)
         RETURN
      END IF
C
C  ROSTER(N,1..6) CONTAIN THE LATTICE CONSTANTS, 7 GIVES THE NUMBER
C  OF HITS AND 8 CONTAINS THE PEDIGREE.
C
      DO 1990 J=1,N4
      DO 1930 I=1,8
 1930 CEL(I)=ROSTER(J,I)
      CALL NORM3
      DO 1940 I=1,8
 1940 ROSTER(J,I)=CEL(I)
      IF (J .EQ. 1) GO TO 1990
      MJ=J-1
C  ---------------------------------------------------------------------
C  FIRST WE MARK 'EQUAL' LATTICES .
C  ---------------------------------------------------------------------
      DO 1980 K=1,MJ
      U=0.0
      DO 1950 I=1,6
 1950 U=U+ABS(ROSTER(K,I)-ROSTER(J,I))
      IF (U .GT. FMINEQ) GO TO 1990
C  ---------------------------------------------------------------------
C  IF TWO LATTICES ARE EQUAL,WE SET ROSTER(N,7)=0 AND GIVE THE
C  PEDIGREE OF THE OTHER IN ROSTER(N,9)
C  ---------------------------------------------------------------------
      IF (ROSTER(K,7) .LT. ROSTER(J,7) ) THEN
         ROSTER(K,7)=0.
         ROSTER(K,9)=ROSTER(J,8)
      ELSE
         ROSTER(J,7)=0.
         ROSTER(J,9)=ROSTER(K,8)
      END IF
 1980 CONTINUE
 1990 CONTINUE
C-----------------------------------------------------------------------
C  THE LATTICES ARE NOW NORMALIZED,AND THE PATHOLOGICAL CASES HAVE
C  BEEN MARKED. THE NEXT STEP IS TO TRY AND REFINE THE LATTICES THAT
C  ARE LEFT.
C  ---------------------------------------------------------------------
2000  PRINT *,'$$$ AT 2000'
      ALL=.FALSE.
      CALL NTEST(NTST,7)
      DO 2100 J=1,N4
     	PRINT *,'$$$ j, N4, ROSTER(J,1&J,7),FMINEQ',
     1 j,N4,ROSTER(J,1),ROSTER(J,7), FMINEQ
      IF (ROSTER(J,1) .LT. FMINEQ) GO TO 2100
      IF (ROSTER(J,7) .LT. FMINEQ) GO TO 2100
      IF (TEST) WRITE(IOUT ,9909) J
C  COUNT THE NUMBER OF TIMES THAT THIS LATTICE GOES THROUGH REFPAR AND
C  NORM3 (CHANGE 150986)
      NRNORM=0
      DO 2030 I=1,10
 2030 CEL(I)=ROSTER(J,I)
      CEL(10)=CEL(8)
 2040 CEL(7)=0.
      ZERTOT=0.
C-----------------------------------------------------------------------
C  A FIRST ROUND OF LEAST SQUARES (IN REFPAR) WILL NOW BE RUN.
C-----------------------------------------------------------------------
     	PRINT *,'$$$ into REFPAR'
      CALL REFPAR(NSOL,ALL)
     	PRINT *,'$$$ AFTER REFPAR'
      IF (CEL(7) .LT. PRNTLN) GO TO 2080
      IF (CEL(8) .LT. PRNTMR) GO TO 2080
      DO 2060 I=1,6
 2060 OLD(I)=CEL(I)
      CALL NORM3
     	PRINT *,'$$$ AFTER NORM3'
      NRNORM=NRNORM+1
      IF (CEL(1) .LT. FMINIM) GO TO 2080
      U=0.
      DO 2070 I=1,6
 2070 U=U+ABS(CEL(I)-OLD(I))
      IF (NRNORM .GT. 3) THEN
         WRITE(IOUT ,9905) J, (CEL(I),I=1,10)
         WRITE(NSUMM,9905) J, (CEL(I),I=1,10)
         GO TO 2080
      END IF
      IF (U .GT. 0.1) GO TO 2040
C  IF NORM3 CHANGES THE RECIPROCAL CELL, WE HAVE TO DETERMINE
C  THE BRAVAIS TYPE AGAIN. (AND REFINE THE NEW PARAMETERS)
C
 2080 DO 2090 K=1,10
 2090 ROSTER(J,K)=CEL(K)
      ROSTER(J,11)=ZERTOT
2100     PRINT *,'END OF 2100 LOOP'
C 2100 CONTINUE
C-----------------------------------------------------------------------
C  PRINT THE FATE OF THE DISCARDED LATTICES AND
C  COLLECT THE GOOD ONES.
C-----------------------------------------------------------------------
      NR=0
      IF (ENL) WRITE(IOUT ,9936) NAME
      DO 2200 J=1,N4
      IF (ROSTER(J,1) .GT. FMINEQ) GO TO 2110
      IF (ENL) WRITE(IOUT ,9910) J,(ROSTER(J,I),I=1,6),ROSTER(J,8)
      GO TO 2200
 2110 IF (ROSTER(J,7) .GT. FMINEQ) GO TO 2120
      IF (ENL) WRITE(IOUT ,9911) J,(ROSTER(J,I),I=1,9)
      GO TO 2200
 2120 IF (ROSTER(J,7) .GT. PRNTLN) GO TO 2150
      IF (ENL) WRITE(IOUT ,9912) J,(ROSTER(J,I),I=1,10)
      GO TO 2200
 2150 IF (ROSTER(J,8) .GT. PRNTMR) GO TO 2160
      IF (ENL) WRITE(IOUT ,9913) J,(ROSTER(J,I),I=1,10)
      GO TO 2200
C  FINALLY,THE LATTICE IS ACCEPTED
 2160 NR=NR+1
      DO 2180 K=1,12
 2180 ROSTER(NR,K)=ROSTER(J,K)
      IF (ENL) WRITE(IOUT ,9914) J,(ROSTER(NR,I),I=1,10)
 2200 CONTINUE
      IF (ENL) THEN
         WRITE(IOUT ,9915)
         WRITE(IOUT ,9916)
      END IF
C-----------------------------------------------------------------------
C  A DOUBLE SORT WILL NOW BE CARRIED OUT, FIRST ON THE NUMBER OF
C  LINES INDEXED(POS. 7)   AND SECONDLY ON THE FIGURE OF MERIT (POS. 8)
C  OF THE LATTICES.
C-----------------------------------------------------------------------
      IF (NR .EQ. 0) THEN
         WRITE(NSUMM,9924)
         WRITE(IOUT ,9990)
          RETURN
      END IF
      WRITE(IOUT ,9950) NAME
      WRITE(IOUT ,9951)
      SNAME=NAME
      WRITE(NSUMM,9980) SNAME
      JJ=100
      DO 3050 L=1,NR
      K=0
      AA=0.
      DO 3030 H=L,NR
      M= ROSTER(H,7)+0.5
      B= ROSTER(H,8)
      IF (M-K) 3030,3020,3010
 3010 K=M
      I=H
      AA=B
      GO TO 3030
 3020 IF(B.LE.AA) GO TO 3030
      AA=B
      I=H
 3030 CONTINUE
      IF (K.GE.JJ) GO TO 3040
      WRITE(IOUT ,9952)
      WRITE(NSUMM,9952)
      JJ=K
 3040 DO 3045 N=1,12
      A= ROSTER(I,N)
      ROSTER(I,N)= ROSTER(L,N)
      ROSTER(L,N)=A
      IF (N .LE. 10) CEL(N)=A
 3045 CONTINUE
C  CALCULATE VOLUME OF UNIT CEL AND PRINT
      CALL UNTCEL(WMOL,DOBS,VOL,0)
      KKK=NINT(ROSTER(L,9))
      TET=TWORAD*ROSTER(L,11)
      WRITE(IOUT ,9953) (ROSTER(L,N),N=1,8),BRAVTP(KKK),
     * ROSTER(L,10), VOL,TET
      WRITE(NSUMM,9983) (ROSTER(L,N),N=1,8),BRAVTP(KKK)
 3050 CONTINUE
C-----------------------------------------------------------------------
C--END OF SORT AND PRINT.
C-----------------------------------------------------------------------
      WRITE(IOUT ,9963)
      CALL NTEST(NTST,8)
C-----------------------------------------------------------------------
C  WE NOW PRINT OUT THE BEST *NSOLMX* SOLUTIONS WITH AN INDEXED LIST
C  OF THE FIRST TWENTY LINES.  IF YOU HAVE SET  'TEST' THEN YOU ALSO
C  GET  THE Q-VALUES PRINTED OUT AS A  LATTICE.
C-----------------------------------------------------------------------
      IF (NR .LT. NSOLMX) NSOLMX=NR
      WRITE(IOUT ,9964) NSOLMX
      N=0
      MIN=-1
      SELF=.TRUE.
      DO 3060 I=1,3
 3060 IF (NSYST(I) .NE.0) N=I
      IF (N .NE. 0) SELF=.FALSE.
      IF (SELF) GO TO 3090
C  SELF=.FALSE. IF WE WANT TO FORCE A CERTAIN CRYSTAL SYSTEM.
      N=0
      M=0
      DO 3070 K=1,3
      IF (NSYST(K) .EQ. 1) N=N+1
      IF (NSYST(K) .EQ.(-1) ) M=M+1
 3070 CONTINUE
      IF (N .GT. 0) WRITE(IOUT ,9965)
      IF (M .GT. 0) WRITE(IOUT ,9966)
      IF (N .GT. 0) MIN=1
      IF (M .GT. 0) MIN=0
      DO 3080 M=1,MAX
      DIF=OBS(M)-OBSMT(M)
      OBSMT(M)=OBSMT(M)-DIF
 3080 OBSPT(M)=OBSPT(M)+DIF
C  THE TOLERANCES HAVE BEEN DOUBLED WHEN YOU WANT TO  FORCE A CERTAIN
C  CRYSTAL SYSTEM.
      WRITE(IOUT ,9967)
 3090 NSOL=0
      DO 3200 I=1,NR
      IF (SELF) GO TO 3120
      NORDE=1
      DO 3110 J=4,6
 3110 IF (ABS(ROSTER(I,J)) .GT. FMINEQ) NORDE=NORDE+1
      IF (NORDE .EQ. 4) NORDE=3
      IF (NSYST(NORDE) .LT. MIN) GO TO 3200
 3120 NSOL=NSOL+1
C  FIND OUT WHETHER THERE ARE IDENTICAL SOLUTIONS AMONG THE BEST ONES,
C  SKIP REFINEMENT + PRINT OUT IF THERE ARE.
      IF (I .EQ. 1) GO TO 3140
      A=0.
      K=I-1
      DO 3130 J=1,6
 3130 A=A+ABS(ROSTER(K,J)-ROSTER(I,J))
      IF (A .GT. FMINEQ) GO TO 3140
C  THE LATTICE IS NEARLY THE SAME AS THE PRECEDING ONE
      WRITE(IOUT ,9949) I,K
      GO TO 3170
 3140 DO 3150 J=1,10
 3150 CEL(J)= ROSTER(I,J)
      IF (NSOL .EQ. 1) WRITE(NSUMM,9970) (CEL(K),K=1,8)
      ZERTOT=ROSTER(I,11)
      CALL REFPAR(NSOL,ALL)
      IF (NSOL .EQ. 1) WRITE(NSUMM,9970) (CEL(K),K=1,8)
      DO 3160 J=1,10
 3160 ROSTER(I,J)=CEL(J)
      ROSTER(I,11)=ZERTOT
 3170 IF (NSOL .EQ. NSOLMX) GO TO 3205
 3200 CONTINUE
     	PRINT *,'$$$ AFTER 3200'
C  AFTER THIS REFINEMENT CYCLE, A NEW SORT IS MADE.
C  THE NEW CRITERION IS: (FIGURE OF MERIT)/ (21-LINES INDEXED)
C  ROSTER:1-6 CELL CONST., 7=LINES INDEXED, 8=M20, 9-BRAVTYP,
C          10=PEDIGREE, 11=ZEROSHIFT,  12=RANKING.
 3205 NSOL=NSOLMX
      IF (NR .LT. NSOLMX) NSOL=NR
      DO 3210 I=1,NSOL
 3210 ROSTER(I,12)=ROSTER(I,8)/(21.0-ROSTER(I,7) )
      DO 3250 I=1,NSOL
      AA=0.
      DO 3230 J=I,NSOL
         IF (ROSTER(J,12) .GT. AA) THEN
            AA=ROSTER(J,12)
            L=J
         END IF
 3230 CONTINUE
      DO 3240 N=1,12
      A=ROSTER(I,N)
      ROSTER(I,N)=ROSTER(L,N)
      ROSTER(L,N)=A
 3240 CONTINUE
 3250 CONTINUE
C  PRINT FINAL TABLE OF LATTICES
      WRITE(IOUT ,9954)  NAME,NSOL
      WRITE(NSUMM,9984) SNAME,NSOL
      DO 3260 I=1,NSOL
      TET=TWORAD*ROSTER(I,11)
      KKK=NINT( ROSTER(I,9) )
      WRITE(IOUT ,9955) ( ROSTER(I,N),N=1,8),  ROSTER(I,10),
     1   BRAVTP(KKK),TET
      WRITE(NSUMM,9983) ( ROSTER(I,N),N=1,8), BRAVTP(KKK)
 3260 CONTINUE
      CALL NTEST(NTST,9)
C-----------------------------------------------------------------------
C--PRINT THE VALUES OF THE PARAMETERS OF THE DIRECT LATTICE.
C-----------------------------------------------------------------------
      WRITE(IOUT ,9956)
      WRITE(NSUMM,9956)
     	WRITE(20,99561)SNAME,NSOL
99561 FORMAT(///1H0,A70/ 35(2H--) //
     1 20X,'THE ',I2, 25H MOST PROBABLE SOLUTIONS ////
     2 55X,'LINES FIGUR LAT' /
     * 60H     A        B        C       ALFA     BETA     GAMMA INDXD,
     4 ' MERIT TYP'  )
      DO 3280 I=1,NSOL
      DO 3270 J=1,10
 3270 CEL(J)= ROSTER(I,J)
      CALL UNTCEL(WMOL,DOBS,VOL,1)
      WRITE(NSUMM,9988) (CEL(J),J=1,7)
      WRITE(20,99880) (CEL(J),J=1,6), (ROSTER(I,N),N=7,8), BRAVTP(KKK)
99880 FORMAT(1H ,6F9.4,F6.0,F6.1,2X,A1)
 3280 CONTINUE
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C  THE BEST SOLUTIONS ARE TESTED ON HIGHER SYMMETRIES,THAT MAY BE
C  CONCEALED IN THE Q-VALUES.
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      DO 3320  ISOL=1,NSOL
      DO 3310 J=1,10
 3310 CEL(J)=ROSTER(ISOL,J)
      CALL SYMTST
 3320 CONTINUE
     	PRINT *,'$$$ AFTER 3320'
      WRITE(IOUT ,9957)
      WRITE(IOUT ,9959)
      WRITE(NSUMM,9959)
      WRITE(IOUT ,9958) NAME
C  ---------------------------------------------------------------------
C  THE BEST SOLUTION IS NOW USED TO TRY AND INDEX ALL LINES ENTERED.
C-----------------------------------------------------------------------
      IF (ROSTER(1,7) .LT. 17.9) RETURN
      IF (ROSTER(1,8) .LT. 10.0) RETURN
      CALL SRTSYM
      IF ((SYMCON(1,3).GT.0.1) .AND. (SYMCON(1,4) .GT. 0.9*MAX))  THEN
         CALL SYMPRT
      ELSE
         ALL=.TRUE.
         DO 3400 I=1,10
 3400    CEL(I)=ROSTER(1,I)
         ZERTOT=ROSTER(1,11)
         WRITE(NSUMM,9970) (CEL(K),K=1,8)
         CALL REFPAR(1,ALL)
         WRITE(NSUMM,9970) (CEL(K),K=1,8)
         WRITE(NSUMM,9971) TWORAD*ZERTOT
         ROSTER(1,11)=ZERTOT
         WRITE(IOUT ,9956)
         WRITE(NSUMM,9956)
         CALL UNTCEL(WMOL,DOBS,VOL,1)
         WRITE(NSUMM,9988) (CEL(J),J=1,7)
         WRITE(NSUMM,'(/,1X,A70)' ) SNAME
      END IF
      NRRUN=2
     	PRINT *,'$$$ AT END'
      RETURN
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
 9905 FORMAT(/' LATTICE',I3,' SWITCHES BACK AND FORTH BETWEEN NORM3 ',
     * 'AND REFINE.' / ' POSSIBLY REFINE CHANGES THE ORDER OF THE ',
     * 'CELL CONSTANTS AND NORM3'/ ' SWITCHES THEM BACK.   THE CELL ',
     * 'IS ACCEPTED AS:' / 6F10.2 / 4F8.0 )
 9909 FORMAT(//,' LATTICE NR.',I3,'  WILL NOW BE REFINED AND',
     * ' NORMALISED',20(2H==)  )
 9910 FORMAT(1H ,I5,6F8.1,24X,F8.0,20H DISCARDED BY NORM3  )
 9911 FORMAT(1H ,I5,6F8.1,F8.0,16X,F8.0,23H IS THE SAME AS LATTICE,
     *  F6.0)
 9912 FORMAT(1H ,I5,6F8.1,4F8.0,25H NOT ENOUGH LINES INDEXED  )
 9913 FORMAT(1H ,I5,6F8.1,4F8.0,28H FIGURE OF MERIT IS TOO LOW.  )
 9914 FORMAT(1H ,I5,6F8.1,4F8.0, 20H LATTICE ACCEPTED.   )
 9915 FORMAT(///47H0***  NOTES  ***    THE SUBROUTINE ** THREED **,
     1       50HCOLLECTS TOO MANY LATTICES (BETTER TOO MANY THAN T,
     2       12HOO LITTLE).    /
     3       50H SOME OF THESE MUST BE DISCARDED BECAUSE THE LATTI,
     4       50HCE IS IMPOSSIBLE (E.G.CONTAINS NEGATIVE Q-VALUES),/
     5       50H OTHERS MUST BE DISCARDED BECAUSE THE ALLEGED SYMM,
     6       50HETRY REDUCES THE RECIPROCAL CELL EDGES BELOW THE  /
     7       50H MINIMUM ALLOWED.                                 )
 9916 FORMAT(50H     AFTER NORMALISING, SOME LATTICES TURN OUT TO ,
     9       50HBE THE SAME. ONLY ONE OF THESE IS RETAINED. THESE /
     A       50H LATTICES ARE INDICATED.                          /
     B       50H     FINALLY, THE PROGRAM TRIES TO INDEX THE PATTE,
     C       52HRN WITH THESE LATTICES AND REFINES THEIR PARAMETERS./
     D       50H IF A LATTICE DOES NOT INDEX ENOUGH LINES OR IF IT,
     E       51HS FIGURE OF MERIT IS BELOW ALL HOPE FOR IMPROVEMENT/
     F       50H THAT LATTICE IS DROPPED TOO.                       )
 9924 FORMAT(/,' NOT ONE ACCEPTABLE LATTICE HAS BEEN FOUND. ' )
 9925 FORMAT(50H NO LATTICES HAVE BEEN FOUND. TRY AGAIN WITH E.G. ,
     1       50HLARGER TOL3 (COL. 16-20).  THIS PROBLEM IS SKIPPED)
 9936 FORMAT(1H1,20X,A80 / 40(3H---)/
     *60H0               LATTICES REDUCED AND PUT INTO STANDARD FORM //
     *56X,24HLINES  FIGURE    BRAVAIS/
     *56H SERIAL   A       B       C        D       E       F     ,
     *34HINDXD  OF MERIT  LATTICE  PEDIGREE   )
 9937 FORMAT(1H ,I5,6F8.1, 2F8.0)
 9949 FORMAT(9H0SOLUTION,I3,35H IS NOT PRINTED, IT IS THE SAME AS ,
     2 9HSOLUTION ,I2)
 9950 FORMAT(1H1,20X,A80/ 40(3H---)  )
 9951 FORMAT(1H0,14X,47HTRIAL LATTICES AFTER INDEXING AND LS REFINEMENT
     *//14X,52HONLY THE FIRST 20 LINES ARE CONSIDERED FOR INDEXING  //
     *59X,24HLINES  FIGURE    BRAVAIS/
     *95H      A        B        C        D         E        F    INDEXE
     *D  OF MERIT  LATTICE  PEDIGREE   ,18HVOLUME  ZEROSHIFT  )
 9952 FORMAT(1H )
 9953 FORMAT(1H ,6F9.1,F9.0,F9.1,8X,A1,F9.0,F8.1,F10.4 )
 9954 FORMAT(1H1,20X,A80/ 40(3H---) // 30X, 4HTHE ,I2,
     1        25H MOST PROBABLE SOLUTIONS //// 58X,
     2 34HLINES  FIGURE OF           BRAVAIS/
     3 60H    Q(A)     Q(B)     Q(C)      Q(D)     Q(E)     Q(F)   IND,
     4 44HEXED   MERIT   PEDIGREE  LATTICE  ZEROSHIFT  )
 9955 FORMAT(1H ,6F9.2,F9.0,F9.1,F9.0,8X,A1,F10.4)
 9956 FORMAT(//1H ,20X,38HTHE DIRECT CONSTANTS OF THESE LATTICES//
     * 60H     A        B        C       ALFA     BETA     GAMMA      ,
     * 6HVOLUME)
 9957 FORMAT(1H0///40X,24H ****   ATTENTION   ****//
     160H 1. ANY *SOLUTION* WITH A FIGURE OF MERIT LESS THAN 4.0 IS W,
     260HORTHLESS AND SHOULD BE REJECTED RIGHT AWAY.                 //
     360H 2. ANY *SOLUTION*THAT LEAVES MORE THAN TWO VERY WEAK LINES ,
     460HUNEXPLAINED  IS WORTHLESS.   HOWEVER, IF THE FIGURE OF MERIT/
     560H    IS GREATER THAN 10.0 IT MIGHT BE WORTHWHILE TO EXAMINE  ,
     660HTHE INPUT DATA  AND THE *SOLUTION* MORE CLOSELY.            //
     760H 3. EVEN *SOLUTIONS* THAT INDEX ALL LINES WITH A FIGURE OF M,
     860HERIT GREATER THAN 10.0 SHOULD NOT BE ACCEPTED UNCRITICALLY  /
     933H    WITHOUT FURTHER INVESTIGATION//)
 9958 FORMAT(1H0,A80)
 9959 FORMAT(// 1X,23(3H***) / 2H *,67X,1H* /,
     1 ' * THIS PROGRAM GIVES A NUMBER OF SUGGESTIONS THAT MAY HELP',
     2 ' YOU TO   *' /' * FIND THE UNIT CELL THAT EXPLAINS YOUR',
     3 ' EXPERIMENTAL DATA. HOWEVER  *'/2H *,67X,'*'/
     4  2H *,5X,19(3HWWW),6H     * /
     4  ,8H *     W,55X,7HW     */ 8H *     W,
     5  '  IT IS ENTIRELY THE USERS RESPONSIBILITY TO DECIDE    W     *'
     6 ,/ 8H *     W,'  WHETHER ANY OF THE SUGGESTED UNIT CELLS IS THE',
     7  7X,'W     *' / ' *     W   CORRECT CELL. ',38X,'W     *' /
     8  8H *     W, 55X, 7HW     * / 7H *      ,
     9 19(3HWWW), 6H     * / 2H *,67X,1H*/1H ,23(3H***) )
 9960 FORMAT(1H1,20X,A80 / 40(3H---)  //
     160H POSSIBLE LATTICES BASED ON THE ABOVE COMBINATION OF ZONES  //
     258H SERIAL    A       B       C        D       E       F       ,
     315HHITS  PEDIGREE   )
 9963 FORMAT(7H0NOTES //
     160H  1. *LINES INDEXED* REFERS TO THE NUMBER OF LINES FROM THE ,
     260HFIRST 20 OBSERVED LINES THAT HAVE BEEN ACCEPTED AS INDEXED  /
     360H     BY THIS LATTICE, UP TO A MAXIMUM OF 20, WHICH WOULD REP,
     460HRESENT COMPLETE SUCCESS.                                    //
     560H  2. *FIGURE OF MERIT* IS DEFINED AS QOBS(20)/(2*MEANDIFF*NC,
     660HALC) WHERE DIFFERENCES ARE TAKEN FOR INDEXED LINES ONLY     /
     760H     AND *NCALC* IS THE NUMBER OF LINES GENERATED UP TO Q=1.,
     860H003*Q(20),REGARDLESS OF WHETHER OR NOT THEY CONTRIBUTE TO   /
     960H     THE INDEXINGS, BUT EXCLUDING SYSTEMATIC ABSENCES DUE TO,
     *27H THE BRAVAIS LATTICE TYPE. /
     A60H     IT IS CLOSELY RELATED TO  M20 ,DE WOLFF(J.APPL.CRYST.19,
     B60H68,1972), ALLOWING FOR LESS THAN 20 SUCCESSFUL INDEXINGS.   )
 9964 FORMAT(4H0THE, I3,
     *       50H BEST SOLUTIONS ARE NOW PRINTED MORE EXTENSIVELY.   )
 9965 FORMAT(50H0THE CHANCES FOR SOME LATTICES ARE ENHANCED. ONLY ,
     *       50HTHOSE ARE PRINTED OUT.                            )
 9966 FORMAT(50H0 THE PRINTING OF SOME  LATTICES IS SUPPRESSED    )
 9967 FORMAT(50H0FOR THE LATTICES IN THE CRYSTAL SYSTEMS THAT ARE ,
     *       50HTREATED PREFERENTIALLY, THE TOLERANCE HAS BEEN    /
     2       54H DOUBLED AND TWO ROUNDS OF LEAST SQUARES ARE EXECUTED.,
     3       50H AFTERWARDS, THESE SOLUTIONS ARE SORTED WITH      /
     4       50H THE OTHER SOLUTIONS ON LINES INDEXED AND FIGURE  ,
     5       50HOF MERIT.                                         /
     6       50H IF YOUR PREFERRED SYSTEM DOES NOT COME OUT ON TOP,
     7       54H THERE MIGHT BE SOMETHING AMISS WITH YOUR ASSUMPTIONS.)
 9970 FORMAT(//'    SOLUTION NR. 1  '/ 1X,6F9.2,F6.0,F6.1  )
 9971 FORMAT(' ZEROSHIFT=',F9.4 / )
 9980 FORMAT(////1H ,A70/ 35(2H--)  //
     1       1H ,14X,47HTRIAL LATTICES AFTER INDEXING AND LS REFINEMENT
     2//14X,52HONLY THE FIRST 20 LINES ARE CONSIDERED FOR INDEXING  //
     3 55X,'LINES FIGUR LAT' /
     4 '      A        B        C        D         E        F  INDXD',
     5 ' MERIT TYP'  )
 9983 FORMAT(1H ,6F9.1,F6.0,F6.1,2X,A1)
 9984 FORMAT(///1H0,A70/ 35(2H--) //
     1 20X,'THE ',I2, 25H MOST PROBABLE SOLUTIONS ////
     2 55X,'LINES FIGUR LAT' /
     3 '    Q(A)     Q(B)     Q(C)     Q(D)      Q(E)     Q(F) INDXD',
     4 ' MERIT TYP'  )
 9988 FORMAT(1H ,6F9.4,F12.2)
 9990 FORMAT(/52H NOT ONE ACCEPTABLE LATTICE HAS BEEN FOUND FOR THIS
     2     , 52H PROBLEM. USUALLY THIS MEANS THAT YOUR DATA ARE NOT
     3     / 50H GOOD ENOUGH. YOU MIGHT ALSO TRY OTHER PARAMETERS.  /
     4       50H0THIS PROBLEM IS SKIPPED.                           )
      END

      SUBROUTINE ARRAY(MODE,I,J,N,M,S,D)
C
C  COPIED FROM  THE I.B.M.  SCIENTIFIC SUBROUTINE PACKAGE.
C  TO BE USED IN COMBINATION WITH THE SUBROUTINE MINV (MATRIX INVERSION
C
      DIMENSION S(64),D(64)
      NI=N-I
      IF(MODE .GT. 1) GO TO 20
      IJ=I*J+1
      NM=N*J+1
      DO 10 K=1,J
      NM=NM-NI
      DO 10 L=1,I
      IJ=IJ-1
      NM=NM-1
      D(NM)=S(IJ)
   10 CONTINUE
      GO TO 40
   20 IJ=0
      NM=0
      DO 30 K=1,J
      DO 25 L=1,I
      IJ=IJ+1
      NM=NM+1
      S(IJ)=D(NM)
   25 CONTINUE
      NM=NM+NI
   30 CONTINUE
   40 RETURN
      END

      SUBROUTINE LISTM
C  THE SUBROUTINE PRODUCES A LIST OF ALL OBSERVED AND CALCULATED
C  LINES, SORTED ON TWOTHETA. IF THE SUCCEEDING VALUES DIFFER
C  MORE THAN THE TOLERANCE(TOLG), A BLANK LINE IS INSERTED.
C
C  LAST CHANGE:  19 JULI 1985.  FORMAT 9914
C
      LOGICAL SELF,ENL,TEST
      CHARACTER NAME*80,CARD*80,INTSY*2
      COMMON /LINES / IDOC,INP,IOUT,MAX,OBS(40),CEL(10)
      COMMON /CHARCM/ NAME,CARD,INTSY(50)
      COMMON /PARMS / TOL2,TOL3,WAVEL,FMINEQ,FMINIM,SELF,ENL,TEST,LIST,
     *                NSYST(3),OBSMT(40),OBSPT(40),TWTHET(40),TOLG
      COMMON /LINSTR/ NMAX,NCALC,N20,NOTIND,INSTOR(200,3),QSTOR(200)
      COMMON /CONSTS/ NSUMM,RAD,TWORAD,RRAD,HLFLAM,WAVSQR,STORE(200,9)
      COMMON /THETAS/ ZERTOT,OBSTET(40),SHFTET(40),TETPT(40),TETMT(40)
      COMMON /HISYM / ISOL,SYMCON(12,6),NMATCH(45)
C-----------------------------------------------------------------------
      SWAVEL=WAVEL/200.
C--SORT CALCULATED VALUES.
      DO 100 I=1,NCALC
      Q=1.E+6
      DO  60 J=I,NCALC
      IF (QSTOR(J).GT. Q) GO TO 60
      Q=QSTOR(J)
      N=J
   60 CONTINUE
      U=QSTOR(I)
      QSTOR(I)=QSTOR(N)
      QSTOR(N)=U
      DO  80 K=1,3
      INX  =INSTOR(I,K)
      INSTOR(I,K)=INSTOR(N,K)
   80 INSTOR(N,K)= INX
  100 CONTINUE
      NO=1
      NC=1
      UL=0.
      WRITE(IOUT,9903) NMAX,TOLG
      WRITE(IOUT,9904)
      WRITE(NSUMM,9914) NMAX,TOLG
C--PRINT QCALC AND QOBS IN ORDER OF TWOTHETA.
  110 CONTINUE
      IF (NO .GT. NMAX) GO TO 115
      IF (NC .GT. NCALC) GO TO 200
      IF (QSTOR(NC) .GT. OBS(NO)) GO TO 200
  115 IF (NC .GT. NCALC) GO TO 300
      W=SQRT(QSTOR(NC))
      U=SWAVEL*W
      W=100./W
      U=TWORAD*ATAN(U/SQRT(1.0-U*U))
      IF (U-UL .GT. TOLG) WRITE(IOUT,9905)
      IF (U-UL .GT. TOLG) WRITE(NSUMM,9905)
      WRITE(IOUT ,9901) U,W,(INSTOR(NC,K),K=1,3),QSTOR(NC)
      WRITE(NSUMM,9901) U,W,(INSTOR(NC,K),K=1,3),QSTOR(NC)
      UL=U
      NC=NC+1
      GO TO 110
  200 TET=SHFTET(NO)
      U=TWORAD*TET
      D=WAVEL/(2.*SIN(TET))
      Q=100./D
      Q=Q*Q
      IF (U-UL .GT. TOLG) THEN
         WRITE(IOUT ,9905)
         WRITE(NSUMM,9905)
      END IF
      IF (NMATCH(NO) .EQ. 0) THEN
         WRITE(IOUT ,9902) U,D, Q ,INTSY(NO)
         WRITE(NSUMM,9902) U,D, Q ,INTSY(NO)
      ELSE
         WRITE(IOUT ,9906) U,D, Q ,INTSY(NO)
         WRITE(NSUMM,9906) U,D, Q ,INTSY(NO)
      END IF
      UL=U
      NO=NO+1
      GO TO 110
  300 RETURN
 9901 FORMAT(1H ,F10.2,F10.3,2X,3I4,F8.1)
 9902 FORMAT(1H ,F10.2,F10.3,14H     *OBSERVED ,F8.1,5X,A2,
     * '  XXXXXXXX')
 9903 FORMAT(52H0LIST OF ALL OBSERVED AND CALCULATED LINES UP TO AND,
     *  14H INCLUDING THE, I4,3HST.  /
     1  34H IF THE TWOTHETAS DIFFER MORE THAN  ,F6.3,
     2  37H  DEGREES, A BLANK LINE IS INSERTED.  )
 9904 FORMAT(45H0  TWOTHETA       D       H   K   L     Q          )
 9905 FORMAT(3X)
 9906 FORMAT(1H ,F10.2,F10.3,14H     *OBSERVED ,F8.1,5X,A2)
 9914 FORMAT(/// 5X,16('***')/ 5X,
     1  '*   USE THE LIST TO FIND SYSTEMATIC ABSENCES   *'/
     2  5X,16('***') //
     3  ' LIST OF ALL OBSERVED AND CALCULATED LINES UP TO THE',
     4  I4,3HST.  /
     5  34H IF THE TWOTHETAS DIFFER MORE THAN  ,F6.3,
     6  '  DEGREES, '/ ' A BLANK LINE IS INSERTED. ' //
     7  45H   TWOTHETA       D       H   K   L     Q          )
      END

      SUBROUTINE LSTET(TYP,A,B,FMAAS, ZERSHF,QCALMX, NR)
C-----------------------------------------------------------------------
C  REFINES THE PARAMETERS (A,B,FMAAS) BY LEAST-SQUARES, PRINTS THE ZONE
C  IF TYP=TRUE  AND FILLS THE ARRAY EVAL FOR THE PROCEDURE  CRITERION.
C  THE ELEMENTS OF  EVAL  REMAIN  ZERO  WHEN Q IS TOO LARGE (Z >QCALMX),
C  THEY BECOME +1 WHEN THE CALCULATED Q-VALUE EXISTS (CORRESPONDS TO
C  AN OBSERVED Q) AND BECOME  -1  IF THE Q-VALUE IS NOT IN THE LIST
C  OF OBSERVED VALUES.VALUES ARE NOW + OR - THE PROBABILITY.
C-----------------------------------------------------------------------
      DIMENSION  REF(4,4),HULP(6),FYN(4),WAR(4),IHELP1(8),IHELP2(8)
      INTEGER  H,HG,HMIN, ORDE,ORP1,ORP2
      LOGICAL SELF,ENL,TEST,RECTAN,TYP
      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 /EVALN / NORDER,RECTAN,PROB(40),EVAL(-8:8,0:8),
     *                UPBNDR(40)
      COMMON /TEMRES/ LINCO,NQ2,NINDEX,INDEX(40),ROSTER(150,12)
      COMMON /CONSTS/ NSUMM,RAD,TWORAD,RRAD,HLFLAM,WAVSQR,STORE(200,9)
      COMMON /THETAS/ ZERTOT,OBSTET(40),SHFTET(40),TETPT(40),TETMT(40)
      DO 20 H=1,4
      FYN(H)=0.
      WAR(H)=0.
      DO 10 K=1,4
   10 REF(H,K)=0.
   20 CONTINUE
      DO 30 H=-8,8
      DO 25 K=0,8
   25 EVAL(H,K)=0.
   30 CONTINUE
      DO 36 J=1,40
   36 INDEX(J)=0
      NINDEX=0
C-----------------------------------------------------------------------
      HG=NINT( SQRT(QCALMX/A) )
      KG=NINT( SQRT(QCALMX/B) )
      IF (HG .GT. 8) HG=8
      IF (KG .GT. 8) KG=8
      IF (FMAAS .LT. 0.1) THEN
         FMAAS=0.0
         HMIN= 0
         RECTAN=.TRUE.
         ORDE=2
      ELSE
         HMIN=-HG
         RECTAN= .FALSE.
         ORDE=3
      END IF
      IF (TYP)  WRITE(IOUT,1)NR
C--WRITES HEADING,  THE CALCULATION OF Q-VALUES STARTS HERE.
      ORP1=ORDE+1
      ORP2=ORDE+2
      DO 70 H=HMIN,HG
         ISPAT=0
         X=A*FLOAT(H*H)
         Y=FMAAS*FLOAT(H)
         DO 60 K=0,KG
            ISPAT=ISPAT+1
            IF (K .EQ. 0  .AND.  H .LE. 0 ) GO TO 60
            Q=X+B*FLOAT(K*K)+Y*FLOAT(K)
            IF (Q .GT.QCALMX) GO TO 65
            IF (TYP) CALL WRISP(ISPAT, Q,0)
            DO 53 I=1,MAX
               IF (Q .GT. UPBNDR(I) ) GO TO 53
               EVAL(H,K)=-PROB(I)
               GO TO 55
   53       CONTINUE
C           WE HAVE DETERMINED THE PROBABILITY.
C           NOW SEE WHETHER THERE IS A FITTING OBSERVED LINE.
   55       Z=TETCAL(Q)
            DO 59 J=1,MAX
               IF (TETPT(J) .LT. Z) GO TO 59
               IF (TETMT(J) .GT. Z) GO TO 60
C           THERE IS A FITTING LINE. USE IT FOR THE REFINEMENT.
               INDEX(J)=1
               EVAL(H,K)=-EVAL(H,K)
               IF (TYP)   CALL WRISP(ISPAT,Q,1)
               C=WAVSQR/SIN(2.*Z)
               HULP(1)=1.
               HULP(2)=H*H*C
               HULP(3)=K*K*C
               HULP(4)=H*K*C
               HULP(ORP2)=OBSTET(J)-Z
               DO 58 N2=1,ORP1
                  WAR(N2)=WAR(N2)+HULP(N2)*HULP(ORP2)
                  DO  57 N1=1,ORP1
   57             REF(N1,N2)=REF(N1,N2)+HULP(N1)*HULP(N2)
   58          CONTINUE
               GO TO 60
   59       CONTINUE
C  59       IS THE END OF THE LOOP ON THE SEARCH FOR A FITTING  Q
   60    CONTINUE
C  60    IS THE END OF THE LOOP ON K
   65    IF (TYP) CALL WRISP(ISPAT, Z,2)
   70 CONTINUE
C  END OF LOOP ON H (I,NH)
      IF (TEST) THEN
         L=-8
         IF (RECTAN) L=0
         WRITE(IOUT,6)  ((EVAL(I,J),J=0,8),I=8,L,-1)
      END IF
      DO 72 J=1,NQ2
   72 NINDEX=NINDEX+INDEX(J)
C-----------------------------------------------------------------------
C  COUNTS THE NUMBER OF LINES THAT ARE INDEXED  AMONGST THE FIRST
C  NQ2 LINES.
C  REFINEMENT-CORRECTIONS ON ZONE CONSTANTS WILL BE CALCULATED FROM
C  HERE TO 170
C-----------------------------------------------------------------------
      CALL ARRAY(2,ORP1,ORP1,4,4,REF,REF)
C  RLS
      CALL MINV(REF,ORP1,DETER,IHELP1,IHELP2)
      CALL ARRAY(1,ORP1,ORP1,4,4,REF,REF)
      IF (ABS(DETER) .LT. 1.0E-6) THEN
         IF (TYP) WRITE(IOUT, 3) NR
         GO TO 200
      END IF
      DO 170 I=1,ORP1
         DO 160 J=1,ORP1
  160    FYN(I)=FYN(I)+REF(I,J)*WAR(J)
  170 CONTINUE
      ZERSHF=FYN(1)
      A=A+FYN(2)*1.E+4
      B=B+FYN(3)*1.E+4
      ZERO=ZERSHF*TWORAD
      IF (ORDE .EQ. 3) FMAAS=FMAAS+FYN(4)*1.E+4
      IF (TYP) WRITE(IOUT,4) NR,A,B,FMAAS,ZERO
      IF (TEST)WRITE(IOUT,5) NINDEX,NQ2
  200 RETURN
    1 FORMAT(10H0 ZONE NR ,I2,27H PRINTED OUT FOR INSPECTION /
     *42H   *A*  UP THE PAGE,  *B* ACROSS THE  PAGE  )
C   2 FORMAT(1H0)
    3 FORMAT(' ZONE ',I3,'  NOT REFINED')
    4 FORMAT(5H ZONE,I3,18H AFTER REFINEMENT , 3F10.2,
     1 '    ZERO SHIFT=',F8.3,' DEGREES TWOTHETA'  )
    5 FORMAT(1H0,  I3,
     * 21H LINES FROM THE FIRST,I3,26H ARE INDEXED BY THIS ZONE       )
    6 FORMAT(11H0ARRAY EVAL / (1H ,9E12.3)  )
      END

      SUBROUTINE MINV(A,N,D,L,M)
C
C  MATRIX INVERSION. COPIED FROM THE I.B.M. SCIENTIFIC SUBROUTINE PACKAG
C
      DIMENSION A(64),L(8),M(8)
      D=1.0
      NK=-N
      DO 80 K=1,N
      NK=NK+N
      L(K)=K
      M(K)=K
      KK=NK+K
      BIGA=A(KK)
      DO 20 J=K,N
      IZ=N*(J-1)
      DO 20 I=K,N
      IJ=IZ+I
   10 IF(ABS(BIGA).GE.ABS(A(IJ))) GO TO 20
      BIGA=A(IJ)
      L(K)=I
      M(K)=J
   20 CONTINUE
      J=L(K)
      IF(J .LE. K) GO TO 35
      KI=K-N
      DO 30 I=1,N
      KI=KI+N
      HOLD=-A(KI)
      JI=KI-K+J
      A(KI)=A(JI)
      A(JI)=HOLD
   30 CONTINUE
   35 I=M(K)
      IF(I .LE. K) GO TO 45
      JP=N*(I-1)
      DO 40 J=1,N
      JK=NK+J
      JI=JP+J
      HOLD=-A(JK)
      A(JK)=A(JI)
      A(JI)=HOLD
   40 CONTINUE
   45 IF (ABS(BIGA) .GT. 0.0001) GO TO 48
      D=0.0
      RETURN
   48 DO 55 I=1,N
      IF(I .EQ. K) GO TO 55
      IK=NK+I
      A(IK)=A(IK)/(-BIGA)
   55 CONTINUE
      DO 65 I=1,N
      IK=NK+I
      HOLD=A(IK)
      IJ=I-N
      DO 65 J=1,N
      IJ=IJ+N
      IF(I .EQ. K) GO TO 65
      IF(J .EQ. K) GO TO 65
      KJ=IJ-I+K
      A(IJ)=HOLD*A(KJ)+A(IJ)
   65 CONTINUE
      KJ=K-N
      DO 75 J=1,N
      KJ=KJ+N
      IF(J .EQ. K) GO TO 75
      A(KJ)=A(KJ)/BIGA
   75 CONTINUE
      D=D*BIGA
      A(KK)=1.0/BIGA
   80 CONTINUE
      K=N
  100 K=(K-1)
      IF(K .LE. 0) GO TO 150
      I=L(K)
      IF(I .LE. K) GO TO 120
      JQ=N*(K-1)
      JR=N*(I-1)
      DO 110 J=1,N
      JK=JQ+J
      HOLD=A(JK)
      JI=JR+J
      A(JK)=-A(JI)
      A(JI)=HOLD
  110 CONTINUE
  120 J=M(K)
      IF(J .LE. K) GO TO 100
      KI=K-N
      DO 130 I=1,N
      KI=KI+N
      HOLD=A(KI)
      JI=KI-K+J
      A(KI)=-A(JI)
      A(JI)=HOLD
  130 CONTINUE
      GO TO 100
  150 RETURN
      END

      SUBROUTINE SYMPRT
C  SYMPRT WILL PRINT THE LIST OF ALL LINES OF A SYMMETRIC LATTICE,
C  E.G. A TETRAGONAL,HEXAGONAL OR CUBIC LATTICE.
      CHARACTER NAME*80,CARD*80,SNAME*70,INTSY*2
      COMMON /LINES / IDOC,INP,IOUT,MAX,OBS(40),CEL(10)
      COMMON /CHARCM/ NAME,CARD,INTSY(50)
      COMMON /CONSTS/ NSUMM,RAD,TWORAD,RRAD,HLFLAM,WAVSQR,STORE(200,9)
      COMMON /HISYM / ISOL,SYMCON(12,6),NMATCH(45)
      SNAME=NAME
      WRITE(IOUT,9950) NAME
      WRITE(NSUMM,9982) SNAME
      ISOL=1
      CALL REFSYM
      CALL LISTM
      QX=SYMCON(1,1)
      QZ=SYMCON(1,2)
      NSYM=NINT(SYMCON(1,3) )
      A=100./SQRT(QX)
      C=100./SQRT(QZ)
      IF (NSYM .EQ. 3) THEN
         WRITE(IOUT ,9943) QX,A
         WRITE(NSUMM,9943) QX,A
      ELSE IF (NSYM .EQ. 4) THEN
         WRITE(IOUT ,9944) QX,QZ,A,C
         WRITE(NSUMM,9944) QX,QZ,A,C
      ELSE IF (NSYM .EQ. 5) THEN
        A=200./SQRT(3.*QX)
        WRITE(IOUT ,9945) QX,QZ,A,C
        WRITE(NSUMM,9945) QX,QZ,A,C
      ELSE IF (NSYM .EQ. 6) THEN
        A=200./SQRT(3.*QX)
        WRITE(IOUT ,9946) QX,QZ,A,C
        WRITE(NSUMM,9946) QX,QZ,A,C
      END IF
      WRITE(IOUT ,9947) SYMCON(1,6), SYMCON(1,5)
      WRITE(NSUMM,9947) SYMCON(1,6), SYMCON(1,5)
      WRITE(NSUMM,9960) SNAME
      WRITE(IOUT ,9961) NAME
      WRITE(IOUT,9948)
      WRITE(NSUMM,9948)
      RETURN
 9943 FORMAT(2X/' CUBIC  INDEXING '/
     1       ' THE RECIPROCAL CONSTANT (QX) IS  ',F10.2 /
     2  ' THE DIRECT CONSTANT IS    A=',F10.4 )
 9944 FORMAT(2X/' TETRAGONAL INDEXING '/
     1       ' RECIPROCAL CONSTANTS (QX,QZ) ARE ',2F10.2 /
     2  ' THE DIRECT CONSTANTS ARE  A=',F10.4,'  C=',F10.4 )
 9945 FORMAT(2X/' RHOMBOHEDRAL INDEXING ON HEXAGONAL AXES '/
     1       ' RECIPROCAL CONSTANTS (QX,QZ) ARE ',2F10.2 /
     2  ' THE DIRECT CONSTANTS ARE  A=',F10.4,'  C=',F10.4 )
 9946 FORMAT(2X/' HEXAGONAL INDEXING  '/
     1       ' RECIPROCAL CONSTANTS (QX,QZ) ARE ',2F10.2 /
     2  ' THE DIRECT CONSTANTS ARE  A=',F10.4,'  C=',F10.4 )
 9947 FORMAT(' VOLUME = ',F9.1 ,'   ZEROSHIFT=',F8.4 )
 9948 FORMAT(// 5X,12(4HWWWW)/
     *  '     W   WARNING:  THE TRUE CELL CAN BE A           W' /
     1  '     W           SIMPLE SUB- OR SUPERLATTICE        W' /
     2  '     W             OF THE CELL GIVEN HERE.          W' /
     3  5X,12(4HWWWW)  )
 9950 FORMAT(1H1,20X,A80/ 40(3H---) )
 9960 FORMAT(//,1X,A70)
 9961 FORMAT(//,1X,A80)
 9982 FORMAT(2X///1H ,A70/35(2H--)/' INDEXED ON HIGH-SYMMETRY LATTICE')
      END

      SUBROUTINE SRTSYM
*  SORTS THE POSSIBLE HIGH-SYMMETRICAL LATTICES ON THE NUMBER
*  OF INDEXED LINES.
*  IF THE NUMBER OF INDEXED LINES IS THE SAME, CHOOSE THE CELL
*  WITH THE SMALLEST VOLUME.
*  PUTS THE BEST LATTICE IN POSITION 1.
*
      COMMON /HISYM / ISOL,SYMCON(12,6),NMATCH(45)
      NINDXD=0
      N=0
      DO 100 I=1,12
         IF (SYMCON(I,3) .LT. 0.1 ) GO TO 100
         N=N+1
         IND=NINT( SYMCON(I,4) )
         IF (IND .LT. NINDXD) THEN
            GO TO 100
         ELSE IF (IND .GT. NINDXD) THEN
            NINDXD=IND
            KBEST=I
         ELSE IF (IND .EQ. NINDXD) THEN
            IF (SYMCON(I,6) .LT. SYMCON(KBEST,6) ) KBEST=I
         END IF
  100 CONTINUE
      IF (N .GT. 0) THEN
         DO 120 J=1,6
         SYMCON(1,J)=SYMCON(KBEST,J)
  120 CONTINUE
      END IF
      RETURN
      END

      SUBROUTINE REFSYM
C-----------------------------------------------------------------------
C  REFINES THE PARAMETERS (QX,QZ) BY LEAST-SQUARES,
C-----------------------------------------------------------------------
C          .... CHANGED:  210585  OUTPUT FOR NSUMM
C          LAST CHANGED:  151085  REFINEMENT OF CUBIC,HEX,TETR.
C-----------------------------------------------------------------------
      DIMENSION REF(4,4),VAR(6),FYN(4),DIFR(4),MATR1(4),MATR2(4),TEMP(4)
      INTEGER  H,HG, ORDE,ORP1,ORP2
      LOGICAL SELF,ENL,TEST,TYP
      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 /LINSTR/ NMAX,NCALC,N20,NOTIND,INSTOR(200,3),QSTOR(200)
      COMMON /CONSTS/ NSUMM,RAD,TWORAD,RRAD,HLFLAM,WAVSQR,STORE(200,9)
      COMMON /THETAS/ ZERTOT,OBSTET(40),SHFTET(40),TETPT(40),TETMT(40)
      COMMON /HISYM / ISOL,SYMCON(12,6),NMATCH(45)
      TYP=.FALSE.
      NMAX=MAX
      QMAX=OBSPT(MAX)+10.
      QX=SYMCON(ISOL,1)
      QZ=SYMCON(ISOL,2)
      NSYM=NINT( SYMCON(ISOL,3) )
      ZERTOT=ROSTER(ISOL,11)
C  LATTICE IS:  CUBIC  TETRAGONAL RHOMBOHEDRAL  HEXAGONAL
C  NSYM =         3         4          5            6
      HG=NINT( SQRT(QMAX/QX) )
      IF (NSYM .EQ. 3) THEN
         LG=HG
         QZ=QX
      ELSE
         LG=NINT( SQRT(QMAX/QZ) )
      END IF
      LMAX=LG
      LSTEP=1
      IF (NSYM .EQ. 5) THEN
         LMAX=(LG+2)/3
         LMAX=3*LMAX
         LSTEP=3
      END IF
      ZERSHF=0.0
      TOLR=0.5*RRAD*TOLG
C  THE THETA- VALUES ARE STORED AS THETA AND IN RADIANS.
      DO 1000 NCYCLE=1,5
      NCALC=0
      DO 20 H=1,4
         FYN(H)=0.
         DIFR(H)=0.
         DO 10 K=1,4
   10       REF(H,K)=0.
   20 CONTINUE
      DO 30 J=1,40
         NMATCH(J)=0
   30 INDEX(J)=0
      F=1.+(5-NCYCLE)/2.0
      TOL=F*TOLR
      DO 40 M=1,MAX
         SHFTET(M)=OBSTET(M)-ZERTOT
         TETPT(M)=SHFTET(M)+TOL
         TETMT(M)=SHFTET(M)-TOL
   40 CONTINUE
      ORDE=2
      IF (NSYM .EQ. 3) ORDE=1
      ORP1=ORDE+1
      ORP2=ORDE+2
      KHMIX=1
      IF (NSYM .LT. 5) KHMIX=0
C-----------------------------------------------------------------------
      DO 80 H=0,HG
         DO 70 K=H,HG
            NN=H*H+K*K+H*K*KHMIX
            A=NN*QX
            IF (A .GT.QMAX) GO TO 80
            IF (NSYM .EQ. 3) THEN
               LMIN=K
            ELSE IF (NSYM .EQ. 5) THEN
               M=MOD(K-H,3)
               IF (M .EQ. 0) THEN
                  LMIN=0
               ELSE
                  LMIN=-M-LMAX
               END IF
            ELSE
               LMIN=0
            END IF
            DO 60 L=LMIN,LMAX,LSTEP
               Q=A+L*L*QZ
               IF (Q .LT. 1.0 ) GO TO 60
               IF (Q .GT. QMAX) GO TO 60
               NCALC=NCALC+1
               IF (NCALC .GT. 200) THEN
                  TYP=.TRUE.
                  GO TO 90
               END IF
               INSTOR(NCALC,1)=H
               INSTOR(NCALC,2)=K
               INSTOR(NCALC,3)=L
               QSTOR(NCALC)=Q
               Z=TETCAL(Q)
               DO 50 M=1,MAX
                  IF (TETPT(M) .LT. Z) GO TO 50
                  IF (TETMT(M) .GT. Z) GO TO 55
C           THERE IS A FITTING LINE. USE IT FOR THE REFINEMENT.
C                 PRINT *,' A FITTING LINE IS NR ',M,TWORAD*OBSTET(M)
                  I=NMATCH(M)+1
                  NMATCH(M)=I
C                 WRITE(IOUT,360) H,K,L,Q,TWORAD*Z,M
C 360 FORMAT(1X,3I4,F10.2,F10.3,'  LINE NR ', I2 ,'  FITS' )
                  CONS=WAVSQR/SIN(2.*Z)
                  W=1.
                  DIF=SHFTET(M)-Z
                  VAR(1)=NN*CONS
                  VAR(2)=L*L*CONS
                  IF (NSYM .EQ. 3) VAR(1)=VAR(1)+VAR(2)
                  VAR(ORP1)=1.
                  DO 42 J=1,ORP1
                     TEMP(J)=W*VAR(J)
                     DIFR(J)=DIFR(J)+W*DIF*VAR(J)
   42             CONTINUE
                  DO 48 J=1,ORP1
                     DO 45 I=1,ORP1
   45                REF(I,J)=REF(I,J)+VAR(I)*TEMP(J)
   48             CONTINUE
   50       CONTINUE
C  50       IS THE END OF THE LOOP ON THE SEARCH FOR A FITTING  Q
   55       CONTINUE
C  55       WRITE(IOUT,365) H,K,L,Q,TWORAD*Z
C 365       FORMAT(1X,3I4,F10.2,F10.3,'  NO FIT ' )
   60    CONTINUE
C  60    IS THE END OF THE LOOP ON L
   70 CONTINUE
C  END OF LOOP ON K
   80 CONTINUE
C  80 IS END OF LOOP ON H
C-----------------------------------------------------------------------
   90 CALL ARRAY(2,ORP1,ORP1,4,4,REF,REF)
      CALL MINV(REF,ORP1,DETER,MATR1,MATR2)
      CALL ARRAY(1,ORP1,ORP1,4,4,REF,REF)
      IF (ABS(DETER) .LT. 1.0E-30) THEN
         WRITE(IOUT,'(A38)' ) ' NOT REFINED, DETER LESS THAN 1.0E-30 '
         GO TO 200
      END IF
      NINDEX=0
      DO 120 I=1,MAX
         IF (NMATCH(I) .GT. 0) NINDEX=NINDEX+1
  120 CONTINUE
      SYMCON(ISOL,4)=NINDEX
      DO 170 I=1,ORP1
         DO 160 J=1,ORP1
  160    FYN(I)=FYN(I)+REF(I,J)*DIFR(J)
  170 CONTINUE
      ZERSHF=FYN(ORP1)
      ZERTOT=ZERTOT+ZERSHF
C  QX AND QZ ARE Q-VALUES,THAT WERE MULTIPLIED BY 1.E+4
      QX=QX+FYN(1)*1.E+4
      IF (ORDE .GT. 1) THEN
         QZ=QZ+FYN(2)*1.E+4
      ELSE
         QZ=QX
      END IF
C     IF (TEST) WRITE(IOUT,4) ISOL,QX,QZ,ZERSHF*TWORAD
  200 CONTINUE
 1000 CONTINUE
      IF (TYP) WRITE(IOUT,7)
      IF (NSYM .LT. 5) THEN
         VOL=1.E+6 / SQRT(QX*QX*QZ)
      ELSE
         VOL=2.E+6 / SQRT(3.*QX*QX*QZ)
      END IF
      IF (NSYM .EQ. 5 ) THEN
         WRITE(IOUT ,5) QX,QZ,NINDEX,MAX,ZERTOT*TWORAD,VOL
         WRITE(NSUMM,5) QX,QZ,NINDEX,MAX,ZERTOT*TWORAD,VOL
      ELSE
         WRITE(IOUT ,6) QX,QZ,NINDEX,MAX,ZERTOT*TWORAD,VOL
         WRITE(NSUMM,6) QX,QZ,NINDEX,MAX,ZERTOT*TWORAD,VOL
      END IF
      SYMCON(ISOL,1)=QX
      SYMCON(ISOL,2)=QZ
      SYMCON(ISOL,5)=ZERTOT*TWORAD
      SYMCON(ISOL,6)=VOL
      RETURN
C   4 FORMAT(' LATTICE ',I3,18H AFTER REFINEMENT , 2F10.3,
C    1 '    ZERO SHIFT=',F8.3,' DEGREES TWOTHETA ,  TOTAL:', F8.3)
    5 FORMAT(' THE RHOMBOHEDRAL LATTICE WITH THE (HEXAGONAL) CONSTANTS'
     1 /1X,2F12.3,'  CAN INDEX',I3,' FROM THE',I3,' LINES ' /
     2 ' ZEROSHIFT= ',F7.3,'  VOLUME=',F10.2 //)
    6 FORMAT(' THE SYMMETRICAL LATTICE WITH THE  CONSTANTS'
     1 /1X,2F12.3,'  CAN INDEX',I3,' FROM THE',I3,' LINES ' /
     2 ' ZEROSHIFT= ',F7.3,'  VOLUME=',F10.2 //)
    7 FORMAT(// 1X,30(2HYY)/ ' SUBR. REFSYM.  ARRAY INSTORE WAS FULL' /
     1  ' GENERATION OF Q-VALUES WAS DISCONTINUED.'//
     2  ' THIS MAKES IT RATHER IMPROBABLE THAT THE HIGH-SYMMETRY ',
     3  'LATTICE IS CORRECT. '/ 1X,30(2HYY) // )
      END

      SUBROUTINE REFINE(IREFCY,ALL,DRUK,PRIN,LIST2)
C-----------------------------------------------------------------------
C
C  CHANGED TO CONFORM WITH FORTRAN 77.  SEPT. 83
C
C  CHANGED TO REFINE DELTA THETA  INSTEAD OF DELTA Q.SEPT. 83
C
C  REFINES CELL CONSTANTS + ZERO SHIFT.  OCT.  83
C-----------------------------------------------------------------------
C  .... CHANGE:  140585,  CHANGED /HISYM / FOR OUTPUT SUMMARY.
C  LAST CHANGE:  151085,  CHANGED FORMATS  FOR OUTPUT SUMMARY.
C-----------------------------------------------------------------------
C
C--REFPAR REFINES  THE PARAMETERS, PRINTS OUT THE LATTICE IF PRIN=.TRUE.
C--AND PRINTS OUT A LIST OF   TWOTHETA,D,Q,H,K,L.
C  IN THE FIRST CALLS OF REFINE, ALL PARAMETERS ARE .FALSE.
C  THE CELL CONSTANTS ARE THEN 'SILENTLY' REFINED.
C  AT THE END, THE BRAVAIS TYPE IS ESTABLISHED.
C  IN THE FOLLOWING CALLS, THE BOOLEANS DO THE FOLLOWING
C  PRIN  COMMANDS THE PRINTING OF A Q-SCHEME
C  DRUK  PRINTS THE TWOTHETA-HKL-Q LISTS
C  LIST2 PRINTS A LIST OF THE OBSERVED AND ALL CALCULATED LINES
C  ALL  TRIES TO INDEX ALL GIVEN LINES(<41).
C-----------------------------------------------------------------------
      DIMENSION HULP(8),VAR(8),FYN(8),WE(10),WAR(8),
     1 MATR1(8),MATR2(8),REF(8,8)
      DIMENSION KSTOR(40,5),LSTOR(40,5),QSTOR(40,5),IBR(41,5)
      INTEGER HSTOR(40,5)
      INTEGER P,  H,H1,HS,HMIN,HMAX,HMAX1,COUNT,ORDE,ORP1,ORP2
      LOGICAL SELF,ENL,TEST,ALL,DUBL
      LOGICAL HEVEN,KEVEN,LEVEN,MIS,DRUK,PRIN,LIST2
      CHARACTER NAME*80,CARD*80,INTSY*2
      COMMON /LINES / IDOC,INP,IOUT,MAX,OBS(40),CEL(10)
      COMMON /CHARCM/ NAME,CARD,INTSY(50)
      COMMON /PARMS / TOL2,TOL3,WAVEL,FMINEQ,FMINIM,SELF,ENL,TEST,LIST,
     *                NSYST(3),OBSMT(40),OBSPT(40),TWTHET(40),TOLG
      COMMON /LINSTR/ NMAX,NCALC,N20,NOTIND,INSTOR(200,3),QQSTOR(200)
      COMMON /CONSTS/ NSUMM,RAD,TWORAD,RRAD,HLFLAM,WAVSQR,STORE(200,9)
      COMMON /THETAS/ ZERTOT,OBSTET(40),SHFTET(40),TETPT(40),TETMT(40)
      COMMON /HISYM / ISOL,SYMCON(12,6),NMATCH(45)
      EQUIVALENCE (INSTOR(1,1),HSTOR(1,1)),(INSTOR(1,2),KSTOR(1,1)),
     1  (INSTOR(1,3),LSTOR(1,1)),(QQSTOR(1),QSTOR(1,1) )
      DATA NIN/0/
C-----------------------------------------------------------------------
      DO 20 M=1,MAX
   20 SHFTET(M)=OBSTET(M)-ZERTOT
      IF (TEST) WRITE(IOUT,9904) CEL
      IBRAV=NINT( CEL(9) )
      NMAX=20
      IF (ALL)  NMAX=MAX
      J20=20
      Q20=OBSPT(20)
      QMAX=OBSPT(NMAX)+1.
      IF (TEST) WRITE(IOUT,9923)
C  BEGIN OF  THE REFINEMENT CYCLES.
C
      DO 1250 NREFCY=1,IREFCY
C  FIRST, WE DETERMINE THE NUMBER OF PARAMETERS TO BE REFINED (ORDE),
C  ORDE IS THE ORDER OF THE MATRIX.
      ORDE=3
      DO 110 J=4,6
      IF (ABS(CEL(J)) .GT. FMINEQ) THEN
         ORDE=ORDE+1
      ELSE
         CEL(J)=.0
      END IF
  110 CONTINUE
      ORP1=ORDE+1
      ORP2=ORDE+2
C  SET THE BRAVAIS-ARRAY TO ZERO
      DO 211 I=1,41
      DO 210 J=1,5
  210 IBR(I,J)=0
      IF (.NOT. LIST2) NMATCH(I)=0
  211 CONTINUE
C  SET UP THE (ENLARGED) TOLERANCES
      TOLRAD=0.5*TOLG*RRAD
      TOL=(IREFCY-NREFCY+1)*TOLRAD
      DO 220 I=1,NMAX
      TETPT(I)=SHFTET(I)+TOL
      TETMT(I)=SHFTET(I)-TOL
  220 CONTINUE
C
C                 START GENERATING  Q-VALUES
C
      NCALC=0
      N20=0
      ISPAT=1
      HMAX=INT(SQRT(QMAX/CEL(1)))+1
      KMAX=INT(SQRT(QMAX/CEL(2)))+1
      LMAX=INT(SQRT(QMAX/CEL(3)))+1
      KEVEN=.FALSE.
      LMIN=-LMAX
      IF(ORDE .EQ. 3) LMIN=0
      HMIN=0
      IF(ORDE .GE. 5) HMIN=-HMAX
      IF (PRIN) WRITE(IOUT,9917)
      KMAX1=KMAX+1
C  K-LOOP STARTS HERE
      DO 1050 K1=1,KMAX1
      K=K1-1
      RK=K
      IF(PRIN) THEN
         WRITE(IOUT,9918)  K
         ISPAT=1
      END IF
      C1=RK*RK*CEL(2)
      C2=RK*CEL(4)
      C3=RK*CEL(6)
      KEVEN=.NOT. KEVEN
      HS=HMIN
      IF(K .EQ.  0) HS=0
      HEVEN=.TRUE.
      IF(2*(HS/2) .EQ. HS) HEVEN=.FALSE.
      HMAX1=HMAX-HS+1
C  H-LOOP STARTS HERE
      DO 1040 H1=1,HMAX1
      H=H1+HS-1
      RH=H
      C4=C1+RH*(C3+RH*CEL(1))
      C5=C2+RH*CEL(5)
      IF (PRIN) THEN
         WRITE(IOUT,9923)
         ISPAT=1
         IF((K+IABS(H)) .EQ. 0 .AND. ORDE .GT. 4) ISPAT=ISPAT+LMAX
         IF(ORDE .EQ. 4 . AND. H .EQ. 0) ISPAT=ISPAT+LMAX
      END IF
      HEVEN=.NOT. HEVEN
      LB=LMIN
      IF (H .EQ. 0) THEN
         IF(K .EQ. 0.OR. ORDE .EQ. 4) LB=0
      END IF
      LEVEN=.TRUE.
      IF(2*(LB/2) .EQ. LB) LEVEN=.FALSE.
      LMAX1=LMAX+1-LB
C  L-LOOP STARTS HERE
      DO 1030 L1=1,LMAX1
      L=L1-1+LB
      LEVEN=.NOT. LEVEN
      Q=C4+FLOAT(L)*(C5+FLOAT(L)*CEL(3))
      IF (Q .LT. 1.0) THEN
         IF (PRIN) CALL WRISP(ISPAT,0. ,0)
         GO TO 1020
      END IF
      IF(Q .GE. QMAX) GO TO 1020
      IF (CEL(7) .LT. 1.) GO TO 700
      GO TO (640,650,660,670,680,700,950),IBRAV
C REPRESENTING A, B , C , I , F , P    LATTICES
  640 IF(KEVEN .EQV. LEVEN) GO TO 700
      GO TO 1020
  650 IF(HEVEN .EQV. LEVEN) GO TO 700
      GO TO 1020
  660 IF(HEVEN .EQV. KEVEN) GO TO 700
      GO TO 1020
  670 IF (2*((H+K+L)/2) .EQ. H+K+L) GO TO 700
      GO TO 1020
  680 IF ((HEVEN .EQV. KEVEN) .AND. (KEVEN .EQV. LEVEN)) GO TO 700
      GO TO 1020
  700 IF (Q .LT. Q20) N20=N20+1
      NCALC=NCALC+1
      IF (LIST2) THEN
         INSTOR(NCALC,1)=H
         INSTOR(NCALC,2)=K
         INSTOR(NCALC,3)=L
         QQSTOR(NCALC)=Q
         GO TO 1030
      END IF
      IF (PRIN) CALL WRISP(ISPAT,Q,0)
      THET=TETCAL(Q)
      MIS=.TRUE.
      NRINDX=0
      DO 800 J=1,NMAX
      IF (TETPT(J) .LT. THET ) GO TO 800
      IF (TETMT(J) .GT. THET ) GO TO 950
      NRINDX=NRINDX+1
C  COUNTS THE NUMBER OF TIMES THAT ONE H,K,L-COMBINATION IS USED.
      P=J
      IF (PRIN) THEN
         CALL WRISP(ISPAT,Q,1)
         MIS=.FALSE.
      END IF
      IF (CEL(7) .GT. .9) GO TO 780
      IF (KEVEN .EQV. LEVEN) IBR(P,1)=IBR(P,1)+1
      IF (HEVEN .EQV. LEVEN) IBR(P,2)=IBR(P,2)+1
      IF (HEVEN .EQV. KEVEN) IBR(P,3)=IBR(P,3)+1
      IF (2*((H+K+L)/2) .EQ. H+K+L) IBR(P,4)=IBR(P,4)+1
      IF ((HEVEN.EQV.KEVEN).AND.(KEVEN.EQV.LEVEN)) IBR(P,5)=IBR(P,5)+1
C     PRINT *, H,K,L,Q, TWORAD*THET, (IBR(P,I),I=1,5)
  780 I=NMATCH(P)+1
      IF (I .GT.  5) GO TO 800
      NMATCH(P)=I
      IF (NRINDX .GE. 2) L=L+1000
C  IF A SET OFH,K,L'S IS USED MORE THAN ONCE, IT IS MARKED.
      HSTOR(P,I)=H
      KSTOR(P,I)=K
      LSTOR(P,I)=L
      QSTOR(P,I)=Q
  800 CONTINUE
  950 CONTINUE
 1020 IF (PRIN) ISPAT=ISPAT+1
C  END OF K,H,L-LOOP ON 1050,1040,1030 RESPECTIVELY
 1030 CONTINUE
      IF (PRIN) CALL WRISP(ISPAT,Q,2)
 1040 CONTINUE
 1050 CONTINUE
C
C             END GENERATING Q
C
      IF (LIST2) THEN
         CALL LISTM
         RETURN
      END IF
      DO 1110 H=1,8
      DO 1105 K=1,8
 1105 REF(H,K)=0.0
 1110 CONTINUE
      DO 1120 H=1,8
      FYN(H)=0.0
      WAR(H)=0.0
 1120 CONTINUE
      T=0.0
      RMS=0.0
      SUMEPS=1.E-3
C  WHEN THE NUMBER OF LINES INDEXED  EQUALS THE NUMBER OF PARAMETERS,
C  SUMEPS WILL REMAIN AT ITS STARTING VALUE.  THE VALUE *ZERO*
C  CAUSES ACCIDENTS.
      COUNT=0
      IF (DRUK)  WRITE(IOUT,9910)
      DO 1180 P=1,NMAX
      NR=NMATCH(P)
      W=SIN(SHFTET(P))/HLFLAM
      QOBS=W*W
      IF (DRUK) THEN
         D=100./W
         U=TWORAD*SHFTET(P)
         IF (NR .EQ. 0) THEN
            WRITE(IOUT,9911) U,D,QOBS,INTSY(P)
         ELSE
            WRITE(IOUT,9912) U,D,QOBS,INTSY(P)
         END IF
      END IF
      IF (NR .EQ. 0) GO TO 1180
      EPS=1000.0
      SUMW=0.0
      DO 1130 N=1,NR
         D=ABS(QOBS-QSTOR(P,N))
         IF(D .LT. EPS)EPS=D
         IF (D .GT. 9.0) D=9.0
C--NECESSARY TO PREVENT UNDERFLOW ON  PDP-10 (AND UNIVAC?)
         W=EXP(-D*D)
         WE(N)=W
         SUMW=SUMW+W
 1130 CONTINUE
      IF (P .LE. J20) SUMEPS=SUMEPS+EPS
      DO 1170 N=1,NR
         DUBL=.FALSE.
         H=HSTOR(P,N)
         K=KSTOR(P,N)
         L=LSTOR(P,N)
         Q=QSTOR(P,N)
         TETC=TETCAL(Q)
         IF (L .GT. 900) THEN
            L=L-1000
            DUBL=.TRUE.
         END IF
C  FLAG DOUBLE USAGE OF 1 SET OF HKL
         IF (DRUK) THEN
            D=100./SQRT(Q)
            U=TWORAD*TETC
            IF (DUBL) THEN
               WRITE(IOUT,9914) U,D,Q,H,K,L
            ELSE
               WRITE(IOUT,9915) U,D,Q,H,K,L
            END IF
         END IF
         DIF=SHFTET(P)-TETC
         COUNT=COUNT+1
         W=WE(N)/SUMW
         RMS=RMS+W*DIF*DIF
         T=T+W
         C=WAVSQR/SIN(2.*TETC)
         VAR(1)=H*H*C
         VAR(2)=K*K*C
         VAR(3)=L*L*C
         VAR(4)=H*L*C
         VAR(5)=K*L*C
         VAR(6)=H*K*C
         VAR(ORP1)=1.
         DO 1150 J=1,ORP1
            HULP(J)=W*VAR(J)
            WAR(J)=WAR(J)+W*DIF*VAR(J)
 1150    CONTINUE
         DO 1160 J=1,ORP1
            DO 1155 I=1,ORP1
 1155       REF(I,J)=REF(I,J)+VAR(I)*HULP(J)
 1160    CONTINUE
 1170 CONTINUE
 1180 CONTINUE
      CALL ARRAY(2,ORP1,ORP1,8,8,REF,REF)
      CALL MINV(REF,ORP1,DETER,MATR1,MATR2)
      CALL ARRAY(1,ORP1,ORP1,8,8,REF,REF)
      IF (ABS(DETER) .LT. 1.E-30) THEN
         IF (TEST) WRITE(IOUT,9931)
         GO TO 1250
      END IF
      IF (T .GT. FLOAT(ORP1+1)) RMS=RMS/(T-FLOAT(ORP1))
      DO 1195 I=1,ORP1
         DO 1190 J=1,ORP1
 1190    FYN(I)=FYN(I)+REF(I,J)*WAR(J)
 1195 CONTINUE
      DO 1200 N=1,ORDE
 1200 FYN(N)=1.E+4*FYN(N)
      DO 1210 N=1,3
 1210 CEL(N)=CEL(N)+FYN(N)
      IF (ORDE .GT. 3) CEL(5)=CEL(5)+FYN(4)
      IF (ORDE .GT. 4) CEL(4)=CEL(4)+FYN(5)
      IF (ORDE .GT. 5) CEL(6)=CEL(6)+FYN(6)
      ZERSHF=FYN(ORP1)
      ZERTOT=ZERTOT+ZERSHF
      IF (TEST) WRITE(IOUT,9925) TWORAD*ZERSHF,CEL
      DO 1240 M=1,MAX
 1240 SHFTET(M)=SHFTET(M)-ZERSHF
C
C
 1250 CONTINUE
      IF (DRUK) WRITE(IOUT,9927) -TWORAD*ZERTOT
C  END OF REFINEMENT CYCLE
C-----------------------------------------------------------------------
C
C  NOW DETERMINE BRAVAIS TYPE AND FIGURE OF MERIT
C
      N=0
      DO 1215 I=1,NMAX
      IF (NMATCH(I) .GE. 1) N=N+1
      J20=I
      IF (N .EQ. 20) GO TO 1220
 1215 CONTINUE
C  N=NUMBER OF INDEXED LINES,  J20 IS THE NUMBER OF EXAMINED LINES,
C  IN ORDER TO FIND 20 INDEXED LINES.
 1220 NOTIND=J20-N
      INDEXD=N
C  DETERMINE THE BRAVAIS TYPE ONLY AT THE FIRST CALL OF *REFPAR
      IF (CEL(7) .GT. 0.9) GO TO 1350
C
      DO 1320 J=1,5
      N=0
      DO 1310 M=1,20
 1310 IF (IBR(M,J) .GT. 0 ) N=N+1
      IBR(21,J)=N
 1320 CONTINUE
      IBRAV=6
      DO 1330 J=1,5
      IF ( IBR(21,J) .GE. INDEXD) IBRAV=J
 1330 CONTINUE
      CEL(9)=IBRAV
C-----------------------------------------------------------------------
C
 1350 CEL(7)=INDEXD
      CEL(8)=FLOAT(INDEXD)*OBS(J20)/(FLOAT(2*N20)*SUMEPS)
      IF (TEST) WRITE(IOUT,9926) CEL
      RETURN
 9904 FORMAT(2X/' AT ENTERING REFINE, CEL= ',6F10.3,4F8.2 )
 9910 FORMAT(2X/41H TWOTHETA     D       Q       H    K    L)
 9911 FORMAT(2X/1H , F8.2, F8.3, F9.1, 18H    *OBSERVED**    ,A2,
     *       '  NOT INDEXED' )
 9912 FORMAT(2X/1H , F8.2, F8.3, F9.1, 18H    *OBSERVED**    ,A2)
 9914 FORMAT(1H ,F8.2,F8.3,F9.1, 3I5,
     *      '   HKL-COMBINATION WAS USED BEFORE' )
 9915 FORMAT(1H ,F8.2,F8.3,F9.1, 3I5)
 9917 FORMAT(36H0LATTICE PRINTED OUT FOR INSPECTION //
     1 56H DISPLAYED IN LAYERS UP --K-- ,    --L-- ACROSS THE PAGE,
     2 22H --H-- DOWN THE PAGE.   //)
 9918 FORMAT( 8H ***  H ,I2,8H  L  ***   // )
 9923 FORMAT(2X / 2X)
 9925 FORMAT(1H ,'ZEROSHIFT IS',F10.4,'  CELL IS NOW:', 6F9.2,4F5.0)
 9926 FORMAT(1H ,' AT END OF SUBROUTINE REFINE, CEL =  '/
     *           6F9.2, F9.0,F9.1, 2F9.0 )
 9927 FORMAT(2X/' ZEROSHIFT IS ',F10.4,'  DEGREES TWOTHETA')
 9931 FORMAT(1H ,14X,19HMATRIX IS SINGULAR )
      END

      SUBROUTINE SYMTST
C  SUBROUTINE SYMTEST TRIES TO FIND OUT WHETHER THE LATTICES FOUND
C  ARE LATTICES OF A HIGHER SYMMETRY IN DISGUISE.
C  LAST CHANGE:  15 OCT. 1985,  REFINING HIGH-SYMM. LATTICES
C
      CHARACTER*12  TEXT(2)
      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 /CONSTS/ NSUMM,RAD,TWORAD,RRAD,HLFLAM,WAVSQR,STORE(200,9)
      COMMON /HISYM / ISOL,SYMCON(12,6),NMATCH(45)
      NPASS=0
      SYMCON(ISOL,3)=0
      A=CEL(1)
      B=CEL(2)
      C=CEL(3)
      E=CEL(5)
      DIFMIN=FMINEQ
      IF (CEL(4)+ABS(CEL(6)) .GT. FMINEQ) GO TO 3010
C-----------------------------------------------------------------------
C  THE MONOCLINIC LATTICES ARE TESTED BETWEEN HERE AND LABEL 1000,
C  THE ORTHORHOMBIC LATTICES ARE TESTED BETWEEN LABELS 1000 AND 3000,
C  THE TRICLINIC LATTICES ARE TESTED AFTER LABEL 3000.
C-----------------------------------------------------------------------
   5  DIFMIN=3.0*DIFMIN
      NPASS=NPASS+1
      IF ( E .LT. FMINEQ) GO TO 1000
C  THE NEXT PART   FINDS OUT WHETHER A MONOCLINIC LATTICE IS A DISGUISED
C  RHOMBOHEDRIC LATTICE.
      IF (ABS(E-2.*B) .GT. DIFMIN) GO TO 20
      IF (ABS(C-3.*B) .GT. DIFMIN) GO TO 10
      X=(B+C+E)/18.
      Z=A-X
      GO TO 200
   10 IF (ABS(A-3.*B) .GT. DIFMIN) GO TO 20
      X=(A+B+E)/18.
      Z=C-X
      GO TO 200
   20 IF (ABS(A+B-C) .GT. DIFMIN) GO TO 30
      IF (ABS(2.*ABS(B-A)-E).GT.DIFMIN) GO TO 30
      X=(C-A)/6.+B/6.
      Z=(C+2.*A-2.*B)/3.
      GO TO 200
   30 IF (ABS(B+C-4.*A).GT. DIFMIN) GO TO 40
      IF (ABS(E-ABS(B-C)) .GT. DIFMIN) GO TO 40
      Z=(C-A)/3.
      X=B/6.+(4.*A-C)/6.
      GO TO 200
   40 IF (ABS(2.*C-3.*E).GT.DIFMIN) GO TO 50
      IF (ABS(2.*B+E-6.*A).GT. DIFMIN*FLOAT(NPASS) ) GO TO 50
      Z=(C+E)/15.
      X=(A-Z)/2+B/6.
      GO TO 200
   50 IF (ABS(2.*A-3.*E) .GT. DIFMIN) GO TO 60
      IF (ABS(2.*B+E-6.*C).GT. DIFMIN*FLOAT(NPASS) ) GO TO 60
      Z=(A+E)/15.
      X=(C-Z)/2. +B/6.
      GO TO 200
   60 IF (NPASS .EQ. 1) GO TO 5
      GO TO 3010
  200 IF (Z .LT. 1.) RETURN
      WRITE(IOUT,900)  (CEL(J),J=1,6)
      WRITE(NSUMM,930) (CEL(J),J=1,6)
      A=200./SQRT(B)
      C=100./SQRT(Z)
      WRITE(IOUT ,901) TEXT(NPASS), X,Z, A, C
      WRITE(NSUMM,931) TEXT(NPASS), X,Z, A, C
      IF (CEL(9) .LT. 6) THEN
         X=4.*X
         A=0.5*A
         WRITE(IOUT ,902) X,A
         WRITE(NSUMM,902) X,A
      END IF
      SYMCON(ISOL,1)=X
      SYMCON(ISOL,2)=Z
      SYMCON(ISOL,3)=5
C  THE VALUES OF SYMCON(I,3) ARE RESERVED AS FOLLOWS: 3 FOR CUBIC,
C  4 FOR TETRAGONAL, 5 FOR RHOMBOHEDRAL AND 6 FOR HEXAGONAL LATTICES.
      CALL REFSYM
      RETURN
C-----------------------------------------------------------------------
C                ORTHORHOMBIC LATTICES FROM HERE TO 3000
C-----------------------------------------------------------------------
 1000 DO 1200 I=1,2
      J=3-I
      G=CEL(J)
      H=INT(A/G+0.5 )
      K=INT(B/G+0.5 )
      L=INT(C/G+0.5 )
      IH=0
      IK=0
      IL=0
      IF (ABS(A-G*FLOAT(H)) .LT. DIFMIN) IH=1
      IF (ABS(B-G*FLOAT(K) ) .LT. DIFMIN) IK=1
      IF (ABS(C-G*FLOAT(L) ) .LT. DIFMIN) IL=1
      II=IH+IK+IL -2
      IF (II) 1200,1020,1010
 1010 WRITE(IOUT,900)  (CEL(J),J=1,6)
      WRITE(NSUMM,930)  (CEL(J),J=1,6)
      X=100./SQRT(G)
      WRITE(IOUT, 911) TEXT(NPASS),G,X
      WRITE(NSUMM,941) TEXT(NPASS),G,X
C  LATTICE IS CUBIC
      SYMCON(ISOL,1)=G
      SYMCON(ISOL,2)=G
      SYMCON(ISOL,3)=3
      CALL REFSYM
      RETURN
 1020 IF (IK*IL .EQ. 0) GO TO 1050
      N=L/K
      X=B
      Z=A
      ICEL9=NINT(CEL(9))
      IF(ICEL9.EQ.1) X=4.0*X
      GO TO 1150
 1050 IF (IH*IL .EQ. 0) GO TO 1100
      N=L/H
      X=A
      Z=B
      ICEL9=NINT(CEL(9))
      IF(ICEL9.EQ.2) X=4.0*X
      GO TO 1150
 1100 N=H/K
      X=B
      ICEL9=NINT(CEL(9))
      IF(ICEL9.EQ.3) X=4.0*X
      Z=C
 1150 WRITE(IOUT,900)  (CEL(J),J=1,6)
      WRITE(NSUMM,930)  (CEL(J),J=1,6)
      A=100./SQRT(X)
      C=100./SQRT(Z)
      IF (N-3) 1160,1170,1180
 1160 WRITE(IOUT ,912) TEXT(NPASS),X,Z, A,C
      WRITE(NSUMM,942) TEXT(NPASS),X,Z, A,C
C  LATTICE IS TETRAGONAL
      SYMCON(ISOL,1)=X
      SYMCON(ISOL,2)=Z
      SYMCON(ISOL,3)=4
      CALL REFSYM
      RETURN
 1170 A=2.0*A/SQRT(3.0)
      WRITE(IOUT ,913) TEXT(NPASS),X,Z,A,C
      WRITE(NSUMM,943) TEXT(NPASS),X,Z,A,C
C  LATTICE IS HEXAGONAL
      SYMCON(ISOL,1)=X
      SYMCON(ISOL,2)=Z
      SYMCON(ISOL,3)=6
      CALL REFSYM
      RETURN
 1180 WRITE(IOUT ,914) TEXT(NPASS),N
      WRITE(NSUMM,944) TEXT(NPASS),N
C  LATTICE IS A FREAK
      RETURN
 1200 CONTINUE
C---------------------------------------------------------------------
C       TRICLINIC LATTICES ARE CHECKED FOR EQUALITY OF CONSTANTS
C               OR SIMPLE INTEGER RATIOS BETWEEN THEM.
C---------------------------------------------------------------------
 3010 AA=OBS(MAX)
      Z=0.
      DIFMIN=FMINEQ
      DO 3020 J=1,6
      X=ABS(CEL(J))
      IF (X .LT. FMINEQ) GO TO 3020
      IF (X .LT. AA) AA=X
C  FINDS THE SMALLEST NON-ZERO VALUE OF THE CONSTANT, AA.
      IF (X .GT. Z ) Z=X
C  AND THE LARGEST VALUE,  Z.
 3020 CONTINUE
      DO 3050 I=1,3
      X=AA/FLOAT(I)
      IF (X .LT. 0.5*FMINIM) RETURN
      L=INT(Z/X +0.5)
      X=Z/FLOAT(L)
C  REFINES THE VALUE OF THE CONSTANT BY USING THE LARGEST MULTIPLE.
      N=0
      DO 3030 J=1,6
         G=ABS(CEL(J))
         L=NINT(G/X)
         IF (ABS(G-X*FLOAT(L)) .LE. DIFMIN) N=N+1
 3030 CONTINUE
      IF (N .EQ. 6) GO TO 3100
 3050 CONTINUE
      RETURN
 3100 WRITE(IOUT,900) (CEL(J),J=1,6)
      WRITE(NSUMM,930) (CEL(J),J=1,6)
      Z=100./SQRT(X)
      WRITE(IOUT ,911) TEXT(NPASS),X,Z
      WRITE(NSUMM,941) TEXT(NPASS),X,Z
C  LATTICE IS CUBIC
      SYMCON(ISOL,1)=X
      SYMCON(ISOL,2)=X
      SYMCON(ISOL,3)=3
      CALL REFSYM
      RETURN
  900 FORMAT(2X/52H THE LATTICE WITH THE RECIPROCAL LATTICE CONSTANTS  ,
     * 6F9.2 )
  901 FORMAT(1H ,A12, 24H RHOMBOHEDRAL, WITH THE   ,
     *45H (HEXAGONAL) RECIPROCAL LATTICE CONSTANTS  X=,F9.1,4H  Z=,F9.1/
     * 29H AND THE DIRECT CONSTANTS  A=, F7.2,4H  C=,F7.2 )
  902 FORMAT(21H PROBABLY X SHOULD BE, F10.1,8H  AND A   ,F8.2)
  911 FORMAT(1H ,A12,44H CUBIC WITH A RECIPROCAL LATTICE CONSTANT OF  ,
     *        F10.1,33H AND A DIRECT LATTICE CONSTANT OF , F10.2)
  912 FORMAT(1H ,A12,42H TETRAGONAL WITH THE RECIPROCAL CONSTANTS  ,
     * 2F10.1,29H AND THE DIRECT CONSTANTS  A=, F7.2,4H  C=,F7.2 )
  913 FORMAT(1H ,A12,41H HEXAGONAL WITH THE RECIPROCAL CONSTANTS   ,
     * 2F10.1,29H AND THE DIRECT CONSTANTS  A=, F7.2,4H  C=,F7.2 )
  914 FORMAT(1H ,A12,43H A FREAK. THERE IS A RATIO IN THE CONSTANTS   ,
     * 32H OF THE RECIPROCAL  LATTICE OF   ,I3)
  930 FORMAT(2X/52H THE LATTICE WITH THE RECIPROCAL LATTICE CONSTANTS  /
     * 6F9.2 )
  931 FORMAT(1H ,A12,' RHOMBOHEDRAL, WITH THE  (HEXAGONAL) RECIPROCAL'/
     * ' LATTICE CONSTANTS  X=',F9.1,4H  Z=,F9.1/
     * 29H AND THE DIRECT CONSTANTS  A=, F7.2,4H  C=,F7.2 )
  941 FORMAT(1H ,A12,44H CUBIC WITH A RECIPROCAL LATTICE CONSTANT OF  /
     *        F10.1,33H AND A DIRECT LATTICE CONSTANT OF , F10.2)
  942 FORMAT(1H ,A12,42H TETRAGONAL WITH THE RECIPROCAL CONSTANTS  /
     * 2F10.1,29H AND THE DIRECT CONSTANTS  A=, F7.2,4H  C=,F7.2 )
  943 FORMAT(1H ,A12,41H HEXAGONAL WITH THE RECIPROCAL CONSTANTS   /
     * 2F10.1,29H AND THE DIRECT CONSTANTS  A=, F7.2,4H  C=,F7.2 )
  944 FORMAT(1H ,A12,43H A FREAK. THERE IS A RATIO IN THE CONSTANTS   /
     * 32H OF THE RECIPROCAL  LATTICE OF   ,I3)
      END

      SUBROUTINE REFPAR(NSOL,ALL)
C  REFPAR TAKES CARE THAT THE CELL PARAMETERS ARE REFINED.
C  THE COMPUTING PART IS DONE IN *REFINE*, BUT THE
C  STEERING PARAMETERS ARE ESTABLISHED HERE.
C---------------------------------------------------------------
C  .... CHANGED:  250386  ADDED NSUMM, OUTPUT 9942
C  LAST CHANGED:  080486  ADDED NSUMM, OUTPUT 9910 + 9912
C---------------------------------------------------------------
      LOGICAL SELF,ENL,TEST,ALL,DRUK,PRIN,LIST2
      CHARACTER NAME*80,CARD*80,INTSY*2
      COMMON /LINES / IDOC,INP,IOUT,MAX,OBS(40),CEL(10)
      COMMON /CHARCM/ NAME,CARD,INTSY(50)
      COMMON /CONSTS/ NSUMM,RAD,TWORAD,RRAD,HLFLAM,WAVSQR,STORE(200,9)
      COMMON /PARMS / TOL2,TOL3,WAVEL,FMINEQ,FMINIM,SELF,ENL,TEST,LIST,
     *                NSYST(3),OBSMT(40),OBSPT(40),TWTHET(40),TOLG
      COMMON /LINSTR/ NMAX,NCALC,N20,NOTIND,INSTOR(200,3),QQSTOR(200)
C-----------------------------------------------------------------------
      DRUK =.FALSE.
      PRIN =.FALSE.
      LIST2=.FALSE.
C  THE THREE PLACES IN THE MAIN PROGRAM WHERE REFPAR IS CALLED
C  CAN BE DISTINGUISHED BY THE BEHAVIOUR OF CEL(7) AND  ALL.
C  THE FIRST TIME,  CEL(7)=0., ALL=.FALSE.
C  THE SECOND TIME, CEL(7)>0 AND  ALL=.FALSE.  THE LAST TIME ALL=.TRUE.
C
      IF (CEL(7) .GT. 1.0) GO TO 1000
      IREFCY=3
      CALL REFINE(IREFCY,ALL,DRUK,PRIN,LIST2)
      RETURN
C
C
 1000 IF (ALL) GO TO 2000
      IREFCY=1
      CALL REFINE(IREFCY,ALL,DRUK,PRIN,LIST2)
      IF (TEST) PRIN=.TRUE.
      DRUK=.TRUE.
      WRITE(IOUT,9916) NAME,NSOL
      WRITE(IOUT,9926) CEL
      CALL REFINE(IREFCY,ALL,DRUK,PRIN,LIST2)
      WRITE(IOUT,9941) CEL,NAME,NSOL
      IF (NSOL .GT. LIST) GO TO 1100
      IF (N20 .GT. 100) THEN
         WRITE(IOUT,9942) N20
         WRITE(NSUMM,9942) N20
         GO TO 1100
      END IF
      LIST2=.TRUE.
      PRIN=.FALSE.
      CALL REFINE(IREFCY,ALL,DRUK,PRIN,LIST2)
      WRITE(IOUT,9941) CEL,NAME,NSOL
 1100 RETURN
C
C
 2000 WRITE(IOUT,9902) NAME
      IREFCY=1
      CALL REFINE(IREFCY,ALL,DRUK,PRIN,LIST2)
      WRITE(IOUT,9943) CEL(8),NOTIND,N20,NCALC
      DRUK=.TRUE.
      IF (NCALC .LE. 4*NMAX) THEN
         LIST2=.TRUE.
      ELSE
         WRITE(NSUMM,9910) NCALC
      END IF
      CALL REFINE(IREFCY,ALL,DRUK,PRIN,LIST2)
      WRITE(NSUMM,9912) NOTIND
      RETURN
 9902 FORMAT(1H1,20X,A80 / 40(3H---)/
     *       50H0THE BEST SOLUTION IS NOW USED TO TRY AND INDEX AL,
     *        9HL LINES.  )
 9910 FORMAT(/,' THE NUMBER OF CALCULATED LINES IS ',I3 /
     *  ' NO LISTING IS PRINTED IF THIS NUMBER IS GREATER THAN'/
     *  ' 4 * THE NUMBER OF OBSERVED LINES. ')
 9912 FORMAT(/,' THERE ARE ',I2,' UNINDEXED LINES AMONGST THE',
     *  ' FIRST 20 INDEXED LINES. ')
 9916 FORMAT(1H1,10X,A80,20X,'SOLUTION NR',I3/ 40(3H---) )
 9926 FORMAT(1H ,6F9.2, F9.0,F9.1, 2F9.0 )
 9941 FORMAT(1H0//46H RECIPROCAL LATTICE CONSTANTS AFTER REFINEMENT/
     1 1H ,6F9.2/14H LINES INDEXED,F4.0,18H   FIGURE OF MERIT,F5.1,
     2 16H    LATTICE TYPE,F4.0, 12H    PEDIGREE,F5.0 /
     3 1H ,A80,14H   SOLUTION NR, I3)
 9942 FORMAT(/ ' THE NUMBER OF CALCULATED LINES (N20) IS ',I4,/
     *   ' NO LISTING IS PRINTED WHEN N20 IS MORE THAN 100 ' )
 9943 FORMAT(50H0TWO CYCLES OF REFINEMENT WITH ALL LINES HAVE BEEN,
     1       50H CARRIED OUT.                                     /
     2       25H THE FIGURE OF MERIT IS  , F7.1 /
     3       50H THE NUMBER OF UNINDEXED LINES AMONG THE FIRST 20 ,
     4       18H INDEXED LINES IS ,I3,4H=X20  /
     5       50H THE NUMBER OF LINES CALCULATED UP TO THE 20ST. IN,
     6       15HDEXED LINE IS  , I3 ,4H=N20  /
     7       42H THE TOTAL NUMBER OF LINES CALCULATED IS  , I4   )
      END

      SUBROUTINE UNTCEL(WMOL,DOBS,VOLUM,ISW)
C-----------------------------------------------------------------------
C  THIS SUBROUTINE CALCULATES THE DIRECT CONSTANTS FROM THE CONSTANTS
C  OF THE RECIPROCAL LATTICE.
C  IF THE MOLECULAR WEIGHT AND THE DENSITY HAVE BEEN GIVEN,
C  IT ALSO CALCULATES THE X-RAY DENSITY,DX,  AND THE NUMBER OF MOLECULES
C  PER UNIT CELL,  N.
C-----------------------------------------------------------------------
      DIMENSION          REC(3,6),REEL(3,6),WAAR(6)
      COMMON /LINES / IDOC,INP,IOUT,MAX,OBS(40),CEL(10)
C
      DO 1100 J=1,6
 1100 WAAR(J)= CEL(J)
      DO 1220 J=1,3
      REC(1,J)=SQRT(WAAR(J))
      REC(1,J+3)=REC(1,J)
 1220 CONTINUE
      DO 1230 J=1,3
      W=WAAR(J+3)/(REC(1,J+1)*REC(1,J+2)*2.0)
      REC(2,J+3)=W
      REC(2,J)=W
      REC(3,J+3)=SQRT(1.0-W*W)
      REC(3,J)=REC(3,J+3)
 1230 CONTINUE
      DO 1240 J=1,3
      W=(REC(2,J+1)*REC(2,J+2)-REC(2,J))/(REC(3,J+1)*REC(3,J+2))
      REEL(2,J+3)=W
      REEL(2,J)=W
      REEL(3,J)=SQRT(1.0-W*W)
      REEL(3,J+3)=REEL(3,J)
 1240 CONTINUE
      FNNN=REEL(3,1)*REEL(3,2)*REC(3,3)
      DO 1250 J=1,3
      CEL(J)=REEL(3,J)/(FNNN*REC(1,J))*100.0
 1250 CONTINUE
      DO 1260 J=4,6
      CEL(J)=90.0-57.29578*ATAN(REEL(2,J)/REEL(3,J))
 1260 CONTINUE
      VOLUM=FNNN*CEL(1)*CEL(2)*CEL(3)
      IF (ISW .EQ. 0) RETURN
      CEL(7)=VOLUM
      WRITE(IOUT,9930) (CEL(J),J=1,6) ,VOLUM
      IF (WMOL .LT. 0.1) RETURN
      PARTD=WMOL/(0.6025*VOLUM)
      IF (DOBS .LT. 0.1) GO TO 1270
      J=IFIX(DOBS/PARTD)
      DMAX=1000.
C  J IS A GOOD APPROXIMATION OF THE NUMBER OF MOLECULES PER UNIT CELL.
C  TRY A FEW VALUES AROUND J  AND RETAIN THE ONE WITH THE CALCULATED
C  DENSITY THAT FITS BEST.
      DO 1265 I=1,4
      DX=FLOAT(I+J-2)*PARTD
      DELTA=ABS(DOBS-DX)
      IF (DELTA .GT. DMAX) GO TO 1265
      DMAX=DELTA
      N=I+J-2
 1265 CONTINUE
      DX=PARTD*FLOAT(N)
      WRITE(IOUT,9931)  N,DX,DOBS
      RETURN
 1270 WRITE(IOUT,9932) PARTD
      RETURN
 9930 FORMAT(1H ,6F9.3,  F12.2)
 9931 FORMAT(1H+,68X, 5H   Z=,I2,6H   DX=, F7.4, 8H   DOBS=, F7.4 )
 9932 FORMAT(1H+,68X,11H   DCALC=Z*, F8.4)
      END

      SUBROUTINE WISSEL(I,J)
      COMMON /LINES / IDOC,INP,IOUT,MAX,OBS(40),CEL(10)
      A=CEL(I)
      CEL(I)=CEL(J)
       CEL(J)=A
      A=CEL(I+3)
      CEL(I+3)=CEL(J+3)
      CEL(J+3)=A
      RETURN
      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

      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
