CDC   PROGRAM ITO(INPUT,OUTPUT,WRIT,TAPE5=INPUT,TAPE6=OUTPUT,
CDC  2 TAPE4=WRIT)
C  ---------------------------------------------------------------------
C
C             THIS PROGRAM WAS WRITTEN BY
C                        J.W. VISSER
C             OF THE TECHNISCH PHYSISCHE DIENST TNO-TH
C             P.O.BOX 155 DELFT NETHERLANDS.
C             IT CONTAINS MANY IDEAS AND SUGGESTIONS FROM
C                       PROF.DR. P.M. DE WOLFF
C             OF THE TECHNICAL UNIVERSITY (T.H) DELFT.
C
C
C             A DESCRIPTION OF THE PROGRAM WAS PUBLISHED IN
C             JOURNAL OF APPLIED CRYSTALLOGRAPHY.(1969) 2,89.
C
C
C             THIS WORK WAS LARGELY SPONSORED BY THE
C             JOINT COMMITTEE ON POWDER  DIFFRACTION STANDARDS
C             NOW AT  1601 PARK LANE, SWARTHMORE, PA. 19081
C             AND THE CENTRAL ORGANIZATION FOR  TNO, THE HAGUE,
C             NETHERLANDS.
C
C
C     THE PURPOSE OF THE PROGRAM IS TO FIND THE UNIT CELL FROM
C     THE POWDER PATTERN  OF COMPOUNDS OF LOW SYMMETRY.
C     ONE- AND TWO- PARAMETER PROBLEMS CAN BE SOLVED BY SIMPLER
C     METHODS.
C
C
C     THANKS ARE DUE TO   ROBIN SHIRLEY    OF THE UNIVERSITEY OF SURREY
C     WHO CONTRIBUTED MANY EXCELLENT IDEAS  AND ALSO MADE THE OUTPUT
C     MUCH BETTER UNDERSTANDABLE,
C                   AND
C     DR. GERALD G. JOHNSON JR. OF THE PENNSYLVANIA STATE UNIVERSITY
C     FOR MANY EXCELLENT IDEAS, SOME AS A RESULT FROM FEEDBACK BY USERS
C     IN THE U.S.A. TO WHOM HE DISTRIBUTED  THIS PROGRAM,   AND FOR
C     FORCING ME TO WRITE THIS PROGRAM IN THE PIDGIN FORTRAN CALLED ANSI
C
C
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
C        THIS PROGRAM HAS PASSED THROUGH SEVERAL VERSIONS AND HAS
C     FOUND THE CORRECT UNIT CELL IN A LARGE NUMBER OF CASES. IT IS
C     POSSIBLE, HOWEVER, THAT UNDER SOME CIRCUMSTANCES, WITH A PARTICULA
C     SET OF DATA, OR WHEN USED WITH A PARTICULAR COMPUTER SYSTEM,
C     ERRONEOUS RESULTS WILL BE OBTAINED.
C     THE AUTHOR, THE  TECHNISCH PHYSISCHE DIENST TNO-TH, OR THE
C     JOINT COMMITTEE ON POWDER DIFFRACTION STANDARDS, THEREFORE, TAKE
C     NO RESPONSIBILITIES WHATSOEVER FOR  THE CORRECTNESS OF THE
C     PROGRAM OR THE RESULTS.
C        IT WOULD BE APPRECIATED IF ANY DISCREPANCIES IN THE PROGRAM
C     ARE BROUGHT TO THE ATTENTION OF THE AUTHOR, FOR IMPROVEMENTS IN
C     LATER VERSIONS.
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C                       VERSION 11, AFTER   MAY   1984
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C        VERSION 1
C             TRIALS TO MAKE THE  FORTRAN  PROGRAM EQUIVALENT TO
C             THE PL/I PROGRAM.
C        VERSION 2
C             PROVISIONS FOR MOLECULAR WEIGHT AND DENSITY.
C             PRINT-OUT OF INPUT CARDS.
C        VERSION 3
C             ADDITION OF THE PROCEDURE SYMTST (TESTING ON
C             HIGHER SYMMETRIES) AND LESSER IMPROVEMENTS.
C        VERSION 4.
C             ADAPTATIONS TO AVOID DIFFICULTIES ENCOUNTERED ON
C             OTHER MACHINES (INCLUDING OTHER IBM-MACHINES).
C             INTRODUCTION OF SUBROUTINE  ECHO IN CONNECTION
C             WITH BACKSPACE DIFFICULTIES.
C             MINOR CHANGES IN THE OUTPUT.
C             SUBROUTINES IN ALFABETICAL  ORDER.
C        VERSION 5
C             INSTRUCTION MANUAL
C             IMPROVED OUTPUT TEXT.
C             MORE TESTS ON INPUT.
C             TESTING ON HIGHER SYMMETRIES  IMPROVED.
C             AN EFFORT HAS BEEN MADE TO COMPLY WITH THE REQUIREMENTS
C             OF  *ANS* FORTRAN.  THIS WAS NOT EVERYWHERE POSSIBLE.
C        VERSION 6
C             MANY WISHES FROM USERS FULFILLED.
C             DISCLAIMERS, BOTH IN COMMENT ON MAIN AND IN OUTPUT.
C             MIXED-MODE ARITHMETIC REMOVED.
C             SUBROUTINE WRISP REWRITTEN FOR FASTER PRINTING.
C             POSSIBILITY TO PRINT AN INSTRUCTION MANUAL ON THE
C             SPECIAL FEATURES OF THE PROGRAM.
C             POSSIBILITY TO ENTER  INTENSITIES WITH THE LINES.
C             POSSIBILITY TO PRODUCE A LIST OF ALL CALCULATED AND
C             OBSERVED LINES.
C             POSSIBILITY TO ENHANCE CERTAIN CRYSTAL SYSTEMS. DUE TO
C             THE GENERAL CONSTRUCTION OF THE PROGRAM, THIS IS ONLY
C             POSSIBLE FOR ORTHORHOMBIC,MONOCLINIC OR TRICLINIC
C             LATTICES.
C             PROBABILITY TO  *HIT*  A LINE IS SUBDIVIDED OVER THE
C             Q RANGE.   CHANGES IN  CRITER.
C             ERROR LIMITS ON OBSERVED Q'S CALCULATED ONCE.  DISPLAYED
C             WHEN TEST=.TRUE.
C             OVERLAYS  TESTED.
C        VERSION 7
C             ALL  A-FORMATS REDUCED TO A1, EXCEPT THOSE THAT PRINT
C             INTENSITIES IN CHARACTER FORM (A3)-  THESE OCCUR  ONLY
C             IN *READER*, *LISTM* AND *REFINE* .
C             TRIES TO INDEX ALL LINES WITH THE UNIT CELL OF
C             THE BEST SOLUTION.
C             INSTRUCTIONS REWRITTEN, OTHER EXPLANATIONS EXTENDED.
C             THE USE OF ONE SET OF HKL'S TO EXPLAIN MORE THAN
C             ONE OBSERVED LINE  IS FLAGGED.
C             NUMBER OF SOLUTIONS THAT  ARE PRINTED OUT EXTENSIVELY
C             CAN BE CHANGED.
C        VERSION 8
C             INPUT FORMAT CHANGED-- ADDITIONAL PARAMETER CARD.
C             ZERO SHIFT (SPECIMEN DISPLACEMENT+TRANSPARENCY
C             CORRECTION ON DIFFRACTOMETER)--THE NEW PARAMETER CARD
C             STARTS WITH THE ZERO CORRECTION (ZERCOR), TO BE
C             APPLIED TO ALL ANGLES.
C             ZERCHK-- THE PROGRAM APPLIES 20 DIFFERENT ZERO SHIFTS
C             AND GIVES FOR EACH SHIFT THE NUMBER OF HIGHER ORDER
C             LINES THAT FIT--   IN THIS WAY IT SHOULD BE POSSIBLE
C             TO DECIDE WHETHER A ZERO SHIFT IS NECESSARY.
C             ONLY SOLUTIONS WITH A FIGURE OF MERIT THAT IS BETTER
C             THAN 4.0 ARE PRINTED--YOU MAY SET YOUR OWN MINIMUM
C             WITH PRNTMR.
C             WITH PRNTLN YOU CAN SET YOUR OWN MINIMUM OF INDEXED
C             LINES THAT IS NEEDED IN ORDER THAT THE LATTICE
C             WILL BE PRINTED.
C             FOR TRIAL, IT IS NOW POSSIBLE TO ENTER A DIRECT CELL.
C             COMBINATIONS OF ZONES WITH ITSELF NOW TRIES BOTH AXES
C             AS 'COMMON' AXES.
C             THE NUMBER OF CYCLES OF LATTICE REFINEMENT IS INCREASED.
C             *REFPAR* IS SPLIT INTO TWO SUBROUTINES,ONE OF WHICH
C             (REFPAR) SETS THE PARAMETERS FOR THE OTHER (REFINE)--
C             THE LATTER SUBROUTINE DOES THE ACTUAL CALCULATIONS.
C             TOLERANCES  ON FITTING LINES NOW EXPRESSED IN TWO THETA.
C             A(3) FORMATS CHANGED INTO A(2),1X.
C        VERSION 8.5
C                  THE MAIN PROGRAM IS NOW SPLIT INTO 3 PARTS,AABITO,
C             AAITO1 AND ANITO2--  THE NEW MAIN PROGRAM (AABITO)
C             COMMENTS, INITIALISES SOME VALUES AND CALLS READER, AAITO1
C             AND ANITO2.
C                  ZERCHK IS NOW OPTIONAL (COL. 33 IN PARAMETER CARD).
C                  THE SUBROUTINE *NTEST* LIMITS THE TEST OUTPUT TO
C             A SMALL PART OF THE PROGRAM, OR PRINTS ALL TEST OUTPUT
C             AFTER A CERTAIN POINT-- THE COMMENT IN THE SUBROUTINE
C             GUIDES YOU TO THE ACTUAL CALLS--  THE SUBROUTINE IS ONLY
C             USEFUL WHEN YOU ARE CHANGING PARTS OF THE PROGRAM.
C                  LATTICES THAT ARE DISCARDED IN *NORM3* ARE NOW
C             EXPLICITLY NAMED--  ALSO NAMED ARE THE LATTICES THAT ARE
C             DISCARDED BECAUSE THEY ARE  *EQUAL* TO OTHERS.
C        VERSION 9
C             SUBROUTINE READER WAS SPLIT IN TWO PARTS, READER  AND
C             REDLIN. THIS LAST SUBROUTINE READS THE INPUT LINES. IT
C             ALSO CALCULATES THE MAXIMUM ERROR (IN TWO THETA) THAT
C             YOUR INPUT DATA COULD HAVE IN ORDER TO BE ABLE TO PASS
C             THE TESTS (FIGURE OF MERIT .GT. 10)--
C             TWICE THIS MAXIMUM ERROR IS USED AS THE DEFAULT TOLG,
C             BUT ONLY IF TOLG REMAINS LESS THAN 0.060 DEGREES. --
C             ENTERING YOUR OWN TOLG OVERRIDES THE DEFAULTS.--
C             IN THE LIST OF ALL CALCULATED AND OBSERVED LINES,
C             THE UNINDEXED LINES ARE NOW MARKED CONSPICUOUSLY (IDEA
C             FROM CAM HUBBARD)--
C             THE FATE OF THE POSSIBLE LATTICES IS PRINTED IN FULL,
C             BUT ONLY IF YOU PRINT INTERMEDIATE RESULTS (ENL=TRUE)--
C             THE MINIMUM VALUES OF THE RECIPROCAL LATTICE CONSTANTS
C             NOW ALLOWS THE FIRST LINE TO BE A FOURTH ORDER. IN AN
C             EXPERIMENTAL PATTERN, WHICH WAS ORTHORHOMBIC  WITH
C             SPACEGROUP FDD2, THIS SITUATION ACTUALLY OCCURRED ---
C        VERSION 10  (A)
C                          FORTRAN 77
C             MODIFIED TO CONCUR WITH FULL FORTRAN 77.THIS ALLOWED A
C             SIMPLER VERSION OF SUBROUTINE ECHO.--
C             EXPLANATION OF THE VARIABLES IN COMMON BLOCKS IS
C             GIVEN IN THIS PROGRAM (BELOW). --
C             INDIVIDUAL PROBABILITIES FOR EACH LINE (FOR FINZON).--
C             CRITER CHANGED TO CONCUR WITH EARLIER IDEAS.--
C             SORTING OF ZONES IN THE NEW SUBROUTINE  *SRTZON* --
C               THE OUTPUT OF EVALUATE IS SORTED ON THE CRITERION
C             CRIT,  BUT IF NO CORRECT LATTICE IS OBTAINED, IT IS USED
C             AGAIN BUT SORTED ON THE NUMBER OF OBSERVED LINES (NO)
C             THAT MATCH THE ZONE. A SUBSORT IS MADE ON THE NUMBER OF
C             CALCULATED LINES (NC). THIS SOMETIMES LEADS TO SUCCES.--
C               THE PAPER BY MIGHELL AND STALICK IN 'ACCURACY IN POWDER
C             DIFFRACTION' HAS CONVINCED ME THAT IT IS NECESSARY TO
C             ENSURE THAT ALL AVAILABLE DATA (UP TO 40 LINES) IS
C             ENTERED INTO THE PROGRAM.
C             IF YOU NOW WANT TO ENTER LESS THAN 32 LINES, YOU GET A
C             WARNING AND YOU HAVE TO ENTER THE CORRECT NUMBER OF LINES
C             ON THE SECOND PARAMETER CARD.--
C                  A TRIAL WAS MADE IN REFINING THE ZONES ON THETA
C             INSTEAD OF Q. THIS DID NOT WORK AS WELL AS REFINING ON Q.
C             THEREFORE, ONLY THE BEST ZONES ARE (AT THE END OF THE
C             EVALUATION PROCESS) REFINED ON THETA IN ORDER TO
C             ESTABLISH  A ZERO-SHIFT FOR EACH ZONE.
C             SINCE I DO NOT KNOW YET HOW USEFUL THIS CAN BE, NOTHING
C             MUCH IS BEING DONE WITH IT.  IF YOU WANT TO, YOU CAN
C             USE THE INFORMATION BY ENTERING A CONVINCING
C             TWOTHETA-SHIFT INTO THE PROGRAM (REVERSE THE SIGN STATED
C             AND ENTER ON SECOND PARAMETER CARD IN COL 1-9). --
C                 LATTICES ARE ALSO REFINED ON TWOTHETA, TOGETHER
C             WITH A ZERO SHIFT. THE ZERO SHIFT IS STORED IN
C             ROSTER(I,11) AND PRINTED.
C
C        VERSION 11 (B)
C             ECHO CHANGED FOR THE USE OF *IOSTAT*. --
C             PRINTING OF Q-SCHEMES FOR THE LATTICES ONLY
C             WHEN TEST=TRUE. --
C
C
C
C
C
C
C
C
C
C
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C-----------------------------------------------------------------------
C  THE INPUT OF THE PROGRAM IS ON UNIT 5,
C  THE OUTPUT IS ON UNIT  6.
C  THE SUBROUTINE  *ECHO* WRITES THE CARD IMAGES OF THE INPUT ON
C  AN EXTRA OUTPUT FILE,  NR. 4 (IDOC=4).
C  THIS EXTRA FILE IS ONLY USED IN *ECHO*
C
C     THE INPUT OF THIS PROGRAM SHOULD CONSIST OF:
C     1. A TITLE CARD,CONTAINING UP TO 80 CHARACTERS
C     2. TWO PARAMETER CARDS,THAT CAN BE COMPLETELY BLANK
C     3. A NUMBER OF LINES CARDS, CONTAINING THE POSITIONS OF THE LINES
C        THESE CAN BE ENTERED AS D-,Q- OR TWOTHETA VALUES. ( 8 F10.5 )
C        YOU MAY LEAVE ANY NUMBER OF POSITIONS PER CARD  BLANK.
C        THIS MEANS THAT YOU CAN ALSO ENTER ONLY ONE LINE PER CARD,
C        IN ANY ONE OF THE POSITIONS OF 8F10.5
C     4. A BLANK CARD, TO FINISH THE INPUT OF ONE PROBLEM.
C
C        REPEAT THE ABOVE  AS MANY TIMES AS PROBLEMS,   -- THEN --
C     5. A CARD CONTAINING THE WORD  --  END  -- IN THE THREE FIRST
C        COLUMNS,  TO FINISH THE INPUT.
C
C  THE PROGRAM WILL TAKE ITS OWN (DEFAULT) PARAMETERS.
C  FOR THE INITIATED THE SUBROUTINE READER GIVES FULL DETAILS ABOUT
C  ALL THE  POSIBILITIES OF THE PROGRAM.
C
C
C-----------------------------------------------------------------------
C              THE SUBROUTINES  ARE  LISTED  IN  ALFABETICAL  ORDER.
C-----------------------------------------------------------------------
C
C                         GENERAL FLOW OF SUBROUTINES
C          (THE ACTUAL FLOW OF SUBROUTINES IS A BIT MORE COMPLICATED)
C-----------------------------------------------------------------------
C      SUBROUTINES TIMER AND NTEST ARE FREQUENTLY CALLED BY THE
C      MAIN PROGRAM. NTEST CALLS ARE MARKED BY 1* , 2* , ETC.
C.......................................................................
C
C
C        *******            *******           **********
C        *     *            *     * <-------> *  ECHO  *        ********
C        *  M  *            *  R  *           **********   !--> *SRTRON*
C        *     *            *     * <-------> * REDLIN * <-!    ********
C      1**  A  * <--------> *  E  *           **********   !--> *TETCAL*
C        *     *            *     * <-------> * RECCEL *        ********
C        *  I  *            *  A  *           **********
C        *     *            *     * <-------> * MANUAL *
C        *  N  *            *  D  *           **********
C   !--->*     *            *     * <-------> * INSTRS *
C   !    *******            *  E  *           **********
C   !       !               *     * <-------> * ZERCHK *
C   !       !               *  R  *           **********
C   !       !               *     * <-------> *        *
C   !       !               *******           **********
C   !       !
C   !  2********            **********
C   !    *     *  <-------> * FINZON *
C   !    *     *            **********
C   !    *  A  *
C   !    *     *            *******           ******         *********
C   !    *     *            *     *  <------> * LS *  <----> * WRISP *
C   !    *     *            *  E  *           ******         *********
C   !    *  A  *            *  V  *           ********
C   !    *     *            *  A  *  <------> * NORM *
C   !    *     *  <-------> *  L  *           ********
C   !    *     *            *  U  *           **********
C   !  3**  I  *            *  8  *  <------> * CRITER *
C   !    *     *            *     *           **********
C   !    *     *            *     *           *********
C   !    *     *            *******  <------> * TEXT1 *
C   !    *     *                              *********
C   !    *  T  *            **********
C   !    *     *  <-------> * SRTZON *
C   !    *     *            **********
C   !  4**     *
C   !    *  O  *            **********
C   !    *     *  <-------> * ONEZON *
C   !    *     *            **********
C   !    *     *
C   !    *  1  *
C   !    *     *            ********
C   !  5**     *  <-------> * REAR *
C   !    *******            ********
C   !
C   !
C   !
C   !    *******            **********
C   !--->*     *  <-------> * THREED *
C        *     *            **********
C      6**     *
C        *     *            *********          **********
C        *     *  <-------> * NORM3 * <------> * WISSEL *
C        *  A  *            *********          **********
C        *     *
C      7**     *
C        *     *
C      8**     *        *******        *******           **********
C        *  N  *        *     *        *     *  <------> *        *
C        *     *        *  R  *        *  R  *           **********
C        *     *        *     *        *     *           *********
C        *     *        *  E  *        *  E  *  <------> * WRISP *
C        *     *        *     *        *     *           *********
C        *  I  * <----> *  F  * <----> *  F  *           *********
C        *     *        *     *        *     *  <------> * ARRAY *
C        *     *        *  P  *        *  I  *           *********
C        *     *        *     *        *     *           ********
C        *     *        *  A  *        *  N  *  <------> * MINV *
C        *  T  *        *     *        *     *           ********
C        *     *        *  R  *        *  E  *           *********
C        *     *        *     *        *     *  <------> * LISTM *
C        *     *        *******        *******           *********
C      9**     *
C        *  O  *
C        *     *            **********
C        *     *  <-------> * UNTCEL *
C        *     *            **********
C        *     *
C        *  2  *
C        *     *            **********
C        *     *  <-------> * SYMTST *
C        *     *            **********
C        *     *
C        *     *
C        *******
C
C-----------------------------------------------------------------------
C                       MINIMISING STORAGE REQUIREMENTS.
C-----------------------------------------------------------------------
C                       ONLY VALID FOR VERSION 6.
C                FOR NEWER VERSIONS IT CAN SERVE AS A GUIDE.
C.......................................................................
C
C     1. OVERLAY.
C
C             THE SUGGESTED OVERLAY STRUCTURE NEEDS ONLY ABOUT 60%
C             OF THE STORAGE NEEDED FOR THE FULL PROGRAM,
C             WITHOUT IMPAIRING ITS EFFICIENCY.
C             THE SUBROUTINES  *FINZON*  AND *THREED*  HAVE THE
C             LARGEST REQUIREMENTS FOR MAIN STORAGE.    THEREFORE, IT IS
C             NOT NECESSARY TO GO BEYOND ONE LEVEL OF OVERLAY.
C             I SUGGEST THE FOLLOWING SCHEME
C
C             ENTRY  MAIN      COMMENT* ROOT SEGMENT, SEGMENT 1
C             OVERLAY ONE
C             INCLUDE SUB(READEM)     COMMENT*  SEGMENT 2
C             OVERLAY ONE
C             INCLUDE SUB(FINZON)     COMMENT* SEGMENT 3
C             OVERLAY ONE
C             ETC.  ETC.  ETC.
C
C             SEGMENT     CONTENTS
C                1        MAIN + COMMONS, EXCEPT *COUNTM*
C                2        READER+ECHO+MANUAL+INSTRS+SRTRON
C                3        FINZON
C                4        EVALU8+LS+WRISP+NORM+CRITER+TEXT1
C                5        ONEZON+REAR
C                6        THREED
C                7        NORM3
C                8        REFPAR+WISSEL+ARRAY+MINV+LISTM+ *PRINS*
C                9        UNTCEL+SYMTST
C
C             THE FOLLOWING CHANGES ARE NECESSARY
C             A.DELETE THE COMMON *COUNTM* THROUGHOUT THE PROGRAM,
C               USE AN ARRAY  *ITL(5000)* INSTEAD OF THIS COMMON
C               IN THE SUBROUTINES *FINZON* AND *THREED*
C             B.DELETE THE SUBROUTINE CALL FOR *NORM3* FROM *REFPAR*
C               INSERT A CARD      95 CONTINUE
C             C.CHANGE THE NAME OF THE SUBROUTINE *WRISP* IN *REFPAR*
C               INTO SOMETHING DIFFERENT, LIKE *PRINS*  AND CHANGE THE
C               FOUR SUBROUTINE CALLS FOR *WRISP* IN *REFPAR*
C               ACCORDINGLY
C             THE CALL FOR *NORM3* IN *REFPAR* IS ONLY SELDOM
C             REALLY NECESSARY,  HOWEVER,
C             IF THE RESULTS OF A PARTICULAR POWDER PATTERN SEEM STRANGE
C             TO YOU, YOU CAN ALWAYS FORCE ANOTHER APPLICATION OF *NORM3
C             BY RE-ENTERING YOUR DATA AND THE STRANGE LATTICES,
C             SEE  MANUAL AND INSTRUCTIONS, USE COL 44 OF PARAMETER CARD
C
C
C     1B OVERLAY 2.
C
C             ON AN IBM 370/158 THE PROGRAM RAN IN  84K INSTEAD OF 140
C             AND WITHOUT MEASURABLE LOSS OF SPEED WITH THE FOLLOWING
C             SCHEME
C             A. COMPILE THE PROGRAM AS GIVEN WITH  PARM.LKED=OVLY
C             B.GIVE  LKED.SYSIN  DD  *
C                ENTRY MAIN
C                INSERT MAIN
C                OVERLAY  ALPHA
C                  INSERT READER,MANUAL,ECHO,INSTRS,SRTRON,IHCSSCN
C                OVERLAY  ALPHA
C                  INSERT FINZON,THREED,COUNTM
C                OVERLAY  ALPHA
C                  INSERT EVALU8,LS,NORM,CRITER,TEXT1,EVALN,LINSTR,IHCSE
C                  INSERT REFPAR,WISSEL,WRISP,ARRAY,MINV,NORM3,LISTM
C                OVERLAY  ALPHA
C                  INSERT ONEZON,REAR,UNTCEL,SYMTST
C             /*
C
C
C     2. CHOPPING UP THE PROGRAM.
C
C             IF YOU HAVE A REALLY SMALL COMPUTER, IT IS STILL
C             POSSIBLE TO RUN THIS PROGRAM.
C             AS THE PROGRAM USES ITS MAIN SUBROUTINES SEQUENTIALLY,
C             IT IS POSSIBLE TO IMITATE AN OVERLAY AND DIVIDE THE MAIN
C             PROGRAM IN A NUMBER OF PIECES THAT WILL WORK AS
C             SEPARATE UNITS.  AFTER EVERY PART IT WILL BE NECESSARY
C             TO SAVE THE VALUES OF MOST COMMON AREAS AND OF THE
C             ARRAY  *STORE*.    THE DIVISION OF THE PROGRAM SHOULD
C             FOLLOW THE SAME MAIN LINE  AS THE  OVERLAY STRUCTURE.
C             I HAVE NOT TRIED THIS MYSELF, AS I HAVE A VERY LARGE
C             COMPUTER AVAILABLE.  I AM QUITE WILLING TO HELP ANYBODY
C             WHO TRIES THE ABOVE IDEA.   PLEASE CONTACT ME .
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C                EXPLANATION OF VARIABLES IN  COMMON.
C
C     1.LINES
C       IDOC     OUTPUT FILE FOR DOCUMENTATION PURPOSES. LISTS INPUT
C                AS GIVEN
C       INP      INPUT UNIT,FROM WHICH ALL INPUT IS READ.
C       IOUT     STANDARD OUTPUT UNIT (LINE PRINTER)
C       MAX      THE NUMBER OF OBSERVED LINES THAT ARE USED (.LE. 40)
C       OBS(40)  CONTAINS THE Q-VALUES OF THE OBSERVED LINES.
C       CEL(10)  USED FOR SEVERAL PURPOSES, CAN CONTAIN UNIT CELLS
C                OR LINE COMBINATIONS OR ZONES.
C    2. CHARCM
C       NAME     NAME OF PROBLEM.
C       CARD     CONTAINS THE IMAGE OF THE 'CARD' THAT WAS LAST READ.
C       INTSY(50) CONTAINS THE INTENSITIES (GIVEN AS  VS , S , ETC.)
C                OF THE OBSERVED LINES.
C    3. PARMS
C       TOL2     TOLERANCE ON 2-DIMENSIONAL SEARCH FOR THE ANGLE BETWEEN
C                THE DIRECTIONS OF THE POINTS REPRESENTED BY Q-VALUES.
C       TOL3     SAME FOR THREE-DIMENSIONAL SEARCH (FOR D)
C       WAVEL    WAVELENGTH  (1.5406 BY DEFAULT)(1.5405891)
C       FMINEQ   MAXIMUM VALUE FOR THE DIFFERENCE BETWEEN 'EQUAL' VALUES
C       FMINIM   MINIMUM VALUE FOR THE FIRST THREE LATTICE CONSTANTS.
C       SELF     LOGICAL
C       ENL      LOGICAL FOR PRINTING INTERMEDIATE RESULTS.
C       TEST     LOGICAL FOR PRINTING ALL TEST RESULTS.
C       LIST     THE NUMBER OF LATTICES THAT ARE TO BE PRINTED OUT
C       NSYST(3) SEE INSTRUCTIONS. POSSIBILITY TO ENHANCE SOME CRYSTAL
C                SYSTEMS
C       OBSMT(40) Q-VALUES OF THE OBSERVED LINES MINUS TOLERANCE.
C       OBSPT(40) Q-VALUES OF THE OBSERVED LINES PLUS THE  TOLERANCE.
C       TWTHET(40) TWO-THETA'S OF THE OBSERVED LINES.
C       TOLG     TOLERANCE (IN HUNDREDTH OF A DEGREE  TWO THETA)
C
C    4. EVALN
C       NORDER   NUMBER OF ORDERS OF AXES' Q-VALUES, USUALLY 4.
C                NORDER DETERMINES THE INDEX-FIELD OF THE ZONES AT
C                EVALUATION.
C       RECTAN   LOGICAL THAT IS .TRUE. WHEN THE ZONE IS RECTANGULAR.
C       PROB(40) PROBABILITY TO 'HIT' AN OBSERVED LINE BY CHANCE.
C       EVAL(-8:8,0:8) ARRAY CONTAINING THE PROBABILITIES OF THE
C                CALCULATED LINES THAT ARE ALSO OBSERVED.
C       UPBNDR(40) UPPER BOUNDARY OF THE PROBABILITY FIELD OF EACH LINE.
C
C    5. TEMRES
C       LINCO    THE NUMBER OF LINE COMBINATIONS
C       NQ2      THE NUMBER OF THE LAST LINE TO BE COMBINED WITH OTHERS.
C       NINDEX   THE NUMBER OF LINES INDEXED AMONG THE FIRST NQ2 LINES
C       INDEX(40) ARRAY TO COUNT WHETHER A LINE WAS INDEXED BY A
C                PARTICULAR ZONE.
C       ROSTER(150,12) ARRAY IN WHICH THE USEFUL LATTICES ARE STORED.
C
C    6. CONSTS
C       RAD,TWORAD FACTORS FOR CALCULATING DEGREES FROM RADIANS.
C       RRAD     FACTOR TO CALCULATE RADIANS FROM DEGREES (RECIPROCAL)
C       HLFLAM   WAVELENGTH/(2*100). FACTOR 100 STEMS FROM THE FACT
C                THAT Q HAS BEEN MULTIPLIED BY 100**2
C       WAVSQR   CONSTANT NEEDED FOR CALCULATING DERIVATIVES.
C       STORE(200,9) ARRAY TO STORE THE ZONES.
C
C    7. THETAS
C       ZERTOT   TOTAL ZERO SHIFT OF THE DATA (RADIANS).
C       OBSTET   OBSERVED THETA.
C       SHFTET   OBSERVED THETAS, SHIFTED BY THE CALCULATED ZEROSHIFT.
C       TETPT    OBSERVED THETA PLUS TOLERANCE
C       TETMT    OBSERVED THETA MINUS TOLERANCE.
C
C    8. COUNTM
C       ITL(0:5000) ARRAY TO FIND THE VALUES OF  *F*  AND  *D* .
C                IT IS NOT NECESSARY TO PUT THE ARRAY(S) IN A COMMON,
C                BUT IT SAVES SPACE.
C
C
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
	LOGICAL SELF,ENL,TEST,SOL,RECTAN
	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 /EVALN / NORDER,RECTAN,PROB(40),EVAL(-8:8,0:8),
     *                UPBNDR(40)
	COMMON /TEMRES/ LINCO,NQ2,NINDEX,INDEX(40),ROSTER(150,12)
	COMMON /COUNTM/ ITL(0:5000)
	IDOC=4
	INP=5
	IOUT=6
	NORDER=4
	NRMAX=200
1	CONTINUE
	NRRUN=0
	DO 10 I=1,150
	   DO 5 J=1,12
5	   ROSTER(I,J)=0.0
10	CONTINUE
	CALL TIMER(0)
C
	CALL READER(NZ1,NZ2,NR,NSOLMX,WMOL,DOBS,PRNTMR,PRNTLN,NTST)
C
	IF (INP .EQ. 12345)  GO TO 9000
	CALL TIMER(1)
100	SOL=.TRUE.
	IF (NR .EQ. 0) THEN
	   NRRUN=NRRUN+1
	   SOL=.FALSE.
	   CALL AAITO1(NZ1,NZ2,NR,NRMAX,NTST,NRRUN)
	END IF
	CALL ANITO2(SOL,NR,NSOLMX,WMOL,DOBS,PRNTMR,PRNTLN,NTST,NRRUN)
	IF (NRRUN .EQ. 1) THEN
	   NR=0
	   GO TO 100
	ELSE
	   GO TO 1
	END IF
9000	STOP
	END

	SUBROUTINE AAITO1(NZ1,NZ2,NR,NRMAX,NTST,NRRUN)
C
	DIMENSION TEMP(0:41)
	LOGICAL SELF,ENL,TEST,RECTAN
	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 /EVALN / NORDER,RECTAN,PROB(40),EVAL(-8:8,0:8),
     *                UPBNDR(40)
	COMMON /TEMRES/ LINCO,NQ2,NINDEX,INDEX(40),ROSTER(150,12)
	COMMON /CONSTS/ RAD,TWORAD,RRAD,HLFLAM,WAVSQR,STORE(200,9)
	IF (NRRUN .EQ. 2) THEN
	   DO 9 J=1,12
	       DO 8 I=1,150
8	       ROSTER(I,J)=.0
9	   CONTINUE
	   DO 10 I=1,200
	   IF (STORE(I,1) .LT. 0.1) GO TO 15
10	   CONTINUE
15	   NR=I-1
	   NSRT=7
	   NZ1=NZ1+2
	   NZ2=NZ2+2
	   WRITE(IOUT,9922)
9922	   FORMAT(1H0,63(1H*)/' *  PROBLEM IS RERUN WITH A DIFFERENT',
     1    ' SORTING OF THE ZONES  ****'/ 1H ,63(1H*) )
	   GO TO 500
	END IF
	DO 30 I=1,200
	DO 20 J=1,9
20	STORE(I,J)=0.0
30	CONTINUE
	DO 40 I=1,MAX
40	TEMP(I)=OBS(I)
	TEMP(0)=0.5*OBS(1)
	TEMP(MAX+1)=2.0*OBS(MAX)-OBS(MAX-1)
	DO 50 I=1,MAX
	UPBNDR(I)=0.5*( TEMP(I)+TEMP(I+1) )
	PROB(I)=2.0*(OBSPT(I)-OBSMT(I))/(TEMP(I+1)-TEMP(I-1))
	IF (PROB(I) .GT. 0.5) PROB(I)=0.5
50	CONTINUE
	IF (TEST) WRITE(IOUT,9915) UPBNDR
	IF (TEST) WRITE(IOUT,9916) PROB
	CALL NTEST(NTST,2)
	NN=0
C-----------------------------------------------------------------------
C--FINDING ZONES STARTS HERE.
C-----------------------------------------------------------------------
2	IF (ENL)  THEN
	   WRITE(IOUT,9909)
	   IF (NN .EQ. 0) WRITE(IOUT,9901)
	   WRITE(IOUT,9905)
	END IF
	NN=NN+1
	NR=0
	DO 350 J=1,LINCO
	IF (ROSTER(J,3) .LT. 0.05  .AND.  ROSTER(J,4) .LT. 0.05)GO TO 320
	NR=NR+1
	DO 310 I=1,4
310	STORE(  NR,I)=ROSTER(J,I)
	STORE(NR,5)=NR
	IF (ENL) WRITE(IOUT,9903) NR  , (STORE(NR  ,K),K=1,4)
	GO TO 350
320	CALL FINZON(ROSTER(J,1),ROSTER(J,2) )
	DO 330 I=1,5
	N=CEL(I+5)+0.5
	IF (N .EQ. 0) GO TO 350
	IF (NR.LT.NRMAX) GOTO 260
	WRITE(IOUT,9920) NRMAX
	GOTO 360
260	CONTINUE
	NR=NR+1
	STORE( NR ,1)=ROSTER(J,1)
	STORE( NR ,2)=ROSTER(J,2)
	STORE( NR ,3)=CEL(I)
	STORE(NR,4)=.0
	STORE(NR,5)=NR
	IF (ENL) WRITE(IOUT,9904)   NR ,(STORE(  NR ,K),K=1,3),N
330	CONTINUE
350	CONTINUE
C  ---------------------------------------------------------------------
C-- END OF FINDING ZONES.
C--THE EVALUATION OF THE ZONES FOUND STARTS HERE
C  ---------------------------------------------------------------------
360	CONTINUE
	CALL TIMER(2)
	MILNR=NR
	IF (TEST) WRITE(IOUT,9917)
	DO 450 I=1,NR
	IF (I .EQ. 6) TEST=.FALSE.
C--THIS GIVES 5 ZONES (10 PAGES) WITH TEST OUTPUT.
C--IF YOU WANT MORE WASTEPAPER, CHANGE  6  INTO  NR.
C  IF YOU WANT TO SEE THE TEST OUTPUT OF A PARTICULAR ZONE,
C  YOU CAN ENTER THIS ZONE (AND MAYBE SOME OTHERS) WITH LINCO
C  (SEE MANUAL AND INSTRUCTIONS).
	IF (STORE(I,4) .GT. 1.0) GO TO 450
	CALL EVALU8(I,1)
	IF (STORE(I,1) .LT. FMINIM) STORE(I,4)=-10.
	IF (TEST)WRITE(IOUT,9930) I,(STORE(I,K),K=1,8),STORE(I,9)*TWORAD
C  ---------------------------------------------------------------------
C  THE FOLLOWING 7 PROGRAM CARDS CAN BE USED TO STORE THE ORIGINAL ZONE
C  TOO, WHEN ONE OF ITS AXES OR BOTH HAVE BEEN HALVED IN THE EVALUATION
C  PROCESS.   .
C  THE FEATURE IS ONLY SELDOM USEFUL AND ADDS ABOUT 10-20% COMPUTING TIM
C  OF THE PROBLEM. I DO NOT RECOMMEND IT.
C     U=0.0
C     DO  410 J=1,3
C 410 U=U+ABS(STORE(I,J)-CEL(J))
C     IF ( U .LT. 2.0*FMINEQ) GO TO 450
C     MILNR=MILNR+1
C     DO 420 J=1,6
C 420 STORE(MILNR,J)=CEL(J)
C  ---------------------------------------------------------------------
C  ---------------------------------------------------------------------
450	CONTINUE
	NR=MILNR
	CALL NTEST(NTST,2)
C--RESTORES TEST FOR OTHER TEST OUTPUT.
C  ---------------------------------------------------------------------
C--ZONES HAVE BEEN  EVALUATED AND WILL NOW BE SORTED ON QUALITY.
C-----------------------------------------------------------------------
	NSRT=4
500	CALL SRTZON(NSRT,NR,NZ2)
	IF (NR .GT. NZ2+2) NR=NZ2+2
C-----------------------------------------------------------------------
C  EVALU8 GIVES A FINAL ROUND OF LS AND PRINTS THE LATTICE WHEN ENL=.TRU
C  AT THE SAME TIME A CHECK IS MADE TO SEE WHETHER ALL OR NEARLY ALL
C  LINES USED FOR ZONE-FINDING BELONG TO ONE ZONE
C  ATTENTION:  THIS LAST CHECK IS NOT (REPEAT   NOT)  MADE WHEN YOU
C  HAVE ENTERED YOUR OWN LINE COMBINATIONS.
C-----------------------------------------------------------------------
	CALL NTEST(NTST,3)
	N=0
	DO 700 I=1,NR
	CALL EVALU8(I,2)
	IF (STORE(I,1) .LT. FMINIM) GO TO 700
	N=N+1
	DO 640 J=1,9
640	STORE(N,J)=STORE(I,J)
C--CHECK WHETHER THE ZONE INDEXES ALL BUT ONE OF THE FIRST NQ2
C--LINES USED FOR ZONEFINDING.
	IF (NINDEX .LT. NQ2-1) GO TO 700
	IF (NN .GE. 2) GO TO 700
C--FIRST TIME THROUGH (NN). MAKE NEW LINE COMBINATIONS.
	IF (ENL) WRITE(IOUT,9908)  NINDEX,NQ2
	DO 680 J=1,N
	DO 660 K=1,4
660	ROSTER(J,K)=STORE(J,K)
680	CONTINUE
	LINCO=N
	CALL ONEZON(STORE(I,1),STORE(I,2),STORE(I,3))
	WRITE(IOUT,9902)  I,NINDEX,NQ2
	GO TO 2
700	CONTINUE
	NR=N
C---------------------------------------------------------------------
C  REFINE ON TWOTHETA IN ORDER TO FIND A  ZERO SHIFT.
C---------------------------------------------------------------------
	DO 750 I=1,NR
	   IF (STORE(I,7) .GT. 3. ) CALL EVALU8(I,3)
750	CONTINUE
	WRITE(IOUT,9910)
	WRITE(IOUT,9911) (I,(STORE(I,J),J=1,8),TWORAD*STORE(I,9),I=1,NR)
	WRITE(IOUT,9912)
	CALL TIMER(3)
C  ---------------------------------------------------------------------
C  THE SUBROUTINE  REAR  REARRANGES ZONES AND TRIES TO FIND A COMMON
C  LINE IN THE RECIPROCAL LATTICE.
C  A COMMON Q-VALUE CLOSE TO THE ORIGIN IS ALL THAT IS NEEDED.
C  ---------------------------------------------------------------------
800	IF (NZ1 .GT. NR) NZ1=NR
	CALL NTEST(NTST,4)
	IF (NZ2 .GT. NR) NZ2=NR
	N=0
	IF (ENL) WRITE(IOUT,9906)
	DO 1850 J=1,NZ1
	M=0
	DO 1805 I=1,3
1805	CEL(I)=STORE(J,I)
	DO 1840 K=J,NZ2
	IF (K .NE. J) GO TO 1808
	CEL(1)=STORE(J,2)
	CEL(2)=STORE(J,1)
1808	DO 1810 I=1,3
1810	CEL(I+3)=STORE(K,I)
	CALL REAR
	IF (CEL(4).LT. 0.1) GO TO 1840
	N=N+1
	N4=4*N
	L=10*J+K
	DO 1820 I=1,6
1820	ROSTER(N4,I)=CEL(I)
	ROSTER(N4,8)=L
	ROSTER(N4,9)=0.5*(STORE(J,9)+STORE(K,9))
	IF (ENL) WRITE(IOUT,9907) J,K,(CEL(I),I=1,3),CEL(5),CEL(6),L,
     * TWORAD*ROSTER(N4,9)
C  CARE MUST BE TAKEN THAT THERE IS NO OVERFLOW IN *ROSTER*
C  ROSTER CAN CONTAIN 4*36+3 POSSIBLE LATTICES (150 IS MAXIMUM)
	IF (N .EQ. 36) THEN
	   WRITE(IOUT,9925) N,J,K
	   GO TO 1860
	END IF
	DO 1830 I=1,3
	M=M+1
1830	CEL(I)=STORE(J,I)
	IF (M .GT. 4) GO TO 1840
	IF (K .EQ. J) GO TO 1808
1840	CONTINUE
1850	CONTINUE
	IF (N .GE. 6) GO TO 1860
	IF (NZ1 .EQ. NR) GO TO 1860
	NZ1=2*NZ1
	NZ2=2*NZ2
	GO TO 800
1860	LINCO=4*N
	CALL TIMER(4)
	RETURN
9901	FORMAT(32H0THE SEARCH FOR ZONES COMMENCES    )
9902	FORMAT(1H0 , 70(1H*)  /
     2 11H *     ZONE,I3,8H INDEXES,I3,15H LINES FROM THE,I3,
     3 28H USED FOR ZONEFINDING      * /
     4 52H *     ZONE-FINDING IS RESTARTED WITH OTHER LINES.  ,
     5 18X,1H* / 1H ,70(1H*)  ///)
9903	FORMAT(1H ,I3,4F8.1)
9904	FORMAT(1H ,I3, 3F8.1, I4)
9905	FORMAT(35H0 NR     A       B     FMAAS  HITS            )
9906	FORMAT(1H0,
     *60H POSSIBLE COMBINATIONS OF ZONES (ON WHICH TO BASE A SEARCH  ,
     *15H FOR LATTICES)  /
     *62H COMBINATION A       B       C       D       E       F   NEWNR,
     * '  ZEROSHIFT' )
9907	FORMAT(1H ,2I3,3F8.1,8H     XXX,2F8.1,I6, F8.4)
9908	FORMAT(18H0THIS ZONE INDEXES ,I3,15H LINES FROM THE  ,I3,
     * 22H USED FOR ZONE-FINDING    )
9909	FORMAT(1H1)
9910	FORMAT('1ZONES AFTER REFINEMENT ON TWOTHETA, ONLY THE',
     1  ' ZEROSHIFT IS USED' //
     2  ' NEWNR    A       B     FMAAS       QUALITY   OLDNR CENTER ',
     3  ' NOBS  NCALC  ZEROSHIFT' / )
9911	FORMAT(1H , I4, 3F8.1,2X,F14.1, 4F6.0, F10.4)
9912	FORMAT( ///' THE ZEROSHIFTS ARE NOT YET REALLY USED IN THE',
     * ' PROGRAM.'/ ' IF YOU SEE A CONVINCING ZERO SHIFT, YOU COULD',
     2 ' ENTER THE ZEROSHIFT (WITH THE OPPOSITE SIGN) ' /
     3 ' IN COL. 1-10 OF THE SECOND PARAMETER CARD, AND RERUN THE '
     4 ' PROBLEM'  )
9915	FORMAT('0UPPER BOUNDARIES FOR THE REGIONS OF THE PROBABILITIES',
     * ' GIVEN BELOW' / (1H ,10F12.2)  )
9916	FORMAT(51H0PROBABILITIES OF FINDING A CALCULATED LINE IN THE ,
     2 38H PARTS OF THE Q-RANGE AS SHOWN  ABOVE   / (1H ,10E12.3) )
9917	FORMAT(50H0TEST OUTPUT, MAINLY FROM SUBROUTINE  ** EVALU8 **    )
9920	FORMAT(/42H SEARCH FOR ZONES TERMINATED - STORAGE FOR,I4,14H ZONES
     1 MAXIMUM/ 70H0TO SEARCH FURTHER COMBINATIONS YOU WILL HAVE TO SUBM
     2IT A SEPARATE RUN/60H WITH THE EXTRA LINE PAIRS SUPPLIED USING THE
     3 LINCO FACILITY)
9925	FORMAT(1H0,78(1H*)/' MAXIMUM CAPACITY OF ROSTER IS REACHED'/
     1 ' NUMBER OF COMBINATIONS IS ',I3,
     2 '   THE LAST COMBINATION WAS  ',I2,'-',I2 / 1H0,78(1H*) )
9930	FORMAT(' ZONE ',I2,5X,8F12.1, F8.4)
	END

	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
	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,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/ RAD,TWORAD,RRAD,HLFLAM,WAVSQR,STORE(200,9)
	COMMON /THETAS/ ZERTOT,OBSTET(40),SHFTET(40),TETPT(40),TETMT(40)
	DATA  BRAVTP(1)/'A'/,BRAVTP(2)/'B'/,BRAVTP(3)/'C'/
     *   ,  BRAVTP(4)/'I'/,BRAVTP(5)/'F'/,BRAVTP(6)/'P'/
	IF (ALL) THEN
	   N4=NR
	   GO TO 2000
	END IF
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)
	CALL 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
	CALL TIMER(5)
	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(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
2000	CALL TIMER(6)
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  ---------------------------------------------------------------------
	ALL=.FALSE.
	CALL NTEST(NTST,7)
	DO 2100 J=1,N4
	IF (ROSTER(J,1) .LT. FMINEQ) GO TO 2100
	IF (ROSTER(J,7) .LT. FMINEQ) GO TO 2100
	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  ---------------------------------------------------------------------
C
	CALL REFPAR(NSOL,ALL)
	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
	IF (CEL(1) .LT. FMINIM) GO TO 2080
	U=0.
	DO 2070 I=1,6
2070	U=U+ABS(CEL(I)-OLD(I))
	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	CONTINUE
C-----------------------------------------------------------------------
C  PRINT THE FATE OF THE DISCARDED LATTICES AND
C  COLLECT THE GOOD ONES.
C-----------------------------------------------------------------------
	NR=0
	IF (ENL) WRITE(IOUT,9936)
	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
	CALL TIMER(7)
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(IOUT,9990)
	   RETURN
	END IF
	WRITE(IOUT,9950) NAME
	WRITE(IOUT,9951)
	WRITE(20,9950) NAME
	WRITE(20,9951)
	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(20,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(20,9953) (ROSTER(L,N),N=1,8),BRAVTP(KKK),
     * ROSTER(L,10), VOL,TET
3050	CONTINUE
C-----------------------------------------------------------------------
C--END OF SORT AND PRINT.
C-----------------------------------------------------------------------
	WRITE(IOUT,9963)
	WRITE(20,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
	WRITE(20,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) WRITE(20,9965)
	IF (M .GT. 0) WRITE(20,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)
	WRITE(20,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
	WRITE(20,9949) I,K
	GO TO 3170
3140	DO 3150 J=1,10
3150	CEL(J)= ROSTER(I,J)
	ZERTOT=ROSTER(I,11)
	CALL REFPAR(NSOL,ALL)
	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
C  AFTER THIS REFINEMENT CYCLE, A NEW SORT IS MADE.
C  THE NEW CRITERION IS: (FIGURE OF MERIT)/ (21-LINES INDEXED)
C  ROSTER: 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
	CALL TIMER(8)
	WRITE(IOUT,9954)  NAME,NSOL
	WRITE(20,9954)  NAME,NSOL
	DO 3260 I=1,NSOL
	TET=TWORAD*ROSTER(I,11)
	KKK=NINT( ROSTER(I,9) )
	WRITE(20,9955) ( ROSTER(I,N),N=1,8),  ROSTER(I,10),BRAVTP(KKK)
     1  ,TET
3260	WRITE(IOUT,9955) ( ROSTER(I,N),N=1,8),  ROSTER(I,10),BRAVTP(KKK)
     1  ,TET
	CALL NTEST(NTST,9)
C-----------------------------------------------------------------------
C--PRINT THE VALUES OF THE PARAMETERS OF THE DIRECT LATTICE.
C-----------------------------------------------------------------------
	WRITE(IOUT,9956)
	WRITE(20,9956)
	DO 3280 I=1,NSOL
	DO 3270 J=1,10
3270	CEL(J)= ROSTER(I,J)
	CALL UNTCEL(WMOL,DOBS,VOL,1)
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  I=1,NSOL
	DO 3310 J=1,10
3310	CEL(J)=ROSTER(I,J)
	CALL SYMTST
3320	CONTINUE
	WRITE(IOUT,9957)
	WRITE(IOUT,9959)
	WRITE(IOUT,9958) NAME
	WRITE(20,9957)
	WRITE(20,9959)
	WRITE(20,9958) NAME
	CALL TIMER(9)
C  ---------------------------------------------------------------------
C  THE BEST SOLUTION IS NOW USED TO TRY AND INDEX ALL LINES ENTERED.
C  IF THIS IS TOO MUCH FOR YOU, YOU MAY DELETE THE FOLLOWING CARDS
C  UNTIL  *RETURN*.
C-----------------------------------------------------------------------
	IF (ROSTER(1,7) .LT. 17.9) RETURN
	IF (ROSTER(1,8) .LT. 10.0) RETURN
	NRRUN=2
	ALL=.TRUE.
	DO 3400 I=1,10
3400	CEL(I)=ROSTER(1,I)
	ZERTOT=ROSTER(1,11)
	CALL REFPAR(1,ALL)
	ROSTER(1,11)=ZERTOT
	WRITE(IOUT,9956)
	CALL UNTCEL(WMOL,DOBS,VOL,1)
	CALL TIMER(9)
	RETURN
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
9910	FORMAT(1H ,I5,6F8.1,24X,F8.0,20H DISCARDED BY NORM3  )
9911	FORMAT(1H ,I5,6F8.1,F2.0,22X,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.                       )
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(
     *60H1               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,A80)
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,A80  /// 30X, 4HTHE ,I2,
     1        25H MOST PROBABLE SOLUTIONS //// 58X,
     234HLINES  FIGURE OF           BRAVAIS/
     360H    Q(A)     Q(B)     Q(C)      Q(D)     Q(E)     Q(F)   IND,
     444HEXED   MERIT   PEDIGREE  LATTICE  ZEROSHIFT  )
9955	FORMAT(1H ,6F9.2,F9.0,F9.1,F9.0,8X,A1,F10.4)
9956	FORMAT(1H0,20X,38HTHE DIRECT CONSTANTS OF THESE LATTICES//
     *60H     A        B        C       ALFA     BETA     GAMMA      ,
     * 6HVOLUME)
9957	FORMAT(1H0///40X,' ****   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('0*****************************************************', 
     1'*****************'/2H *,68X,1H*/2H *,'  THIS PROGRAM GIVES A ', 
     2'NUMBER OF SUGGESTIONS THAT MAY HELP YOU TO   *'/ ' *  FIND ', 
     3'THE UNIT CELL THAT EXPLAINS YOUR EXPERIMENTAL DATA. HOWEVER  *' 
     4/2H *,68X,1H*/ 2H *,5X,'WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW', 
     5'WWWWWWWWWWWWWWWWWWWWW',6H     */8H *     W,56X,7HW     */
     6 8H *     W,'   IT IS ENTIRELY THE USERS RESPONSIBILITY TO ', 
     7'DECIDE    W*'/' *     W   WHETHER ANY OF THE SUGGESTED UNIT ',
     8'CELLS IS THE    W     *  '/' *     W   CORRECT CELL.',
     9'                W     * '/ 8H *     W, 56X, 7HW     * / 
     A7H *      , 58(1HW), 6H     * / 2H *,68X,1H*/1H ,70(1H*) )
9960	FORMAT(1H1,A80  //
     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.)
9990	FORMAT(52H0NOT 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     CRITER(K1,L1,CENTER,CRIT,NO,NC)
C  THIS SUBROUTINE CALCULATES THE PROBABILITY THAT A ZONE IS FOUND BY
C  PURE CHANCE. A LINE WITH ITS HIGHER ORDERS IS COUNTED AS ONE
C  OBSERVATION.
C  INPUT PARAMETERS:   K1 IS THE STEPSIZE FOR K
C                      L1 IS THE STEPSIZE FOR L (H)
C                      CENTER OBVIOUSLY GOVERNS THE TEST ON CENTERING
C                      PROB IS THE PROBABILTY OF FINDING A LINE BY CHANC
C  OUTPUT PARAMETER:   CRIT IS THE INVERSE OF THE PROBABILITY THAT THIS
C                      ZONE IS FOUND BY PURE CHANCE
C
C  EVAL IS
C             ZERO ON THE POSITIONS WHERE THE CALCULATED Q WAS TOO LARGE
C
C             NEGATIVE WHEN THE Q-VALUE WAS ONLY CALCULATED BUT NOT
C             OBSERVED,  AND
C
C             POSITIVE WHEN THE CALCULATED Q WAS ALSO OBSERVED.
C
C  THE VALUE OF EVAL IS A SORT OF PROBABILITY THAT A CALCULATED LINE
C  HITS AN OBSERVED LINE BY CHANCE.
C
C
	INTEGER  H,HG
	DIMENSION EVA(-8:8,0:8)
	LOGICAL     CALC,PRES,CENTER,RECTAN
	COMMON /EVALN / NORDER,RECTAN,PROB(40),EVAL(-8:8,0:8),
     1                UPBNDR(40)
	ZERO=1.0E-30
	NCENTR=1
	IF ( CENTER) NCENTR=2*K1
	DO 12 I=-8,8
	DO 10 J=0,8
10	EVA(I,J)=EVAL(I,J)
12	CONTINUE
	KG=NORDER*K1
	HG=NORDER*L1
	NC=0
	NO=0
	CRIT=1.
	LSTEP=L1
	IF (CENTER) LSTEP=NCENTR
	DO 110 K=0,KG, K1
	   LMIN=0
	   IF (CENTER) LMIN=MOD(K,NCENTR)
	   IF (K .EQ. 0) LMIN=LSTEP
	   DO 100 L=LMIN,HG, LSTEP
	   H=L
30	   IF (ABS(EVA(H,K)) .LT. ZERO) GO TO 100
	   PRES=.FALSE.
	   DO 90 N=1,8
	      NH=N*H
	      NK=N*K
	      IF (IABS(NH).GT.HG .OR. NK.GT.KG) GO TO 95
	      IF (PRES) GO TO 85
	      IF (CENTER .AND. (MOD((NH+NK),NCENTR) .NE. 0) ) GO TO 90
	      IF (EVA(NH,NK) .GT. ZERO) THEN
	         CRIT=CRIT*EVA(NH,NK)
	         PRES=.TRUE.
	      END IF
85	      EVA(NH,NK)=0.
90	   CONTINUE
95	   IF (PRES) NO=NO+1
	   NC=NC+1
	   IF (RECTAN) GO TO 100
	   IF (H .LE. 0) GO TO 100
	   H=-L
	   GO TO 30
100	   CONTINUE
110	CONTINUE
	MIN=3
	IF (RECTAN) MIN=2
	NO=NO-MIN
	NC=NC-MIN
C  WE REQUEST  THAT THE NUMBER OF OBSERVED LINES IS AT LEAST
C  25 PERCENT OF THE NUMBER OF CALCULATED LINES.
	IF (4*NO .GE. NC) THEN
	   J=NO
	   IF (2*NO .GT. NC) J=NC-NO
	   IF (J .EQ. 0) GO TO 150
	   DO 120 II=1,J
	      I=II-1
	      CRIT=CRIT*FLOAT(NC-I)/FLOAT(II)
120	   CONTINUE
	ELSE
	   CRIT=100.
	END IF
150	CRIT=1.0/CRIT
	IF (CRIT .GT. 1.0E+10) CRIT=9.999E+9
	RETURN
	END

	SUBROUTINE ECHO
C-----------------------------------------------------------------------
C  THE PURPOSE OF THIS SUBROUTINE IS TO PRINT THE CARD IMAGES OF THE
C  INPUT 'CARDS' ON A SEPARATE PAGE.
C  THIS IS VERY USEFUL FOR TEST PURPOSES AND FOR DOCUMENTATION.
C
C  THE INPUT DATA ARE READ AS CHARACTERS FROM THE 'INP' UNIT AND
C  WRITTEN ON THE *DOCUMENTATION FILE* (4) (IDOC=4).
C  THE OTHER SUBROUTINES READ FROM THE 'INTERNAL FILE' *CARD*
C
C  ECHO ALSO TESTS FOR THE END OF THE DATA. ('END' IN THE FIRST 4 COLS)
C-----------------------------------------------------------------------
	CHARACTER NAME*80,CARD*80,INTSY*2
	COMMON /LINES / IDOC,INP,IOUT,MAX,OBS(40),CEL(10)
	COMMON /CHARCM/ NAME,CARD,INTSY(50)
	DATA NCRD/0/
	NCRD=NCRD+1
C      IF (NCRD .EQ. 1) WRITE(IDOC,9900)
	READ(INP,9901,IOSTAT=IOS)  CARD
	IF (IOS .EQ. 0) THEN
C         WRITE(IDOC ,9902) NCRD, CARD
	   I=INDEX(CARD, 'END' )
	   IF ( (I .EQ. 1) .OR. (I .EQ. 2) ) THEN
	      WRITE(IOUT,9998)
	      STOP
	   END IF
	ELSE  IF (IOS .LT. 0 ) THEN
	   WRITE(IOUT,9995)
	ELSE
	   WRITE(IOUT,9996)
	   STOP
	END IF
	RETURN
9900	FORMAT(48HON THIS PAGE THE CARD IMAGES OF THE INPUT CARDS ,
     1 14H ARE PRINTED.  //  20H  NR     CONTENTS       //
     2 7X,2H 0,8X,2H10,8X,2H20,8X,2H30,8X,2H40,8X,2H50,8X,2H60,
     3 8X,2H70,8X,2H80/ 8X,8(10HI....!....),1HI)
9901	FORMAT(A80)
9902	FORMAT(1H ,I3,5X,A80)
9995	FORMAT('0UNEXPECTED END OF INPUT, PROGRAM TRIES TO CONTINUE')
9996	FORMAT('0ERROR IN INPUT, PROGRAM STOPS.')
9998	FORMAT('0END OF INPUT,  NORMAL END OF PROGRAM')
	END

	SUBROUTINE EVALU8(NR,ISW)
C  THIS SUBROUTINE EVALUATES THE LATTICES FOUND IN FINZON,
C  IT TAKES CARE THAT THEY ARE NORMALISED, FINDS OUT WHETHER ANY OF
C  THE AXES SHOULD BE HALVED AND REFINES THE PARAMETERS.
C-----------------------------------------------------------------------
	INTEGER H
	DIMENSION  HULP(6) ,NOAR(6),NCAR(6)
	LOGICAL RECTAN,REPEAT
	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 /EVALN / NORDER,RECTAN,PROB(40),EVAL(-8:8,0:8),
     *                UPBNDR(40)
	COMMON /CONSTS/ RAD,TWORAD,RRAD,HLFLAM,WAVSQR,STORE(200,9)
	REPEAT=.FALSE.
	A=STORE(NR,1)
	B=STORE(NR,2)
	FMAAS=STORE(NR,3)
	CR=STORE(NR,4)
	CENTR=STORE(NR,6)
20	TYP=.FALSE.
	REPEAT=.NOT. REPEAT
	IF (TEST) TYP=.TRUE.
	IF (ENL .AND. (ISW .EQ. 2) ) TYP=.TRUE.
	IF (TYP )  WRITE(IOUT, 1 ) NR,A,B,FMAAS,TWORAD*ZERSHF
C--WRITES HEADING, INITIAL VALUES.
	CALL NORM(A,B,FMAAS)
	IF (A .LT. FMINIM) GO TO 130
	D=(A+B+FMAAS)*16.5
	IF (D .GT. OBSPT(MAX) ) D=OBSPT(MAX)
C-----------------------------------------------------------------------
	IF (ISW .EQ. 2)  THEN
	
	   GO TO 125
	ELSE IF (ISW .EQ. 3) THEN
	   CALL LSTET(TYP,A,B,FMAAS,ZERSHF,D,NR)
	   STORE(NR,9)=ZERSHF
	   RETURN
	END IF
	CALL LS(TYP,A,B,FMAAS,ZERSHF,D,NR)
	CALL NORM(A,B,FMAAS)
	IF (TEST) WRITE(IOUT,2) A,B,FMAAS,TWORAD*ZERSHF
	IF (A .LT. FMINIM) GO TO 130
	A=0.25*A
	B=0.25*B
	FMAAS=0.25*FMAAS
	CALL  LS(TYP    ,A,B,FMAAS,ZERSHF,D,NR)
C--FROM HERE UNTIL LABEL 50  WE EVALUATE DIFFERENT FORMULATIONS
C--OF THE ZONE.
	NN=0
	DO 40 I=1,2
	DO 30 J=1,2
	H=I
	K=J
	NN=NN+1
	CALL CRITER(H,K,.FALSE.,CRIT,NO,NC)
	NOAR(NN)=NO
	NCAR(NN)=NC
30	HULP(NN)=CRIT
40	CONTINUE
	DO 50 J=1,2
	I=J
	NN=NN+1
	CALL CRITER(I,I,.TRUE.,CRIT,NO,NC)
	NOAR(NN)=NO
	NCAR(NN)=NC
50	HULP(NN)=CRIT
C     CRTMIN=1.2*FLOAT(NOAR(4))/FLOAT(NCAR(4))
C  THE NEW LATTICE SHOULD BE DECIDEDLY  BETTER THAN THE OLD ONE.
	AA=10.0*HULP(4)
	NN=4
	DO 60 N=1,6
C     CRIT=FLOAT(NOAR(N))/FLOAT(NCAR(N))
C     IF (CRIT .LE. CRTMIN) GO TO  60
C     CRTMIN=CRIT
	IF (HULP(N) .LE. AA) GO TO 60
	AA=HULP(N)
	NN=N
60	CONTINUE
	CENTR=0.
	IF (TEST) CALL TEXT1(NN,HULP,NOAR,NCAR)
	GO TO (101,102,103,104,105,106),NN
C--CHOOSE THE BEST.
102	A=4.*A
	FMAAS=2.*FMAAS
	GO TO  101
103	B=4.*B
	FMAAS=2.*FMAAS
	GO TO  101
104	A=4.*A
	B=4.*B
	FMAAS=4.*FMAAS
	GO TO 101
106	CENTR =1.
	GO TO 104
105	CENTR =1.
101	IF (FMAAS .LT. 0.01) GO TO 120
	IF (CENTR .LT. 0.5 ) GO TO 120
C  CENTERED SKEW ZONES ARE REFORMULATED HERE AND RECYCLED.
	E=2.*(B-A)
	A=A+B-FMAAS
	B=A+2.*FMAAS
	FMAAS=E
	IF (TEST) WRITE(IOUT, 3) A,B,FMAAS,NR
	IF (REPEAT) GO TO 20
C  AT THE FIRST PASS, REPEAT IS TRUE
120	CALL  LS(TYP    ,A,B,FMAAS,ZERSHF,D,NR)
	CR= HULP(NN)
	STORE(NR,7)=NOAR(NN)
	STORE(NR,8)=NCAR(NN)
	CALL NORM(A,B,FMAAS)
	IF (TEST) WRITE(IOUT,4) CR
	GO TO 130
125	CALL LS(TYP,A,B,FMAAS,ZERSHF,D,NR)
130	STORE(NR,1)=A
	STORE(NR,2)=B
	STORE(NR,3)=FMAAS
	STORE(NR,4)=CR
	STORE(NR,6)=CENTR
	STORE(NR,9)=ZERSHF
	RETURN
1	FORMAT(10H1 ZONE NR , I4,21H     ORIGINAL  VALUES ,3F10.2,F8.4)
2	FORMAT(21H ZONE REFORMULATED   , 3F10.2,F8.4)
3	FORMAT(43H0A CENTERED SKEW ZONE HAS BEEN REFORMULATED,3F10.2,I4)
4	FORMAT(1H+,90X,12HPROBABILITY  ,E12.3)
	END

	SUBROUTINE FINZON(A,B)
C-----------------------------------------------------------------------
C  SEARCH FOR POSSIBLE ZONES, BASED ON LINES A AND B.
C-----------------------------------------------------------------------
	INTEGER  H
	DIMENSION  BE(10),IBE(10)
	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 /COUNTM/ ITL(0:5000)
C-----------------------------------------------------------------------
	DO 10 I=0,5000
10	ITL(I)=0
	DO 20 I=1,10
	   IBE(I)=0
	   CEL(I)=.0
20	   BE(I)=.0
	F2=6.0/TOL2
	RF2=1.0/F2
	FMAAX=1.95*SQRT(A*B)
	LIMT=NINT(F2*FMAAX)
	IF(LIMT.GT.5000) THEN
	   IF (ENL) WRITE(IOUT,992) LIMT
992	   FORMAT(37H0ARRAY ITL TOO SMALL, SIZE SHOULD BE , I4)
	   LIMT=5000
	END IF
	LASTH=2
	LASTK=2
	MIN=2
	QMAX=OBS(MAX)
	IF (4.*(A+B+FMAAX) .GT. QMAX) THEN
C  RANGE OF H AND K IS NOW NON-STANDARD
	   LASTH=INT(SQRT(QMAX/A))
	   IF (LASTH .GT. 4) LASTH=4
	   LASTK=INT(SQRT(QMAX/B))
	   MIN=(4*LASTH*LASTK)/5
	   IF (LASTH .EQ. 2  .AND.  LASTK .EQ. 2) GO TO 35
	   IF (TEST) WRITE(IOUT,991) A,B,LASTH,LASTK,MIN
	END IF
C  DETERMINATION OF POSSIBLE VALUES OF FMAAS BEGINS HERE.
35	DO 50 H=1,LASTH
	   DO 45 K=1,LASTK
	      Y=A*FLOAT(H*H)+B*FLOAT(K*K)
	      Z=F2/FLOAT(H*K)
	      DO 40 M=1,MAX
	         L=NINT( Z*ABS(OBS(M)-Y)  )
	         IF(L.LE.LIMT) ITL(L)=ITL(L)+1
40	      CONTINUE
45	   CONTINUE
50	CONTINUE
C  ALL POSSIBLE VALUES ARE STORED IN ITL
C     IF (TEST) WRITE(IOUT,994) (ITL(J),J=1,LIMT)
C 994 FORMAT('0ARRAY ITL: '/  (40I3)  )
	LIMT=LIMT-6
	KG=NINT(F2*FMINEQ)
	IF (KG .EQ. 0) KG=1
	DO 55 I=1,KG
55	ITL(0)=ITL(0)+ITL(I)
	ITL(0)=2*ITL(0)
C  THE VALUE ZERO CAN OCCUR ONLY HALF AS MUCH AS NON-ZERO VALUES
	IBE(1)=ITL(0)
	BE(1)=0.
	MIN=1
	NR=1
C  WE NOW TRY TO FIND COINCIDING VALUES (MORE THAN 'MIN' VALUES IN
C  6 SUCCESIVE POSITIONS OF ITL).
C  MANY VALUES (ESPECIALLY THE CORRECT ONES) WILL BE FOUND MORE THAN
C  THAN ONCE. THEREFORE WE MUST STORE AT LEAST 10 VALUES.
C  START OF THE SEARCH IN ITL
	DO 90 I=KG,LIMT,2
	   K=0
	   DO 65 L=I,I+5
65	   K=K+ITL(L)
	   IF (K .LE. MIN) GO TO 90
	   IU=0
	   DO 70 L=I,I+5
70	   IU=IU+ITL(L)*L
C........CALCULATION OF THE WEIGHTED MEAN.
	   DO 80 L=1,NR+1
	      IF (K .GT. IBE(L) ) THEN
	         DO 75 N=NR,L,-1
	            IBE(N+1)=IBE(N)
	            BE(N+1)=BE(N)
75	         CONTINUE
C  ALL 'LOWER' VALUES HAVE NOW BEEN SHIFTED ONE POSITION
	         IBE(L)=K
	         BE(L)=RF2*FLOAT(IU)/FLOAT(K)
	         IF (IBE(9) .GT. MIN) MIN=IBE(9)
	         NR=NR+1
	         IF (NR .GT. 9) NR=9
	         GO TO 90
	      END IF
80	   CONTINUE
90	CONTINUE
C     IF (TEST) WRITE(IOUT,993) BE, IBE
C 993 FORMAT('0ARRAYS BE AND IBE:' / ' ',10F8.1,/ ' ',10I8 )
C  NOW SELECT THE 10 MOST FREQUENT VALUES OF FMAAS AND STORE THEM IN BE.
	II=0
	DO 110 I=1,NR-1
	   FMAAS=BE(I)
	   DO 100 L=I+1,NR
	      IF(ABS(BE(L)-FMAAS).GE.FMINEQ) GO TO 100
C  CHECK ON 'EQUAL' VALUES.  USE THE MEAN.
	      FMAAS=0.5*(FMAAS+BE(L))
	      BE(L)=100000.
100	   CONTINUE
101	   IF(FMAAS.GE.10000.) GO TO 105
	   II=II+1
	   CEL(II)=FMAAS
	   CEL(II+5)=IBE(I)
C  RETURN THE VALUES FOUND VIA CEL
	   IF (II .EQ. 5) RETURN
105	   CONTINUE
110	CONTINUE
	RETURN
991	FORMAT(1H+, 40X,26H NON-STANDARD ZONE SEARCH   ,2F7.1,
     *12H /H/ FROM 1-,I1,10H K FROM 1-,I1,23H MINIMUM NUMBER OF HITS,I3)
	END

	SUBROUTINE INSTRS
C  PRINTS INSTRUCTIONS FOR THE SPECIAL FEATURES.
	COMMON /LINES / IDOC,INP,IOUT,MAX,OBS(40),CEL(10)
	WRITE(IOUT,10)
	WRITE(IOUT,20)
	WRITE(IOUT,30)
	WRITE(IOUT,40)
	WRITE(IOUT,50)
	WRITE(IOUT,60)
	WRITE(IOUT,70)
	WRITE(IOUT,80)
	WRITE(IOUT,90)
	WRITE(IOUT,100)
	WRITE(IOUT,110)
	WRITE(IOUT,120)
	WRITE(IOUT,130)
	WRITE(IOUT,135)
	WRITE(IOUT,140)
	WRITE(IOUT,150)
	WRITE(IOUT,200)
	RETURN
10	FORMAT(35H1 INSTRUCTIONS FOR SPECIAL FEATURES   ,
     1       ' ****  VERSION  11  ***** ' //
     2       1H0,10X,42HCOL  3  **INTENS**                        /
     3       50H  IF YOU WANT TO ENTER INTENSITIES WITH YOUR LINE ,
     4       50HPOSITIONS YOU SHOULD SET INTENS(COL 3) AT 1.      /
     5       50H THE LINE POSITIONS AND INTENSITIES ARE THEN READ ,
     6       50HIN THE FORMAT 8(F7.3,A2,1X),WHICH MEANS THAT THE  /
     7       50H INTENSITIES ARE READ AS CHARACTERS. YOU MAY ENTER,
     8       50H THEM AS A 2-DIGIT FIGURE BUT ALSO AS   VS, W , VW/
     9       50H ETC. YOU COULD ALSO USE THIS TO MARK ONLY THE VER,
     A       50HY WEAK LINES OR SUSPECTED LINES.                  /
     B       50H NO CALCULATIONS ARE CARRIED OUT WITH THE INTENSIT,
     C        4HIES. )
20	FORMAT(1H0,10X,42HCOL  4  **NSOLMX**                        /
     *       50H NSOLMX IS THE MAXIMUM NUMBER OF SOLUTIONS FOR WHI,
     1       50HCH A LIST OF (INDEXED) REFLECTIONS IS PRINTED.    /
     2       50H THE CORRESPONDING Q-SCHEME IS ONLY PRINTED WHEN  ,
     3       50HINTERMEDIATE RESULTS ARE GIVEN (COL. 34)  .AND.   /
     4       50H THE FIGURE OF MERIT IS LARGER THAN 4.0.          )
30	FORMAT(1H0,10X,42HCOL 5-6, 7-8, 9-10   **NSYST**            /
     *       50H WHEN YOU ARE REASONABLY SURE THAT THE SOLUTION MU,
     1       50HST BELONG TO A CERTAIN CRYSTAL SYSTEM, YOU CAN    /
     2       50H INDICATE THIS IN COLS 5-10 (3I2). THE COLUMNS SER,
     3       50HVE THE ORTHORHOMBIC, MONOCLINIC AND TRICLINIC     /
     4       50H SYSTEM RESPECTIVELY.  A ZERO MEANS THAT YOU ARE  ,
     5       50HINDIFFERENT (DEFAULT CASE), A +1 INDICATES THAT   /
     6       50H YOU ARE PRETTY SURE THAT THE CORRECT SOLUTION BEL,
     7       53HONGS TO THIS CRYSTAL SYSTEM, A -1 MEANS THE OPPOSITE./
     8       50H THIS OPTION WORKS ONLY ON THE LAST PART OF THE PR,
     9       50HOGRAM AND DOUBLES THE TOLERANCES THERE..          /
     O       50H THIS IS ONLY POSSIBLE FOR ORTHORHOMBIC,MONOCLINIC,
     A       50H AND TRICLINIC LATTICES.                           )
40	FORMAT(1H0,10X,22H COL 11-15   **TOL2**    /
     *           11X,22H COL 16-20   **TOL3**    /
     *       50H THE TOLERANCE ON THE VALUE OF R (THE TRIAL VALUE ,
     1       50HOF F USED IN THE SEARCH FOR ZONES, SEE MY PAPER,  /
     2       50H EQ. 4) IS GIVEN BY *TOL2*. WHEN YOU HAVE VERY GOO,
     3       50HD DATA (0.01 DEGREE TWOTHETA) ON A FAIRLY LARGE   /
     4       50H UNIT CELL,YOU MIGHT LOWER IT TO 2.5 OR 2.0 .WITH ,
     5       50HA SMALL UNIT CELL YOU MIGHT INCREASE IT SOMEWHAT  /
     6       50H TO SAY 4.0. TOL3 IS THE SAME KIND OF TOLERANCE, B,
     7       50HUT NOW ON THE TRIAL VALUE OF D IN THE SEARCH FOR  /
     8       50H LATTICES (SEE PAPER, PAGE 91).ITS VALUE CAN BE CH,
     9       50HANGED IN THE SAME WAY AS THAT OF TOL2.            )
50	FORMAT(1H0,10X,42HCOL 31-31  **LINCO**                      /
     *       50H LINCO ENABLES YOU TO ENTER YOUR OWN LINE COMBINAT,
     1       50HIONS OR YOUR OWN ZONES.                           /
     2       44H WHEN LINCO IS NOT ZERO, THE PROGRAM EXPECTS      ,
     3       50H CARDS (AFTER THE LINES CARDS+THE BLANK CARD)     /
     4       50H AND READS THESE IN THE FORMAT  3F10.2,F10.1      /
     4       50H YOU MAY ENTER ANY NUMBER OF ZONES (LESS THAN 40) ,
     5       50HBUT YOU MUST END WITH A BLANK CARD.               /
     5       50H WHEN YOU ENTER ONLY 2 VALUES, THE PROGRAM INTERPR,
     6       50HETS THIS AS A COMBINATION OF LINES, TO BE TRIED   /
     7       50H FOR FINDING ZONES.                               )
60	FORMAT(50H WHEN YOU ENTER 3 VALUES, THE PROGRAM WILL TREAT  ,
     1       50HTHESE AS  A ZONE AND WILL EVALUATE IT. IF YOU WANT/
     2       50H TO ENTER A RECTANGULAR ZONE, YOU MUST SET ITS THI,
     3       50HRD CONSTANT AT 0.08 , OTHERWISE THE PROGRAM WILL  /
     4       50H TREAT IT AS A LINE COMBINATION.                  )
70	FORMAT(50H IF YOU WANT TO ENSURE THAT A SPECIFIED ZONE WILL ,
     1       50HBE INCLUDED IN THE SEARCH FOR COMPLETE LATTICES,  /
     2       50H YOU SHOULD ALSO ENTER A HIGH QUALITY FIGURE (FOUR,
     3       50HTH VALUE ON THE CARD), PUNCHING A SINGLE DIGIT IN /
     4       50H COL. 31 ASSURES YOU OF A HIGH QUALITY FIGURE.    /
     5       50H IF THE PROGRAM DOES NOT FIND THE SOLUTION AUTOMAT,
     6       50HICALLY, IT CAN OFTEN BE HELPED BY ADDING ZONES IN /
     7       50H THIS WAY. ZONES WHICH ARE CORRECT, ALTHOUGH HAVIN,
     8       51HG A LOW QUALITY FIGURE, CAN SOMETIMES BE RECOGNISED/
     9       50H IN THE LIST OF INTERMEDIATE RESULTS (K, COL.34)  ,
     *       50HBY THEIR SMALL AREA (I.E. LARGE *A* AND *B*) AND  /
     A       50H GOOD COVERAGE NEAR THE ORIGIN (I.E. MOST OF THE  ,
     B       50HLOW ANGLE CALCULATED LINES HAVE BEEN MATCHED WITH /
     C       50H OBSERVED LINES). ZONES FROM ORIENTED SPECIMENS OR,
     D       50H FIBRES, OR FROM ELECTRON DIFFRACTION STUDIES, CAN/
     E       20H ALSO BE USED.         )
80	FORMAT(1H0,10X,42HCOL 33   ***LZERCHK***                    /
     1       50H  YOU CAN ALSO ASK THE PROGRAM TO PERFORM A CHECK ,
     2       50HON ZERO ERRORS (COL. 33). THE SUBROUTINE GIVES    /
     3       20H FULL PARTICULARS.   )
90	FORMAT(1H1)
100	FORMAT(1H0,10X,42HCOL 43-44   **NR**                        /
     *       50H    WHEN YOU WANT TO ENTER A NUMBER OF TRIAL LATTI,
     1       50HCES OF YOUR OWN, YOU CAN DO THIS BY MAKING NR POS./
     2       50H THE PROGRAM WILL NOW EXPECT CARDS WITH TRIAL LATT,
     4       50HICES AFTER THE INPUT OF THE LINES (+BLANK CARD).  /
     5       50H THERE SHOULD BE 6 DATA PER CARD (FORMAT 6F10.2), ,
     7       50HTHE COEFFICIENTS A,B,C,D,E AND F OF               /
     8       50H  Q=H*H*A+K*K*B+L*L*C+K*L*D+L*H*E+H*K*F           /
     9       50H WHERE Q=10000/(D*D)  (D=INTERPLANAR SPACING IN A)/
     O       50H0*****YOU MAY ALSO ENTER THE TRIAL LATTICE AS A DI,
     A       50HRECT LATTICE WITH 6 CONSTANTS IN ANGSTROMS AND (DE,
     B       15HCIMAL) DEGREES. /
     C       50H THE TRIAL LATTICES ARE TERMINATED BY A CARD WITH ,
     D       50HZERO OR BLANK IN COLS 1-10.                       )
110	FORMAT(1H0,10X,42HCOL 47-48  **LIST**                       /
     *       50H LIST. FOR THE BEST SOLUTION, THE PROGRAM WILL PRO,
     1       50HDUCE A LIST OF ALL OBSERVED AND CALCULATED LINES, /
     2       50H SORTED ON TWOTHETA, UNLESS IT WOULD BE TOO LONG (,
     3       50HMORE THAN 100 CALCULATED LINES). IF YOU WANT THIS /
     4       50H LIST FOR MORE SOLUTIONS, YOU SHOULD ENTER THE NUM,
     5       52HBER THAT YOU WANT IN COL. 48.  YOU CAN ALSO SUPPRESS/
     6       50H THIS OPTION BY MAKING LIST NEGATIVE. DEFAULT=1   )
120	FORMAT(1H0,10X,42HCOL 71-78   **TOLG**                      /
     *       50H TOLG IS THE TOLERANCE ON THE FIT BETWEEN OBSERVED,
     1       52H AND CALCULATED LINES. IT IS EXPRESSED IN HUNDREDTHS/
     2       52H OF A DEGREE TWOTHETA, E.G. THE DEFAULT VALUE OF 6.0,
     3       50H MEANS THAT A DIFFERENCE OF 0.060 DEGREES TWOTHETA/
     4       50H BETWEEN CALCULATED AND OBSERVED LINES IS JUST ACC,
     5       10HEPTABLE.  /
     6       50H TOLG APPLIES TO BOTH THE ZONE AND LATTICE SEARCHE,
     6       50HS, AND TO THE LATTICE REFINEMENT.                 /
     7       47H FOR VERY ACCURATE DATA IT CAN BE ADVANTAGEOUS ,
     8       47HTO LOWER TOLG TO 4.0 OR 3.5.                   /
     9       50H FOR INACCURATE DATA IT CAN BE INCREASED, BUT SUCH,
     *       50H DATA IS NOT USUALLY WORTH PROCESSING UNLESS THE  /
     A       50H CELL IS VERY SMALL.                              )
130	FORMAT(1H0,10X,42HCOL  79-80  *NTST**                       /
     *       50H THE TEST OUTPUT CAN BE LIMITED TO THE SUBROUTINE ,
     1       50H WHERE YOU REALLY NEED IT, BY ENTERING A PARAMETER/
     2       50H FOR *NTEST* IN COLUMN 80 OF THE FIRST PARAMETER  ,
     3       50HCARD. (SEE SUBROUTINE FOR DIRECTIONS)             /
     4       50H0IF YOU WANT THE TEST OUTPUT TO CONTINUE UNTIL THE,
     5       50H END OF THE PROGRAM, YOU MUST MULTIPLY THE DESIRED/
     6       50H NUMBER BY 11 AND USE COLUMNS 79 AND 80 OF THE FIR,
     7       50HST PARAMETER CARD. PUNCH E.G. 88 IF THE LAST PART /
     8       50H OF THE PROGRAM GIVES YOU TROUBLE.                //)
135	FORMAT(1H0,30X,20(1H*) /31X,20H*  ON SECOND CARD  * /
     1       31X, 20(1H*)   )
140	FORMAT(1H0,10X,42HCOL 11-20  **PRNTMR**                     /
     *       50H **PRNTMR**THE PROGRAM DOES NOT PRINT ANY LATTICES,
     1       50H THAT HAVE A FIGURE OF MERIT THAT IS BELOW PRNTMR./
     2       51H THE DEFAULT VALUE IS 4.0, WHICH MUST BE CONSIDERED,
     3       50H AS AN ABSOLUTE MINIMUM FOR A LATTICE THAT HAS A  /
     4       50H CHANCE TO IMPROVE ITS PERFORMANCE.               )
150	FORMAT(1H0,10X,42HCOL 21-30  **PRNTLN**                     /
     *       50H **PRNTLN**THE PROGRAM DOES NOT PRINT ANY LATTICES,
     1       50H THAT HAVE LESS THAN *PRNTLN* INDEXED LINES.      /
     2       50H LATTICES WITH SUCH A LOW NUMBER OF INDEXED LINES ,
     3       50HARE SOMETIMES A SUBLATTICE OF THE TRUE LATTICE.   /
     4       50H A LATTICE CAN ONLY BE CONSIDERED AS SATISFACTORY ,
     5       50HWHEN IT EXPLAINS ALL LINES.                       )
200	FORMAT(///47H0**NOTE**  A USEFUL STRATEGY IS TO IMPROVE THE ,
     1       50HOBSERVED DATA BY USING LINES FOR WHICH 2 OR MORE  /
     2       50H ORDERS HAVE BEEN OBSERVED. THE LOWEST ORDERS CAN ,
     3       50HBE REPLACED BY A MORE ACCURATE VALUE CALCULATED   /
     4       50H FROM THE HIGHER ORDERS, AND SOMETIMES UNOBSERVED ,
     5       50HFIRST OR SECOND ORDER LINES CAN BE INFERRED AND   /
     6       50H ADDED TO THE  *OBSERVED* PATTERN.  LOW ORDER ABSE,
     7       50HNCES IN WELL-ESTABLISHED ZONES CAN BE FILLED IN A /
     8       50H SIMILAR WAY. MULTIPLE ORDERS CAN ALSO BE USED TO ,
     9       50HCHECK FOR ZERO ERRORS IN TWOTHETA.                /
     O       50H SEE ALSO SUBROUTINE ZERCHK.                      )
	END

	SUBROUTINE  LISTM(NMATCH)
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.
	DIMENSION NMATCH(45)
	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/ RAD,TWORAD,RRAD,HLFLAM,WAVSQR,STORE(200,9)
	COMMON /THETAS/ ZERTOT,OBSTET(40),SHFTET(40),TETPT(40),TETMT(40)
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)
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)
	WRITE(IOUT,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) WRITE(IOUT,9905)
	WRITE(IOUT,9902) U,D, Q ,INTSY(NO)
	IF (NMATCH(NO) .EQ. 0) WRITE(IOUT,9906)
	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)
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+,50X,12H XXXXXXXXXXX )
	END

	SUBROUTINE   LS(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,6),HULP(6),FYN(3)
	INTEGER  H,HG,HMIN, ORDE,ORP1
	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)
	DO 20 H=1,4
	DO 10 K=1,6
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
	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	      CONTINUE
	      DO 59 J=1,MAX
	         IF (OBSPT(J) .LT. Q) GO TO 59
	         IF (OBSMT(J) .GT. Q) 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)
	         HULP(1)=H*H
	         HULP(2)=K*K
	         HULP(3)=H*K
	         HULP(ORP1)=OBS(J)-Q
	         DO 58 N2=1,ORP1
	            DO  57 N1=1,ORDE
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, Q,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-----------------------------------------------------------------------
	IF (RECTAN) THEN
	   Z=REF(1,1)*REF(2,2)-REF(2,1)*REF(1,2)
	   IF (ABS(Z) .LT. 0.1 ) THEN
	      IF (TEST) WRITE(IOUT,3) NR
	      RETURN
	   END IF
	   FYN(1)=+(REF(1,3)*REF(2,2)-REF(2,3)*REF(1,2))/Z
	   FYN(2)=-(REF(1,1)*REF(2,3)-REF(2,1)*REF(1,3))/Z
	   GO TO 140
	END IF
	DO 80 H=1,2
	   DO 75 K=1,3
75	   REF(K,H+4)=REF(K,H)
80	CONTINUE
	DO 90 H=1,4
	   K=H+1
	   L=H+2
	   REF(4,H)=REF(1,H)*(REF(2,K)*REF(3,L)-REF(3,K)*REF(2,L))
     1           -REF(2,H)*(REF(1,K)*REF(3,L)-REF(3,K)*REF(1,L))
     2           +REF(3,H)*(REF(1,K)*REF(2,L)-REF(2,K)*REF(1,L))
90	CONTINUE
	IF (ABS(REF(4,1)) .LT. 1.0) THEN
	   IF (TYP) WRITE(IOUT, 3) NR
	   RETURN
	END IF
	DO 100 J=1,3
100	FYN(J)=REF(4,J+1)/REF(4,1)
	FMAAS=FMAAS+FYN(3)
140	A=A+FYN(1)
	B=B-FYN(2)
	IF (TYP) WRITE(IOUT,4) NR,A,B,FMAAS
	IF (TEST)WRITE(IOUT,5) NINDEX,NQ2
	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 )
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   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),HELP1(4),HELP2(4)
	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/ 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)
	CALL MINV(REF,ORP1,DETER,HELP1,HELP2)
	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  MANUAL
C  THIS SUBROUTINE PRINTS AN INSTRUCTION MANUAL HOW TO ARRANGE
C  THE INPUT PARAMETERS.
	COMMON /LINES / IDOC,INP,IOUT,MAX,OBS(40),CEL(10)
	WRITE(IOUT,9950)
	WRITE(IOUT,9952)
	WRITE(IOUT,9954)
	WRITE(IOUT,9956)
	WRITE(IOUT,9958)
	WRITE(IOUT,9960)
	WRITE(IOUT,9962)
	WRITE(IOUT,9964)
	WRITE(IOUT,9966)
	WRITE(IOUT,9968)
	WRITE(IOUT,9970)
	WRITE(IOUT,9972)
	WRITE(IOUT,9974)
	WRITE(IOUT,9976)
	WRITE(IOUT,9978)
C     WRITE(IOUT,9980)
	RETURN
9950	FORMAT( 1H1,10X,29H***  INSTRUCTION  MANUAL  ***, 20X,
     * 20H---  VERSION 11  --- //
     2 35H IN CASE OF PROBLEMS, WRITE OR CALL /
     3114H     J.W.VISSER,  TECHNISCH PHYSISCHE DIENST, P.O.BOX 155,
     4DELFT, NETHERLANDS.   PHONE  015 - 569300, EXT. 249 / 40X,2HOR /
     5 115H     G.G.JOHNSON JR., MATERIALS RESEARCH LAB., PENN. STATE
     6UNIVERSITY, UNIVERSITY PARK, PENN. 16802, 814-865-1637/40X,2HOR/
     7 120H     R.SHIRLEY, CHEMICAL PHYSICS DEPT., UNIV. OF SURREY, GUIL
     8DFORD, SURREY GU2 5XH, ENGLAND.   PHONE 0483-71281 X449/427)
9952	FORMAT(50H0THIS PROGRAM TRIES TO FIND THE UNIT CELL FROM THE,
     1       50H POWDER PATTERN.  THE PRINCIPLES ARE DESCRIBED IN /
     2       50H A PAPER BY VISSER IN THE JOURNAL OF APPLIED CRYST,
     3       50HALLOGRAPHY (1969) 2, 89.                          //
     4       50H MOST OF THE CALCULATIONS ARE CARRIED OUT IN TERMS,
     5       50H OF Q-VALUES, THESE ARE (100/D)**2 ,WHERE D IS THE/
     6       50H INTERPLANAR DISTANCE IN ANGSTROMS.  ALL VALUES IN,
     7       50H THE PROGRAM ARE EXPRESSED IN QU (Q-UNITS), UNLESS/
     8       30H INDICATED OTHERWISE.            / )
9954	FORMAT(50H THE INPUT OF THE PROGRAM IS ON UNIT 5, THE REGULA,
     A       50HR OUTPUT IS ON UNIT 6. THE CARD-IMAGES OF THE     /
     B       50H INPUT CARDS ARE PRINTED ON UNIT 4.               / )
9956	FORMAT(50H THE USUAL INPUT OF THIS PROGRAM SHOULD CONSIST OF/
     A       50H 1.A TITLE CARD, CONTAINING UP TO 80 CHARACTERS   /
     B       50H 2.A PARAMETER CARD, WHICH MAY BE COMPLETELY BLANK/
     X       50H 3.ANOTHER PARAMETER CARD. CAN ALSO BE BLANK.      /
     C       50H 4.A NUMBER OF LINES CARDS.                       /
     D       50H   THE OPTIMUM NUMBER OF LINES TO BE ENTERED IS   ,
     *       50HSOMEWHERE BETWEEN 30 AND 35.                      /
     1       50H   THE ABSOLUTE MINIMUM IS 20. THERE IS NO MAXIMUM,
     2       50H, BUT THE PROGRAM WILL USE ONLY THE FIRST 40 LINES/
     3       50H   THE LINES SHOULD BE ENTERED IN ASCENDING ORDER ,
     4       50HOF TWOTHETA ( =DESCENDING ORDER OF D).            /
     E       50H   THEY CAN BE ENTERED AS 2THETA, Q- OR D- VALUES ,
     F       50H IN THE FORMAT  8F10.5                              )
9958	FORMAT(50H   YOU MAY LEAVE ANY NUMBER OF POSITIONS ON THE CA,
     G       50HRD BLANK. THIS MEANS THAT YOU CAN ALSO ENTER ONLY /
     I       50H   ONE LINE PER CARD,IN ANY OF THE POSITIONS OF 8F,
     J       50H10.5  MOREOVER, ONLY THE FIRST SIX LINES MUST BE  /
     K       50H   IN THEIR CORRECT ORDER. AFTER DETERMINING WHETH,
     L       50HER THESE ARE IN ASCENDING(2THETA,Q) OR DESCENDING /
     M       50H   (D) ORDER, THE PROGRAM WILL SORT ALL LINES.      )
9960	FORMAT(50H   IF YOU WANT TO DROP A QUESTIONABLE INPUT LINE, ,
     1       50H CHANGE ONLY ITS SIGN (PUNCH A MINUS).            /
     2       50H   THE VALUE WILL REMAIN READABLE, BOTH ON THE CAR,
     3       50HD AND IN THE TEST OUTPUT.                         /
     4       50H 5.A BLANK CARD (OR A ZERO IN COL. 1 IF YOU WISH  ,
     4       50HTO BE MORE EXPLICIT) TO FINISH THE INPUT OF LINES /
     4       50H   OF ONE PROBLEM.                                )
9962	FORMAT(50H0  REPEAT THE POINTS 1-4 AS MANY TIMES AS PROBLEMS,
     1       50H THEN                                             /
     2       50H 6.A CARD CONTAINING THE WORD -- END -- IN THE FIR,
     3       50HST THREE COLUMNS, TO FINISH THE INPUT.            /  )
9964	FORMAT(1H1)
9966	FORMAT(50H THE PARAMETER CARD CAN CONTAIN THE FOLLOWING PARA,
     1       50HMETERS--                                          //
     2       50H COLUMN  FORMAT  NAME   DEFAULT   MEANING         /
     3       50H   1      A1     MAN      0       PRINTS THESE INS,
     4       50HTRUCTIONS, UNLESS MAN=9                           )
9968	FORMAT(50H   2      A1     INSTR    0       PRINTS INSTRUCTI,
     1       50HONS FOR SPECIAL FEATURES IF INSTR=1               /
     2       50H   3      A1     INTENS   0       READS INTENSITIE,
     3       50HS BESIDES LINE POSITIONS (8(F7.3,A3)) IF INTENS=1 /
     A       50H   4      I1     NSOLMX   4       NUMBER OF SOLUTI,
     B       50HONS PRINTED OUT EXTENSIVELY.                      /
     4       50H  5- 6    I2     NSYST(1) 0       ORTHORHOMBIC SYS,
     5       50HTEM,  +1=YES, -1=NO, 0=INDIFFERENT.               /
     6       50H  7- 8    I2     NSYST(2) 0       MONOCLINIC   SYS,
     7       50HTEM,  +1=YES, -1=NO, 0=INDIFFERENT.               /
     8       50H  9-10    I2     NSYST(3) 0       TRICLINIC    SYS,
     9       50HTEM,  +1=YES, -1=NO, 0=INDIFFERENT.               )
9970	FORMAT(50H 11-15   F5.2    TOL2    3.0 (QU) TOLERANCE(RANGE),
     6       50H ON 2-DIMENSIONAL SEARCH.                         /
     7       50H 16-20   F5.2    TOL3    4.5 (QU) TOLERANCE(RANGE),
     8       50H ON 3-DIMENSIONAL SEARCH.                         /
     9       50H 21-30   F10.5   WAVEL   1.5406(A)WAVELENGTH      /
     O       50H 31-32    I2     LINCO    0       IF LINCO .GT. 0,,
     O       51H ENTER LINE COMBINATIONS OR ZONES(4F10.2), END WITH/
     A 35X,  55H A BLANK CARD,OR A CARD WITH A ZERO IN THE FIRST COLUMN/
     B       50H    33    I1     LZERCK   0       GIVES A CHECK ON,
     C       50H ZERO-ERRORS WHEN LZERCK > 0                      /
     D       50H    34    I1     K (ENL)  0       GIVES OUTPUT OF ,
     E       50HINTERMEDIATE RESULTS, IF K .GE. 1                   )
9972	FORMAT(50H 35-36    I2     NQ1      3       NQ1 AND NQ2 DETE,
     1       50HRMINE THE NUMBERS OF LINES TO BE USED IN          /
     2       50H 37-38    I2     NQ2      6           COMBINATIONS,
     3       50H FOR FINDING ZONES.                               /
     4       50H 39-40    I2     NZ1      6       NZ1 AND NZ2 DO T,
     5       50HHE SAME FOR FINDING LATTICES FROM ZONES           /
     6       50H 41-42    I2     NZ2      6                       /
     7       50H 43-44    I2     NR       0       REFINEMENT AND E,
     8       52HVALUATION ONLY IF NR .GT. 0. TRIAL LATTICES REQUIRED/
     9       50H 45-46    I2     INDAT    5       NUMBER OF THE UN,
     O       50HIT ON WHICH TO READ THE LINE POSITIONS            )
9974	FORMAT(50H 47-48    I2     LIST     1       PRINTS A LIST OF,
     *       50H ALL CALC AND OBS LINES FOR *LIST* LATTICES       /
     A       50H 49-50    I2     J(TEST)  0       TEST OUTPUT IS  ,
     B       50HPRINTED WHEN  J.GE.1                              /
     C       50H 51-60   F10.5   WMOL    0.0      MOLECULAR WEIGHT/
     D       50H 61-70   F10.5   DOBS    0.0      OBSERVED DENSITY/
     E       50H 71-78    F8.4   TOLG    6.0      TOLERANCE ON THE,
     F       50H MATCH BETWEEN CALCULATED AND OBSERVED TWO THETAS./
     G       51H 79-80    I2     NTST     0       NUMBER,INDICATING,
     H       53H WHERE IN THE MAIN PROGRAM,  TEST OUTPUT SHOULD START)
9976	FORMAT(50H0ON SECOND CARD                                   /
     1       50H0 1-10   F10.4   ZERCOR  0.0      ZERO CORRECTION ,
     2       50H(DEGREES TWOTHETA), TO BE APPLIED TO ALL LINES    /
     3       50H 11-20   F10.1   PRNTMR  4.0      THE MINIMUM FIGU,
     4       50HRE OF MERIT FOR A LATTICE TO BE PRINTED           )
9978	FORMAT(50H 21-30   F10.0   PRNTLN  14.      THE MINIMUM NUMB,
     1       50HER OF INDEXED LINES FOR A LATTICE TO BE PRINTED.  /
     2       50H 31-40   F10.0   RMAX     0.      THE NUMBER OF LI,
     3       52HNES ENTERED. ONLY NECESSARY IF RMAX IS LESS THAN 33 )
	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 NORM(A,B,FMAAS)
C-----------------------------------------------------------------------
C  NORM PUTS ZONES INTO THEIR  STANDARD (NORMALISED) FORM..
C  IT LOOKS FOR SYMMETRY AND MAKES  SKEW ZONES  AS NEAR TO RECTANGULAR
C  AS POSSIBLE.    A AND B ARE INTERCHANGED IF  A>B.
C-----------------------------------------------------------------------
	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
	IF (TEST) WRITE(IOUT,1) A,B,FMAAS
10	IF (A .LT. FMINIM) GO TO 65
	IF ( B-A .GE.0.0) GO TO 20
	AA=A
	A=B
	B=AA
20	IF (ABS(B-FMAAS).GE.FMINEQ) GO TO 30
	AA=A
	A=0.25*B
	B=AA-A
	FMAAS=0.0
	GO TO 10
30	IF  (   FMAAS.LE. FMINEQ) GO TO 40
	IF  (ABS(A-B).GT. FMINEQ) GO TO 40
	A=0.25*(A+B-FMAAS)
	B=A+0.5*FMAAS
	FMAAS=0.0
	GO TO 10
40	IF (ABS(A-FMAAS).GT.FMINEQ) GO TO 50
	A=0.25*A
	B=B-A
	FMAAS=0.0
	GO TO 10
50	IF ((A-FMAAS).GT. 0.0) GO TO 60
	B=A+B-FMAAS
	FMAAS=ABS(FMAAS-2.*A)
	GO TO 10
60	IF (FMAAS .LT. FMINEQ) FMAAS=0.0
65	IF (TEST) WRITE(IOUT,2) A,B,FMAAS
	RETURN
1	FORMAT(15H INPUT OF NORM  , 3F10.2 )
2	FORMAT(1H+, 50X, 8H OUTPUT  , 3F10.2 )
	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=10000.0
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 NTEST(NTST,I)
C  NTEST CAN LIMIT THE TEST OUTPUT TO SMALL PARTS OF THE PROGRAM.
C  IT IS PROBABLY ONLY USEFUL WHEN YOU ARE MODIFYING THE PROGRAM.
C  THE WHOLE PROGRAM IS SUBDIVIDED IN 10 PARTS, EACH WITH ITS OWN
C  TEST OUTPUT.
C
C                                  CALL IS BETWEEN
C    NTST    TEST OUTPUT FOR       LABELS        IN SUBROUTINE
C      1     READER+BEGIN MAIN       11 -  110      READER
C      2     FINZON + EVALU8        450 -  460      AAITO1
C      3     EVALU8 + ONEZON        550 -  670      AAITO1
C      4     REAR                   570 - 1805      AAITO1
C      5     THREED                1860 - 1905      ANITO2
C      6     NORM3                 1910 - 1920      ANITO2
C      7     REFPAR+UNTCEL         2055 - 2060      ANITO2
C      8     REFPAR (2)            3050 - 3052      ANITO2
C      9     UNTCEL + SYMTST       3110 - 3115      ANITO2
C
C  IF YOU WANT YOUR TEST OUTPUT TO CONTINUE AFTER IT HAS STARTED,
C  YOU SHOULD MULTIPLY THE NUMBER GIVEN ABOVE WITH 11,E.G. YOU
C  ENTER 88 IF YOU SUSPECT THAT THE TROUBLE IS IN THE LAST PART OF
C  THE PROGRAM.
C
C              *****************************
C  N.B.        *  NTEST  OVERRULES   TEST  *
C              *****************************
C
C
C  IF YOU HAD SET  *TEST* (COL 50 IN FIRST PARAMETER CARD),THAN THIS
C  PARAMETER IS OVERRULED WHEN YOU SET *NTEST*. THE PHILOSOPHY IS THAT
C  WHEN YOU USE  *NTEST* YOU KNOW EXACTLY WHAT TEST OUTPUT YOU WANT.
C  BY SETTING *TEST* AND NOT USING *NTEST* YOU GET ALL TEST OUTPUT.
C
	LOGICAL TEMP,SELF,ENL,TEST
	CHARACTER*8 TEXT(9)
	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
	DATA TEXT/' READER ',' FINZON ',' EVALU8 ','  REAR  ',' THREED ',
     2  '  NORM3 ',' REFPAR ',' REFPAR2',' UNTCEL '/
	IF (I .EQ. 1) TEMP=ENL
C--KEEPS *ENL* FOR THE REST OF THE PROGRAM.
	TEST=.FALSE.
	ENL=TEMP
	IF (NTST .EQ. I) GO TO 50
	IF (NTST .EQ. 99) GO TO 50
	IF (NTST .LT. 20) RETURN
	J=NTST/11
	IF (I .GE. J) GO TO 50
C  THIS MAKES IT POSSIBLE TO OBTAIN ALL TEST OUTPUT AFTER
C  A CERTAIN POINT.
	RETURN
50	TEST=.TRUE.
	ENL=.TRUE.
	WRITE(IOUT, 1) I,TEXT(I)
C
C     SAVE TEMP
C
	RETURN
1	FORMAT(15H0TEST OUTPUT **,I3,32H **   MAINLY FROM SUBROUTINE  **,
     2  A8, 2H** )
	END

	SUBROUTINE ONEZON(A,B,FMAAS)
C-----------------------------------------------------------------------
C  WHEN TOO MANY LINES AMONG THE FIRST SIX BELONG TO ONE ZONE,
C  THIS SUBROUTINE PRODUCES COMBINATIONS OF ONE LINE OF THE ZONE
C  WITH   A  LINE THAT DOES NOT BELONG TO THE ZONE.
C  THESE COMBINATIONS ARE USED TO ENTER  FINZON  AGAIN.
C-----------------------------------------------------------------------
	DIMENSION ZON(4)
	COMMON /LINES / IDOC,INP,IOUT,MAX,OBS(40),CEL(10)
	COMMON /TEMRES/ LINCO,NQ2,NINDEX,INDEX(40),ROSTER(150,12)
C     WRITE(IOUT, 1) MAX,LINCO,(INDEX(I),I=1,40)
	ZON(1)=A
	ZON(2)=B
	ZON(3)=A+B+FMAAS
	ZON(4)=A+B-FMAAS
	NEND=4
	IF (ZON(3)-ZON(4) .LT. 0.01) NEND=3
	DO 50 I=1,NEND
	N=0
	DO 40 J=1,MAX
	IF (INDEX(J) .EQ. 1) GO TO 40
	N=N+1
	LINCO=LINCO+1
	ROSTER(LINCO,1)=ZON(I)
	ROSTER(LINCO,2)=OBS(J)
	ROSTER(LINCO,3)=.0
	ROSTER(LINCO,4)=.0
C     WRITE(IOUT, 2) LINCO,(ROSTER(LINCO,K),K=1,4)
	IF (N .EQ. 4) GO TO 50
40	CONTINUE
50	CONTINUE
	RETURN
C  1  FORMAT(1H ,2I3,5X,40I2)
C  2  FORMAT(1H ,I3, 4F8.2)
	END

	SUBROUTINE READER(NZ1,NZ2,NR,NSOLMX,WMOL,DOBS,PRNTMR,PRNTLN,
     1  NTST)
C  ---------------------------------------------------------------------
C----THE FIRST CARD MUST CONTAIN THE NAME OF THE PROBLEM (TITLE).
C    THE NEXT CARDS ARE USUALLY COMPLETELY BLANK. THE PROGRAM
C    THEN TAKES THE DEFAULT VALUES.. IF YOU WANT TO SPECIFY YOUR
C    OWN PARAMETERS, YOU SHOULD DO THIS AS FOLLOWS:
C  COLUMN  FORMAT  NAME   DEFAULT   MEANING
C    1      A1     MAN      0       PRINTS INSTRUCTIONS,UNLESS MAN=9
C    2      A1     INSTR    0       PRINTS INSTRUCTIONS FOR SPECIAL
C                                   FEATURES IF INSTR=1
C    3      A1     INTENS   0       READS INTENSITIES BESIDES LINE
C                                   POSITIONS IF INTENS=1
C    4      I1     NSOLMX   4       MAXIMUM NR OF SOLUTIONS TO BE PRINTE
C   5- 6    I2     NSYST(1) 0       ORTHORHOMBIC SOLUTION. +1=YES, -1=NO
C   7- 8    I2     NSYST(2) 0       MONOCLINIC SYSTEM.   0 MEANS  INDIFF
C   9-10    I2     NSYST(3) 0       TRICLINIC SYSTEM
C  11-15   F5.2    TOL2    3.0      TOLERANCE ON 2-DIMENSIONAL SEARCH
C  16-20   F5.2    TOL3    4.5      TOLERANCE ON 3-DIMENSIONAL SEARCH
C  21-30   F10.5   WAVEL   1.54060  WAVELENGTH
C  31-32    I2     LINCO    0       NUMBER OF GIVEN LINE COMBINATIONS
C                                   FOR THE 2-DIMENSIONAL SEARCH
C                                   IF LINCO>0, THE LINECOMBINATIONS ARE
C                                   READ (4F10.2) AFTER THE LINES DECK.
C                                   TERMINATION BY A BLANK CARD (OR ZERO
C     33    I1     LZERCK   0       CHECK ON ZERO-ERROR IF LZERCK>0
C     34    I1     K (ENL)  0       ENL: OUTPUT OF INTERMEDIATE RESULTS(
C  35-36    I2     NQ1      3       NQ1 AND NQ2 DETERMINE THE NUMBERS OF
C  37-38    I2     NQ2      6       LINES TO BE USED IN COMBINATIONS
C                                   FOR FINDING ZONES .
C  39-40    I2     NZ1      6       NZ1 AND NZ2 DO THE SAME FOR
C  41-42    I2     NZ2      6       FINDING LATTICES FROM ZONES.
C  43-44    I2     NR       0       REFINEMENT AND EVALUATION RUN ONLY
C                                   IF NR>0. THE DECK OF COMPLETE
C                                   LATTICES TO BE REFINED (6F10.2)
C                                   MUST FOLLOW THE LINES DECK, ALSO
C                                   TERMINATED BY A CARD WITH ZERO OR
C                                   BLANKS IN COLS 1-10.
C  45-46    I2     INDAT    5       NUMBER OF THE UNIT ON WHICH TO READ
C                                   THE LINE POSITIONS.
C  47-48    I2     LIST     1       LIST OF CALCULATED LINES FOR THE FIR
C                                   *LIST* LATTICES.  SUPPRESS: LIST=-1
C  49-50    I2     J (TEST) 0       TEST OUTPUT IF J .GT. 0
C  51-60   F10.5   WMOL    0.0      MOLECULAR WEIGHT.
C  61-70   F10.5   DOBS    0.0      OBSERVED DENSITY.
C  71-78    F8.5   TOLG    6.0      TOLERANCE ON MATCH BETWEEN CALCULATE
C                                   AND OBSERVED LINES IN HUNDREDTH OF.
C                                   A DEGREE TWOTHETA.
C  79-80    I2     NTST     0       THE NUMBER IN THE MAIN PROGRAM
C                                   WHERE TEST OUTPUT IS GIVEN.
C  ---------------------------------------------------------------------
C    ON THE NEXT CARD, THE FOLLOWING PARAMETERS SHOULD BE ENTERED
C   1-10   F10.5   ZERCOR  0.0      ZEROSHIFT, TO APPLY TO ALL LINES.
C  11-20   F10.5   PRNTMR  4.0      THE MINIMUM VALUE OF THE FIGURE
C                                   OF MERIT OF A LATTICE TO BE PRINTED.
C  21-30   F10.5   PRNTLN  14.      THE MINIMUM NUMBER OF INDEXED LINES
C                                   FOR A LATTICE TO BE PRINTED.
C  31-40   F10.5   RMAX     0.      THE NUMBER OF LINES, ENTER AS REAL.
C                                   ONLY NECESSARY WHEN 20.<RMAX<33. --
C  ---------------------------------------------------------------------
	LOGICAL SELF,ENL,TEST
	CHARACTER  NAME*80,CARD*80,INTSY*2
	CHARACTER*1 MAN,INSTR,INTENS,NINE,IONE,IBLA,PLUS*2
	COMMON /LINES / IDOC,INP,IOUT,MMAX,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/ RAD,TWORAD,RRAD,HLFLAM,WAVSQR,STORE(200,9)
	COMMON /THETAS/ ZERTOT,OBSTET(40),SHFTET(40),TETPT(40),TETMT(40)
	DATA NINE/'9'/,IONE/'1'/,IBLA/' '/,PLUS/'++'/
	INPFIX=INP
	NROSTR=150
C  ---------------------------------------------------------------------
C  ---------------------------------------------------------------------
C  AFTER THE PARAMETER CARD WE MUST ENTER THE LINES,  EITHER AS
C  TWOTHETA, D- OR Q- VALUES.
C  THE PROGRAM WILL FIND OUT WHAT KIND OF VALUES IT HAS.
C  THE NUMBER OF LINES MUST BE AT LEAST 20  AND MUST NOT EXCEED 200.
C  A MAXIMUM OF 40 LINES IS ACTUALLY USED, THE REST IS HARMLESS.
C  ONLY THE FIRST SIX LINES SHOULD BE IN THEIR CORRECT ORDER, THE REST
C  CAN BE ENTERED IN A RANDOM WAY.  EVEN LINES THAT SHOULD BE AMONG
C  THE FIRST SIX CAN BE ENTERED LATER.
C  THE LINES CARDS SHOULD BE TERMINATED BY A BLANK CARD.
C  THE PROGRAM TESTS ON THIS  BLANK CARD AND THEN COUNTS THE NUMBER
C  OF REFLECTIONS.
C  ---------------------------------------------------------------------
9990	CONTINUE
C  ---------------------------------------------------------------------
C--IF YOU SET LINCO  UNEQUAL  ZERO YOU MUST NOW ADD LINCO CARDS.
C  THERE ARE SEVERAL POSSIBILITIES:
C  A. YOU ENTER TWO Q-VALUES (2F10.5) PER CARD.
C    IN THIS CASE THE PROGRAM TREATS THEM AS ONE OF HIS OWN
C    LINE COMBINATIONS AND TRIES TO FIND THE ANGLE IN FINZON  ETC.
C  B. YOU ENTER THREE CONSTANTS (3F10.5).
C     THE PROGRAM RECOGNIZES THIS AS A COMPLETE ZONE, BUT WILL
C     DETERMINE THE QUALITY OF THIS ZONE. THIS MAY LEAD TO AN
C     UNEXPECTED DOUBLING OF AXES.
C     NOTE: A RECTANGULARZONE MUST BE ENTERED AS  A,B, 0.1  FOR READING
C     IN CORRECTLY.  THE PROGRAM WILL VERY SOON  REPLACE 0.1 BY 0.0
C  C. YOU ENTER 4 VALUES  (4F10.5).  THE FOURTH VALUE IS THE
C     QUALITY OF THE ZONE.  IF YOU SET THIS HIGH ENOUGH, THE ZONE
C     WILL CERTAINLY BE USED IN LATTICE-FINDING.
C     ENTERING A SINGLE DIGIT IN COLUMN 31 ASSURES YOU OF A HIGH QUALITY
C  ---------------------------------------------------------------------
C-----------------------------------------------------------------------
C
	CALL ECHO
	READ(CARD,9901) NAME
5	DO 8 I=1,40
8	INTSY(I)=IBLA
	ENL=.TRUE.
	TEST=.TRUE.
	SELF=.FALSE.
	CALL ECHO
C--PARAMETER CARD
	READ(CARD,9902) MAN,INSTR,INTENS,NSOLMX,NSYST,TOL2,TOL3,WAVEL,
     * LINCO,LZERCK,K,NQ1,NQ2,NZ1,NZ2,NR,INDAT,LIST,J,WMOL,DOBS,TOLG,
     3 NTST
	IF (INTSY(50) .EQ. PLUS) GO TO  11
	IF (MAN .NE. NINE)  CALL MANUAL
	IF (INSTR .EQ. IONE) CALL INSTRS
11	CONTINUE
	WRITE(IOUT,9920) NAME
C     WRITE(IOUT,9956)
	WRITE(IOUT,9950)
	IF (NSOLMX .EQ. 0) NSOLMX=4
	IF (TOL2.LT. 0.1) TOL2=3.0
	IF (TOL3.LT. 0.1) TOL3=4.5
	IF (WAVEL.LT.0.1) WAVEL=1.54060
	IF ( K   .EQ. 0) ENL =.FALSE.
	IF (NQ1 .EQ. 0) NQ1=3
	IF (NQ2 .EQ. 0) NQ2=6
	IF (NZ1 .EQ. 0) NZ1=6
	IF (NZ2 .EQ. 0) NZ2=6
	IF (INDAT .EQ. 0) INDAT=INP
	IF (LIST .EQ. 0) LIST=1
	IF ( J   .EQ. 0) TEST=.FALSE.
	TOLKEP=TOLG
	IF (TOLG .LT. 0.1) TOLG=6.0
	IF ( (NTST .EQ. 0) .AND. TEST ) NTST=99
	SELF=ENL
	CALL NTEST(NTST,1)
	INTSY(50)=PLUS
	TOLG=0.01*TOLG
	TOLKEP=0.01*TOLKEP
	TWORAD=90./ATAN(1.)
	RRAD=ATAN(1.)/45.
	RAD=1.0/RRAD
	HLFLAM=0.005*WAVEL
	WAVSQR=(0.5*WAVEL)**2
C  INP WILL BE RESORED AT THE END OF THE ROUTINE.
C  (SEE FIRST AND LAST STATEMENT)
	INP=INDAT
C
C  READ THE NEXT PARAMETER CARD AND ALL LINES
C
	CALL REDLIN(PRNTMR,PRNTLN,INDAT,INTENS,TOLKEP)
C
C
C     PRINT *,' MMAX IS ',MMAX
	IF (MMAX .LT. 20) GO TO 9990
	IF (FMINIM .LT. 2.0) GO TO 9990
	TOLRAD=0.5*TOLG*RRAD
	DO 20 M=1,MMAX
	   Z=ASIN(HLFLAM*SQRT(OBS(M)))
	   OBSTET(M)=Z
	   TETPT(M)=Z+TOLRAD
	   TETMT(M)=Z-TOLRAD
C        PRINT *,' THETA ARRAYS',OBSTET(M),TETPT(M),TETMT(M)
20	CONTINUE
	IF (LZERCK .GT. 0) CALL ZERCHK
	IF (NQ1 .GT. NQ2) THEN
	   I=NQ2
	   NQ2=NQ1
	   NQ1=I
	END IF
	IF (NZ1 .GT. NZ2) THEN
	   I=NZ2
	   NZ2=NZ1
	   NZ1=I
	END IF
	IF (TEST) ENL=.TRUE.
	IF (.NOT. ENL) WRITE(IOUT,9927)
	IF (ENL .AND. (.NOT. TEST)  ) WRITE(IOUT,9923)
	IF (TEST) WRITE(IOUT, 9924)
	WRITE(IOUT,9922) TOL2,TOL3,TOLG
	IF (ABS(TOLG-0.06) .GT. 1.0E-4) WRITE(IOUT,9928) TOLG
	IF (NR .EQ. 0) GO TO 315
C---READ THE COMPLETE LATTICES (REFINEMENT RUN ONLY)
C-----------------------------------------------------------------------
	WRITE(IOUT,9941)
	NNR=NROSTR+1
	I=0
290	CALL ECHO
	READ(CARD,9940) (CEL(J),J=1,6)
C     PRINT 91, (CEL(J),J=1,6)
C  91 FORMAT(' CEL IS',6F10.2)
	IF (CEL(1) .LT. 0.01) GO TO 314
	I=I+1
	IF (NNR-I) 290,300,310
300	WRITE(IOUT,9943) NROSTR
	GO TO 290
310	DO 312 J=1,6
312	ROSTER(I,J)=CEL(J)
C     PRINT 92,I, (ROSTER(I,J),J=1,6)
C  92 FORMAT(' I=',I3,'  ROSTER(I,J)=', 6F10.2 )
	ROSTER(I,7)=99
	ROSTER(I,8)=I
	INDEX(40)=I
	CALL RECCEL
	WRITE(IOUT, 9942) I,(ROSTER(I,J),J=1,6)
	GOTO 290
314	NR=I
	IF (NR.GT.NROSTR) NR=NROSTR
	WRITE(IOUT,9944) NR
	IF (LIST .EQ. 1) LIST=NR
C  WITH THE IDEA, THAT YOU WANT TO SEE AS MUCH AS POSSIBLE
C  FROM THE LATTICES THAT YOU ENTERED.
	GO TO 500
315	IF (LINCO .EQ. 0) GO TO 350
C-----------------------------------------------------------------------
C  THE COMBINATIONS OF LINES WILL BE READ HERE
C-----------------------------------------------------------------------
	NQ1=MMAX
	NQ2=MMAX
	DO 320  I=1,40
	CALL ECHO
	READ(CARD, 9919)   (ROSTER(I,J),J=1,4 )
	IF (ROSTER(I,1) .LT. 0.1) GO TO 330
320	CONTINUE
330	LINCO=I-1
	WRITE(IOUT,9931) LINCO
	WRITE(IOUT,9917) (I,(ROSTER(I,J),J=1,4),I=1,LINCO)
	GO TO 400
C  THE NUMBER OF COMBINATIONS OF NQ1 AND NQ2 IS  NQ1*(2*NQ2-NQ1+1)/2
C  THIS NUMBER SHOULD NOT BE LARGER THAN 40, AS THERE ARE 5LINES PER
C  COMBINATION, AND ONLY ROOM FOR 200 POSSIBILITIES.
C-----------------------------------------------------------------------
C  THE COMBINATIONS OF LINES WILL BE MADE NOW
C-----------------------------------------------------------------------
350	WRITE(IOUT,9932) NQ1,NQ2
	DO 370 I=1,NQ1
	DO 360 J=I,NQ2
	IF ( ABS(4.0*OBS(I)-OBS(J)) .LT. 6.0) GO TO 360
	LINCO=LINCO+1
C     PRINT *,' LINCO=',LINCO,'  I RESP J =', I, J
	ROSTER(LINCO,1)=OBS(I)
	ROSTER(LINCO,2)=OBS(J)
360	CONTINUE
370	CONTINUE
400	CONTINUE
	WRITE(IOUT,9933) NZ1,NZ2
500	CONTINUE
	IF (NSYST(1) .EQ. 1) WRITE(IOUT,9952)
	IF (NSYST(2) .EQ. 1) WRITE(IOUT,9953)
	IF (NSYST(3) .EQ. 1) WRITE(IOUT,9954)
	WRITE(IOUT,9934) WMOL
	WRITE(IOUT,9935) DOBS
	INP=INPFIX
	RETURN
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C  THE FOLLOWING CARDS ARE USEFUL FOR TEST PURPOSES:
C     WRITE(IOUT,9945) TOL2,TOL3,TOLG,WAVEL,LINCO,SELF,ENL,TEST,
C    * NQ1,NQ2,NZ1,NZ2 ,NR,INP,WMOL,DOBS
C9945 FORMAT(53H   TOL2  TOL3  TOLG  WAVEL   LINCO   SELF   ENL  TEST,
C    *33H NQ1 NQ2 NZ1 NZ2 NR INP WMOL DOBS/
C    *  3F6.1, F8.5, I8, 3L6,    6I4 ,2F8.5)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
9901	FORMAT(A80)
9902	FORMAT( 3A1,I1,3I2,2F5.2,F10.5, I2,2I1,8I2, 2F10.5,F8.4,I2)
9917	FORMAT(15H INPUT ZONE NR., I3, 4F8.1)
9919	FORMAT(3F10.2, F10.1)
9920	FORMAT(50H1THIS PROGRAM TRIES TO FIND THE UNITCELL FROM THE ,
     *18H POWDER PATTERN OF ///20X,A80///26H WITH A METHOD PROPOSED BY,
     *61H  DE WOLFF AND PROGRAMMED BY VISSER (J.APPL.CRYST.2,89(1969))/
     *40H            FORTRAN  VERSION    11        )
9922	FORMAT(50H0(COL 11-15) TOL2, THE TOLERANCE ON THE TRIAL F IS,
     1  F4.1//
     2  50H (COL 16-20) TOL3, THE TOLERANCE ON THE TRIAL D IS , F4.1//
     3  48H (COL 71-78) TOLG, THE TOLERANCE ON TWO THETA IS  ,F6.3,
     4  10H  DEGREES.)
9923	FORMAT(38H0INTERMEDIATE RESULTS ARE PRINTED OUT./
     * 58H PUNCH A 1 IN COLUMN 50 OF THE PARAMETER CARD IF YOU WANT  ,
     * 18H THE TEST OUTPUT.   )
9924	FORMAT(39H0TESTRUN, ALL TEST RESULTS ARE PRINTED.)
9927	FORMAT(55H0ONLY A MINIMUM OF INTERMEDIATE RESULTS WILL BE PRINTED/
     * 59H PUNCH A 1 IN COLUMN 34 OF THE PARAMETER CARD IF YOU WANT  ,
     * 14H MORE OUTPUT.  )
9928	FORMAT(48H ***********************************************/
     1       40H * THE OVERALL TOLERANCE IS CHANGED INTO,F6.2,2H */
     2       48H *********************************************** )
9931	FORMAT(21H0THE PROGRAM WILL USE ,I3,24H GIVEN LINE COMBINATIONS)
9932	FORMAT(35H0THE PROGRAM WILL COMBINE THE FIRST ,I3,
     * 21H LINES WITH THE FIRST  ,I3)
9933	FORMAT(42H0THE PROGRAM WILL TRY TO COMBINE THE FIRST  ,I3,
     * 21H ZONES WITH THE FIRST ,I3)
9934	FORMAT(30H0THE GIVEN MOLECULAR WEIGHT IS,  F10.3)
9935	FORMAT(30H0THE GIVEN OBSERVED DENSITY IS,  F10.4)
9940	FORMAT( 8F10.2 )
9941	FORMAT(57H0INPUT OF COMPLETE LATTICES FOLLOWS (REFINEMENT RUN ONLY
     1)/)
9942	FORMAT(18H INPUT LATTICE NR., I3, 8F10.2 )
9943	FORMAT(41H0** EXCEEDS STORAGE ALLOCATED, ONLY FIRST,I4,
     1  21H WILL BE PROCESSED **/  )
9944	FORMAT(30H0LATTICE INPUT FINISHED, FIRST,I4,18H WILL BE PROCESSED)
9950	FORMAT(65H0*******************************************************
     1*********/  2H *,62X,1H*/ 65H *   ARE YOUR INPUT DATA REALLY THE B
     2EST THAT CAN BE OBTAINED   */2H *,62X,1H*/   65H *****************
     3***********************************************    ///
     480H FINDING THE UNIT CELL  DEPENDS FOR 95 PERCENT ON THE QUALITY O
     5F THE INPUT DATA./118H A RANDOM ERROR OF 0.03 DEGREES TWOTHETA CAN
     6 USUALLY BE TOLERATED, BUT A SYSTEMATIC (ZERO POINT) ERROR OF 0.02
     7 DEGREES/77H IS PROBABLY DISASTROUS.  CHECK YOUR INPUT LINES AGAIN
     8ST THEIR HIGHER ORDERS.//  83H DO NOT USE A DEBIJE-SCHERRER CAMERA
     9  UNLESS THE UNIT CELL TO BE EXPECTED IS SMALL./76H WOULD YOU LIKE
     O TO SOLVE A JIGSAW PUZZLE WHEN HALF THE PIECES ARE MISSING---)
9952	FORMAT(50H ORTHORHOMBIC LATTICES ARE GIVEN AN EXTRA CHANCE   )
9953	FORMAT(50H MONOCLINIC   LATTICES ARE GIVEN AN EXTRA CHANCE   )
9954	FORMAT(50H TRICLINIC    LATTICES ARE GIVEN AN EXTRA CHANCE   )
C9956 FORMAT(50H0ATTENTION,  THE INPUT FORMAT HAS BEEN CHANGED
C    1  ,    50H XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  )
C9998 FORMAT(16H END OF PROGRAM   )
	END

	SUBROUTINE REAR
	DIMENSION ZZ(7,2),IJ(2)
C  THIS SUBROUTINE LOOKS FOR THE COMMON ROW IN TWO ZONES . THE EXISTENCE
C  OF TWO EQUAL Q-VALUES IS ENOUGH.
	COMMON /LINES / IDOC,INP,IOUT,MAX,OBS(40),CEL(10)
	J=0
	I=1
10	ZZ(1,I)=CEL(J+1)
	ZZ(2,I)=CEL(J+2)
	ZZ(3,I)=CEL(J+1)+CEL(J+2)+CEL(J+3)
	ZZ(4,I)=CEL(J+1)+CEL(J+2)-CEL(J+3)
	ZZ(5,I)=CEL(J+3)
	ZZ(6,I)=CEL(J+1)*2.0+CEL(J+3)
	ZZ(7,I)=CEL(J+1)*2.0-CEL(J+3)
	IF (J.GT.1) GO TO 20
	I=2
	J=3
	GO TO 10
20	DO 50 I=1,4
	DO 40 J=1,4
	IJ(1)=I
	IJ(2)=J
	IF (ABS(ZZ(I,1)-ZZ(J,2)).LT. 1.0 ) GO TO 100
40	CONTINUE
50	CONTINUE
	CEL(4)=0.0
	RETURN
100	CEL(1)=0.0
	DO 150 N=1,2
	L=7-N
	K=IJ(N)
	CEL(1)=0.5*ZZ(K,N)+CEL(1)
	IF (K.GT.1) GO TO 120
	CEL(N+1)=ZZ(2,N)
	CEL(L)   =ZZ(5,N)
	GO TO 150
120	CEL(N+1)=ZZ(1,N)
	CEL(L)   =ZZ(K+3,N)
150	CONTINUE
	RETURN
	END

	SUBROUTINE RECCEL
C  CHECKS WHETHER THE INPUT LATTICE IS IN ANGSTROMS AND DEGREES
C  AND CALCULATES A,B,C,D,E,F IF IT IS.
	DIMENSION WAAR(6),WAR(6),CL(6)
	COMMON /LINES / IDOC,INP,IOUT,MAX,OBS(40),CEL(10)
	COMMON /TEMRES/ LINCO,NQ2,NINDEX,INDEX(40),ROSTER(150,12)
	RAD=57.29578
	NR=INDEX(40)
	DO 10 J=1,6
10	WAAR(J)=ROSTER(NR,J)
	IF(WAAR(4) .EQ. 0.0 .AND. WAAR(6) .EQ. 0.0) GO TO 150
	IF(ABS(WAAR(5)) .LT. 2.0*SQRT(WAAR(1)*WAAR(3))) GO TO 150
	WRITE(IOUT,9909)(WAAR(J),J=1,6)
	A=1.0
	B=2.0
	DO 130 J=1,3
	C=COS(WAAR(J+3)/RAD)
	CL(J)=C
	CL(J+3)=C
	A=A-C*C
	B=B*C
	WAR(J)=WAAR(J)
	WAR(J+3)=WAR(J)
	CEL(J)=WAAR(J+3)/RAD
	CEL(J+3)=CEL(J)
130	CONTINUE
	FVOL=10000.0/(A+B)
	DO 140 J=1,3
	WAAR(J)=FVOL*(SIN(CEL(J))/WAR(J))**2
	WAAR(J+3)=2.0*FVOL*(CL(J+1)*CL(J+2)-CL(J))/(WAR(J+1)*WAR(J+2))
140	CONTINUE
	DO 145 J=1,6
145	ROSTER(NR,J)=WAAR(J)
150	RETURN
9909	FORMAT(30H INPUT LATTICE OF DIRECT CELL   ,6F10.4)
	END

	SUBROUTINE REDLIN(PRNTMR,PRNTLN,INDAT,INTENS,TOLKEP)
	DIMENSION QIN(200)
	LOGICAL SELF,ENL,TEST,REPEAT
	CHARACTER NAME*80,CARD*80,INTSY*2,INTENS*1,IBLA*1
	CHARACTER*2 IIN(200),INTY(8),BLANK
	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 /CONSTS/ RAD,TWORAD,RRAD,HLFLAM,WAVSQR,STORE(200,9)
	DATA IBLA/' '/,BLANK/'  '/
	INPFIX=INP
	REPEAT=.FALSE.
	DO 6 I=1,8
	CEL(I)=0.
6	INTY(I)=BLANK
	DO 7 I=1,200
7	IIN(I)=BLANK
C
C  READ SECOND PARAMETER CARD. IF NOT PRESENT,REPAIR SITUATION.
C
	CALL ECHO
	READ(CARD,9903)  ZERCOR,PRNTMR,PRNTLN,RMAX, (CEL(J),J=5,8)
	MMAX=NINT(RMAX)
	L=0
	IF (ZERCOR .LT. 3.0) GO TO 9
C  MAKES REPAIRS IF THE  NEW SECOND PARAMETER CARD WAS NOT INSERTED.
	CEL(4)=RMAX
	CEL(3)=PRNTLN
	CEL(2)=PRNTMR
	CEL(1)=ZERCOR
	ZERCOR=0.
	PRNTMR=4.0
	PRNTLN=14.
	MMAX=0
	WRITE(IOUT,9960)
	GO TO 14
C  SECOND PARAMETER CARD WAS PRESENT.
9	IF ( PRNTMR .LT. 0.1) PRNTMR=4.0
	IF (PRNTLN .LT. 0.1) PRNTLN=14.
	IF (INDAT .NE. INP) INP=INDAT
C  NOTE THAT INPFIX=INP IN THE FIRST STATEMENT OF THIS ROUTINE
C  AND THAT  INP=INPFIX IN THE LAST STATEMENT.
C  INDAT IS THUS ONLY VALID  WITHIN THIS SUBROUTINE.
10	CONTINUE
	CALL ECHO
	IF (INTENS .EQ. IBLA) THEN
	   READ(CARD,9903)  (CEL(J),J=1,8)
	ELSE
	   READ(CARD,9912)  (CEL(J),INTY(J),J=1,8)
	END IF
14	U=0.0
	DO 15 I=1,8
15	U=U+ABS(CEL(I) )
C  IT IS POSSIBLE TO GET RID OF A LINE BY SIMPLY PUNCHING A MINUS  SIGN
C  IN FRONT OF THE VALUE(PROVIDED THERE WAS ENOUGH SPACE IN THE FORMAT).
C  THIS TRICK NOW ALSO WORKS WITH ONLY ONE LINE PER CARD.
	IF (U .LT. 0.1) GO TO 30
	DO  20 I=1,8
	IF (CEL(I) .LT. 0.1) GO TO 20
	L=L+1
	QIN(L)  =CEL(I)
	IIN(L)=INTY(I)
20	CONTINUE
	IF (L .GT. 190) L=190
C  PREVENTS OVERFLOW OF THE INPUT  ARRAYS.  NOW THERE IS NO LIMIT ON
C  THE NUMBER OF LINES THAN CAN BE READ.
	GO TO 10
C-----------------------------------------------------------------------
C--END OF READING LINES CARDS
C-----------------------------------------------------------------------
30	L=L+1
	DO 40 I=1,8
	MAX=L-I
	IF (QIN(MAX) .GT. 0.1) GO TO 90
40	CONTINUE
90	WRITE(IOUT,9911)  WAVEL
	B=1.0E+4/WAVSQR
	WRITE(IOUT, 9918)  MAX
	WRITE(IOUT, 9915) ( QIN(M),IIN(M) ,M=1,MAX)
C-----------------------------------------------------------------------
C--WE NOW KNOW THE NUMBER OF LINES. NEXT WE TRY TO FIND OUT,
C--WHAT KIND OF VALUES WE HAVE  AND PRINT THEM OUT.
C-----------------------------------------------------------------------
	IF (MAX .GE. 40) MAX=40
	DO 140 M=1,MAX
	OBS(M)=QIN(M)
	INTSY(M)=IIN(M)
140	CONTINUE
	DO 160 M=1,5
	IF(OBS(M+1) .LT. OBS(M)) GO TO 195
160	CONTINUE
C-----------------------------------------------------------------------
C--IF INPUT IS IN DESCENDING ORDER, WE HAVE  D-VALUES
C-----------------------------------------------------------------------
	ISTYG=1
	CALL SRTRON(ISTYG)
	IF(OBS(MAX) .GE. 180.0) GO TO 190
C-----------------------------------------------------------------------
C--IF THE LAST INPUT VALUE IS GREATER THAN 180 WE DO NOT HAVE TWOTHETAS.
C-----------------------------------------------------------------------
	WRITE(IOUT,9906)
	WRITE(IOUT,9909)
	DO 180 M=1,MAX
	W=OBS(M)+ZERCOR
	TWTHET(M)=W
	D=SIN(0.5*RRAD*W)
	OBS(M)=B*D*D
180	CONTINUE
	WRITE(IOUT,9911)  WAVEL
	GO TO 220
190	CONTINUE
C  INPUT IS  Q-VALUES
	WRITE(IOUT,9921)  MAX
	WRITE(IOUT,9906)
	WRITE(IOUT,9908)
	DO 192 M=1,MAX
	TWTHET(M)=TWORAD*TETCAL(OBS(M))+ZERCOR
	IF (ABS(ZERCOR) .LT. 1.E-8) GO TO 192
	D=SIN(0.5*RRAD*TWTHET(M))
	OBS(M)=B*D*D
192	CONTINUE
	GO TO 230
C  INPUT IS D-VALUES.
195	ISTYG=0
	CALL SRTRON(ISTYG)
	WRITE(IOUT,9906)
	WRITE(IOUT,9907)
	DO 210 M=1,MAX
	SIT=0.5*WAVEL/OBS(M)
	TWTHET(M)=TWORAD*ASIN(SIT)+ZERCOR
	D=SIN(0.5*RRAD*TWTHET(M))
	OBS(M)=B*D*D
210	CONTINUE
220	WRITE(IOUT,9921)  MAX
	IF (ABS(ZERCOR) .GT. 1.0E-4) WRITE(IOUT,9945) ZERCOR
	WRITE(IOUT,9910)
230	WRITE(IOUT,9916)  ( OBS(M),INTSY(M) ,M=1,MAX)
	IF (MAX .GT. 32 ) THEN
	   GO TO 240
	ELSE IF (MAX .LT. 20) THEN
	   WRITE(IOUT, 9930) MAX
	ELSE IF (MMAX .EQ. 0) THEN
	   WRITE(IOUT,9931)
	ELSE IF (MMAX .NE. MAX) THEN
	   WRITE(IOUT,9932) MAX,MMAX
	ELSE
	   WRITE(IOUT,9934) MAX
	   GO TO 240
	END IF
	MAX=10
	GO TO 9990
240	CONTINUE
C-----------------------------------------------------------------------
	
C--THE ARRAY OF Q-VALUES HAS NOW BEEN ESTABLISHED.
C-----------------------------------------------------------------------
	IF (ABS(ZERCOR) .LT. 1.E-4) GO TO 265
	IF (.NOT. TEST) GO TO 265
	WRITE(IOUT,9946) ZERCOR, (TWTHET(M),M=1,MAX), (OBS(M),M=1,MAX)
265	CONTINUE
	IF (TEST) WRITE(IOUT,9936)
	U=0.
	DO 270 M=1,MAX
	SIT=SIN(0.5*RRAD*(TWTHET(M)+TOLG) )
	W=SIT/HLFLAM
	OBSPT(M)=W*W
	OBSMT(M)=2.*OBS(M)-OBSPT(M)
	IF (M .LE. 20) U=U+OBSPT(M)-OBSMT(M)
	IF (TEST) WRITE(IOUT,9937) OBS(M),OBSMT(M),OBSPT(M)
270	CONTINUE
C--THE RANGE FOR MATCHING Q-VALUES HAS BEEN SET
	U=U/40.
	IF (TEST) WRITE(IOUT,9938) TOLG,U
	IF (REPEAT) GO TO 280
	Q20=OBS(20)
	N20=50
	M20=10
	DELQMX=Q20/FLOAT(2*N20*M20)
C  DELQMAX IS THE MAXIMUM VALUE OF THE MEAN DELTA Q THAT IS
C  STILL CONSISTENT WITH AN ELEMENTARY CELL THAT HAS A
C  FIGURE OF MERIT THAT IS GREATER THAN 10.
C    THE MEAN TOLERATED DELTA Q IS U, CORRESPONDING TO A
C  MEAN DELTA TWO THETA OF *TOLG*
C  THE MAXIMUM TOLERABLE TOLG IS THUS (DELQMX/ U )*TOLG
	TOLMAX=TOLG*DELQMX/U
	WRITE(IOUT,9964)
	WRITE(IOUT,9965) TOLMAX
	IF ( ABS(TOLG-TOLKEP) .LT. 1.E-4) GO TO 280
C  IF A TOLERANCE WAS GIVEN ORIGINALLY, TOLG=TOLKEP.
C  IN THAT CASE, EVERYTHING REMAINS AS IT WAS.
C  IF NO TOLERANCE WAS GIVEN, WE SET IT AT 2*TOLMAX.
	TOLG=2.*TOLMAX
	IF (TOLG .GT. 0.06) TOLG=0.06
	TOLKEP=TOLG
	REPEAT=.TRUE.
	GO TO 265
280	FMINEQ=OBS(20)*5.0E-4
	IF (FMINEQ .LT. 0.8)  FMINEQ=0.8
	A=8.0 *FMINEQ
	C=1.6*FMINEQ
	DO 290 M=2,MAX
	B=OBS(M)-OBS(M-1)
	IF  (B .GT. A) GO TO 290
	A=B
	IF ( B .GE. C) GO TO 290
	I=M-1
	D=RAD*WAVEL*B/( 200.0*SQRT(OBS(M)))
	WRITE(IOUT,9914) I,B,D
290	CONTINUE
	FMINEQ=0.125*A
	IF (FMINEQ .LT. 0.5)  FMINEQ=0.5
	WRITE(IOUT,9926) FMINEQ
	FMINIM=OBS(1)/17.
	A=100./SQRT(FMINIM)
	WRITE (IOUT,9925) FMINIM,A
	IF (FMINIM  .LE. 2.0) WRITE(IOUT,9929)
9990	INP=INPFIX
	RETURN
9903	FORMAT( 8F10.5 )
9906	FORMAT(23H0THE INPUT CONSISTED OF   )
9907	FORMAT(1H+, 23X,10H D-VALUES.)
9908	FORMAT(1H+, 23X,10H Q-VALUES.)
9909	FORMAT(1H+, 23X,12H TWO THETAS. )
9910	FORMAT(45H0THE INPUT LEADS TO THE FOLLOWING Q-VALUES   )
9911	FORMAT(20H0THE WAVELENGTH IS  , F8.4)
9912	FORMAT( 8(F7.3,A2,1X) )
9914	FORMAT(40H0XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX      /
     * 41H THE DIFFERENCE IN Q-VALUES BETWEEN LINE ,  I3,
     * 22H AND THE NEXT IS ONLY  , F6.2 /
     * 14H THIS IS ONLY  ,F5.2,  18H DEGREES TWO THETA     /
     * 40H            PLEASE  CHECK  YOUR  INPUT     /
     * 56H THE PROGRAM WILL CONTINUE, BUT MAY NOT GIVE THE CORRECT,
     * 59H ANSWERS SINCE THE LIMIT FOR  *EQUAL* VALUES CAN BE TOO LOW /
     * 40H XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX      )
9915	FORMAT( 10(1X,F8.3,1X,A2,1X )  )
9916	FORMAT( 10(1X,F7.1,2X,A2,1X )  )
9918	FORMAT(37H0THE INPUT CONSISTED OF THE FOLLOWING,I4,7H LINES   )
9921	FORMAT(41H0THE NUMBER OF LINES THAT WILL BE USED IS ,I4)
9925	FORMAT(44H0THE LOWER LIMIT FOR THE VALUES OF A,B,C IS ,  F8.2/
     * 40H THE MAXIMUM EDGE OF THE UNIT CELL IS    ,F6.1,10H ANGSTROM )
9926	FORMAT(67H0Q-VALUES WILL BE CONSIDERED TO BE EQUAL WHEN THEY DIFFE
     *R LESS THAN , F4.1 )
9929	FORMAT(50H0FROM YOUR DATA ONE EXPECTS A VERY LARGE UNIT CELL,
     1       50H.  THIS REQUIRES AN ACCURACY THAT IS BETTER THAN  /
     2       50H 0.005 DEGREES TWO THETA IN YOUR OBSERVATIONS.    /
     4       50H     READ THE PAPER BY P.E.WERNER IN J.APPL.CRYST.,
     5       50H(1976) 9, P.216.                                  /
     6       50H THE PROGRAM CANNOT HANDLE THIS CASE WITHOUT SOME ,
     7       50HMODIFICATION. THIS PROBLEM IS SKIPPED.            )
9930	FORMAT(32H0THE NUMBER OF USEFUL LINES  IS ,I4 /
     *60H SINCE THIS IS LESS THAN THE MINIMUM NUMBER OF LINES THAT TH ,
     *60HE PROGRAM NEEDS (20),YOU CANNOT EXPECT RELIABLE  RESULTS.   /
     *24H THIS PROBLEM IS SKIPPED  //
     4 ' IF THIS IS REALLY THE MAXIMUM NUMBER OF LINES THAT CAN' ,
     5 ' BE MEASURED,' / ' THE PROBLEM CAN ONLY BE SOLVED WHEN THE ',
     6 'LATTICE HAS AT LEAST ORTHORHOMBIC SYMMETRY.' /
     7 ' THERE ARE MANY PROGRAMS THAT ARE MUCH BETTER SUITED TO THOSE',
     8 ' SYMMETRIES '  )
9931	FORMAT('0  THE NUMBER OF LINES ENTERED FOR THIS PROBLEM IS RATHER
     1SMALL.'/' WITH ONLY 20 LINES YOUR CHANCES ARE 50:50 THAT YOU ',
     2 'WILL FIND A SUBLATTICE' /' OR A DERIVED LATTICE.  WITH 40 LINES
     3 YOUR CHANCES OF FINDING THE CORRECT' / ' LATTICE IMMEDIATELY',
     4 ' ARE MUCH BETTER.  PLEASE ENTER MORE LINES.'/
     5 ' IF YOU CAN FIND ONLY 32 LINES OR LESS, ENTER THE EXACT ' /
     6 ' NUMBER OF LINES IN COLUMNS 31-40 OF THE SECOND PARAMETER CARD'/
     7 ' AS A REAL NUMBER (E.G.  27.0 ) AND RUN THE PROBLEM AGAIN.'  )
9932	FORMAT(/' THE NUMBER OF LINES ACTUALLY READ IS ',I2,/
     2 ' WHILE THE NUMBER OF LINES GIVEN WAS ',I2, /
     3 ' OBVIOUSLY SOME ERROR WAS MADE.  THIS PROBLEM IS SKIPPED.'  )
9934	FORMAT(1H0,39X,'******',/ 39X,'*THE   *'/38X,'*NUMBER  *'
     1 ,/37X,'*OF LINES  *',/ 36X,
     2 '*ENTERED IS  *'/35X,'*ONLY ',I2,'. THIS *'/34X,
     3 '*MEANS THAT YOU  *'/33X,'*MAY FIND A LATTICE*'/32X,
     4 '*THAT IS A DERIVATIVE*'/31X,'*OF THE CORRECT LATTICE*'/
     5 30X, 26(1H*)  )
9936	FORMAT(63H0THE Q-VALUES WITH THEIR LIMITS OF ERROR,Q +/- TOLERANCE
     1 (TOLG)/ 32H       Q        Q-TOL     Q+TOL    )
9937	FORMAT(1H ,3F10.1)
9938	FORMAT(20H0WITH A TOLERANCE OF, F6.3,18H DEGREES TWOTHETA,   ,
     2 25H THE MEAN DELTA Q BECOMES, F7.2)
9945	FORMAT(21H0A ZERO CORRECTION OF   , F6.3,
     1       44H  DEGREES TWOTHETA WAS APPLIED TO ALL LINES   )
9946	FORMAT(56H0TWOTHETAS AND Q S AFTER APPLYING A ZERO CORRECTION OF
     2  , F10.4/ 4(1H ,10F12.2/)  // 4(1H ,10F12.1/)  )
C9950 FORMAT(1H0///80(1H*)/2H *,78X,1H*/2H *,5X,16HTHE MAXIMUM MEAN
C    2 22H ERROR IN YOUR DATA IS, F6.3,18H DEGREES TWO THETA, 11X,
C    3 1H*/2H *,78X,1H*/1H ,80(1H*) )
9960	FORMAT(50H0*************************************************/
     2       50H *    YOU FORGOT TO INSERT THE NEW PARAMETER CARD*/
     3       50H *THE PROGRAM TOOK DEFAULT VALUES FOR  ZERCOR    */
     4       50H * AND PRNTMR     AND  CONTINUES.......          */
     5       50H *************************************************)
9964	FORMAT(48H0THIS PROGRAM CALCULATES A NUMBER OF UNIT CELLS.
     1 /     51H THESE CELLS ARE ALL TESTED BY DE WOLFF'S CRITERION,
     2        8H (READ:   /
     3       50H A SIMPLIFIED CRITERION FOR THE RELIABILITY OF A  ,
     4       50HPOWDER PATTERN INDEXING, BY P.M. DE WOLFF ,       /
     5       50H JOURNAL APPLIED CRYSTALLOGRAPHY (1968) 1, 108  )./
     6       50H IN ORDER THAT THE UNIT CELL BE ABLE TO PASS THE  ,
     7       50HTEST (M20 .GT. 10), AND ASSUMING THAT YOU DID NOT /
     8       50H MISS TOO MANY LINES (N20=50), THE CRITERION DEMAN,
     9       50HDS THAT THE MEAN ERROR IN YOUR DATA IS            )
9965	FORMAT(1H0,20X,50(1H*)/1H ,20X,14H*   LESS THAN , F8.3,5X,
     1       17HDEGREES TWO THETA,5X,1H* / 1H ,20X,50(1H*) )
	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
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),QPTOL(40),QMTOL(40)
	DIMENSION KSTOR(40,5),LSTOR(40,5),QSTOR(40,5)
     1  ,NMATCH(45),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/ RAD,TWORAD,RRAD,HLFLAM,WAVSQR,STORE(200,9)
	COMMON /THETAS/ ZERTOT,OBSTET(40),SHFTET(40),TETPT(40),TETMT(40)
	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(NMATCH)
	   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)
	   WRITE(IOUT,9911) U,D,QOBS,INTSY(P)
	END IF
	IF (NR .EQ. 0) THEN
	   IF (DRUK) WRITE(IOUT,9912)
	   GO TO 1180
	END IF
	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
	      WRITE(IOUT,9914) U,D,Q,H,K,L
	      IF (DUBL) WRITE(IOUT,9915)
	   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('0AT ENTERING REFINE, CEL= ',6F10.3,4F8.2 )
9910	FORMAT(41H0TWOTHETA     D       Q       H    K    L)
9911	FORMAT(1H0, F8.2, F8.3, F9.1, 18H    *OBSERVED**    ,A2)
9912	FORMAT(1H+,45X ,  12H NOT INDEXED)
9914	FORMAT(1H ,F8.2,F8.3,F9.1, 3I5)
9915	FORMAT(1H+,40X,32H HKL-COMBINATION WAS USED BEFORE  )
9917	FORMAT(36H0LATTICE PRINTED OUT FOR INSPECTION //
     156H DISPLAYED IN LAYERS UP --K-- ,    --L-- ACROSS THE PAGE,
     222H --H-- DOWN THE PAGE.   //)
9918	FORMAT( 8H ***  H ,I2,8H  L  ***   // )
9923	FORMAT(1H0)
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(1H0,' ZEROSHIFT IS ',F10.4,'  DEGREES TWOTHETA')
9931	FORMAT(1H ,14X,19HMATRIX IS SINGULAR )
	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---------------------------------------------------------------
	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 /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-100) 1050,1050,1060
1050	LIST2=.TRUE.
	PRIN=.FALSE.
	CALL REFINE(IREFCY,ALL,DRUK,PRIN,LIST2)
	WRITE(IOUT,9941) CEL,NAME,NSOL
	GO TO 1100
1060	WRITE(IOUT,9942) N20
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) LIST2=.TRUE.
	CALL REFINE(IREFCY,ALL,DRUK,PRIN,LIST2)
	RETURN
9902	FORMAT(1H1,20X,A80 /
     *       50H0THE BEST SOLUTION IS NOW USED TO TRY AND INDEX AL,
     *        9HL LINES.  )
9916	FORMAT(1H1,A80,14H   SOLUTION NR,I3)
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(35H0THE NUMBER OF CALCULATED LINES IS  ,I4,
     *50H. 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 SRTRON(ISTYG)
C
C
C     THIS SUBPROGRAM SORTS THE TWOTHETA-,D- OR Q-VALUES (POSSIBLY
C     TOGETHER WITH INTENSITIES).THE TWOTHETA- OR Q-VALUES ARE PUT IN
C     INCREASING ORDER,THE D-VALUES IN DECREASING ORDER.
C
C     ISTYG=1 - DATA ARE TWOTHETA- OR Q-VALUES (TO PUT IN INCREASING
C               ORDER).
C     ISTYG=0 - DATA ARE D-VALUES (TO PUT IN DECREASING ORDER).
	CHARACTER NAME*80,CARD*80,INTSY*2,IINT*2
	COMMON /LINES / IDOC,INP,IOUT,MAX,OBS(40),CEL(10)
	COMMON /CHARCM/ NAME,CARD,INTSY(50)
C
	FAKTOR=1.0
	IF(ISTYG.EQ.0)FAKTOR=-1.0
20	IEIND=MAX-1
C
	DO 70 I=1,IEIND
	IF(FAKTOR*(OBS(I+1)-OBS(I)))30,70,70
30	OBSTR=OBS(I+1)
	IINT=INTSY(I+1)
	DO 40 J=1,I
C
	IF(FAKTOR*(OBSTR-OBS(J)))50,40,40
40	CONTINUE
C
50	MEIND=I+1-J
C
	DO 60 M=1,MEIND
	K=I+1-M
	OBS(K+1)=OBS(K)
	INTSY(K+1)=INTSY(K)
60	CONTINUE
	OBS(J)=OBSTR
	INTSY(J)=IINT
70	CONTINUE
	DO 80 I=1,IEIND
	IF (ABS(OBS(I+1)-OBS(I) ) .GT. 0.0001) GO TO 80
	WRITE(IOUT,99) OBS(I)
	OBS(I+1)=OBS(MAX)
	INTSY(I+1)=INTSY(MAX)
	MAX=MAX-1
	GO TO 20
80	CONTINUE
C
	RETURN
99	FORMAT(41H0XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX   /
     * 17H THE INPUT VALUE  ,F10.5, 15H APPEARS TWICE   /
     * 27H           CHECK YOUR INPUT    /
     * 41H XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX   )
	END

	SUBROUTINE  SRTZON(NSRT,NR,NZ2)
C--THE SUBROUTINE SORTS ZONES THAT ARE STORED IN ARRAY  *STORE*
C--SORTING WILL BE DONE ON ELEMENT *NSRT*, WHICH CAN BE EITHER 4
C--OR 7. IN THE LATTER CASE WE MAKE A DOUBLE SORT ON ELEMENT 8.
	LOGICAL SELF,ENL,TEST,EXCH,ABSNT
	COMMON /PARMS / TOL2,TOL3,WAVEL,FMINEQ,FMINIM,SELF,ENL,TEST,LIST,
     *                NSYST(3),OBSMT(40),OBSPT(40),TWTHET(40),TOLG
	COMMON /LINES / IDOC,INP,IOUT,MAX,OBS(40),CEL(10)
	COMMON /CONSTS/ RAD,TWORAD,RRAD,HLFLAM,WAVSQR,STORE(200,9)
	IF (ENL) WRITE(IOUT,9910)
C--FIRST FIND LARGEST ELEMENT
	N=1
460	U=-100.
	DO 470 M=N,NR
	   IF (STORE(M,NSRT) .GT. U) THEN
	      U=STORE(M,NSRT)
	      NRM=M
	   END IF
470	CONTINUE
C--PUT LARGEST ELEMENT IN PROPER POSITION
	DO 480 M=1,9
	   U=STORE(N,M)
	   STORE(N,M)=STORE(NRM,M)
480	STORE(NRM,M)=U
	IF (ENL) WRITE(IOUT,9911) N,(STORE(N,I),I=1,8),TWORAD*STORE(N,9)
	N=N+1
	IF (N .LE. NR) GO TO 460
C-----------------------------------------------------------------------
C--END OF SORT ON QUALITY. NOW SEARCH FOR *EQUAL* ZONES.
C--STORE ONLY THE *UNIQUE* ZONES IN  *STORE*
C  ---------------------------------------------------------------------
	IF (NSRT .EQ. 7) GO TO 580
C--THE *EQUAL* ZONES NEED TO BE THROWN AWAY ONLY ONCE.
	NUNIQ=1
	DO 560 I=2,NR
C--CHECK WHETHER THIS ZONE IS ALREAY PRESENT IN THE  *UNIQUE*  SET.
	   ABSNT=.TRUE.
	   DO 520 J=1,NUNIQ
	      U=.0
	      DO 510 K=1,3
510	      U=U+ABS(STORE(I,K)-STORE(J,K))
	      IF (U .LT. FMINEQ) ABSNT=.FALSE.
520	   CONTINUE
C--IF ABSENT THEN ADD THIS ZONE
	   IF (ABSNT) THEN
	      NUNIQ=NUNIQ+1
	      DO 530 K=1,9
530	      STORE(NUNIQ,K)=STORE(I,K)
	   END IF
560	CONTINUE
	STORE(NUNIQ+1,1)=.0
	NR=MIN(NUNIQ,(NZ2+2) )
	GO TO 750
C--MAKE THE DOUBLE SORT ON NC (NO HAS BEEN SORTED ON)
580	NSRT1=NSRT+1
600	EXCH=.FALSE.
	DO 700 I=1,NR
	IF ( (STORE(I,NSRT)-STORE(I+1,NSRT) ) .GT. 0.1 ) GO TO 700
	IF ( (STORE(I,NSRT1)-STORE(I+1,NSRT1) ) .GT. 0.1 ) THEN
	   DO 610 K=1,9
	      U=STORE(I+1,K)
	      STORE(I+1,K)=STORE(I,K)
610	      STORE(I,K)=U
	   EXCH=.TRUE.
	END IF
700	CONTINUE
	IF (EXCH) GO TO 600
	NR=MIN( NR,(NZ2+2) )
750	WRITE(IOUT,9912)
	WRITE(IOUT,9911) (I,(STORE(I,J),J=1,8),TWORAD*STORE(I,9),I=1,NR)
	WRITE(IOUT,9913)
	RETURN
9910	FORMAT('1ZONES AFTER EVALUATION, IN DESCENDING ORDER OF ' ,
     1  ' RELIABILITY'///
     2  ' NEWNR    A       B     FMAAS       QUALITY   OLDNR CENTER ',
     3  ' NOBS  NCALC  ZEROSHIFT' / )
9911	FORMAT(1H , I4, 3F8.1,2X,F14.1, 4F6.0, F10.4)
9912	FORMAT(38H1ZONES SELECTED FOR FURTHER REFINEMENT  ///
     * ' NEWNR    A       B     FMAAS       QUALITY   OLDNR CENTER ',
     2 ' NOBS  NCALC  ZEROSHIFT ')
9913	FORMAT(////47H NOTE--SOME OF THESE ZONES ARE PROBABLY CORRECT,
     *56H EVEN IF THE PROGRAM FAILS TO BUILD A PLAUSIBLE LATTICE ,
     *11H FROM THEM.)
	END

	SUBROUTINE  SYMTST
C  SUBROUTINE SYMTEST TRIES TO FIND OUT WHETHER THE LATTICES FOUND
C  ARE LATTICES OF A HIGHER SYMMETRY IN DISGUISE.
	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
	NPASS=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
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)
	IF (NPASS .EQ. 1) WRITE(IOUT,921)
	IF (NPASS .EQ. 2) WRITE(IOUT,922)
	WRITE(20,900)  (CEL(J),J=1,6)
	IF (NPASS .EQ. 1) WRITE(20,921)
	IF (NPASS .EQ. 2) WRITE(20,922)
	A=200./SQRT(B)
	C=100./SQRT(Z)
	WRITE(IOUT,901) X,Z, A, C
	WRITE(20,901) X,Z, A, C
	X=4.*X
	A=0.5*A
	IF (CEL(9).LT. 6.) WRITE(IOUT,902) X,A
	IF (CEL(9).LT. 6.) WRITE(20,902) X,A
	RETURN
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)
	IF (NPASS .EQ. 1) WRITE(IOUT,921)
	IF (NPASS .EQ. 2) WRITE(IOUT,922)
	WRITE(20,900)  (CEL(J),J=1,6)
	IF (NPASS .EQ. 1) WRITE(20,921)
	IF (NPASS .EQ. 2) WRITE(20,922)
	X=100./SQRT(G)
	WRITE(IOUT, 911) G,X
	WRITE(20, 911) G,X
	RETURN
1020	IF (IK*IL .EQ. 0) GO TO 1050
	N=L/K
	X=B
	Z=A
	GO TO 1150
1050	IF (IH*IL .EQ. 0) GO TO 1100
	N=L/H
	X=A
	Z=B
	GO TO 1150
1100	N=H/K
	X=B
	Z=C
1150	WRITE(IOUT,900)  (CEL(J),J=1,6)
	IF (NPASS .EQ. 1) WRITE(IOUT,921)
	IF (NPASS .EQ. 2) WRITE(IOUT,922)
	WRITE(20,900)  (CEL(J),J=1,6)
	IF (NPASS .EQ. 1) WRITE(20,921)
	IF (NPASS .EQ. 2) WRITE(20,922)
	A=100./SQRT(X)
	C=100./SQRT(Z)
	IF (N-3) 1025,1030,1040
1025	WRITE(IOUT,912) X,Z, A,C
	WRITE(20,912) X,Z, A,C
	RETURN
1030	A=2.0*A/SQRT(3.0)
	WRITE(IOUT,913) X,Z,A,C
	WRITE(20,913) X,Z,A,C
	RETURN
1040	WRITE(IOUT,914) N
	WRITE(20,914) N
1200	CONTINUE
C     IF (NPASS .EQ. 1) GO TO 5
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=INT(G/X+0.5)
	IF (ABS(G-X*FLOAT(L)) .LE. DIFMIN) N=N+1
3030	CONTINUE
	IF (N .LT. 6) GO TO 3050
	GO TO 3100
3050	CONTINUE
	RETURN
3100	WRITE(IOUT,900) (CEL(J),J=1,6)
	IF (NPASS .LE. 1) WRITE(IOUT,921)
	IF (NPASS .GE. 2) WRITE(IOUT,922)
	WRITE(20,900) (CEL(J),J=1,6)
	IF (NPASS .LE. 1) WRITE(20,921)
	IF (NPASS .GE. 2) WRITE(20,922)
	Z=100./SQRT(X)
	WRITE(IOUT,911) X,Z
	WRITE(20,911) X,Z
	RETURN
900	FORMAT(52H0THE LATTICE WITH THE RECIPROCAL LATTICE CONSTANTS  ,
     * 6F9.2 )
901	FORMAT(1H+, 12X, 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+,12X,44H CUBIC WITH A RECIPROCAL LATTICE CONSTANT OF  ,
     *        F10.1,33H AND A DIRECT LATTICE CONSTANT OF , F10.2)
912	FORMAT(1H+,12X,42H TETRAGONAL WITH THE RECIPROCAL CONSTANTS  ,
     * 2F10.1,29H AND THE DIRECT CONSTANTS  A=, F7.2,4H  C=,F7.2 )
913	FORMAT(1H+,12X,41H HEXAGONAL WITH THE RECIPROCAL CONSTANTS   ,
     * 2F10.1,29H AND THE DIRECT CONSTANTS  A=, F7.2,4H  C=,F7.2 )
914	FORMAT(1H+,12X,43H A FREAK. THERE IS A RATIO IN THE CONSTANTS   ,
     * 32H OF THE RECIPROCAL  LATTICE OF   ,I3)
921	FORMAT(12H IS PROBABLY )
922	FORMAT(12H   MIGHT BE  )
	END

	FUNCTION TETCAL(Q)
C--THIS FUNCTION CALCULATES THETA FROM Q
	COMMON /CONSTS/ RAD,TWORAD,RRAD,HLFLAM,WAVSQR,STORE(200,9)
	ARG=HLFLAM*SQRT(Q)
	IF (ARG .GT. 1.0) THEN
	   ARG=1.
	   PRINT 1,Q
	END IF
	TETCAL=ASIN(ARG)
	RETURN
1	FORMAT(' Q-VALUE TOO BIG FOR TETCAL, Q=',F10.2)
	END

	SUBROUTINE TEXT1(NN,HULP,NOAR,NCAR)
C  OUTPUT OF THE PROBABILITY FOR DIFFERENT FORMULATIONS OF THE ZONES
C  ONLY USED AS TEST OUTPUT
	DIMENSION  HULP(6),NOAR(6),NCAR(6)
	COMMON /LINES / IDOC,INP,IOUT,MAX,OBS(40),CEL(10)
	WRITE(IOUT,1)
	DO 50 J=1,6
	GO TO (11,12,13,14,15,16), J
11	WRITE(IOUT,21) NOAR(J),NCAR(J),HULP(J)
	GO TO 40
12	WRITE(IOUT,22) NOAR(J),NCAR(J),HULP(J)
	GO TO 40
13	WRITE(IOUT,23) NOAR(J),NCAR(J),HULP(J)
	GO TO 40
14	WRITE(IOUT,24) NOAR(J),NCAR(J),HULP(J)
	GO TO 40
15	WRITE(IOUT,25) NOAR(J),NCAR(J),HULP(J)
	GO TO 40
16	WRITE(IOUT,26) NOAR(J),NCAR(J),HULP(J)
40	IF (J .EQ.NN) WRITE(IOUT,2)
50	CONTINUE
	RETURN
1	FORMAT(51H0 NO  NC    FORMULATION OF ZONE            QUALITY /)
2	FORMAT(1H+,53X,13H  ** USED **  )
21	FORMAT(1H ,I3,I4,32H     A  AND  B  DOUBLED          , F12.2)
22	FORMAT(1H ,I3,I4,32H     B  DOUBLED                  , F12.2)
23	FORMAT(1H ,I3,I4,32H     A  DOUBLED                  , F12.2)
24	FORMAT(1H ,I3,I4,32H     UNCHANGED ZONE              , F12.2)
25	FORMAT(1H ,I3,I4,32H     A AND B DOUBLED, CENTERED   , F12.2)
26	FORMAT(1H ,I3,I4,32H     UNCHANGED ZONE,  CENTERED   , F12.2)
	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
	LIMH=1
	NR=1
	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
C      KG=KG-5
	KST=-KG
	KG=KG-5
C      NR=1
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 TIMER(N)
C  TIMER SUBROUTINES ARE INSTALLATION DEPENDENT. YOU MAY ENTER YOUR OWN
C  SUBROUTINE HERE.
	CALL ZTIMER(N)
	RETURN
	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
	WRITE(IOUT,9930) (CEL(J),J=1,6) ,VOLUM
	WRITE(20,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
	WRITE(20,9931)  N,DX,DOBS
	RETURN
1270	WRITE(IOUT,9932) PARTD
	WRITE(20,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 WRISP(NPOS,VAR,ICODE)
C  THE FUNCTION OF  WRISP  IS TO WRITE THE Q-VALUES OF THE CALCULATED
C  LINES IN A NEAT SCHEME, RECOGNISABLE FOR THOSE WHO  KNOW HOW TO
C  HANDLE THESE SCHEMES.
C     IT WRITES  VAR  (THE Q-VALUE) IN THE FORMAT F8.1, PRECEDED BY
C  (NPOS-1)*9 SPACES.   IF THE CALCULATED VALUE OF Q CORRESPONDS TO
C  AN OBSERVED VALUE, AN X IS ADDED LATER.
C
C  THE PARAMETERS OF WRISP ARE:
C  1. NPOS     = POSITION OF Q-VALUE
C  2. VAR      = VALUE OF Q
C  3. ICODE    = CODE TO TELL WRISP WHAT  TO DO
C  ICODE=0  PUT VAR IN POSITION
C       =1  PUT X IN POSITION.
C       =2  PRINT THE LINE WITH Q-VALUES.
C
C  A LINEAR ARRAY  *ITXT(133)*, REPRESENTING ONE LINE ON THE PRINTER,
C  IS FILLED WITH CHARACTERS.  WHEN THE LINE IS COMPLETELY READY,
C  IT IS PRINTED.
C
C  THIS PARTICULAR FORM OF WRISP TAKES A LITTLE BIT OF COMPUTING TIME
C  BUT IT IS KIND TO THE PRINTER.  IT SHOULD BE USED WHEN YOU HAVE A
C  RELATIVELY FAST COMPUTER AND A SLOW PRINTER (E.G. A SLOW TERMINAL
C  ON A FAST COMPUTER).
C
	CHARACTER*1 IBLANK,IX,IPNT,ISTAR
	CHARACTER*1 IFIGRS(10),ITXT(133)
	COMMON /LINES / IDOC,INP,IOUT,MAX,OBS(40),CEL(10)
	DATA IBLANK/' '/,IX/'X'/,IPNT/'.'/,IFIGRS(1)/'0'/,IFIGRS(2)/'1'/,
     * IFIGRS(3)/'2'/,IFIGRS(4)/'3'/,IFIGRS(5)/'4'/,IFIGRS(6)/'5'/,
     2 IFIGRS(7)/'6'/,IFIGRS(8)/'7'/,IFIGRS(9)/'8'/,IFIGRS(10)/'9'/
	DATA NUSED/0/,ISTAR/'*'/
	IF (NUSED .GT. 0) GO TO 20
	GO TO 400
20	NUSED=NUSED+1
	IF (ICODE-1) 100,200,300
C
100	IF (NPOS .LE. 14) GO TO 150
C  ONLY 132 PRINT POSITIONS AVAILABLE,  15*9=135]
	WRITE(IOUT,9990)ITXT
	NPOS=1
	DO 110 K=1,133
110	ITXT(K)=IBLANK
150	N=9*(NPOS-1)
	IF (VAR .GT. 0.1) GO TO 156
	DO 152 I=3,9
	K=N+I
152	ITXT( K )=ISTAR
	ITXT(N+6)=IFIGRS(1)
	GO TO 180
156	NR=IFIX(VAR*10.0+0.5)
	IREST=NR/10
	NFIG=NR-10*IREST
	ITXT(N+8)=IFIGRS(NFIG+1)
	ITXT(N+7)=IPNT
	NR=IREST
	DO 160 I=1,6
	INR=7-I
	IREST=NR/10
	NFIG=NR-10*IREST
	NR=IREST
	K=N+INR
	ITXT(  K  )=IFIGRS(NFIG+1)
	IF (NR .EQ. 0) GO TO 180
160	CONTINUE
180	RETURN
C
200	N=9*NPOS
	ITXT(N)=IX
	RETURN
C
300	WRITE(IOUT,9990)ITXT
	NPOS=1
C
400	ITXT(1)=IFIGRS(1)
	DO 410 K=2,133
410	ITXT(K)=IBLANK
	IF (NUSED .EQ. 0) GO TO 20
	RETURN
9990	FORMAT(133A1)
	END

	SUBROUTINE ZERCHK
C  THIS SUBROUTINE TRIES TO FIND SPECIMEN DISPLACEMENT ERRORS OR
C  ERRORS IN THE CALIBRATION OF THE ZERO POINT OF THE DATA.
C
	DIMENSION TTT(40),QPT(40),QMT(40)
	DIMENSION Q(40),NFT(21),DEL(21),NR(21)
	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
	WRITE(IOUT,1)
	WRITE(IOUT,2)
	HRAD=ATAN(1.)/90.
	SWAVEL=WAVEL/200.
	DO 20 M=1,MAX
20	TTT(M)=TWTHET(M)-0.11
	DO 100  I=1,21
	K=I-11
C     WRITE(IOUT, 5) K
	NR(I)=K
	DO 30 M=1,MAX
	TTT(M)=TTT(M)+0.01
	SIT=SIN(HRAD*TTT(M))/SWAVEL
	W=SIT*SIT
	Q(M)=W
	TTP=TTT(M)+TOLG
	SIT=SIN(HRAD*TTP)/SWAVEL
	QPT(M)=SIT*SIT
	QMT(M)=2.*W-QPT(M)
30	CONTINUE
C     WRITE(IOUT,6) (Q(M), Q(M+10), Q(M+20), Q(M+30), M=1,10)
	QMAX= QPT(MAX)+10.
	NFIT=0
	SUMDEL=0.
	DO 90 M=1,6
	IF (4.*Q(M) .GT. QMAX) GO TO 95
	DO 60 J=1,5
	II=J+1
	FIII=II*II
	QQ= FIII*Q(M)
	IF (QQ .GT. QMAX) GO TO 90
	DO 50 K=1,MAX
	IF (QPT(K) .LT. QQ) GO TO 50
	IF (QMT(K) .GT. QQ) GO TO 60
C  NOW QQ FITS A Q-VALUE.
	NFIT=NFIT+1
	SUMDEL=SUMDEL+ABS(QQ-Q(K))/FIII
C     WRITE(IOUT, 7) Q(M), II,  Q(K) ,SUMDEL
50	CONTINUE
60	CONTINUE
90	CONTINUE
95	DEL(I)=SUMDEL
	NFT(I)=NFIT
100	CONTINUE
	WRITE(IOUT, 8)
	WRITE(IOUT,3) NR,NFT,DEL
	WRITE(IOUT, 9)
	WRITE(IOUT, 9)
	N=0
	DO 200 I=1,21
	IF (NFT(I) .LE. N) GO TO 200
	M=I
	N=NFT(I)
	DIF=DEL(I)
200	CONTINUE
C     WRITE(IOUT,4)M,N, DIF
	RETURN
1	FORMAT(50H0  THE PURPOSE OF SUBROUTINE ZERCHK IS TO FIND OUT,
     2       50H WHETHER THERE ARE SPECIMEN DISPLACEMENT ERRORS OR/
     3       50H ZERO ERRORS IN THE DATA.  WE TRY TO CHECK THESE B,
     4       51HY LOOKING FOR HIGHER ORDERS OF THE LOW-ANGLE LINES./
     5       50H   THEREFORE WE APPLY A ZERO-CORRECTION TO ALL LIN,
     6       50HES AND FIND THE NUMBER OF HIGHER ORDERS.          /
     7       50H THE NUMBER OF HIGHER ORDERS IS PRINTED TOGETHER  ,
     8       50HWITH THEIR MEAN DEVIATION.                        /
     9       50H THE BEST ZERO CORRECTION IS PROBABLY THE ONE THAT,
     O       51H GIVES THE LARGEST NUMBER OF HIGHER ORDERS WITH THE/
     A       50H SMALLEST MEAN DEVIATION.  THE PROGRAM DOES NOT AP,
     B       50HPLY THIS CORRECTION AS IT IS NOT ALWAYS MEANINGFUL/
     C       50H BECAUSE THE AMOUNT OF USEFUL DATA IS SO SMALL.   )
2	FORMAT(50H0  IT IS ENTIRELY YOUR RESPONSABILITY TO DECIDE WH,
     2       50HETHER THE BEST ZERO CORRECTION IS A REASONABLE ONE/
     3       50H FOR YOUR EXPERIMENTAL CONDITIONS.                /
     4       50H YOU CAN APPLY YOUR CHOSEN CORRECTION BY ENTERING ,
     5       50HIT AS THE FIRST PARAMETER ON THE SECOND PARAMETER /
     6       50H CARD AND RERUNNING THE PROGRAM.                  )
3	FORMAT(15H0SHIFTS (X0.01), 21I5/ 15H NR OF FITS    , 21I5/
     2       15H MEAN DELTA Q  , 21F5.1 //
     3       56H0IF YOU  WANT TO APPLY A ZERO CORRECTION TO ALL TWOTHETA
     4 ,     55H VALUES, YOU CAN DO THIS ON THE SECOND PARAMETER CARD.
     5  /    56H THIS IS NOT NECESSARY IF YOU USED AN INTERNAL STANDARD
     6  ,    11HCORRECTLY.  )
C   4 FORMAT(1H ,2I3, F8.3)
C   5 FORMAT(23H0SHIFT IN HUNDREDTHS=  , I3)
C   6 FORMAT(    1H , 4F10.1  )
C   7 FORMAT(1H ,F10.1, I5,2F10.1)
8     FORMAT(56H0ALL INPUT TWO THETA VALUES WERE SHIFTED BY A NUMBER OF  
     2 ,     55HHUNDREDTHS OF A DEGREE TWO THETA, AS INDICATED BELOW.   
     3 /     56H THE NEW VALUES WERE TESTED ON THE OCCURRENCE OF SECOND- 
     4  ,    56H AND HIGHER ORDER LINES.  THE NUMBER OF FITS AND THE    
     5  /    56H MEAN DEVIATION ARE GIVEN BELOW.                        
     6  /    56H THE RESULTS MAY GIVE AN INDICATION OF A SPECIMEN DISPLA 
     7  ,    50HCEMENT CORRECTION OR A ZERO ERROR.                   )
9	FORMAT(56H0WARNING.  MESSING AROUND WITH MEDIOCRE DATA IS NO  (REP
     2  ,    56HEAT NO) SUBSTITUTE FOR GOOD EXPERIMENTAL METHODS.       
     3  /    56H ------------------------------------------------------- 
     4  ,    55H-------------------------------------------------------)
	END

	SUBROUTINE ZTIMER(N)
C  IBM TIMER
	COMMON /LINES / IDOC,INP,IOUT,MAX,OBS(40),CEL(10)
	DIMENSION TIMSP(10)
	CHARACTER*8  TEXT(10)
	DATA TEXT/' READER ',' FINZON ','EVALUATE','REAR ETC',' THREED ',
     * ' NORM3  ',' REFPAR ',' REFPAR2','   TOTAL',' RESERVE'/
	IF (N .GT. 0) GO TO 10
C  SETTIM STARTS  THE TIMER ON OUR MACHINE
	LTIM=0
	DO 5 I=1,10
5	TIMSP(I)=0
        LTIM = LIB$STAT_TIMER(2,LTIM)
	RETURN
10	CONTINUE
C     CALL ASKTIM(ITIME)
C  ASKTIM GIVES THE TIME IN HUNDREDTHS OF SECONDS SINCE THE LAST CALL
C  OF SETTIM
CVAX  ISTAT = LIB$STAT_TIMER(2,ITIME)
CVAX  TIME = FLOAT( ITIME ) / 100.0
CVAX  CPUSEC = TIME - START
        ITIME = LIB$STAT_TIMER(2,ITIME)
	TIME=0.01*ITIME
	TIM=0.01*(ITIME-LTIM)
	IF (N .EQ. 9) TIM=TIME
	WRITE(IOUT,99) TEXT(N),TIM,TIME
	LTIM=ITIME
	TIMSP(N)=TIM
	IF (N .EQ. 9) WRITE(IOUT,98) (TEXT(I),I=1,9),(TIMSP(I),I=1,9)
	RETURN
98	FORMAT(16H0TIME SPENT IN  ,9(A8,2X) / 10X, 9F10.2)
99	FORMAT(14H0TIME USED IN  ,1A8,F8.2,8H SECONDS,5X,6H TOTAL,F8.2)
	END
