Title:  		LATTPARM    A Personal Computer Program for finding
			            Lattice Parameters from Powder Data
KeyWords:  		Powder Diffraction, Lattice Parameters
Computer:        	PC, XT, AT
Operating System:     	MS(PC)-DOS
Programming Language:  	TURBO Pascal
Hardware:              	Math Co-Processor desirable
Author:                  Roy Garvey
Correspondence Address:  Department of Chemistry
                         North Dakota State University
                         Fargo, North Dakota
                         58105 - 5516     U S A

      Using the algorithm of Visser (1969), the problem of finding the
lattice parameters from powder diffraction data becomes one of finding the
six constants, A, B, C, D, E, and F where
Q(hkl)  =  10**4 / d**2  =  hh A + kk B + ll C + kl D + hl D + hk F       < 3 >
The unknown unit cell so derived will be triclinic.
      A perfectly general approach to the problem is to first find zones of
the lattice (net-planes of the reciprocal lattice containing the origin).
Any two such zones will have a line of intersection.  If the angle between
two such zones is found, a reciprocal lattice is determined.  Sometimes this
lattice can be reduced (described on symmetry axes in the reciprocal lattice).

Title:  		LATTPARM    A Personal Computer Program for finding
			            Lattice Parameters from Powder Data
KeyWords:  		Powder Diffraction, Lattice Parameters
Computer:        	PC, XT, AT
Operating System:     	MS(PC)-DOS
Programming Language:  	TURBO Pascal
Hardware:              	Math Co-Processor desirable
Author:                  Roy Garvey
Correspondence Address:  Department of Chemistry
                         North Dakota State University
                         Fargo, North Dakota
                         58105 - 5516     U S A
Documentation:
-----------------------------------------------------------------------------

      Using the algorithm of Visser (1969), the problem of finding the
lattice parameters from powder diffraction data is formulated with reference
to the reciprocal lattice:

(2sin(theta))2     1         ( h  k  l )  ( a*.a*  a*.b*  a*.c* ) ( h )
   --------    =  ---   =                 ( a*.b*  b*.b*  b*.c* ) ( k )  < 1 >
(   lambda  )2     d2                     ( a*.c*  b*.c*  c*.c* ) ( l )

which reduces to

                      1
                     ---  =  h2 a*2  +  k2 b*2  +  l2 c*2
                      d2
                                 +  2 k l b* c* cos(alpha*)
                                 +  2 h l a* c* cos(beta*)
                                 +  2 h k a* b* cos(gamma*)              < 2 >

Defining

  Q(hkl)  =  104 / d2  =  h2 A + k2 B + l2 C + kl D + hl D + hk F       < 3 >

(retaining the notation of Visser) the problem becomes one of finding the six
constants, A, B, C, D, E, and F.  The very fact that this relationship is
used to determine the lattice parameters means that the unknown unit cell so
derived will be triclinic.

      A perfectly general approach to the problem is to first find zones of
the lattice (net-planes of the reciprocal lattice containing the origin).
Any two such zones will have a line of intersection.  If the angle between
two such zones is found, a reciprocal lattice is determined.  Sometimes this
lattice can be reduced (described on symmetry axes in the reciprocal
lattice).

      The program LATTPARM performs the following steps:

      1)      Data entry.  The option to input data as d, Q, or 2-theta values
                is retained.  A mixed mode of input is not supported.
      2)      Test whether the base vectors should be halved.  Refine the
                parameters with a least-squares method and calculate the
                probability that the zone is found by pure chance (quality
                value).
      3)      Find pairs of zones with a common row and determine the angle
                between these zones.
      4)      Reduce the lattices found and transform if necessary so that
                the lattice is described in a standard way.
      5)      Index the first twenty lines of the pattern and repeat after a
                least-squares refinement of the parameters, note the number of
                lines actually indexed and calculate a "figure of merit".

                         F i n d i n g      Z o n e s

      Any two points in the reciprocal lattice, together with the origin,
define a plane in the reciprocal lattice corresponding to a crystallographic
zone.   Since

      ( 1 / d2hk0 )  =  h2 a*2  +  k2 b*2  +  2 h k a* b* cos(gamma*)   < 4 >

it follows that the equation of the plane has the form

                       Qhk  =  h2 Q'  +  k2 Q"  +  h k F                < 5 >

where                    F  =  2  / Q' Q"  cos(gamma*)                  < 6 >

gamma* being the angle between the directions from the origins to the points
with values Q' and Q".  Rearranging <5>

                             ( Q(hk)  -  h2 Q'  -  k2 Q" )
                     F  =   -------------------------------             < 7 >
                                         h k

Combinations of the first three with the first six lines of the Q list are
tried as Q' and Q" unless specific combinations are given in the input.

      By inserting for Q(hk) all observed Q values up to a reasonable limit,
inserting for h and k a few positive integer values, and storing the absolute
value of F, a great number of : F :  values are obtained.  Some of these
values are equal within the limits of error.

      In the program, F is multiplied by a constant factor, rounded to the
nearest integer value, and 1 is added to a corresponding ITL (tally of the
distribution) array element. If the integer exceeds the bounds of the ITL
array, an error message is placed on the screen and the point ignored. Lipson
and Steeple (1970) present a graphical interpretation of this procedure.

      When finished, the ITL array is processed in small overlapping steps,
summing the contents of six successive elements. If the sum exceeds a preset
number (e.g. 4 if Q' <> Q") then a weighted mean is calculated and stored. 
The tolerance in F has a fixed absolute value determined by the constant by
which F is multiplied, TOL1. 

      When the multiplicity-factor of a powder diffraction line is greater
than two, its Q value represents more than one significant point in the
reciprocal lattice.  Therefore it is useful to try combinations of Q values
of the type Q' = Q" (as illustrated in the expanded listing of intermediate
results "the Search for Zones commences").  In an orthorhombic lattice, for
example, each Q(hkl) with h,k,l <> 0 represents 8 points in the reciprocal
lattice.  These points define three different zones.  In the program up to
four significant values of F are accepted for the same pair Q', Q".  A
transformation of the zone axes is carried out for  F > Q' and if there are
equalities between Q', Q", and F.

      When R > Q' there is at least one reciprocal lattice point (i.e. 1, -1)
with Q = Q' + Q" - F less than that of one or both of the axes on which the
lattice was described.  In this case a transformation is carried out to
describe the zone on shorter axes.  This process is repeated as necessary.

      When there are equalities between Q', Q", and F there are pairs of Q's
which are systematically equal such that there is a mirror line in the zone.
The lattice is redefined on orthogonal axes.


                    Improvement  and  Evaluation  of  Zones

      The zones found to this point are recaluclated up to 0.8 Qmax with
Q'/4, Q"/4, and F/4 as parameters.  If as an examble the number of observed Q
values which match calculated Q values with odd h (and even k) is sufficient
then the program assumes Q'/4 is the correct parameter.  The same sort of
test is carried out for Q" (conditions:  odd k, even h) and for centering of
the old zone (conditions: h = 2n, k = 2n, h + k = 4n) and the new zone (h + k
= 2n).  With the parameters thus obtained a least-squares refinement is
carried out.  At the same time the "Zones after Evaluation" and the lines of
a zone "original values" can be printed.  A summary list of zones "Selected
for Refinement" is printed followed by a similar listing of the zone "after
refinement on two-theta".

      During this refinement process, the probability that a zone has been
found by pure chance is calculated.  Both the parameters of the zones and the
observed Q values will contain deviations from the ideal value.  When the 
difference between Qcalc and Qobs is less than TOL1 in Q = 104/d2, this near
coincidence is regarded as support for the zone.  Associated with each Qobs
in the list there is a certain range, /_ Q, of significant Q values.  Any
Qcalc falling within this range is considered to "fit".  If Qmax is the
highest observed Q value in the zone, then the probability, p, that an
arbitrary Q value fits is obtained as

                                     /_ Q
                             p  =   ------                             < 8 >
                                     Qmax

If there are Nc calculated Q values in the zone, No of which give a "fit"
then the propability that this should occur by pure chance is

                             Nc!
           chance   =  --------------  pNo ( 1 - p )Nc-No              < 9 > 
                       No! (Nc - No)!

where No and Nc should refer to significant points.

      A line, together with its higher orders is therefore counted as a
single observation, both for No and Nc.  Also the number of parameters
involved is deduced from the number of significant points to obtain No and
Nc.  The reciprocal value of c is assigned as a "quality index" of the xone.
The zones are sorted on this quality index.  Only the best six zones are used
in subsequent stages of processing by the program.
   Combination  of  Zones  and  Determination  of the  Angle  between  them

      In order to find a lattice, all possible combinations (including the
combinations of each zone with itself) of the six best zones (derived as
described above) are tried.  For each pair of zones the program seeks to find
the line of intersection and subsequently to determine the anble between
them.

      The line of intersection of two zones is a row of common points; the Q
value of the first point of this row must be a low Q value occurring in both
zones.  In order to find this point the Q values of the four points nearest
to the origin in each zone ( Q', Q", Q' + Q" + F, Q' + Q" - F) are
computed and compared with the points of the other zone.  If a common Q value
is found the zones are redefined (if necessary) with the common value as one
of the axes and the smallest remaining Q value as the other.  Provisionally
labeling the parameters of the first zone A, B and F and those of the second
zone A, C and E, equation < 3 > is used to determine the remaining parameter,
D:
                   ( Q - h2 A - k2 B - l2 C - hl E - hk F )
            D  =   ----------------------------------------            < 10 >
                                   k l

All observed Q values in the list are used with
           -2 < h < 2;   k = -2, -1, 1, 2;  l = 1, 2   (k l <> 0 !)
The actual calculation is carried out in much the same way as the selection
of values of F described in "Finding Zones" above.  Apart from the most
frequently occurring value of D, three other values are retained, but only if
they are not much less frequent that the best values.  The results are
described in lists of "Possible Combinations of Zones on which to base search"
and "Possible lattices based on the above Combination of Zones". 


                 R e d u c t i o n     of     L a t t i c e s

      Some of the sets of constants, A ... F, derived above might be
different descriptions of the same lattice.  Other lattices could be
described in a simpler way.  Therefore the lattices are reduced and brought
to a standard description.

      The reduction of the lattices is carried out in two steps.  In the
first step the shorest non-coplanar translations in the reciprocal lattice
are determined by calculating the Q values of half the number of reciprocal
lattice points close to the origin (the inversion-equivalent points are left
out), selecting the smallest Q values and storing their indices (hn, kn, ln,
n= 1,2,3).  The translations are coplanar if the determinant of the matrix

                                [ h1  k1  l1 ]
                                [ h2  k2  l2 ]    =  0
                                [ h3  k3  l3 ]

In this case the next smallest Q value is used as the third translation, etx.
until a non-coplanar set is found.  The second step in the reduction is a
test for symmetry.
      Only two kinds of symmetry relations are tested:

      a)      Two of the translational constants (A, B, C) are equal while
           the angle between their directions differs from 90o, e.g. A = B and
           F <> 0. In this case the reduced constants (primed) are:

              A' = ( A + B - F ) / 4             B' = ( A + B + F ) / 4
              C' =  C                            D' = ( D - E ) / 2
              E' = ( D + E ) / 2                 F' =  0

      b)      An angle factor is equal to one of the translational constants
           which defines one limit of the angle, e.g. F = A or F = B.
           If F = A the reduced constants are:

              A' = A / 4                         B' = B - A / 4
              C' = C                             D' = D - E / 2
              E' = E / 2                         F' = 0

      The third possible relation, e.g. D = E only leads to reduction if A = B
which is treated above, or if F = 2 A or F = 2 B, which is impossible when
the lattice has been defined by the shortest translations.  By cyclic
permutation of the constants they are all made to pass the few simple test
described above.  In this way all orthorhombic, and most of the monoclinic
lattices are recognized as such.  Only one extra test is applied in order to
pick out a few special cases of monoclinic centered lattices.  Both symmetry
relations point to a probable centering of the lattice.  This centering is,
however, disregarded at this point of the processing.  The reduction of the
lattice will establish the Bravais type which follows from the final
indexing.  The "Lattices Reduced and put into Standard Form" are listed.


   F i n a l     T e s t     o n     t h e     L a t t i c e s     F o u n d

      With the lattices found, the program tries to index the first 20 lines.
At the same time a least-squares refinement of the parameters is performed.
By counting the number of lines that can be indexed on the basis of a
centered lattice ( A, B, C, I, of F ) the Bravais-type is established.  This
procedure is repeated once, after which the number of non-indexable lines is
noted and a "figure of Merit" (de Wolff, 1968) is calculated.  This figure of
merit is in fact the ratio of the expected average discrepancy between the
observed and calculated Q values for an arbitrary reciprocal lattice of
the same size and the actual average discrepancy.  Parameters for the "best
lattices" are printed.   If requested the lattice may be printed in layer by
layer detail where the calculated Q values that match an observed value are
marked with an x.  It is suggested that this scheme makes any additional
extinction more easily identifiable by inspection.  In a final pass through
the refinement, the "best solution is used to index all lines".
                          I M P L E M E N T A T I O N

      The executable program for the derivation of Lattice Parameters from
Powder Diffraction data is distributed in two versions.  LATTPARM.com and
LATTPARM.000 are the main machine code program file and an Overlay file
respectively requiring the presence of the "math coprocessor" for proper 
execution.  UNITCELL.com and UNITCELL.000 are the corresponding files used
when no 8087 (or 80287) is avaialble.  Not only does UNITCELL require a
longer time for execution, but the smaller dynamic range of numeric values may
also affect actual data processing. 

      Because the user of a personal computer is often impatient when longer
computations are taking place, a visually "active screen" has been
incorporated in the latest version of the program.  In addition to rather
conventional "Program Identification", a running display of "Progress" in
processing of the data is presented along the lower line of the display
screen.  In addition, certain phases of the program present additional
information on typical data being processed.  Most screens display A, B, F,
etc, values with no identifying labels.  Should problems leading to "Run Time
Errors" cause the program to halt, a "PrtSc" dump of the screen to the
printer can be helpful in tracing the source of the error.  The following
out-line of the Visser algorithm as implemented might prove useful in
following the "Progress Line".

            LATTICE  PARAMETERS   from   POWDER  DIFFRACTION  DATA

Initiate Visser Algorithm.
      Name of System.
      Printer Output File Selection.
      Printer Output File Initialization.

Read in Control and Diffraction Data.
      Open Input Data File.
      Read Control Parameters.
      Display Screen and Print Report Headers.
      Check and Adjust Input Control Data.
      Read Diffraction Lines Data.
              Print the Observed Values.
              Compute  Q  Values.
              Print  Q  Values.
              Compute Tolerances for Input Data.
      Read Line and Lattice Combinations (if any).
      Check for systematic Zero-Shift (if requested).

Crystallographic Zone Finding.
      Initiate Probabilities.
      Search for Zones.
      Find Zones.
      Evaluate Zones.
              Normalize Zone.
              Least-Squares for Theta.
              Least_Squares for Lattice Vectors.
      Combination of Zones.

Crystallographic Lattice Finding
      Lattice Search.
              Intersect Zones.
                   TriClinic Lattice.
                   MonoClinic Lattice.
              Find Most Frequent Values.
              Best Values.
      Normalize Lattice.
      Refine Parameters.
              Generate Q Values.
              Bravais Lattice Type.
      Discard Lattices.
      Most Probable Solution.


                              Program  Initiation

      The program receives information on the InPut data file (and optionally
on the printer OutPut file) through entries on the command line.  Following
the DOS system prompt > the entry of a name matching a .COM or .EXE program
file will cause the program to be loaded into memory and execution will
begin.  If no matching program file is found, the DOS system looks for a .BAT
file of commands.  Under predetermined conditions, information typed
following these "Commands" as a "tail" can be read by the program and used
for a variety of purposes.  The execution of Lattice Parameter searches
requires that the necessary control and diffraction data be included in a
TEXT file stored on disk.  No interactive data entry is expected.

      To process the diffraction data for a file named  YOUR.dat  residing on
a diskette in drive  B  you would type the information following > assuming
that your program is in the default directory.

              default directory > LATTPARM  B:YOUR.dat

The printed output would be sent to the "device" referenced to LST:, LPT1:,
or PRN: depending upon the whim of your DOS environment.  The rather extensive
printed output can slow processing because of a slow printer.  If you have a
"spooler" or ram buffer this makes little difference.  However, as an
alternative, the program includes the option of sending the entire "printed
output, including control characters" to a disk file for later printing
and/or input to a word processor for incorporation into a report.  This is
accomplished by the following "command tail".

              default directory > LATTPARM  B:YOUR.dat  C:\REPORT\YOUR.out

C:\REPORT\  identifies that YOUR.out file will be placed in a REPORT
subdirectory found on drive C:.


                             Printer Control Codes

      An additional feature was introduced to lessen the impact wrought by
the multifarious range of printer options and control codes.  A very limited
subset of functions was selected for implementation:
              Compressed mode printing           15 - 17 cpi
              Double width printing              2 x current width
              Normal mode printing               10 or 12 cpi
The program looks on the "default drive" for a  PRINTER.$pr  file of three
lines of up to five decimal number equivalents to the CHR( n ) ASCII codes.
If found, the file is read for the codes to set your printer in the three
modes described one line each in the order listed.  On some printers this is
sufficient to include line spacing codes as well.  For example, the codes
corresponding to 15 27 48 when sent to an IBM-Epson printer will "change the
printer to compressed print mode" and "change the paper feed to 1/8 inch".

      If the PRINTER.$pr file is not found, the program forces you to select
between a standard Okidata microline and an Epson type printer.  The
corresponding printer codes are generated internal to the program.

                      D A T A      I N P U T      F i l e

      LATTPARM (or UNITCELL) reads an ASCII TEXT file of input records from
disk.  The usual input should consist of:

1.    A Title record containing up to 80 characters.

2.    A parameter record including the following.  All variables must be given
        appropriate values.  The numbers are all "place holders" and are read
        format free except for the requirement that at least one space
        separate adjacent numerical values.  Zeros must be entered.

      MaxSoln  --  The maximum number of solutions for which a list of
         (indexed) reflections is printed.  The corresponding Q-scheme is only
         printed when intermediate results are requested ( ENL below) and the
         figure of merit is larger than 4.0.   4 is a reasonable value.

      Nsyst[1]  --  Orthorhombic      Used to indicate that the solution must
      Nsyst[2]  --  Monoclinic           belong to a certain crystal system.
      Nsust[3]  --  Triclinic      =  0  means your are indifferent.
                                   = +1  indicates you are confident that the
                                            solution belongs to this system.
                                   = -1  means the opposite.
            This option works only on the last part of the data processing and
           doubles the tolerances there.  Only the three crystal systems are
           included.  Generally the 0 is used.

      Tol2  --  Tolerance range for 2-dimensional search.  3.0 Q-units is
                  generally used.  When you have very good data (0.01 degree
                  two-theta) on a fairly large unit cell, you might lower it
                  to 2.5 or 2.0. With a small unit cell you might increase it
                  somewhat to say 4.0.

      Tol3  --  Tolerance range for 3-dimensional search.  4.5 Q-units is
                  generally satisfactory.  Tol3 is the same kind of tolerance,
                  but now on the trial value of  D  in the search for
                  lattices. Its value can be changed in the same manner as
                  that described for Tol2.

      Lambda  --  X-ray wavelength used in recording the diffraction pattern.
                    Values are given in Angstrom units.

      LinCombo  --  Generally = 0.  For value greater than zero, the value
                    represents the number of line-combination records (or
                    zones) you will enter. When non-zero, the program expects
                    records (after the lines records + "-1 record").  See
                    below for more information.

      Lzerck  --  > 0  the program checks on "zero-errors" .

      Enl  --  Program prints "intermediate" results for a non-zero value.

      Nq1  --  Determine the numbers of lines to be used in combinations
      Nq2  --    for finding zones.  Values of 3 and 6 respectively are good.

      Nz1  --  Determine the numbers of zones to be used in combinations to
      Nz2  --    find lattices.  Values of 6 and 6 respectively are good.

      Nr  --  Generally zero.  For positive non-zero value, the program
                 expects records for Nr trial lattices (entered after the
                 lines records + "-1 record").  See below for further details.

      List  --  The number of "best solutions" to be printed.  Among other
                 information, the program will produce a list of all observed
                 and calcualted lines, sorted on two-theta, unless it would be
                 too long (more than 100 calculated lines).  Generally List is
                 set at 1.

      MolMass  --  If the formula mol mass is entered, calculated density and
                     number of units in the cell are printed.  Enter zero if
                     no experimental value is available.

      Density  --  Observed density.  Used with volume of calculated unit
                     cell and MolMass to estimate  Z.  If no value is
                     available, enter zero.

      Tolg  --  Tolerance on the match between the calculated and the
                  observed two-theta values.  A value of 6 is generally
                  acceptable. Values are expressed in hundredths of a degree
                  two-theta and apply to both the zone and lattice searches
                  and the lattice refinement. For very accurate data it can be
                  advantageous to use a value of 4.0 to 3.5. For inaccurate
                  data the value can be increased, but such data is not
                  usually worth processing unless the cell is very small.


3.    A second parameter record with the numeric values for:

      ZeroCorr  --  Zero correction (as degrees two-theta) to be applied to
                      all lines.

      PrintMerit  --  The minimum figure of merit for a lattice to be
                        printed.  The value of 4.0 must be considered as an
                        absolute minimum for a lattice that has a chance to
                        imporve its performance.

      PrintLine  --  The minimum number of indexed lines for a lattice to be
                       printed.  Lattices with a lower number of indexed lines
                       are sometimes a sub-lattice of the true lattice.  A
                       lattice can only be considered as satisfactory when it
                       explains all lines. A value of 14 is considred minimal.
                       Note that this parameter is entered with no decimal
                       point.


4.     Diffraction Lines.

            Observed diffraction lines, one to a line.
              May be input as "D-vlaues", "Q-values" or as "two-theta" values.

5.    -1       used to indicate end of diffraction lines data.

6.    Trial Lattices can be input at this point.  There must be six fields of
        data for each record.  The first field of the last record must be -1.
        The six data fields of a record correspond to the coefficients A, B,
        C, D, E, and F of

            Q  =  1000 / ( d * d )
               =  H*H*A  +  K*K*B  +  L*L*C  +  K*L*D  +  L*H*E  +  H*K*F

             where  d  is the interplanar spacing in Angstroms.

7.    Line combinations or Zone input.  Each record must have four entries.
        You may enter any number of zones (less than 40), but must end the set
        with a -1 for the first of the four values.  When you enter only two
        non-zero values (separated from adjacent numbers by at least one blank
        space) the program interprets this as a combination of lines to be
        tried for finding zones.  When you enter three non zero values, the
        program will treat these as a zone and will evaluate it.  If you want
        to enter a rectangular zone, you must set this third constant at 0.08,
        otherwise the program will treat this as a line combination.  If you
        want to ensure that a specified zone will be included in the search
        for complete lattices, you should also enter a high quality figure
        (fourth value of the record).  If the program does not find the
        solution automatically, it can often be helped by adding zones in this
        way.  Zones which are correct, although having a low quality figure,
        can sometimes be recognized in the listing of intermediate results
        (ENL <> 0) by their small area ( i e, large *A* and *B*) and good
        coverage near the origin (i e, most of the low angle calculated lines
        have been matched with observed lines). Zones from oriented specimens
        or fibers, or from electron diffraction studies, can also be used.


                               Acknowledgements

      The author wishes to thank JCPDS-ICDD for financial support for the
initiation of the project leading to generation of the Lattice Parameter
programme.  Special acknowledgement is accorded Dr G McCarthy, his students
who used earlier versions of the program, and especially to Paul Maistrelli
whose kind comments and criticisms guided evolution of the "people interface"
to the computer programme and prompted generation of this document.

                              R E F E R E N C E S

Azaroff, L V and M J Buerger (1958) "The Powder Method in X-ray
     Crystallography", McGraw-Hill Book Company, Chapter 10.

Lipsom, H and H Steeple (1971) "Interpretation of X-ray Powder Diffraction
     Patterns", MacMillan, Chapter 5.

Visser, J W (1969) "A fully automated program for finding the unit cell from
     powder data", J Appl Cryst, 2, 89-95.

de Wolff, P M (1968), J appl Cryst, 1, 108.


****   Sample DATA input file:

PDF  29-0733  Ilmenite test data set.
4  0  0  0  3.000  4.500  1.54060  0  0  0  3  6  6  6  0  1  0.0  0.0  6.0
0.000  4.0  14
   3.7370
   2.7540
   2.5440
   2.3490
   2.2370
   2.1772
   2.1032
   1.8683
   1.8309
   1.7261
   1.6535
   1.6354
   1.6206
   1.5057
   1.4686
   1.4342
   1.3757
   1.3421
   1.3337
   1.2834
   1.2719
   1.2453
   1.2279
   1.2101
   1.2040
   1.1871
   1.1744
   1.1547
   1.1185
   1.0884
   1.0758
   1.0663
   1.0516
   1.0156
   1.0042
   0.9872
   0.9813
   0.9719
   0.9617
   0.9421
   -1
END

NOTE:   The  LATTPARM.scr  file must be decreased in size to be compatible
 with the memory segment limitation of Turbo Pascal 2. and 3. (I do not yet
 have a copy of version 4. so I do not know of the limitations).  Using the
 non-document mode of your favorite word processor, <block> and save the last
 60K of the file as  LATTPAR2.scr.  Delete this block from the active 
 workspace.  Next add the line
                        {$I LATTPAR2.scr }
 as the last line of the active workspace and block the remaining first 50K of
 source statements.  Save this BLOCK as the file LATTPARM.pas.  Terminate your
 wordProcessor program [but do not save the changes to disk].

      Initialize TURBO Pascal.  Make the working file  LATTPARM.pas  and
 compile your version of LATTPARM.  Section headers in the source code
 suggest segmentation of the program used in its development.  Each segment
 was stored in an individual file and these were {$I  included  } into a 
 master supervisor file.

-----------------------------------------------------------------------------
End of File
Title:  		LATTPARM    A Personal Computer Program for finding  
		                      Lattice Parameters from Powder Data 
KeyWords:  		Powder Diffraction, Lattice Parameters            
Computer:         	PC, XT, AT                                  
Operating System:     	MS(PC)-DOS                              
Programming Language:  	TURBO Pascal                           
HardWAre:              	Math Co-Processor desirable            
Author:                  Roy Garvey                           
Correspondence Address:  Department of Chemistry              
                         North Dakota State University        
                         Fargo, North Dakota                  
                         58105 - 5516     U S A               

Source Code:
-----------------------------------------------------------------------------

PROGRAM   Unit_Cell_from_Powder_Pattern;
   CONST
      PgmTitle : STRING[54] =
                    'LATTICE  PARAMETERS   from   POWDER DIFFRACTION DATA';
      Version :  STRING[30] = 'EMMPDL      October, 1987';

   VAR                          {  Common to the three  OVERLAY  Programs  }
      Answer                                              : CHAR;
      Enl, Soln, Self, Rectangl                           : BOOLEAN;
      JobTitle, ScreenLine                                : STRING[80];
      Message                                             : STRING[40];
      Location                                            : STRING[20];
      MaxSoln, MaxObsLns, Lincombo, NrMax, Nroster,
        Nz1, Nz2, Nq1, Nq2,          {  Zone run,  Q run  }
        Nrun, Nindex, Norder, Nr, List, Prnt,
        Lzerck, PrintLn, dx, cx, HoldX, HoldY             :  INTEGER;
      Tol2, Tol3, Tolg,              {  tolerance for Zone, Lattice, line   }
        Lambda, MolMass, Density,
        Fminim, Fmineq, Rrad,
        TolRad, QcalMx, Volume, Theta, TolKeep,
        ZeroCorr, ZerTot, PrintMerit, WaveSqr,
        TwoRadian, HalfLambda, Rad                        :  REAL;
      Obs, TwoTheta, ObsTheta,                   {  Observed values          }
        Obspt, Obsmt, ThetaPt, ThetaMt,          {  value + and - tolerance  }
        Prob, ShfTheta, UpBndr      :  ARRAY[1..40]        of REAL;
      Roster                        :  ARRAY[1..150,1..12] of REAL;
      Store                         :  ARRAY[1..200,1.. 9] of REAL;
      Cel                           :  ARRAY[1..10]        of REAL;
      ArrayS                        :  ARRAY[1..64]        of REAL;
      ArrayD                        :  ARRAY[1..8,1..8]    of REAL;
      Eval                          :  ARRAY[-8..8,0..8]   of REAL;
      Index                         :  ARRAY[1..40]     of INTEGER;
      Nsyst                         :  ARRAY[1..3]      of INTEGER;
      ChrWdt                        :  ARRAY[1..3,1..3] of STRING[5];
      Lstr                          :  TEXT;

   TYPE
      DateStr = STRING[10];
      TimeStr = STRING[5];
      RegPack = RECORD
                ax, bx, cx, dx, bp, si, ds, es, flags  : INTEGER;   END;

                {  spaces, time, data, pgmHeader   }

  FUNCTION  Spaces( Nbr : INTEGER ) : CHAR;
     VAR
        I  :  INTEGER;
     BEGIN
          FOR I := 1 to Nbr-1 DO   WRITE( Lstr, ' ' );    Spaces := ' ';
     END;    {   Spaces to  printer  }

  FUNCTION Date: DateStr;
    VAR
       RecPack     :   regpack;                { record for MsDos call }
       Month, Day  :   STRING[2];
       Year        :   STRING[4];

    BEGIN
       WITH recpack DO    ax := $2a shl 8;
       MsDos( recpack );
       WITH recpack DO  BEGIN
          STR( cx, year );          STR( dx mod 256,day );
          STR( dx shr 8,month );    END;
       IF Length(Month) = 1 THEN Month := ' ' + Month;
       IF Length(Day)   = 1 THEN Day   := ' ' + Day;
       date := Year + '/' + month + '/' + day;
     END;    {   D a t e   }

  FUNCTION Time: TimeStr;
     VAR
       recpack          :        regpack;           {record for MsDos call}
       Hours, Minutes:         STRING[2];

     BEGIN
        WITH recpack DO    ax := $2c shl 8;
        MsDos(recpack);                         { call function }
        WITH recpack DO  BEGIN
           STR( cx MOD 256, Minutes);    STR( cx shr 8, Hours);   END;
        IF Length(Hours)   = 1 THEN Hours   := ' ' + Hours;
        IF Length(Minutes) = 1 THEN Minutes := '0' + Minutes;
        Time := Hours + ':' + Minutes;
     END;              {  T i m e  }

  FUNCTION  Response : CHAR;
     VAR
        HoldLine  :  INTEGER;
        Respond   :  CHAR;
     BEGIN
        HoldLine := WhereY;      GoToXY( 1,25);  ClrEol;
        HighVideo;  WRITE( Message );  LowVideo;   GoToXY(50,25);
        READ( Kbd, Respond );      Response := UPCASE( Respond );
        GoToXY(1,HoldLine);
     END;                       {   Response  Character  Input   }

  PROCEDURE ScreenWrite( Code : BYTE );
     VAR
       Line, X  :  INTEGER;
     BEGIN
        Line := Length( ScreenLine );
        CASE  Code  OF
           0, 10 :  X := 1;                              { Left Justify  }
           1, 11 :  X := Round( (79-line)/2 );           { Centered text }
           1, 12 :  X := 79 - Line;       END;           { Right Justify }
        GoToXY( X, WhereY );  WRITELN( ScreenLine );
        IF Code > 9 THEN  WRITELN( Lstr, Spaces( x ), ScreenLine );
     END;                {   Place text line on Screen  with justification }

  PROCEDURE  Programme_Identification_Header;
     BEGIN
        ClrScr;           ScreenLine := PgmTitle;            ScreenWrite( 1 );
        WRITELN;          Screenline := 'Roy G. Garvey';     ScreenWrite( 1 );
        Screenline := 'Department of Chemistry';             ScreenWrite( 1 );
        Screenline := 'North Dakota State University';       ScreenWrite( 1 );
        Screenline := 'Fargo,  ND   58105-5516';             ScreenWrite( 1 );
        WRITELN;        ScreenLine := Version;               ScreenWrite( 1 );
        WRITELN;
     END;                 {  Programme Identification   Screen display  }

  PROCEDURE  No_PRINTER_File;
    VAR
      Answer  :  CHAR;
   BEGIN
     Sound( 525 );    Delay( 750 );    NoSound;
     ScreenLine := 'Initialization file   PRINTER.$pr   not Found';
     ScreenWrite( 1 );    WRITELN;
     WRITELN('Check the Printer !   O   >Okidata Version.');
     WRITELN('                      E   >Epson   Version.');
     Message := '         Your Selection: ';
     ChrWdt[1,1] := Chr( 29 );   ChrWdt[2,1] := Chr( 15 );  { Compressed }
     ChrWdt[1,2] := Chr( 31 );   ChrWdt[2,2] := Chr( 14 );  { Expanded   }
     ChrWdt[1,3] := Chr( 24 );   ChrWdt[2,3] := Chr( 18 );  { Normal     }
     REPEAT Answer := Response UNTIL Answer IN ['O', 'E'];  WRITELN;
     CASE  Answer OF
          'O' :  BEGIN  WRITE(Lstr, #29, #27, #56);  Prnt := 1;  END;
          'E' :  BEGIN  WRITE(Lstr, #15, #27, #48);  Prnt := 2;  END;
     END;  END;       {  No  PRINTER.$pr  file found on default drive  }

  PROCEDURE  PRINTER_File( VAR  PrnFile : TEXT );
     VAR
        I, J, Num  :  INTEGER;
     BEGIN
        FOR I := 1 to 3 DO BEGIN              {   Compressed   }
           ChrWdt[3,i] := '';                 {   Expanded     }
           FOR J := 1 to 5 DO BEGIN
              READ( PrnFile, Num );
              ChrWdt[3,i] := ChrWdt[3,i] + Chr( Num );   END;
           READLN( PrnFile );   END;
        CLOSE( PrnFile );   Prnt := 3;
     END;    {  PRINTER.$pr  file  found on default disk drive  }

  PROCEDURE  Printer_Initialization;
    VAR
       PrnFile :  TEXT;
    BEGIN
       ASSIGN( PrnFile, 'PRINTER.$pr' );
       {$I-}  RESET( PrnFile );   {$I+}
       IF IOResult <> 0  THEN  No_PRINTER_File
                         ELSE     PRINTER_File( PrnFile );
    END;                   {  Initialization of Printer  }

 PROCEDURE  Printed_Output_File_Selection;
    BEGIN
       Message := ParamStr( 2 );
       IF Message = '' THEN  Message := 'LST:';
       ASSIGN( Lstr, Message );
       IF Message = 'LST:' THEN BEGIN
          WRITELN('  OutPut records sent directly to the Printer Device.');
          RESET( Lstr );   END  ELSE  BEGIN
          WRITELN('  OutPut Records sent to file  ', Message );
          REWRITE( Lstr );   END;
    END;          {    Command Tail processing for Printer  output  }

 PROCEDURE  Progress_Line;
     VAR
        Isp   :  INTEGER;
     BEGIN
        HoldY := WhereY;     HoldX := HoldX + 1;            WRITELN;
        GoToXY( 1,24);       HoldX := HoldX MOD 35;         ClrEol;
        WRITE('Time   '); HIGHVIDEO;  WRITE(Time,'  '); LOWVIDEO;
        FOR Isp := 1 to HoldX DO   WRITE('$');        GoToXY(55,24);
        HIGHVIDEO;      WRITE( '@ - ', Location:20);       LOWVIDEO;
        IF HoldY > 20 THEN BEGIN   FOR Isp := 4 to 23 DO BEGIN
             GoToXY(1,isp);   ClrEol;   END;   GoToXY(1,4);  END
           ELSE GoToXY(1,HoldY+1);
     END;                         {   Progress  Line  }

PROCEDURE  Printed_Report_Header;
  BEGIN
   WRITELN(Lstr,'Date -  ',date,'  - - -   Time  [24:00] clock -    ', Time);
   WRITELN(Lstr, ChrWdt[Prnt,2] );     WRITELN(Lstr, PgmTitle );
   WRITELN(Lstr, ChrWdt[Prnt,3] );     WRITELN(Lstr, JobTitle );
   WRITELN(Lstr );     HoldX := 0;
  END;              {    Printed  Report  Header  }

 PROCEDURE  Screen_Display_Header;
    BEGIN
       ScreenLine := 'Date -  ' + date + '  - - -   Time  [24:00] clock -    '
                                + Time;         ClrScr;  ScreenWrite( 1 );
       ScreenLine := PgmTitle;   ScreenWrite( 1 );
       ScreenLine := JobTitle;   ScreenWrite( 1 );   WRITELN;
    END;                {  Screen  Display  Header  }

           {  initialization, wrispec, theta claculate  }

   FUNCTION  Theta_Calculate ( Q : REAL ) : REAL;
     VAR
        Arg  :  REAL;
     BEGIN
        IF    Q > 1E-4 THEN  Arg := HalfLambda * Sqrt( Q ) ELSE  Arg := 0;
        IF  Arg > 1.0  THEN  Arg := 1.0;
        Theta_Calculate := ArcTan( Arg/ Sqrt(1- Sqr(Arg)) );
     END;                    {   Calculate Theta from Q value   }

  PROCEDURE  Name_Of_System;
     BEGIN
        Programme_Identification_Header;        WRITELN;
        ScreenLine := 'Based upon the algorithm of Visser';  ScreenWrite( 1 );
        ScreenLine := 'J Appl Crystal (1969)  2, 89';        ScreenWrite( 1 );
        WRITELN;
     END;                 {  Name of system Screen display  }

 PROCEDURE  WriSp( VAR Position :  INTEGER; Value : REAL;  Code : BYTE );

   PROCEDURE  Print_Line;
      BEGIN
         WRITELN( Lstr );        Position := 1;
      END;        {  Print the Line  }

   PROCEDURE  Put_Variable;
      BEGIN
         IF position > 12 THEN Print_Line;
         WRITE(Lstr, Value:10:1 );
      END;        {   Put the Variable  }

     BEGIN                 {  W r i t e    S p e c i a l  }
        CASE Code OF
             0 :   Put_Variable;
             1 :   WRITE( Lstr,'x' );
             2 :   Print_Line;        END;
     END;                {   Write Special Format of Q Values  }

 PROCEDURE  Initiate_Visser_Algorithm;
    BEGIN
       TEXTMODE( bw80 );
       Nrun  := 0;     Nr := 0;      Nroster := 150;
       NrMax := 200;   Norder := 4;
       FILLCHAR( Roster,    SIZEOF( Roster ),   0 );
       FILLCHAR( Store,     SIZEOF( Store  ),   0 );
       FILLCHAR( Eval,      SIZEOF( Eval   ),   0 );
       FILLCHAR( ArrayS,    SIZEOF( ArrayS ),   0 );
       FILLCHAR( Index,     SIZEOF( Index  ),   0 );
       FILLCHAR( ArrayD,    SIZEOF( ArrayD ),   0 );
       FILLCHAR( Obs,       SIZEOF( Obs    ),   0 );
       FILLCHAR( Obspt,     SIZEOF( Obspt  ),   0 );
       FILLCHAR( Prob,      SIZEOF( Prob   ),   0 );
       FILLCHAR( Obsmt,     SIZEOF( Obsmt  ),   0 );
       FILLCHAR( Upbndr,    SIZEOF( Upbndr ),   0 );
       FILLCHAR( Thetapt,   SIZEOF( Thetapt ),  0 );
       FILLCHAR( Thetamt,   SIZEOF( Thetamt ),  0 );
       FILLCHAR( TwoTheta,  SIZEOF( TwoTheta ), 0 );
       FILLCHAR( ObsTheta,  SIZEOF( ObsTheta ), 0 );
       FILLCHAR( ShfTheta,  SIZEOF( ShfTheta ), 0 );
       Name_Of_System;
       Location := 'Initiation.';     Progress_Line;
       Printed_Output_File_Selection;
       Printer_Initialization;
    END;       {   Initialize Variables In  Main   }

                 {  matrix inversion procedure   }

PROCEDURE  Arrange_Array_Single( I, M : INTEGER );
    VAR
       J, L, Im  :  INTEGER;
    BEGIN
       Im := 0;
       FOR J := 1 to I DO FOR L := 1 to I DO  BEGIN
           Im := Im + 1;   ArrayS[im] := ArrayD[j,l];  END;
    END;     {  Arrange Array in Single Format  }

PROCEDURE  Arrange_Array_Double( I, M : INTEGER );
    VAR
       J, L, Im  :  INTEGER;
    BEGIN
       Im := 0;
       FOR J := 1 to I DO FOR L := 1 to I DO  BEGIN
          Im := Im + 1;   ArrayD[j,l] := ArrayS[im];   END;
    END;     {  Arrange Array in Double Dimension Format  }

 PROCEDURE  MatInvert( N : INTEGER; VAR  D : REAL );
     VAR
        Nk, K, Iz, Ij, j, Ki, Jp, Ji, Jk, Ik, Kj,
          Kk, Jq, Jr, I                          :  INTEGER;
        BigA, Hold                               :  REAL;
        L, M                                     :  ARRAY[1..8] of BYTE;

 PROCEDURE  Switch_Them( I1, I2 : INTEGER );
     BEGIN
         Hold := - ArrayS[i1];  ArrayS[i1] := ArrayS[i2];  ArrayS[i2] := Hold;
     END;    {  Switch them  }

 PROCEDURE  Reorder_Array_Elements;
     BEGIN
        J := L[k];
        IF J > K THEN BEGIN
           Ki := K - N;
           FOR I := 1 to N DO BEGIN
               Ki := Ki + N;   Ji := Ki - K + J;
               Switch_Them( ki, ji );  END;  END;
       I := M[k];
       IF I > K THEN BEGIN
          Jp := N * (I - 1);
          FOR J := 1 to N DO BEGIN
              Jk := Nk + J;   Ji := Jp + J;
              Switch_them( jk,ji );  END;  END;
     END;    {   Reorder  Array  Elements   }

 PROCEDURE  Initialize_Working_Variables;
     BEGIN
         Nk := Nk + N;  L[k] := k;  M[k] := k;  Kk := Nk + K;
         BigA := ArrayS[kk];
         FOR J := K to N DO BEGIN
             Iz := N * (J - 1);
             FOR I := K to N DO BEGIN
                Ij := Iz + I;
                IF ABS( BigA ) < ABS( ArrayS[ij] ) THEN BEGIN
                   BigA := ArrayS[ij];  L[k] := I;  M[k] := j;  END;
         END;  END;
     END;   {  Initialize  Working  Variables  }

 PROCEDURE  Reduce_Array_Elements;
     BEGIN
        FOR I := 1 to N DO BEGIN
            IF I <> K THEN BEGIN
               Ik := Nk + I;  ArrayS[ik] := ArrayS[ik] / - BigA;   END; END;
        FOR I := 1 to N DO BEGIN
            Ik := Nk + I;  Hold := ArrayS[ik];  Ij := I - N;
            FOR J := 1 to N DO BEGIN            Ij := Ij + N;
                IF I <> K THEN IF J <> K THEN BEGIN
                   Kj := Ij - I + K;
                   ArrayS[ij] := Hold * ArrayS[kj] + ArrayS[ij]; END; END; END;
        Kj := K - N;
        FOR J := 1 to N DO BEGIN            Kj := Kj + N;
            IF J <> K THEN ArrayS[kj] := ArrayS[kj] / BigA;   END;
        D := D * BigA;   ArrayS[kk] := 1 / BigA;
     END;    {   Reduce  Array  Elements    }

 PROCEDURE  Switch_them_2 (J1, J2 :  INTEGER );
     BEGIN
         Hold := ArrayS[j1];  ArrayS[j1] := - ArrayS[j2]; ArrayS[j2] := Hold;
     END;    {  Switch Them 2  }

 PROCEDURE  Reorder_Array_2;
     BEGIN
        Jq := N * (K-1);   Jr := N * (I-1);
        FOR J := 1 to N DO BEGIN
            Jk := Jq + j;   Ji := Jr + J;   Switch_Them_2( jk, ji );   END;
     END;    {  Reorder  Array   2  }

 PROCEDURE  Reorder_Array_3;
     BEGIN
        Ki := K - N;
        FOR I := 1 to N DO BEGIN
            Ki := Ki + N;    Ji := Ki - K + j;  Switch_Them_2( ki, ji );  END;
     END;   {   Reorder  Array  3   }

     LABEL  100, Return;
     BEGIN                  {    M a t r i x     I n v e r s i o n   }
        D := 1;  Nk := - N;
        FOR K := 1 to N DO BEGIN
            Initialize_Working_Variables;
            Reorder_Array_Elements;
           IF ABS( BigA ) <= 0.0001 THEN BEGIN  D := 0;  GOTO Return;   END;
           Reduce_Array_Elements;        END;
        K := N;
 100:     K := K - 1;    IF K <= 0 THEN                  GoTo Return;
          I := L[k];
          IF I > K THEN  Reorder_Array_2;
          J := M[k];
          IF J <= K THEN                                 GOTO 100;
          Reorder_Array_3;                               GOTO 100;
 Return:  END;  {  Matrix Inversion from IBM Scientific Subroutine  }

  OVERLAY PROCEDURE  Read_in_the_Control_and_Diffraction_Data;

          {  reads input control and lines data from file   }

     VAR
        DataFile        :  TEXT;
        D, B, Sit, Z    :  REAL;
        Ix, Kr          :  BYTE;

  PROCEDURE  Open_Input_File( VAR InFile : TEXT );
     VAR
        FileName  :  STRING[50];
     BEGIN
        FileName := ParamStr( 1 );
        ASSIGN( InFile, FileName );
        {$I-}     RESET( InFile );    {$I+}
        IF IOResult <> 0 THEN BEGIN
           WRITELN( #7, #7, '     InFile    ->   ',
                    FileName, #7, #7,
                  '    is missing.', #7, #7 );
           WRITELN( #7 );    HALT;     END;
     END;               {  Open InFile for Data Input  }

   PROCEDURE  Read_Program_Control_Parameters;
      BEGIN
         Open_Input_File( DataFile );
         READLN( DataFile, JobTitle );
         READLN( DataFile, MaxSoln, Nsyst[1], Nsyst[2], Nsyst[3],
                 Tol2, Tol3, Lambda, LinCombo, Lzerck, kr, Nq1, Nq2,
                 Nz1, Nz2, Nr, List, MolMass, Density, Tolg );
         READLN( DataFile, ZeroCorr, PrintMerit, PrintLn );
      END;               {   Read Programme Control Parameters  }
         
   PROCEDURE  Check_and_Adjust_Input;
      BEGIN
         Enl := True;                TolKeep := Tolg;
         IF kr = 0 THEN Enl := False;    Self := Enl;
         IF Tolg < 0.1 THEN Tolg := 6.0;
         Tolg := 0.01 * Tolg;  TolKeep := TolKeep * 0.01;
         TwoRadian := 90.0 / ArcTan( 1.0 );
         Rrad := ArcTan( 1.0 )/ 45.0;
         Rad := 1.0 / Rrad;  HalfLambda := 0.005 * Lambda;
         WaveSqr := SQR( 0.5*Lambda);
      END;                        {   Check and Adjust Input  }

  PROCEDURE  Sort_Input_Values( Istyg : BYTE );
     VAR
        Faktor      :  REAL;
        I, Iend     :  INTEGER;

   PROCEDURE  Insert_the_Value( I : INTEGER );
     VAR
        J, K, M, Mend  :  INTEGER;
        Obstr          :  REAL;
     LABEL 50;
     BEGIN
        Obstr := Obs[i+1];
        FOR J := 1 to I DO BEGIN
            IF (Faktor * (Obstr - Obs[j])) < 0 THEN GoTo 50;  END;
 50:    Mend := I + 1 - J;
        FOR M := 1 to Mend DO BEGIN
            K := I + 1 - M;     Obs[k+1] := Obs[k];   END;
        Obs[j] := Obstr;
     END;               {   Insert the Value in ordered list  }

     LABEL  20;
     BEGIN                              {  Sort_Input_Values  }
        Faktor := 1.0;
        IF  Istyg = 0 THEN   Faktor := - 1.0;
 20 :   Iend := MaxObsLns - 1;
        FOR I := 1 to Iend DO
           IF ( Faktor * (Obs[i+1] - Obs[i])) < 0 THEN
                 Insert_the_Value( i );
        FOR I := 1 to Iend DO
           IF ABS( Obs[i+1] - Obs[i] ) < 0.0001 THEN BEGIN
              Obs[i+1]  := Obs[maxobslns];
              MaxObsLns := MaxObsLns - 1;           GoTo 20;   END;
    END;          {    Sort  Two Theta, d-  or  Q  values  input  }

   PROCEDURE  Print_the_Observed_Values;
      VAR
         I, J, L  :  INTEGER;
      BEGIN
         WRITELN(Lstr, MaxObsLns:2,'  lines input.', ChrWdt[Prnt,1] );
         I := 1;
         REPEAT
            J := I + 15;  IF J > MaxObsLns THEN  J := MaxObsLns;
            FOR L := I to J DO WRITE(Lstr, Obs[l]:8:3); WRITELN(Lstr);
            I := J + 1;
         UNTIL J = MaxObsLns;
      END;                    {   Print the Observed input values  }

  PROCEDURE  Compute_Q_Values( Iz  : BYTE );
     VAR
       M  :  INTEGER;
    BEGIN
       WRITE( Lstr, ChrWdt[Prnt,3], Spaces(20), 'as  ');
       CASE Iz OF
         1:  WRITELN(Lstr, 'D -  V A L U E S');
         2:  WRITELN(Lstr, 'Q -  V A L U E S');
         3:  WRITELN(Lstr, '2 -  T H E T A       V A L U E S');
       END;  WRITELN(Lstr);
       IF Iz = 1 THEN   Sort_Input_Values( 0 )
                 ELSE   Sort_Input_Values( 1 );
       FOR M := 1 to MaxObsLns DO BEGIN
          CASE Iz OF
             1:  BEGIN      Sit := 0.5 * Lambda / Obs[m];
                    TwoTheta[m] := TwoRadian* ArcTan( Sit /Sqrt(1-Sqr(Sit)))
                                   + ZeroCorr;        END;
             2:     TwoTheta[m] := TwoRadian * Theta_Calculate( Obs[m] )
                                   + ZeroCorr;
             3:     TwoTheta[m] := Obs[m] + ZeroCorr;    END;
          D := Sin( 0.5 * Rrad * TwoTheta[m] );
          Obs[m] :=  D * D * (1.0E+4 / WaveSqr);         END;
    END;                 {    Compute  Q  Values   }

 PROCEDURE  Print_the_Q_Values;
     VAR
        I, J, M  :  INTEGER;
     BEGIN
        I := 1;
        WRITELN(Lstr, 'The input leads to the following  Q -   V A L U E S');
        WRITELN(Lstr );       WRITE(Lstr,'   Where   Q = (10E+4)/(d*d) =  ');
        WRITELN(Lstr, 'h h A + k k B + l l C + k l D + h l E + h k F');
        WRITELN(Lstr, ChrWdt[Prnt,1] );
        REPEAT
            J := I + 15;  IF J > MaxObsLns THEN J := MaxObsLns;
            FOR M := I to J DO WRITE(Lstr, Obs[m]:8:1); WRITELN(Lstr);
            I := J + 1;
           UNTIL J = MaxObsLns;
        WRITELN(Lstr, ChrWdt[Prnt,3] );
        WRITELN(Lstr,Spaces(30),'Zero Correction of ',ZeroCorr:6:3,' degrees',
                                ' has been made.');
        WRITELN(Lstr);
     END;                          {   Print Observed Q Values  }

   PROCEDURE  Print_Summary_Result;
      VAR
        A, B, C  :  REAL;
        L        :  INTEGER;
      BEGIN
         Fmineq := Obs[20] * 5.0E-4;
         IF  Fmineq < 0.8  THEN  Fmineq := 0.8;
         A := 8.0 * Fmineq;    C := 1.6 * Fmineq;
         FOR L := 2 to MaxObsLns DO BEGIN
             B := Obs[l] - Obs[l-1];
             IF B <  A  THEN   A := B;     END;
         Fmineq := 0.125 * A;
         IF Fmineq < 0.5 THEN  Fmineq := 0.5;
         Fminim := Obs[1] / 17;     A := 100 / SQRT( Fminim );
         WRITELN(Lstr,'Q-values considered equal when differ by less than ',
                           Fmineq:4:1);
         WRITELN(Lstr, 'Lower limit on A, B, C is ', Fminim:8:2);
         WRITELN(Lstr, 'Maximum edge of the unit cell is ',A:6:1,' Angstrom');
         WRITELN(Lstr, 'Tolerances for Zone, Lattice, Theta are:  ',
                       Tol2:10:2, Tol3:10:2, Tolg:10:2 );
         WRITELN(Lstr);
      END;                   {  Print Summary Results   }

  PROCEDURE  Compute_Tolerances;
     VAR
        U, SinTheta   :  REAL;
        L                    :  INTEGER;
        DoItAgain            :  BOOLEAN;

     PROCEDURE  Evaluate_Tolerances;
        VAR
           DelQmx, TolMax  :  REAL;
        BEGIN
           DoItAgain := False;
           DelQmx := Obs[20] / 1000;
           TolMax := Tolg * DelQmx / U;
           IF ABS( Tolg - TolKeep ) > 1.0E-4 THEN BEGIN
              Tolg := 2 * TolMax;
              IF Tolg > 0.06  THEN  Tolg := 0.06;
              TolKeep := Tolg;  DoItAgain := True;    END;
       END;                      {  Evaluate  Tolerances   }

      LABEL  265;
      BEGIN
         DoItAgain := True;
 265:    U := 0;
         IF DoItAgain = True THEN BEGIN
             FOR L := 1 to MaxObsLns DO BEGIN
                 SinTheta := Sin( 0.5 * Rrad * (TwoTheta[l] + Tolg));
                 Obspt[l] := SQR( SinTheta / HalfLambda );
                 Obsmt[l] := 2 * Obs[l] - Obspt[l];
                 IF l <= 20 THEN  U := U + Obspt[l] - Obsmt[l];   END;
             U := U / 40;      Evaluate_Tolerances;    GoTo 265;  END;
      END;               {   Compute the Tolerances  }

  PROCEDURE  Read_Diffraction_Lines;
     VAR
        L  :  INTEGER;
      BEGIN
         Location := 'Diffraction Lines.';      Progress_Line;
         L := 0;
         REPEAT
             L := L + 1;         READLN( Datafile, Obs[l] );
         UNTIL  (Obs[l] < 0.001) OR (L = 41);
         MaxObsLns := L - 1;      Print_the_Observed_Values;
         IF Obs[4] < Obs[3] THEN   Compute_Q_Values( 1 )       {  d-values }
           ELSE BEGIN
             IF Obs[ MaxObsLns ] > 180.0 THEN  Compute_Q_Values( 2 )   { Q }
                ELSE   Compute_Q_Values( 3 );       END;      {  2-thetas  }
         Print_the_Q_Values;         Compute_Tolerances;
         Print_Summary_Result;
     END;                    {  Read Diffraction Lines  }

             {  reads zones and complete lattices  }

 PROCEDURE  Lattice_Line_Combinations;
    VAR
       I, J, K  :  INTEGER;

 PROCEDURE  Recalculate_Cell;
    VAR
       Waar, War, Cl         :  ARRAY[1..6] of REAL;
       A, B, C, Fvol         :  REAL;
    LABEL  Return;
    BEGIN
       Nr := Index[40];
       FOR J := 1 to 6 DO Waar[j] := Roster[nr,j];
       IF (Waar[4] = 0) AND (Waar[6] = 0) THEN               GoTo Return;
       IF ABS( Waar[5] ) < 2*SQRT( Waar[1]*Waar[3] ) THEN    GoTo Return;
       WRITE(Lstr,'Input Lattice of DIRECT CELL  ');
       FOR J := 1 to 6 DO  WRITE(Lstr, Waar[j]:10:4 );
       WRITELN(Lstr);      WRITELN(Lstr);
       A := 1;  B := 2;
       FOR J := 1 to 3 DO BEGIN
          C := COS( Waar[j+3] / Rad );   Cl[j] := C;  Cl[j+3] := C;
          A := A - SQR( C );    B := B * C;   War[j] := Waar[j];
          War[j+3] := War[j];   Cel[j] := Waar[j+3] / Rad;
          Cel[j+3] := Cel[j];    END;
      Fvol := 10000 / (A + B );
      FOR J := 1 to 3 DO BEGIN
         Waar[j] := Fvol * SIN( Cel[j] ) / SQR( War[ j ] );
         Waar[j+3] := 2 * Fvol * (Cl[j+1] * Cl[j+2] - Cl[j])
                       / (War[j+1] * War[j+2]);      END;
      FOR J := 1 to 6 DO Roster[nr,j] := Waar[j];
    Return:  END;               {   Recalculate  Cell  }

  PROCEDURE  Read_Complete_Lattice;
     BEGIN
        I := 0;
        REPEAT
           I := I + 1;  IF I > Nroster THEN  I := Nroster;
           FOR K := 1 to 6 DO  READ( DataFile, Roster[i,k] );
           READLN( DataFile );
           Roster[i,7] := 99;   Roster[i,8] := I;  Index[40] := I;
           IF Roster[i,1] > 0 THEN BEGIN
              Recalculate_Cell;
              WRITE(Lstr, 'Input Lattice Nr. ', I:3 );
              FOR K := 1 to 6 DO WRITE(Lstr, Roster[i,k]:8:2 );
              WRITELN( Lstr );    END;
        UNTIL (I = Nroster) OR (Roster[i,1] < 0);
        Nr := i - 1;     WRITELN(Lstr);
        IF List = 1 THEN    List := Nr;
     END;                      {   Read Complete Lattices  }

   PROCEDURE  Read_Line_Combinations;
     BEGIN
        Nq1 := MaxObsLns;   Nq2 := MaxObsLns;  I := 0;
        REPEAT
           I := I + 1;
           FOR J := 1 to 4 DO READ( DataFile, Roster[i,j] );
        UNTIL Roster[i,1] < 0.01;
        Lincombo := I - 1;
        WRITELN(Lstr,'Program Using ',LinCombo:3,' given Line Combinations');
        FOR I := 1 to LinCombo DO BEGIN
            WRITE(Lstr,'Input Zone Nr. ',I:3);
            FOR J := 1 to 4 DO  WRITE(Lstr, Roster[i,j]:8:1);
            WRITELN(Lstr);   END;     WRITELN(Lstr);
     END;                        {  Read  Line  Combinations   }

   PROCEDURE  Print_Search_Ranges( Icd, i1, i2 : INTEGER );
      BEGIN
        WRITE(Lstr, 'Programme will try first ',I1:3 );
        CASE  Icd  OF
              1 :  WRITE(Lstr, '  LINES' );
              2 :  WRITE(Lstr, '  ZONES' );   END;
        WRITELN(Lstr, '  with first', I2:3);
      END;                   {   Print Search Ranges  }

     BEGIN                                 {  Lattice Line Combinations  }
         Location := 'Line - Lattice.';    Progress_Line;
         IF  Nq1 > Nq2 THEN BEGIN
             ix := Nq2;   Nq2 := Nq1;   Nq1 := Ix;  END;
         IF  Nz1 > Nz2 THEN BEGIN
             Ix := Nz2;   Nz2 := Nz1;   Nz1 := Ix;  END;
         IF Nr <> 0 THEN  Read_Complete_Lattice
         ELSE BEGIN
            IF Lincombo <> 0 THEN  Read_Line_Combinations ELSE BEGIN
               FOR I := 1 to Nq1 DO FOR J := I to Nq2 DO BEGIN
                   IF ABS( 4 * Obs[i] - Obs[j] ) > 6 THEN BEGIN
                      Lincombo := LinCombo + 1;
                      Roster[Lincombo,1] := Obs[i];
                      Roster[Lincombo,2] := Obs[j];  END;  END;  END;
            Print_Search_Ranges( 1, Nq1, Nq2 );
            Print_Search_Ranges( 2, Nz1, Nz2 );
         Message := '  No  ';
         IF Nsyst[1] = 1 THEN   Message := '   ORTHORHOMBIC  ';
         IF Nsyst[2] = 1 THEN   Message := '   MONOCLINIC    ';
         IF Nsyst[3] = 1 THEN   Message := '   TRICLINIC     ';
         WRITELN(Lstr, Message, 'Lattice is given an extra chance.' ); END;
     END;             {   Lattice  Line   Combinations   }

             {  checks lines data for systematic zero-shift  }

 PROCEDURE  Two_Theta_Tolerance;
     VAR
        M  :  INTEGER;
     BEGIN
         TolRad := 0.5 * Tolg * Rrad;
         FOR M := 1 to MaxObsLns DO BEGIN
             Z := HalfLambda * SQRT( Obs[m] );
             ObsTheta[m] := ARCTAN( Z / SQRT( 1-SQR(Z) ));
             ThetaPt[m]  := ObsTheta[m] + TolRad;
             ThetaMt[m]  := ObsTheta[m] - TolRad;    END;
     END;                     {    Two  Theta  Tolerance  }

 PROCEDURE Zero_Check;
     VAR
        M, I, J, K, Nfit, N, II                  :  INTEGER;
        Nrx , Nft                                :  ARRAY[1..21] of Byte;
        Q, TTT, Qmt, Qpt, Del                    :  ARRAY[1..40] of REAL;
        Hrad, Swavel, W, Sit,Qmax, SumDel,
         Dif, Qq, Fiii                           :  REAL;

 PROCEDURE  Compute_Q_Values;
     BEGIN
        K := I - 11;   Nrx[i] := K;
        FOR  M := 1 to MaxObsLns DO BEGIN
             TTT[m] := TTT[m] + 0.01;
             W := SQR( SIN(Hrad * TTT[m])/ Swavel);
             Q[m] := W;
             Qpt[m] := SQR( SIN( Hrad * (TTT[m] + Tolg) )/ Swavel );
             Qmt[m] := 2 * W - Qpt[m];   END;
       Qmax := Qpt[maxobslns] + 10;   Nfit := 0;   SumDel := 0;
     END;                     {   Compute  Q  Values  }

 PROCEDURE  Zero_Check_Results;
     BEGIN
        N := 0;
        FOR I := 1 to 21 DO BEGIN
          IF Nft[i] > N THEN BEGIN
             M := I;   N := Nft[i];  Dif := Del[i];  END;   END;
        Writeln(Lstr, 'Results of Zero Check -->    Index = ', M:3,
        '  Line = ', N:3, '  Deviation = ',Dif:8:3);
     END;                            {  Zero  Check  Results    }

     LABEL 50, 60, 90, 95;
     BEGIN
        Location := 'Zero  Check.';         Progress_Line;
        Hrad := ARCTAN( 1 )/90 ;   Swavel := Lambda / 200;
        FOR M := 1 to MaxObsLns DO   TTT[m] := TwoTheta[m] - 0.11;
        FOR I := 1 to 21 DO BEGIN
           Compute_Q_Values;
           FOR M := 1 to 6 DO BEGIN
               IF 4 * Q[m] > Qmax THEN                              GoTo 95;
               FOR J := 1 to 5 DO BEGIN
                   Ii := J + 1;    Fiii := II * II;   QQ := Fiii * Q[m];
                   IF Qq > Qmax THEN                                GoTo 90;
                   FOR K := 1 to MaxObsLns DO BEGIN
                       IF Qpt[k] < Qq THEN                          GoTo 50;
                       IF Qmt[k] > Qq THEN                          GoTo 60;
                       Nfit := Nfit + 1;
                       SumDel := SumDel + ABS( Qq-Q[k] )/Fiii;
  50:              END;
  60:         END;
  90:     END;
  95:     Del[i] := SumDel;    Nft[i] := Nfit;    END;
       Zero_Check_Results;
    END;                     {  Zero Check  }

     BEGIN                     {  Reader  main loop   }
         Location := 'Read Data';           Progress_Line;
         Read_Program_Control_Parameters;
         Screen_Display_Header;        Printed_Report_Header;
         Check_and_Adjust_Input;      Read_Diffraction_Lines;
         Two_Theta_Tolerance;
         IF LzerCk > 0  THEN     Zero_Check;
         Lattice_Line_Combinations;        CLOSE( DataFile );
     END;             {   Reader of input data file   }

  OVERLAY PROCEDURE   Crystallographic_Zone_Finding;
     VAR
        Nsrt                          :  BYTE;
        Temp                          :  ARRAY[0..41] of REAL;
        Nn, M, N, N4, I, J, K, L      :  INTEGER;
        Itl                           :  ARRAY[0..10000]   of BYTE;

          {  Sort zones, write zone report headers   }

 PROCEDURE  Zone_Head;
     BEGIN
        WRITELN(Lstr,'NewNr    A       B        F          Quality   OldNr ',
                    'Center  Nobs  Ncalc     Delta');
     END;

 PROCEDURE  Write_The_Zone( Iz :  INTEGER );
     VAR
        J   :  INTEGER;
     BEGIN
        WRITE(Lstr, Iz:4 );
        FOR J := 1 to 3 DO WRITE(Lstr, Store[iz,j]:8:1);
        WRITE(Lstr, Store[iz,4]:16:0 );
        FOR J:= 5 to 8 DO WRITE(Lstr, Store[iz,j]:6:0);
        WRITELN(Lstr,'  ',TwoRadian*Store[iz,9]:14:4);
     END;                      {   Write  the  Zone   }

 PROCEDURE  Write_Zones( Ix : BYTE );
     VAR
        I, J  :  INTEGER;
     BEGIN
        WRITELN(Lstr);
        CASE Ix OF
          2: WRITELN(Lstr, Spaces(20),'ZONES  after refinement on two-theta.');
          1: WRITELN(Lstr, Spaces(20),'ZONES  selected for  refinement.');
        END;
        Zone_Head;
        FOR I := 1 to Nr DO Write_The_Zone( I );         WRITELN(Lstr);
     END;                        {   Write Zones   }

 PROCEDURE Sort_Zone( Nsrt  :  INTEGER );
     VAR
        N, M, Nrm, Nuniq, K, Nsrt1, I, J   :  INTEGER;
        U                                  :  REAL;
        Absnt, Exch                        :  BOOLEAN;

PROCEDURE  SortZone1;
     BEGIN
        FOR M := N to Nr DO IF Store[m,nsrt] > U THEN BEGIN
            U := Store[m,nsrt];   Nrm := M;  END;
        FOR M := 1 to 9 DO BEGIN
            U := Store[n,m];   Store[n,m] := Store[nrm,m];
            Store[nrm,m] := U;  END;
        IF Enl THEN   Write_The_Zone( N );
        N := N + 1;
     END;                    {  Sort Zone   }

 PROCEDURE  Nsort_Not_7;
     BEGIN
        Nuniq := 1;
        FOR I := 2 to Nr DO BEGIN
            Absnt := True;
            FOR J := 1 to NUniq DO BEGIN
                U := 0;
                FOR K := 1 to 3 DO  U := U + ABS(Store[i,k]-Store[j,k]);
                IF U < FminEq THEN  Absnt := False;  END;
            IF Absnt  THEN BEGIN
                Nuniq := Nuniq + 1;
                FOR K := 1 to 9 DO  Store[Nuniq,k] := Store[i,k];  END;  END;
        Store[Nuniq+1,1] := 0;
        IF Nuniq < (Nz2+2) THEN Nr := Nuniq ELSE  Nr := Nz2 + 2;
     END;         {    Nsort Not  7   }

 PROCEDURE  Exchange;
     BEGIN
         Nsrt1 := Nsrt + 1;
         REPEAT
             Exch := False;
             FOR I := 1 to Nr DO BEGIN
                 IF (Store[i,nsrt]-Store[i+1,nsrt]) <= 0.1  THEN
                     IF (Store[i,nsrt1]-Store[i+1,nsrt1]) > 0.1 THEN BEGIN
                        FOR K := 1 to 9 DO BEGIN
                           U := Store[i+1,k];  Store[i+1,k] := Store[i,k];
                           Store[i,k] := U;   END;
                        Exch := True;   END;  END;
         UNTIL Exch = False;
         IF Nr > (Nz2+2) THEN Nr := Nz2 + 2;
     END;             {   Exchange   }

     BEGIN                        {  S o r t   Z o n e  }
        N := 1;
        IF Enl   THEN BEGIN
           WRITELN(Lstr,Spaces(20),'  * *    Zones after Evaluation');
           Zone_Head;      END;
        REPEAT     U := - 100;     SortZone1;     UNTIL  N > Nr;
        IF Nsrt <> 7 THEN  Nsort_Not_7
          ELSE  Exchange;
        Write_Zones( 1 );
     END;               {   Sort Zone  }

                 {  Find probable zones    }

 PROCEDURE  Find_Zone( A, B : REAL );
     VAR
        Be                                 :  ARRAY[1..10] of REAL;
        Ibe                                :  ARRAY[1..10] of INTEGER;
        F2, Rf2, Fmaax, Fmaas, Qmax, Y, Z  :  REAL;
        H, K, II, Iu, Limt, M, LastH,
          LastK, Min, Kg, I, J, L, Nrz     :  INTEGER;

  PROCEDURE  Initialize_Variables_Find;
    BEGIN
      FILLCHAR( Itl, SIZEOF( Itl ), 0 );    FILLCHAR( Ibe, SIZEOF( Ibe ), 0 );
      FILLCHAR( Be , SIZEOF( Be  ), 0 );    FILLCHAR( Cel, SIZEOF( Cel ), 0 );
      F2 := 6 / Tol2;   Rf2 := 1 / F2;  Fmaax := 1.95 * Sqrt( A * B );
      Limt := ROUND( F2 * Fmaax );
      IF Limt > 10000 THEN BEGIN   WRITELN;
           WRITE('   ITL size < ',Limt:7);  Limt := 10000;  END;
      LastH := 2;   LastK := 2;   Min := 2;  Qmax := Obs[maxobslns];
      IF (4 * (A+B+Fmaax)) > Qmax THEN BEGIN
           LastH := TRUNC( SQRT( Qmax / A) );
           IF LastH > 4 THEN LastH := 4;
           LastK := TRUNC( SQRT( Qmax / B ) );
           Min   :=  (4 * LastH * LastK) DIV 5;   END;
      FOR H := 1 to LastH DO FOR K := 1 to LastK DO BEGIN
            Y := A * H * H + B * K* K;   Z := F2 /(H * K);
            FOR M := 1 to MaxObsLns DO BEGIN
                L := ROUND( Z * ABS( Obs[m] - Y ) );
                IF L <= Limt THEN Itl[l] := Itl[l] + 1;  END;    END;
      Limt := Limt - 6;   Kg := ROUND( F2 * Fmineq );
      IF Kg = 0 THEN Kg := 1;
      FOR I := 1 to Kg DO Itl[0] := Itl[0] + Itl[i];
      Itl[0] := 2 * itl[0];   Ibe[1] := Itl[0];  Be[1] := 0;  Min := 1;
      Nrz := 1;
    END;           {   Initialize  Variables for  Find Zones  }

     LABEL  90, 105, Return;
     BEGIN                                         {   Find  Zones   }
        GoToXY(20,WhereY);
        WRITE('Finding the Zones of  A, B ! ', A:9:2, B:10:2);
        Initialize_Variables_Find;    I := Kg;
        WHILE I <=  Limt DO BEGIN
            K := 0;
            FOR L := I to I+5  DO K := K + Itl[l];
            IF K <= Min THEN                                    GoTo 90;
            Iu := 0;
            FOR L := I to I+5  DO Iu := Iu + Itl[l] * L;
            FOR L := 1 to Nrz+1 DO   IF K > Ibe[l] THEN BEGIN
                 FOR N := Nrz downto L DO BEGIN
                       Ibe[n+1] := Ibe[n];   Be[n+1] := Be[n];  END;
                 Ibe[l] := K;  Be[l] := Rf2 * Iu / K;
                 IF Ibe[9] > Min THEN Min := Ibe[9];
                 Nrz := Nrz + 1;
                 IF Nrz > 9 THEN Nrz := 9;                  GoTo 90;  END;
 90:       I := I + 2;   END;  { While }
       II := 0;
       FOR I := 1 to Nrz-1 DO BEGIN
           Fmaas := Be[i];
           FOR L := I+1 to Nrz DO IF ABS( Be[l] - Fmaas) < Fmineq THEN BEGIN
                    Fmaas := 0.5 * (Fmaas + Be[l]);  Be[l] := 100000.0;  END;
           IF Fmaas < 10000.0 THEN BEGIN
              II := II + 1;  Cel[ii] := Fmaas;  Cel[ii+5]:= Ibe[i];
              IF II = 5 THEN                      GoTo Return; END;   END;
Return: END;                  {   Find Zones   }

             {  Criticise and Normalize Zones   }

 PROCEDURE  Criterion( k1, l1 : INTEGER; Center : BOOLEAN;
                     VAR Crit : REAL;
                   VAR No, Nc : INTEGER );
     VAR
        H, Hg, I, J, K, Lmin,
         Kg, Lstep, Min, Ii, Ncentr    :  INTEGER;
        Eva                            :  ARRAY[-8..8,0..8] of REAL;
        Zero                           :  REAL;
        Calc, Pres                     :  BOOLEAN;

 PROCEDURE  Initial_Criterion;
     BEGIN
        Zero := 1E-30;  Ncentr := 1;
        IF Center  THEN Ncentr := 2 * k1;
        FOR I := -8 to 8 DO FOR J := 0 to 8 DO Eva[i,j] := Eval[i,j];
        Kg := Norder * K1;     Hg := Norder * L1;   Nc := 0;   No := 0;
        Crit := 1;  Lstep := l1;
        IF Center  THEN Lstep := Ncentr;
     END;     {   Initialize Criter Variables   }

 PROCEDURE  Criterion1;
     VAR
         Nh, Nk, L, N   :  INTEGER;
     LABEL  30, 85, 90, 95, 100;
     BEGIN
         L := Lmin;
         REPEAT
            H := L;
 30:        IF ABS( Eva[h,k] ) < Zero THEN                     GoTo 100;
             Pres := False;
             FOR N := 1 to 8 DO BEGIN
                 Nh := N * H;    Nk := N * K;
                 IF (ABS( Nh ) > Hg) OR (Nk > Kg)  THEN        GoTo  95;
                 IF Pres   THEN                                GoTo  85;
                 IF (Center ) AND (((nh+nk) MOD Ncentr) <> 0)
                           THEN                                GoTo  90;
                 IF Eva[nh,nk] > Zero THEN BEGIN
                       Crit := Crit * Eva[nh,nk];   Pres := True;  END;
 85:             Eva[nh,nk] := 0;
 90:        END;
 95:        IF Pres  THEN No := No + 1;
            Nc := Nc + 1;
            IF RecTangl  THEN                                  Goto 100;
            IF h <= 0 THEN                                     GoTo 100;
            H := - L;                                          GoTo  30;
100:        L := L + Lstep;
       UNTIL L > Hg;
    END;                   {   Criterion  1   }

     BEGIN                                      {  C r i t e r i o n  }
        Initial_Criterion;        K := 0;
        REPEAT
            Lmin := 0;
            IF Center  THEN  Lmin := K MOD Ncentr;
            IF K = 0 THEN Lmin := Lstep;
            Criterion1;   K := K + K1;
        UNTIL  K > Kg;
        Min := 3;
        IF RecTangl  THEN Min := 2;
        No := No - Min;    Nc := Nc - Min;
        IF 4 * No >= Nc THEN BEGIN
             J := No;
             IF 2 * No > Nc THEN J := Nc - No;
             IF J <> 0 THEN     FOR Ii := 1 to J DO BEGIN
                  I := ii - 1;  Crit := Crit * (Nc - I)/ Ii;  END;
        END  ELSE Crit := 100;
        Crit := 1 / Crit;
        IF Crit > 1E+10 THEN Crit := 9.9999999E+9;
     END;                     {     Criterion   }

 PROCEDURE  Normalize_Zone( VAR A, B, Fmxxs  :  REAL );
     VAR
        Aa   :  REAL;
     LABEL  10, 20, 30, 40, 50, 60, 65;
     BEGIN
  10:    IF A < Fminim                THEN                        GoTo 65;
         IF (B - A) >= 0              THEN                        GoTo 20;
         Aa := A;   A := B;  B := Aa;
  20:    IF ABS( B - Fmxxs) >= Fmineq THEN                        GoTo 30;
         Aa := A;   A := 0.25 * B;  B := Aa - A; Fmxxs := 0;      GoTo 10;
  30:    IF Fmxxs <= Fmineq           THEN                        GoTo 40;
         IF ABS( A - B ) > Fmineq     THEN                        GoTo 40;
         A := 0.25 * (A + B - Fmxxs);   B :=  A + 0.5 * Fmxxs;
         Fmxxs := 0;                                              GoTo 10;
  40:    IF ABS( A - Fmxxs) > Fmineq  THEN                        GoTo 50;
            A := 0.25 * A;    B := B - A;   Fmxxs := 0;           GoTo 10;
  50:    IF (A - Fmxxs) > 0           THEN                        Goto 60;
            B := A + B - Fmxxs;  Fmxxs := ABS( Fmxxs - 2 * A );   GoTo 10;
  60:    IF Fmxxs < Fmineq THEN Fmxxs := 0;
  65:  END;                   {     Normalize  Zones   }

               {  Least-Squares refinement of zones  }

 PROCEDURE  Evaluate_Zones( Nr8, Isw8 : INTEGER );
     VAR
        Hulp                                      :  ARRAY[1..6]      of REAL;
        Noar, Ncar                                :  ARRAY[1..6]   of INTEGER;
        Ref                                       :  ARRAY[1..4,1..6] of REAL;
        Fyn, War                                  :  ARRAY[1..4]      of REAL;
        H, Nn, No, Nc, Hg, Kg, Hmin, Orde,
           Orp1, Orp2, I, J, K, Ispat             :  INTEGER;
        Repeet, Typ                               :  BOOLEAN;
        A, B, Fmaas, Aa, Crit, Cr, Centr,
           ZeroShift, Q, Z, X, Y, Zero, Deter     :  REAL;

 PROCEDURE  Initialize_Variables_LS;
     BEGIN
        FILLCHAR( Fyn,   SIZEOF( Fyn ),   0 );
        FILLCHAR( War,   SIZEOF( War ),   0 );
        FILLCHAR( Ref,   SIZEOF( Ref ),   0 );
        FILLCHAR( Eval,  SIZEOF( Eval ),  0 );
        FILLCHAR( Index, SIZEOF( Index ), 0 );
        Nindex := 0;
        Hg := ROUND( SQRT( QcalMx / A ));   Kg := ROUND( SQRT( QcalMx / B ) );
        IF  Hg > 8 THEN Hg := 8;       IF  Kg  > 8 THEN Kg := 8;
        IF Fmaas < 0.1 THEN BEGIN
           Fmaas := 0;   Hmin := 0;  RecTangl := True;   Orde := 2;  END
           ELSE BEGIN
              Hmin := - Hg;  RecTangl := False;  Orde := 3;  END;
        IF Typ THEN
               WRITELN(Lstr,'           *A* up page.     *B* across the page');
        Orp1 := Orde + 1;
     END;           {   Initialize variables for Least-Squares Refinement   }

 PROCEDURE  Loop_The_Variables( Ixz : INTEGER );
     VAR
        N1, N2, Itrm   :  INTEGER;

 PROCEDURE  Ls_Variables;
     BEGIN
         Hulp[1] := SQR ( H );   Hulp[2] := SQR( K );   Hulp[3] := H * K;
         Hulp[Orp1] := Obs[j] - Q;  Itrm := Orde;
     END;                   {    Ls  Variables   }

 PROCEDURE  LsTheta_Variables;
    VAR
       C   :  REAL;
    BEGIN
        C := WaveSqr / SIN( 2 * Z );        Hulp[1] := 1;  Itrm := Orp1;
        Hulp[2] := SQR( H ) * C;            Hulp[3] := SQR( K ) * C;
        Hulp[4] := H * K * C;            Hulp[Orp2] := ObsTheta[j] - Z;
    END;     {    LsTheta  Variables   }

     LABEL 55, 59, 60, 65;
     BEGIN                                          {  Loop the Variables  }
        IsPat := 0;    X := A * SQR( H );    Y := Fmaas * H;
        FOR K := 0 to Kg DO  BEGIN
            Ispat := Ispat + 1;
            IF (K = 0) AND (H <= 0 ) THEN                       GoTo 60;
            Q := X + B * SQR( K ) + Y * K;
            IF Q > QcalMx THEN                                  GoTo 65;
            IF Typ  THEN Wrisp( Ispat, Q, 0 );
            FOR I := 1 to MaxObsLns DO  IF Q < UpBndr[i] THEN BEGIN
                Eval[h,k] := - Prob[i];                         GoTo 55;   END;
 55:        IF ixz = 2 THEN  Z := Theta_Calculate( Q );
            FOR J := 1 to MaxObsLns DO BEGIN
               CASE Ixz OF
                  1  :  BEGIN
                       IF ObsPt[j] < Q THEN                     GoTo 59;
                       IF ObsMt[j] > Q THEN                     GoTo 60;   END;
                  2  :  BEGIN
                       IF ThetaPt[j] < Z THEN                   GoTo 59;
                       IF ThetaMt[j] > Z THEN                   GoTo 60;   END;
               END;
               Index[j] := 1;     Eval[h,k] := - Eval[h,k];
               IF Typ  THEN  WriSp( Ispat, Q, 1 );
               CASE Ixz OF
                   1:  Ls_Variables;
                   2:  LsTheta_Variables;     END;
               FOR N2 := 1 to Orp1 DO BEGIN
                  IF Ixz = 2 THEN  War[n2] := War[n2] + Hulp[n2] * Hulp[Orp2];
                  FOR N1 := 1 to Itrm DO
                    Ref[n1,n2] := Ref[n1,n2] + Hulp[n1] * Hulp[n2];  END;
                  GoTo 60;
  59:      END;
  60:  END;
  65:  CASE Ixz OF
           1 : IF Typ  THEN WriSp( Ispat, Q, 2 );
           2 : IF Typ  THEN Wrisp( Ispat, Z, 2 );   END;
     END;                         {  Loop the Variables   }

 PROCEDURE  Least_Squares_for_Theta;

 PROCEDURE  Update_Variables;
     VAR
        I, J   :  INTEGER;
     BEGIN
        FOR I := 1 to Orp1 DO FOR J := 1 to Orp1 DO
           Fyn[i] := Fyn[i] + Ref[i,j] * War[j];
        ZeroShift := Fyn[1];    A := A + Fyn[2] * 1E+4;
        B := B + Fyn[3] * 1E+4;  Zero := ZeroShift * TwoRadian;
        IF Orde = 3 THEN Fmaas := Fmaas + Fyn[4] * 1E+4;
        IF Typ  THEN
            WRITELN(Lstr, Spaces(12), 'Zone   ',Nr8:2, ' after refinement:   ',
             A:10:1,'  ',B:10:1,'  ',Fmaas:10:1,'   Zero shift = ',Zero:10:4);
     END;              {   Update  Variables   }

     BEGIN                           {   Least Squares on Theta   }
        Initialize_Variables_LS;       Orp2 := Orde + 2;
        FOR H := Hmin to Hg DO  Loop_The_Variables( 2 );
        FOR J := 1 to Nq2 DO   Nindex := Nindex + Index[j];
        FOR I := 1 to 4 DO FOR J:= 1 to 4 DO ArrayD[i,j] := Ref[i,j];
        Arrange_Array_Single( Orp1, 4 );        MatInvert( Orp1, Deter );
        Arrange_Array_Double( Orp1, 4 );
        FOR I := 1 to 4 DO FOR J := 1 to 4 DO Ref[i,j] := ArrayD[i,j];
        IF ABS( Deter ) < 1.0E-6 THEN BEGIN
           IF Typ  THEN  WRITELN(Lst,'Zone ',Nr8:2,'  not refined.');
             END  ELSE  Update_Variables;
     END;                 {  Least  Squares  on  Theta   }

                    {  Evaluates zones found above       }

 PROCEDURE Least_Squares_for_Lattice_Vectors;

  PROCEDURE  Fill_Reference_Array;
     VAR
        H, K, L  :  INTEGER;
     BEGIN
       FOR H := 1 to 2 DO FOR K := 1 to 3 DO Ref[k,h+4] := Ref[k,h];
       FOR H := 1 to 4 DO  BEGIN
           K := H + 1;    L := H + 2;
           Ref[4,h] :=  Ref[1,h] * (Ref[2,k] * Ref[3,l] - Ref[3,k] * Ref[2,l])
                       -Ref[2,h] * (Ref[1,k] * Ref[3,l] - Ref[3,k] * Ref[1,l])
                       +Ref[3,h] * (Ref[1,k] * Ref[2,l] - Ref[2,k] * Ref[1,l]);
            END;
      END;                {   Fill  Reference  Array   }

     LABEL  140, Return;
     BEGIN                                      {  Least Squares  }
        Initialize_Variables_LS;
        FOR H := Hmin to Hg DO  Loop_The_Variables( 1 );
        FOR J := 1 to Nq2 DO   Nindex := Nindex + Index[j];
        IF RecTangl  THEN BEGIN
           Z := Ref[1,1]*Ref[2,2] - Ref[2,1]*Ref[1,2];
           IF ABS( Z )< 0.1 THEN BEGIN    IF Typ  THEN
              WRITELN(Lstr,'Zone ',Nr8:2,'  not refined. -Z-');
                                                            GOTO Return;  END;
           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;
       END  ELSE BEGIN
           Fill_Reference_Array;
           IF ABS( Ref[4,1] ) < 1.0 THEN BEGIN
              IF Typ  THEN WRITELN(Lstr,'Zone ',Nr8:2,'  not refined.  -Ref-');
                                                            GOTO Return;  END;
           FOR J := 1 to 3 DO  Fyn[j] := Ref[4,j+1] / Ref[4,1];
           Fmaas := Fmaas + Fyn[3];   END;
       A := A + Fyn[1];      B := B - Fyn[2];
       IF Typ  THEN WRITELN(Lstr, Spaces( 12 ),
                   'Zone    ',Nr8:2,'  after refinement:  ',
                    A:10:1,'  ', B:10:1,'  ',Fmaas:10:1);
Return :  END;           {  Least Squares for Lattice  Vectors   }

 PROCEDURE  Evaluate_Formulations_Of_Zone;
     VAR
        I, J, K, N   :  INTEGER;

 PROCEDURE What_You_Think( wxz : BOOLEAN );
     BEGIN
         Nn := Nn + 1;             Criterion( I, J, Wxz, Crit, No, Nc);
         NoAr[nn] := No;      NcAr[nn] := Nc;        Hulp[nn] := Crit;
     END;                 {   What  do You  think  }

     BEGIN
         Nn := 0;
         FOR I := 1 to 2 DO For J := 1 to 2 DO  What_You_Think( False );
         FOR J := 1 to 2 DO BEGIN   I := J;     What_You_Think( True  );  END;
         Aa := 10 * Hulp[4];  Nn := 4;
         FOR  N := 1 to 6 DO  IF Hulp[n] > Aa THEN BEGIN
                 Aa := Hulp[n];   Nn := N;  END;
     END;             {    Evaluate Formulations of the Zone   }

 PROCEDURE  A_B_Adjustments;
     BEGIN
        A := 4 * A;     Fmaas := 4 * Fmaas;    B := 4 * B;
     END;                 {     A B Adjustments    }

 PROCEDURE  Centered_Screw_Zones;
     VAR
        E  :  REAL;
     BEGIN
          E := 2 * (B - A);  A := A + B - Fmaas;   B := A + 2 * Fmaas;
          Fmaas := E;
     END;                 {    Centered  Screw   Zones   }

PROCEDURE  Initial_Evaluate_Parameters;
     BEGIN
         GoToXY(5,WhereY);
         WRITE('Evaluate the Zone  ',Nr8:7,'      Stage = ',Isw8:3);
         Repeet := False;            A := Store[Nr8,1];       B := Store[Nr8,2];
         Fmaas  := Store[Nr8,3];    Cr := Store[Nr8,4];   Centr := Store[Nr8,6];
         ZeroShift := 0;           Typ := False;
         IF (Enl) AND (Isw8 = 2) THEN  Typ := True;
     END;                   {   Initial  Evaluate  Parameters   }

     LABEL   120, 130, 125, 200;
     BEGIN                                {   E v a l u a t e   }
        Initial_Evaluate_Parameters;
        REPEAT
           Repeet := NOT Repeet;
           IF Typ  THEN WRITELN(Lstr,'Zone Nr ',Nr8:2,'  Original',
                    ' Values: ',A:12:1,'  ',B:10:1,'  ',Fmaas:10:1,'  ',
                        TwoRadian*ZeroShift:10:4);
           Normalize_Zone( A, B, Fmaas );
           IF A < Fminim THEN                                    GOTO 130;
           QcalMx := 16.5 * (A + B + Fmaas);
           IF QcalMx > ObsPt[maxobslns] THEN  QcalMx := ObsPt[maxobslns];
           IF Isw8 = 2 THEN                                      GOTO 125
             ELSE   IF Isw8 = 3 THEN BEGIN
               Least_Squares_for_Theta;
               Store[Nr8,9] := ZeroShift;            GOTO 200; END;
           Least_Squares_for_Lattice_Vectors;
           Normalize_Zone( A, B, Fmaas );
           IF A < Fminim THEN                                    GOTO 130;
           A := 0.25 * A;   B := 0.25 * B;  Fmaas := 0.25 * Fmaas;
           Least_Squares_for_Lattice_Vectors;
           Evaluate_Formulations_Of_Zone;           Centr := 0;
           CASE Nn OF
                2 :  BEGIN  A := 4 * A;   Fmaas := 2 * Fmaas;  END;
                3 :  BEGIN  B := 4 * B;   Fmaas := 2 * Fmaas;  END;
                4 :  A_B_Adjustments;
                5 :  Centr := 1;
                6 :  BEGIN   Centr := 1;  A_B_Adjustments;  END;
                1 :           END;
           IF Fmaas < 0.01 THEN                              GOTO 120;
           IF Centr < 0.5  THEN                              GOTO 120;
           Centered_Screw_Zones;
         UNTIL Repeet = False;
 120:    Least_Squares_for_Lattice_Vectors;
         Cr := Hulp[nn];   Store[Nr8,7] := NoAr[nn];
         Store[Nr8,8] := NcAr[nn];
         Normalize_Zone( A, B, Fmaas );                       GOTO 130;
 125:    Least_Squares_for_Lattice_Vectors;
 130:    Store[Nr8,1] := A;   Store[Nr8,2] := B;      Store[Nr8,3] := Fmaas;
         Store[Nr8,4] := Cr;  Store[Nr8,6] := Centr;  Store[Nr8,9] := ZeroShift;
 200: END;                  {   Evaluate   }

            {  Zone search algorithm    }

  PROCEDURE  OneZone( A, B, C : REAL );
     VAR
        Zone           :  ARRAY[1..4] of REAL;
        Nend, N, I, J  :  INTEGER;
     LABEL 50;
     BEGIN
        Zone[1] := A;    Zone[2] := B;    Zone[3] := A + B + C;
        Zone[4] := A + B - C;              Nend := 4;
        IF (Zone[3] - Zone[4]) < 0.01 THEN Nend := 3;
        FOR I := 1 to Nend DO BEGIN
            N := 0;
            FOR J := 1 to MaxObsLns DO  IF Index[j] <> 1  THEN BEGIN
                   N := N + 1;       LinCombo := LinCombo + 1;
                   Roster[lincombo,1] := Zone[i]; Roster[lincombo,2] := Obs[j];
                   Roster[lincombo,3] := 0;       Roster[lincombo,4] := 0;
                   IF N = 4 THEN                           GOTO 50;  END;
  50:   END;
     END;   {  One Zone  }

  PROCEDURE  Rearrange_Zones;
     VAR
        Zz             :  ARRAY[1..7,1..2] of REAL;
        Ij             :  ARRAY[1..2]   of INTEGER;
        I, J, L, N, K  :  INTEGER;

 PROCEDURE   Fill_Zz;
     BEGIN
         Zz[1,i] := Cel[j+1];   Zz[2,i] := Cel[j+2];  Zz[5,i] := Cel[j+3];
         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[6,i] := Cel[j+1] * 2 + Cel[j+3];
         Zz[7,i] := Cel[j+1] * 2 - Cel[j+3];
     END;      {   Fill  Zz  Array   }

     LABEL  100, 150;
     BEGIN                                 {  Rearrange    Z o n e s   }
        J := 0;   I := 1;     Fill_Zz;
        J := 3;   I := 2;     Fill_Zz;
        FOR I := 1 to 4 DO  FOR J := 1 to 4 DO BEGIN
            Ij[1] := I;    Ij[2] := J;
            IF ABS( Zz[i,1] - Zz[j,2] ) < 1.0 THEN         GoTo 100;   END;
        Cel[4] := 0;                                       GoTo 150;
 100:   Cel[1] := 0;
        FOR N := 1 to 2 DO BEGIN
            L := 7 - N;    K := Ij[n];   Cel[1] := 0.5 * Zz[k,n] + Cel[1];
            IF K > 1 THEN BEGIN
               Cel[n+1] := Zz[1,n];   Cel[l] := Zz[k+3,n];
            END  ELSE  BEGIN
               Cel[n+1] := Zz[2,n];   Cel[l] := Zz[5,n];   END;  END;
 150: END;  {   Common Row of Two Zones   }

  PROCEDURE  Different_Sorting_Of_Zones;
     VAR
        I, J    :  INTEGER;
     BEGIN
         FILLCHAR ( Roster, SIZEOF( Roster ), 0 );      I := 0;
         REPEAT      I := I + 1;         UNTIL    Store[i,1] < 0.1;
         Nr := I - 1;  Nsrt := 7;  Nz1 := Nz1 + 2;   Nz2 := Nz2 + 2;
     END;  {   Create a Different Sorting of the Zones  }

 PROCEDURE  Initiate_Probabilities;
     VAR
        I, J  :  INTEGER;
     BEGIN
        Location := 'Probabilities.';    Progress_Line;
        FILLCHAR( Store, SIZEOF( Store ), 0 );
        FOR I := 1 to MaxObsLns DO   Temp[i] := Obs[i];
        Temp[0] := 0.5 * Obs[1];
        Temp[Maxobslns+1] := 2 * Obs[maxobslns] - Obs[maxobslns-1];
        FOR I := 1 to MaxObsLns DO BEGIN
            Upbndr[i] := 0.5 * ( Temp[i] + Temp[i+1] );
            Prob[i] := 2 * (Obspt[i] - Obsmt[i]) / (Temp[i+1] - Temp[i-1]);
            IF Prob[i] > 0.5 THEN Prob[i] := 0.5;   END;
        Nn := 0;
     END;              {   Initiate   Probabilities    }

 PROCEDURE  Search_Begins;
     BEGIN
        WRITELN( Lstr );
        IF Nn = 0 THEN
          WRITELN(Lstr, Spaces(12),
                  '* *   The Search for ZONES Commences   * *');
        WRITELN(Lstr, ' NR     A       B       F    Hits');
        WRITELN(Lstr);
     END;                   {   Search  Begins  }

 PROCEDURE  Print_Zone;
     VAR
        Ipz  :  INTEGER;
     BEGIN
         WRITE(Lstr, Nr:3);
         FOR Ipz := 1 to 3 DO WRITE(Lstr, Store[nr,ipz]:8:1);
     END;                  {  Print Zone  }

 PROCEDURE  Print_Possible_Combinations;
     VAR
        I  :  INTEGER;
     BEGIN
         WRITE(Lstr, J:3, K:3);
         FOR I := 1 to 3 DO WRITE(Lstr, Cel[i]:8:1);
         WRITE(Lstr,'     xxx');
         FOR I := 5 to 6 DO WRITE(Lstr, Cel[i]:8:1);
         WRITELN(Lstr, L:6, TwoRadian*Roster[n4,9]:14:4);
     END;             {   Print Possible Combinations   }

             {  Seek zones for diffraction data   }
                                              
  PROCEDURE  Search_for_Zones;
     LABEL  FindZone, LoopEnd, Return;
     BEGIN
        Nn := Nn + 1;  Nr := 0;
        FOR J := 1 to LinCombo DO BEGIN
           IF (Roster[j,3] < 0.05) AND
              (Roster[j,4] < 0.05) THEN              GoTo FindZone;
           Nr := Nr + 1;
           FOR I := 1 to 4 DO  Store[Nr,i] := Roster[j,i];
           Store[nr,5] := Nr;
           IF Enl = True THEN  BEGIN
               Print_Zone;   WRITELN(Lstr, Store[nr,4]:8:1);
               END;                                 GoTo  LoopEnd;
  FindZone:  Find_Zone( Roster[j,1], Roster[j,2] );
           FOR I := 1 to 5 DO BEGIN
              N  := ROUND( Cel[i+5] );
              IF N  = 0 THEN                        GoTo LoopEnd;
              IF Nr >= NrMax THEN                   GoTo Return;
              Nr := Nr + 1;                  Store[nr,5] := Nr;
              Store[nr,1] := Roster[j,1];    Store[nr,2] := Roster[j,2];
              Store[nr,3] := Cel[i];         Store[nr,4] := 0;
              IF Enl = True THEN BEGIN
                 Print_Zone;   WRITELN(Lstr, N:4);   END;   END;
 LoopEnd:   END;
 Return:   WRITELN;    WRITELN;
     END;           {   Search  for  Zones   }

  PROCEDURE  Combination_of_Zones;
     LABEL  LoopIn, Done, ByPas, Rear;
     BEGIN
        Location := 'Combine Zones.';  Progress_Line;
        WRITELN;         Nr := N;
        FOR I := 1 to Nr DO  IF Store[i,7] > 3 THEN Evaluate_Zones( i,3 );
        WRITELN;    WRITELN;     Write_Zones( 2 );
LoopIn:  IF Nz1 > Nr THEN Nz1 := Nr;
        IF Nz2 > Nr THEN Nz2 := Nr;
        N := 0;
        IF Enl  THEN BEGIN
         WRITELN(Lstr,'Possible Combination of ZONES  on which to base search');
         WRITELN(Lstr,'Combinat    A       B       C       D       E       F',
                     '    NewNr     Delta');  END;
        FOR J := 1 to Nz1 DO BEGIN
             M := 0;
             FOR I := 1 to 3 DO Cel[i] := Store[j,i];
             FOR K := J to Nz2 DO BEGIN  IF K = J THEN BEGIN
                    Cel[1] := Store[j,2];   Cel[2] := Store[j,1];  END;
Rear:            FOR I := 1 to 3 DO Cel[i+3] := Store[k,i];   Rearrange_Zones;
                 IF Cel[4] < 0.1 THEN                              GoTo ByPas;
                 N := N + 1;    N4 := 4 * N;    L := 10 * J + K;
                 FOR I := 1 to 6 DO  Roster[n4,i] := Cel[i];
                 Roster[n4,8] := L;
                 Roster[n4,9] := 0.5 * (Store[j,9] + Store[k,9]);
                 IF Enl  THEN    Print_Possible_Combinations;
                 IF N = 36 THEN                                     GoTo Done;
                 FOR I := 1 to 3 DO BEGIN
                     M := M + 1;     Cel[i] := Store[j,i];     END;
                 IF M > 4 THEN                                     GoTo ByPas;
                 IF K = J THEN                                      GoTo Rear;
ByPas: END;   END;
      IF N  >= 6  THEN                                              GoTo Done;
      IF Nz1 = Nr THEN                                              GoTo Done;
      Nz1 := 2 * Nz1;      Nz2 := 2 * Nz2;                        GoTo LoopIn;
Done: LinCombo := 4 * N;
     END;                 {   Combination of Zones   }

     LABEL  Abcd, Skip, Done;
     BEGIN                               {   Aaito1     Finding zones  }
        IF Nrun = 2 THEN BEGIN  Different_Sorting_Of_Zones;    GoTo Skip;  END;
        Initiate_Probabilities;     WRITELN(Lstr, ChrWdt[Prnt,1] );
 Abcd:  IF Enl = True THEN    Search_Begins;
        Location := 'Search for Zones.';    Progress_Line;
        Search_for_Zones;
        Location := 'Evaluate Zones 1.';    Progress_Line;
        FOR I := 1 to Nr DO  IF Store[i,4] <= 1 THEN BEGIN
             Evaluate_Zones( i, 1 );
             IF Store[i,1] < Fminim THEN Store[i,4] := -10;  END;
        Nsrt := 4;  WRITELN;   WRITELN;
 Skip:  Sort_Zone( Nsrt );
        IF Nr > Nz2 + 2 THEN  Nr := Nz2 + 2;
        Location := 'Evaluate Zones 2.';  Progress_Line;    N := 0;
        FOR I := 1 to Nr DO BEGIN
            Evaluate_Zones( i,2 );
            IF Store[i,1] < Fminim THEN                             GoTo Done;
            N := N + 1;
            FOR J := 1 to 9 DO Store[n,j] := Store[i,j];
            IF Nindex < (Nq2 - 1) THEN                              GoTo Done;
            IF Nn >= 2 THEN                                         GoTo Done;
            FOR J := 1 to N DO FOR K := 1 to 4 DO Roster[j,k] := Store[j,k];
            LinCombo := N;    OneZone( Store[i,1], Store[i,2], Store[i,3] );
            WRITELN(Lstr,'ZONE ',I:2,' indexes ',Nindex:2,' lines from ',
                                   Nq2:2,' used for Zone finding');
            WRITELN;  WRITELN;                                     GoTo Abcd;
Done:    END;
        Combination_of_Zones;
     END;                     {  Aaito1  finding zones  }

  OVERLAY PROCEDURE   Crystallographic_Lattice_Finding;
    VAR
       N4, Kkk, Mj, Nsol, Min            :  INTEGER;
       BravType                          :  ARRAY[1..6] of CHAR;

           {  check for symmetry of lattice     }

PROCEDURE  Symmetry_Test;
   VAR
      Npass                                :  BYTE;
      A, B, C, E, Difmin, X, Z, G, Aa      :  REAL;
      N, I, J, H, K, L, Ih, Ik, Il, Ii     :  INTEGER;

PROCEDURE Describe_Lattice;
    BEGIN
      WRITELN(Lstr, 'The LATTICE with Reciprocal Constants');
      FOR I := 1 to 6 DO WRITE(Lstr, Cel[i]:9:2 );  WRITELN( Lstr );
      IF Npass = 1 THEN WRITE(Lstr, ' is probably ');
      IF Npass = 2 THEN WRITE(Lstr, '  might be ');
    END;   {  Describe the Lattice  }

PROCEDURE  Rhombohedral;
    BEGIN
      Describe_Lattice;      A := 200 / SQRT( B );      C := 100 / SQRT( Z );
      WRITELN(Lstr, 'RHOMBOHEDRAL, with (hexagonal) reciprocal lattice ',
                    'constants  X = ',X:9:1,'    Z = ', Z:9:1);
      WRITELN(Lstr,'   and Direct Lattice Constants of     A = ',A:7:2,
                      '     C = ',C:7:2);
      X := 4*x;           A := 0.5 * A;
      IF Cel[9] < 6 THEN WRITELN(Lstr, ' Probably X should be ',X:10:1,
           '  and  A ',A:8:2);
    END;             {   Rhombohedral  Lattice    }

PROCEDURE  Cubic( G  :  REAL );
    BEGIN
        Describe_Lattice;                X := 100 / SQRT( G );
        WRITELN(Lstr, 'CUBIC with reciprocal constant = ', G:9:1,
                       '  and direct constant ', X:9:2);
    END;              {   Cubic Lattice  }

PROCEDURE TetragonalHexagonal;
    BEGIN
        Describe_Lattice;     A := 100 / SQRT( X );   C := 100 / SQRT( Z );
        IF N < 3 THEN BEGIN
            WRITELN(Lstr,'TETRAGONAL with reciprocal constants',X:9:1, Z:9:1,
               ' and direct constants ',A:7:2, C:7:2);    N := -999;   END;
        IF N = 3 THEN BEGIN       A := 2*A / SQRT( 3 );
             WRITELN(Lstr,'HEXAGONAL with Reciprocal Constants ',X:9:1,Z:9:1,
                '  and Direct Constants ',A:7:2, C:7:2);  N := -998;  END;
        IF N > 3 THEN WRITELN(Lstr,'A FREAK.  There is a ratio in reciprocal ',
                'constants of ',N:3);
    END;               {   Tetragonal  and  Hexagonal  Lattices  }

    LABEL 5,10,20,30,40,50,60,200,1000,1050,1100,1150,3010,3020,3100, Return;
    BEGIN
      Location := '-*-  Symmetry Test.';      Progress_Line;
      Npass := 0;   A := Cel[1];  B := Cel[2];  C := Cel[3];  E := Cel[5];
      DifMin := Fmineq;
      IF (Cel[4]+ABS( Cel[6] )) > Fmineq THEN                     GoTo 3010;
   5: DifMin := DifMin * 3;   Npass := Npass + 1;
      IF E < Fmineq THEN                                          GoTo 1000;
      IF ABS( E-2*B ) >          DifMin THEN                      GoTo   20;
      IF ABS( C-3*B ) >          DifMin THEN                      GoTo   10;
      X := (B+C+E)/18;   Z := A-X;                                GoTo  200;
  10: IF ABS( A-3*B) >           DifMin THEN                      GoTo   20;
      X := (A+B+E) / 18;   Z := C - X;                            GoTo  200;
  20: IF ABS( A+B-C ) >          DifMin THEN                      GoTo   30;
      IF ABS( 2*ABS( B-A )-E ) > DifMin THEN                      GoTo   30;
      X := (C-A)/6 + B/6;   Z := (C+2*A-2*B) / 3;                 GoTo  200;
  30: IF ABS( B+C-4*A ) >        DifMin THEN                      GoTo   40;
      IF ABS( E-ABS( B-C ) ) >   DifMin THEN                      GoTo   40;
      Z := (C-A) / 3;   X := B/6 + (4*A-C) / 6;                   GoTo  200;
  40: IF ABS( 2*C-3*E ) >        DifMin THEN                      GoTo   50;
      IF ABS( 2*B+E-6*A ) >      DifMin*Npass THEN                GoTo   50;
      Z := (C+E) / 15;   X := (A-Z) / 2 + B / 6;                  GoTo  200;
  50: IF ABS( 2*A-3*E ) >        DifMin THEN                      GoTo   60;
      IF ABS( 2*B+E-6*C) >       DifMin*Npass THEN                GoTo   60;
      Z := (A+E) / 15;   X := (C-Z) / 2 + B/6;                    GoTo  200;
  60: IF Npass = 1 THEN                                           GoTo    5;
                                                                  GoTo 3010;
  200: IF Z < 1 THEN                                            GoTo Return;
      Rhombohedral;                                             GoTo Return;
 1000:  FOR I := 1 to 2 DO BEGIN
           J := 3 - i;   G := Cel[j];   H := ROUND( A/G );
           K := ROUND( B/G );           L := ROUND( C/G );
           Ih := 0;  Ik := 0;  Il := 0;
           IF ABS( A-G*H ) < DifMin THEN Ih := 1;
           IF ABS( B-G*K ) < DifMin THEN Ik := 1;
           IF ABS( C-G*L ) < DifMin THEN Il := 1;
           Ii := Ih + Ik + Il - 2;
           IF Ii > 0 THEN BEGIN   Cubic( G );              GoTo Return;  END;
           IF Ii = 0 THEN BEGIN
              IF Ik*Il = 0 THEN                                    GoTo 1050;
              N := L DIV K;   X := B;   Z := A;             GoTo 1150;   END;
 1050:     IF Ih*Il = 0 THEN BEGIN
                N :=  H DIV K;   X := B;   Z := C;  END
                 ELSE BEGIN   N :=  L DIV H;   X := A;   Z := B;  END;
 1150:     TetragonalHexagonal;            IF N < 900 THEN GoTo Return;  END;
 3010: Aa := Obs[maxobslns];  Z := 0;   DifMin := Fmineq;
       FOR J := 1 to 6 DO BEGIN              X := ABS( Cel[j] );
          IF X < Fmineq THEN                                      GoTo 3020;
          IF X < Aa THEN Aa := X;
          IF X >  Z THEN Z  := X;
 3020: END;
       FOR I := 1 to 3 DO BEGIN          X := Aa / I;
          IF X < 0.5*Fminim THEN                                  GoTo Return;
          L := TRUNC( Z/X + 0.5 );       X := Z/L;       N := 0;
          FOR J := 1 to 6 DO BEGIN
              G := ABS( Cel[j] );        L := ROUND( G/X );
             IF ABS( G-X*L ) <= DifMin THEN  N := N + 1;   END;
          IF N > 5 THEN                       GoTo 3100;   END;   GoTo Return;
 3100: Cubic( X );
 Return: WRITELN( Lstr );
   END;                       {  Symmetry Test  }

            {  compute unit cell parameters      }

PROCEDURE  UnitCell( Isw : BYTE );
    VAR
       Rec, Reel                                :  ARRAY[1..3,1..6] of REAL;
       Waar                                     :  ARRAY[1..6]      of REAL;
       I, J, N                                  :  INTEGER;
       W, Fnnn, Volum, PartD, Delta, Dmax, Dx   :  REAL;
    LABEL  Return;
    BEGIN
       FOR J := 1 to 6 DO Waar[j] := Cel[j];
       FOR J := 1 to 3 DO BEGIN
          Rec[1,j] := SQRT( Waar[j] );  Rec[1,J+3] := Rec[1,j];   END;
       FOR J := 1 to 3 DO BEGIN
          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 - W * W );  Rec[3,j] := Rec[3,j+3];  END;
       FOR J := 1 to 3 DO BEGIN
          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 - W*W );
          Reel[3,j+3] := Reel[3,j];  END;
       Fnnn := Reel[3,1] * Reel[3,2] * Rec[3,3];
       FOR J := 1 to 3 DO Cel[j] := Reel[3,j] / (Fnnn * Rec[1,j]) * 100.0;
       FOR J := 4 to 6 DO
         Cel[j] := 90.0 - 57.29578 * ARCTAN( Reel[2,j] / Reel[3,j] );
       Volum := Fnnn * Cel[1] * Cel[2] * Cel[3];
       IF Isw = 0 THEN                                        GoTo Return;
       FOR J := 1 to 6 DO WRITE(Lstr, Cel[j]:9:3 );
       WRITELN( Lstr, Volum:12:2 );
       IF MolMass < 0.1 THEN                                  GoTo Return;
       PartD := MolMass /  (0.6025 * Volum);
       IF Density < 0.1 THEN BEGIN
          WRITELN(Lstr, Spaces(20), 'Density calculated = Z * ', PartD:8:4);
                                                          GoTo Return;  END;
       J := TRUNC( Density / PartD );   Dmax := 1000;
       FOR I := 1 to 4 DO BEGIN
          Dx := (i + j-2)* PartD;   Delta := ABS( Density - Dx );
          IF Delta < Dmax THEN BEGIN
             Dmax := Delta;   N := I + J - 2;  END;       END;
       Dx := PartD * N;
       WRITELN(Lstr, Spaces( 20 ), 'Z = ',N:2,'   Dx = ',Dx:7:4,
                     '    D obs = ',Density:7:4);
    Return:  END;                {   Unit Cell Calculation  }

               {  standard settings of lattices     }

PROCEDURE  Wissel( I, J :  INTEGER );
   VAR
      A   :  REAL;
   BEGIN
      A := Cel[i];     Cel[i] := Cel[j];      Cel[j] := A;
      A := Cel[i+3]; Cel[i+3] := Cel[j+3];  Cel[j+3] := A;
   END;   {  Wissel  }

PROCEDURE  Normalize3;
   VAR
      A, B, C, D, E, F, Ped,
        Qmax, U, X, Y, Z, Q      :  REAL;
      Ii, I, J, K, L, N1,
        K1, H, H1, M, N, L1      :  INTEGER;
      Change                     :  BOOLEAN;
      Old                        :  ARRAY[1..6] of REAL;
      Zone                       :  ARRAY[1..3,1..6] of INTEGER;
      Qval                       :  ARRAY[1..5,1..5,1..5] of REAL;

PROCEDURE  Working_Variables;
   BEGIN
      A := Cel[1];  B := Cel[2];  C := Cel[3]; D := Cel[4]; E := Cel[5];
      F := Cel[6];  Ped := Cel[8];  
   END;

PROCEDURE  Return_to_Cell;
   BEGIN
      Cel[1] := A;  Cel[2] := B;  Cel[3] := C;  Cel[4] := D;  Cel[5] := E;
      Cel[6] := F;  Cel[8] := Ped;
   END;

PROCEDURE  Lattice_Has_Higher_Symmetry;
   LABEL  230, 1000;
   BEGIN
      FOR Ii := 1 to 3 DO BEGIN
         X := A; A := B;  B := c;  C := X;   X := D;  D := E;  E := F; F := X;
         IF B < Fminim THEN BEGIN
            Qmax := -999;  A := 0;                            GoTo 1000;  END;
         IF D <= 0 THEN BEGIN   D := - D;  E := - E;  END;
         IF E <= 0 THEN BEGIN   E := - E;  F := - F;  END;
         IF ABS( A-C ) > Fmineq THEN                          GoTo 230;
         IF E < Fmineq          THEN                          GoTo 230;
         A := (A + C + E) / 4;   C := A - E / 2;   F := (D + F) / 2;
         D := ABS( D-F );   E := 0;
 230:    IF ABS( E-A ) <= Fmineq  THEN BEGIN        A := A / 4;
            C := C - A;  E := 0;  F := F / 2;  D := ABS( D-F ); END;
         IF ABS( E-C ) <= Fmineq THEN BEGIN           C := C/4;
            A := A-C;  E := 0; D := D/2;  F := ABS( D-F );  END;
         IF ABS( D ) <= Fmineq THEN       D := 0;
         IF ABS( E ) <= Fmineq THEN       E := 0;
         IF ABS( F ) <= Fmineq THEN       F := 0;
     END;
 1000:  END;    {   Lattice Has Higher Symmetry  }

PROCEDURE  Change_Parameters;
   BEGIN
      U := 0;
      FOR I := 1 to 6 DO BEGIN
         U := U + ABS( Cel[i] - Old[i] );   Old[i]:= Cel[i];  END;
      Change := False;
      IF U > 0.1 THEN Change := True;
   END;   {   Change Parameters  }

PROCEDURE  Find_Smallest_Q;
   LABEL  410;
   BEGIN
      N1 := 0;
 410: N1 := N1 + 1;   U := Qmax - 10;
      FOR L := 3 to 4 DO FOR K := 2 to 4 DO FOR H := 2 to 4 DO
        IF Qval[h,k,l] < U THEN BEGIN
           U := Qval[h,k,l];  Zone[n1,3] := L;  Zone[n1,1] := H;
           Zone[n1,4] := H;   Zone[n1,2] := k;  Zone[n1,5] := K;  END;
      Qval[Zone[n1,1], Zone[n1,2], Zone[n1,3]] := Qmax;
      IF N1 < 3 THEN                                          GoTo 410;
      Ii := 0;
      FOR I := 1 to 3 DO BEGIN
         K := Zone[1,i+1] - 3;  L := Zone[2,i+2] - 3;  M := Zone[2,I+1]-3;
         N := Zone[1,i+2] - 3;  J := K * L - M * N;
         Ii := Ii + J * (Zone[3,i] - 3);  END;
      IF Ii = 0 THEN BEGIN N1 := 2;                           GoTo 410;  END;
   END;   {  Find Smallest Q  }

FUNCTION  Mazen : REAL;
   VAR
      I1, I2, I3   :  INTEGER;
   BEGIN
      I1 := Zone[h,1] + Zone[k,1] - 3;  I2 := Zone[h,2] + Zone[k,2] - 3;
      I3 := Zone[h,3] + Zone[k,3] - 3;   X := Qval[i1, i2, i3];
      I1 := Zone[h,1] - Zone[k,1] + 3;  I2 := Zone[h,2] - Zone[k,2] + 3;
      I3 := Zone[h,3] - Zone[k,3] + 3;   Y := Qval[i1,i2,i3];
      MaZen := (X - Y)/2;
   END;    {   Maxen  }

PROCEDURE Q_Value;
   BEGIN
       Q := A*h*h + B*K*K + C*L*L + D*k*l + E*H*L + F*H*K;
       IF Q < -0.1 THEN BEGIN  A := 0;  Q := -999;  END
         ELSE   Qval[h1,k1,l1] := Q;
   END;    {   Q  Value   }

   LABEL  200, 700, 750, Return;
   BEGIN                                       {  N o r m a l i z e 3  }
      Change := True;
      FOR I := 1 to 6 DO Old[i] := Cel[i];
      Qmax := 10000;
     200: REPEAT
         Working_Variables;
         Lattice_Has_Higher_Symmetry;  IF Qmax = -999 THEN   GoTo Return;
         Return_To_Cell;    Change_Parameters;
      UNTIL Change = False;
      FOR L1 := 3 to 4 DO BEGIN
         L := L1 - 3;
         FOR K1 := 2 to 4 DO BEGIN
            K := K1 - 3;
            FOR H1 := 2 to 4 DO BEGIN
               H := H1 - 3;
               IF (L = 0) AND (K < 1)  THEN Qval[h1,k1,l1] := Qmax
               ELSE BEGIN   Q_Value;  IF Q = -999 THEN   GoTo Return; END;
     END;  END;   END;
     Qval[4,3,3] := A;   Find_Smallest_Q;
     FOR L1 := 1 to 5 DO BEGIN
       L := L1 - 3;
       FOR K1 := 1 to 5 DO BEGIN
          K := K1 -3;
          FOR H1 := 1 to 5 DO BEGIN
             H := H1 - 3;  Q_Value;  IF Q = -999 THEN         GoTo Return;
     END;  END;  END;
     B := Qval[Zone[1,1], Zone[1,2], Zone[1,3]];
     A := Qval[Zone[2,1], Zone[2,2], Zone[2,3]];
     C := Qval[Zone[3,1], Zone[3,2], Zone[3,3]];
     H := 1; K := 3; D := Mazen;   K := 2;   F := Mazen; H := 3; E := Mazen;
     Return_To_Cell;   Change_Parameters;
     IF Change  THEN                                               GoTo 200;
     L := 0;  J := 1;   Return_To_Cell;
     FOR I := 4 to 6 DO BEGIN
        L := L + I;
        IF ABS( Cel[i] ) <= Fmineq THEN BEGIN
           J := J + 1;  K := i;  L := L - I;  Cel[i] := 0;  END;  END;
        Working_Variables;
        IF J = 1  THEN                                             GoTo 750;
        IF J = 2  THEN                                             GoTo 700;
        IF J = 4  THEN                                             GoTo Return;
        IF J = 3  THEN BEGIN  L := L - 3;
            IF L <> 2 THEN Wissel( l,2 );  Working_Variables;
            IF A >  C THEN Wissel( 1,3 );  Working_Variables;   E := ABS( E );
            IF E <= A THEN GoTo Return;    C := C + A - E;
            E := A + A - E;  Return_To_Cell;                   GoTo 200;  END;
 700 : K := K - 3;  IF K <> 3 THEN Wissel( k,3 ); Working_Variables;
 750 : IF D < 0 THEN BEGIN  D := - D;  E := - E;  END;
       IF E < 0 THEN BEGIN  E := - E;  F := - F;  END;
       IF ABS( A+B+F-D-E ) > Fmineq THEN                     GoTo Return;
       E := ABS( A-B ) * 0.5;   A := (A+B+F) / 4;   B := C-A;
       C := A-F/2;  D := 0;  F := 0;  Return_To_Cell;        GoTo 200;
  Return:  Return_to_Cell;
  END;                           {  Normalize3  }

            {  angle between intersecting zones  }
 
PROCEDURE Threed;
   VAR
      Nr3d, Min, J3d, LastNr, NewNr, Nc, LimH,
         L, K, H, It, Kg, Kst, M, N, I          :  INTEGER;
      F3, Rf3, A, B, C, E, F, W, Z, BigD       :  REAL;
      Som, Be                               :  ARRAY[1..10] of REAL;
      Isom, Ibe                             :  ARRAY[1..10] of INTEGER;
      Itl                                   :  ARRAY[-5000..5000] of BYTE;

PROCEDURE  Initialize_ThreeD_Variables;
   BEGIN
      J3d := J3d + 4;       LastNr := NewNr;        A := Roster[j3d,1];
      B := Roster[j3d,2];   C := Roster[j3d,3];     E := Roster[j3d,5];
      F := Roster[j3d,6];   NewNr := ROUND( Roster[j3d,8] );   Nc := 0;
      W := SQRT( (4*A*C - E*E) * (4*a*b - F*F) );
      BigD := (W-E*F) / (2.05*A);             LimH := ROUND( F3*BigD );
      IF LimH > 5000 THEN         LimH := 5000;
      FILLCHAR( Itl, SIZEOF( Itl ), 0 ); FILLCHAR( Isom, SIZEOF( Isom ), 0 );
      FILLCHAR( Som, SIZEOF( Som ), 0 ); FILLCHAR( Be,   SIZEOF( Be  ),  0 );
      FILLCHAR( Ibe, SIZEOF( Ibe ), 0 );
      Min := 4 +  MaxObsLns DIV 10 ;      Nr3d := 1;   WRITE( J3d:4);
    END;    {   Initialize  Three D  Variables   }

  PROCEDURE  MonoClinic_Lattice;
    VAR
       X, Y, Z            :  REAL;
       I, L, K, H, M, It  :  INTEGER;
    BEGIN
       FOR L := 1 to 2 DO  FOR K := 1 to 2 DO BEGIN
          X := B * SQR( K ) + C * SQR( L );    Z := F3 / (K * L);
          FOR H := 0 to 2 DO BEGIN         Y := X + A * SQR( H );
             FOR M := 1 to MaxObsLns DO BEGIN
                It := ROUND( Z * ABS( Obs[m]-Y ) );
                IF It < LimH THEN   Itl[it] := Itl[it] + 1;
             END;   END;  END;
       Kg := LimH - 5;   Kst := 2;   LimH := 1;   Nr3d := 1;
       FOR I := 0 to 2 DO         Itl[0] := Itl[0] + Itl[i];
       Itl[0] := 2 * Itl[0];
       IF Itl[0] >= Min THEN BEGIN
          Isom[Nr3d] := Itl[0];        Som[Nr3d] := 0;  END;
    END;     {   Monoclinic  Lattices   }

  PROCEDURE  Triclinic_Lattice;
    VAR
       X, Y, Z, V       :  REAL;
       L, K, H, It      :  INTEGER;
    BEGIN
       Kg := TRUNC( F3*(W+E*F)/(2.05*A) );
       IF Kg > LimH THEN     Kg := LimH;
       FOR L := 1 to 2 DO FOR K := -2 to 2 DO IF K <> 0 THEN BEGIN
          X := SQR( K ) * B + SQR( L ) * C;   V := K * F + L * E;
          Z := F3 / (K * L);
          FOR H := -2 to 2 DO BEGIN
             Y := X + A * SQR( H ) + V * H;
             FOR M := 1 to MaxObsLns DO BEGIN
                It := ROUND( Z * (Obs[m] - Y) );
                IF ABS( It ) <= Kg THEN   Itl[it] := Itl[It] + 1;
             END;  END;  END;      Kst := -Kg;      Kg := Kg - 5;
    END;                    {   Triclinic  Lattice  }

  PROCEDURE  Find_Most_Frequent_Values;
    VAR
       V   :  REAL;
    LABEL  380;
    BEGIN
       K := 0;             M := H + 5;
       FOR L := H to M DO  K := K + Itl[l];
       IF K > Min THEN BEGIN
          FOR N := 1 to Nr3d DO IF K > Isom[n] THEN BEGIN
             FOR Nc := N to Nr3d DO BEGIN
                I := Nr3d + N - Nc;   Isom[i+1] := Isom[i];
                Som[i+1] := Som[i];   END;
             Isom[n] := k;   Nr3d := Nr3d + 1;
             IF Nr3d > 8 THEN Nr3d := 8;
             V := 0;
             FOR L := H to M DO  V := V + L * Itl[l];
             Som[n] := Rf3 * ( V/K );                         GoTo 380;  END;
 380:    IF Isom[4] > Min THEN Min := Isom[4];  END;
    END;               {   Find the Most Frequent Values of D  }

 PROCEDURE  Best_Values;
    VAR
       Ln, H, K   :  INTEGER;

 PROCEDURE  FillRoster( L : INTEGER );
    VAR
       H, M, K  :  INTEGER;
    BEGIN
      FOR H := 1 to N DO BEGIN
          M := J3d + H - 1;
          FOR K := 1 to 6 DO      Roster[m,k] := Roster[j3d,k];
          Roster[m,4] := Be[h];   Roster[m,7] := Ibe[h];
          Roster[m,8] := L + H;              END;
    END;                    {   Fill   Roster   }

    LABEL  440, 500;
    BEGIN
       FOR H := 1 to Nr3d-1 DO FOR K := H+1 to Nr3d DO BEGIN
          IF (ABS(Som[h] - Som[k])) > Fmineq THEN          GoTo 440;
          IF Isom[k] < 1 THEN                              GoTo 440;
             Som[h] := (Isom[h]*Som[h] + Isom[k]*Som[k])
                            / (Isom[h]+Isom[k]);       Isom[k] := 0;
 440: END;
      N := 0;       Min := Isom[1] - 5;
      FOR H := 1 to Nr3d DO BEGIN
          IF Isom[h] > Min THEN BEGIN
          N := N + 1;    Ibe[n] := Isom[h];     Be[n] := Som[h];
          IF N = 4 THEN                      GoTo 500; END; END;
 500: Ln := 0;
      IF NewNr = LastNr THEN       Ln := 5;
      L := 10 * NewNr + Ln;        FillRoster( L );
    END;               {  Best Values  }

    BEGIN                              {  T h r e e D  }
       F3 := 6 / Tol3;    Rf3 := Tol3 / 6;   NewNr := 0;  J3d := 0;
       Location := '@  Three D.';       Progress_Line;
       REPEAT
          Initialize_ThreeD_Variables;
          IF (ABS( E )+ABS( F )) > 1 THEN Triclinic_Lattice
             ELSE  MonoClinic_Lattice;
          FOR H := Kst to Kg DO Find_Most_Frequent_Values;
          IF Isom[1] >= Min THEN   Best_Values;
       UNTIL  (J3d + 4) > LinCombo;    WRITELN;
    END;                     {   Three D  }

              {  generate  Q  values for lattice   }

PROCEDURE Refine_Parameters;
    VAR
       Druk, Prin, List2               :  BOOLEAN;
       Nrr, N20, Ncalc, NotInd, Nmax   :  INTEGER;

PROCEDURE  Refine( IrefCy  :  INTEGER );
    VAR
       Instor                               :  ARRAY[1..200,1..3] of INTEGER;
       Nmatch                               :  ARRAY[1..45]       of INTEGER;
       Ibr, Hstor, Kstor, Lstor             :  ARRAY[1..41,1..5]  of INTEGER;
       Qstor                                :  ARRAY[1..40,1..5]  of REAL;
       Ref                                  :  ARRAY[1.. 8,1..8]  of REAL;
       QQstor                               :  ARRAY[1..200]      of REAL;
       Hulp, Vax, Fyn, War                  :  ARRAY[1.. 8]       of REAL;
       We                                   :  ARRAY[1..10]       of REAL;
       P, H, Count, Orde, Orp1, Orp2, Ibrav,
         J20, NrefCy, Indexd, Ispat, J, K, L, M,
         NrIndx, N                                 :  INTEGER;
       Thetc, T, Rms, W, SumW, D, Q20, Qmax, Q,
         SumEps, Eps, Qobs, U, Zershf, Dif, C      :  REAL;

PROCEDURE  Initialize_Refine_Variables;
   VAR
      M  :  INTEGER;
   BEGIN
       FILLCHAR( Instor, SIZEOF( Instor ), 0 );
       FILLCHAR( Qqstor, SIZEOF( Qqstor ), 0 );
       FILLCHAR(  qstor, SIZEOF(  qstor ), 0 );
       FILLCHAR(  Hstor, SIZEOF(  Hstor ), 0 );
       FILLCHAR(  Kstor, SIZEOF(  Kstor ), 0 );
       FILLCHAR(  Lstor, SIZEOF(  Lstor ), 0 );
       FOR M := 1 to maxObsLns DO  ShfTheta[m] := ObsTheta[m] - ZerTot;
       Ibrav := ROUND( Cel[9] );
        Nmax := 20;     J20 := 20;              Q20 := ObsPt[20];
       IF Soln  THEN   Nmax := MaxObsLns;      Qmax := ObsPt[nmax] + 1;
   END;              {  Initialize  Refine  Variables  }

PROCEDURE  Generate_Q;
   VAR
       H1, Hs, K1, L1, Lb, Lmax, Lmax1, Lmin, Hmin, Kmin,
         Hmax, Hmax1, Kmax, Kmax1                          :  INTEGER;
       Heven, Keven, Leven, Mis                            :  BOOLEAN;
       C2, C1, C3, C4, C5, Tol                             :  REAL;

PROCEDURE  Set_Up_K_Loop;
   VAR
      I, J  :  INTEGER;
   BEGIN
      Orde := 3;
      FOR J := 4 to 6 DO IF ABS( Cel[j] ) > FminEq THEN Orde := Orde + 1
             ELSE  Cel[j] := 0;
      Orp1  := Orde + 1;            Orp2 := Orde + 2;
      Ncalc := 0;          N20 := 0;      Ispat := 1;
      FILLCHAR( Ibr, SIZEOF( Ibr ), 0 );
      IF not List2 THEN FILLCHAR( Nmatch, SIZEOF( Nmatch ), 0 );
      TolRad := 0.5 * Tolg * Rrad;   Tol := (Irefcy - Nrefcy + 1) * TolRad;
      FOR I := 1 to Nmax DO BEGIN
          Thetapt[i] := Shftheta[i] + Tol;
          Thetamt[i] := Shftheta[i] - Tol;  END;
      Hmax := TRUNC( SQRT( Qmax/ Cel[1] ) ) + 1;
      Kmax := TRUNC( SQRT( Qmax/ Cel[2] ) ) + 1;
      Lmax := TRUNC( SQRT( Qmax/ Cel[3] ) ) + 1;
      Keven := False;   Lmin := - Lmax;  Hmin := 0;  Kmax1 := Kmax + 1;
      IF Orde = 3 THEN Lmin := 0;      IF Orde >= 5 THEN Hmin := -Hmax;
   END;                 {  Set Up K Loop   }

PROCEDURE  Set_Up_H_Loop;
   BEGIN
      K := K1 - 1;         C1 := SQR( K ) * Cel[2];      C2 := K * Cel[4];
     C3 := K * Cel[6];  Keven := NOT Keven;              Hs := Hmin;
     IF K = 0 THEN         Hs := 0;                   Heven := True;
     IF (2*(Hs DIV 2)) = Hs THEN   Heven := False;    Hmax1 := Hmax + 1 - Hs;
   END;                   {   set up H loop   }

PROCEDURE  Set_Up_L_Loop;
   BEGIN
       H := H1 + Hs - 1;         C4 := C1 + H * (C3 + H *Cel[1] );
      C5 := C2 + H *Cel[5];   Heven := NOT Heven;      Lb := Lmin;
       IF H = 0 THEN  IF (K = 0) OR (Orde = 4) THEN    Lb := 0;  Leven := True;
       IF 2*(Lb DIV 2) = Lb  THEN  Leven := False;      Lmax1 := Lmax + 1 - Lb;
   END;                 {   Set up  L  Loop  }

PROCEDURE  Store_Q_Values;
   BEGIN
      IF Ncalc <= 200 THEN BEGIN
          Instor[Ncalc,1] := H;  Instor[Ncalc,2] := K;
          Instor[Ncalc,3] := L;    Qqstor[Ncalc] := Q;   END;
   END;                  {   Store Q Values   }

PROCEDURE  Store_Q_Values_2( P :INTEGER );
   VAR
      I  :  INTEGER;

PROCEDURE  Bravais_Set;
   BEGIN
      IF Keven = Leven THEN                        Ibr[p,1] := Ibr[p,1] + 1;
      IF Heven = Leven THEN                        Ibr[p,2] := Ibr[p,2] + 1;
      IF Heven = Keven THEN                        Ibr[p,3] := Ibr[p,3] + 1;
      IF 2*((H+K+L) DIV 2) = (H+K+L) THEN          Ibr[p,4] := Ibr[p,4] + 1;
      IF (Heven = Keven) AND (Keven = Leven) THEN  Ibr[p,5] := Ibr[p,5] + 1;
   END;                  {  Bravais Set  }

   BEGIN
       NrIndx := NrIndx + 1;
       IF Prin  THEN BEGIN  WriSp( Ispat, Q, 1 );   Mis := False;  END;
       IF Cel[7] < 1 THEN   Bravais_Set;
       I := Nmatch[p] + 1;
       IF I <= 5 THEN  BEGIN     Nmatch[p] := I;
          IF NrIndx >= 2 THEN     L := L + 1000;
          Hstor[p,i] := H;      Kstor[p,i] := K;
          Lstor[p,i] := L;      Qstor[p,i] := Q;  END;
   END;                  {  Store Q Values in alternate array  }

   LABEL  700, 1020, 1030;
   BEGIN                                {  Generation of Q values   }
      GoToXY(1,Wherey);
      FOR k := 1 to 10 DO WRITE(cel[k]:7:0);
      Set_Up_K_Loop;                         holdy := wherey + 1;
      FOR K1 := 1 to Kmax1 DO BEGIN
         Set_Up_H_Loop;
         FOR H1 := 1 to Hmax1 DO BEGIN
            Set_Up_L_Loop;
            FOR L1 := 1 to Lmax1 DO BEGIN
               L := L1 - 1 + Lb;  Leven := NOT Leven;
               Q := C4 + L *( C5 + L * Cel[3] );
               IF Q < 1      THEN BEGIN
                   IF Prin   THEN  WriSp( Ispat, 0, 0 );       GoTo 1020;  END;
               IF Q >= Qmax  THEN                              GoTo 1020;
               IF Cel[7] < 1 THEN                              GoTo  700;
               CASE Ibrav OF
                  1 :  IF Keven = Leven THEN  GoTo  700 ELSE   GoTo 1020;
                  2 :  IF Heven = Leven THEN  GoTo  700 ELSE   GoTo 1020;
                  3 :  IF Heven = Keven THEN  GoTo  700 ELSE   GoTo 1020;
                  4 :  IF 2*((H+K+L) DIV 2)=(H+K+L) THEN
                                                 GoTo 700 ELSE GoTo 1020;
                  5 :  IF (Heven = Keven) AND (Keven = Leven)
                              THEN              GoTo  700 ELSE GoTo 1020;
                  6 :                                          GoTo  700;
                  7 :
                END;                                           GoTo 1020;
  700:        IF Q < Q20 THEN N20 := N20 + 1;         Ncalc := Ncalc + 1;
              IF List2 THEN BEGIN Store_Q_Values;              GoTo 1030; END;
              IF Prin  THEN  WriSp( Ispat, Q, 0 );
              Theta := Theta_Calculate( Q );   Mis := True;   NrIndx := 0;
              FOR J := 1 to NMax DO IF ThetaPt[j] >= Theta THEN BEGIN
                 IF ThetaMt[j] > Theta   THEN                  GoTo 1020;
                 Store_Q_Values_2( J );             END;
 1020:       IF Prin  THEN Ispat := Ispat + 1;
 1030:    END;
          IF Prin  THEN WriSp( Ispat, Q, 2 );
    END;  END;   END;            {   Generation of Q  Values   }

              {  identify matching calc with obs   }

PROCEDURE  Bravais_Lattice;

  PROCEDURE  Lines_Indexed;
     VAR
        N, I   :  INTEGER;
     LABEL  1220;
     BEGIN
        N := 0;
        FOR I := 1 to Nmax DO BEGIN
           IF Nmatch[i] >= 1 THEN    N := N + 1;
           J20 := I;
           IF N = 20 THEN                         GoTo 1220;  END;
 1220:  NotInd := J20 - N;      Indexd := N;
     END;          {    Number of Lines Indexed by this Lattice   }

  PROCEDURE  Bravais_Type;
    VAR
        J, N, M  :  INTEGER;
    BEGIN
        FOR J :=1 to 5 DO BEGIN           N := 0;
           FOR M := 1 to 20 DO IF Ibr[m,j] > 0 THEN N := N + 1;
           Ibr[21,j] := N;  END;
        IBrav := 6;
        FOR J := 1 to 5 DO IF Ibr[21,j] >= Indexd THEN Ibrav := J;
        Cel[9] := Ibrav;
    END;             {   Bravais   Type    }

    BEGIN                 {   Bravais  Lattice   }
        Lines_Indexed;
        IF Cel[7] < 1 THEN     Bravais_Type;
        Cel[7] := Indexd;
        Cel[8] := Indexd * Obs[J20] / (2 * n20 * SumEps);
        GoToXY(1,Wherey);
        FOR k := 1 to 10 DO WRITE(cel[k]:7:0);
    END;              {   Bravais  Lattice   }

 PROCEDURE  List_Matching_Lines;
   VAR
      Swavel, Q,  Ulast, U, W     :  REAL;
      K, No, Nc                   :  INTEGER;

 PROCEDURE   Arrange_For_Print;
   VAR
      I, J, N, K, Inx  :  INTEGER;
   BEGIN
      FOR  I := 1 to Ncalc DO BEGIN
         Q := 1E+6;
         FOR J := I to Ncalc DO
            IF Qqstor[j] < Q THEN BEGIN
               Q := Qqstor[j];   N := J; END;
         U := Qqstor[i];    Qqstor[i] := Qqstor[n];    Qqstor[n] := U;
         FOR K := 1 to 3 DO BEGIN
            Inx := Instor[i,k];  Instor[i,k] := Instor[n,k];
            Instor[n,k] := Inx;         END;        END;
      No := 1;        Nc := 1;      Ulast := 0;     Swavel := Lambda/200;
    WRITELN(Lstr, 'Observed and Calculated lines to',Nmax:4);
    WRITELN(Lstr, 'Blank line inserted when 2-thetas differ more than',
                   Tolg:6:3);
    WRITELN( Lstr ); Writeln(Lstr,'  TwoTheta       D       H   K   L      Q');
  END;                       {   Arrange for  Printing  }

   LABEL 110, 115, 200, Return;
   BEGIN                          {   List Matching Lines   }
      Arrange_For_Print;
 110: IF         No > Nmax     THEN               GoTo    115;
      IF         Nc > Ncalc    THEN               GoTo    200;
      IF Qqstor[nc] > Obs[no]  THEN               GoTo    200;
 115: IF         Nc > Ncalc    THEN               GoTo Return;
      IF Qqstor[nc] > 0 THEN BEGIN
        W := SQRT( Qqstor[nc] );
        U := Swavel * W;        W := 100 / w;
        U := TwoRadian * ARCTAN( U/SQRT( 1-U*U ));
        IF U-Ulast > TolG THEN       WRITELN( Lstr );
        WRITE(Lstr, U:10:2, W:10:3 );
        FOR K := 1 to 3 DO WRITE(Lstr, Instor[nc,k]:5 );
        WRITELN(Lstr, Qqstor[nc]:8:1);       Ulast := U;
        Nc := Nc + 1;     END;                        GoTo 110;
 200: Theta := ShfTheta[no];             U := TwoRadian * Theta;
     D := Lambda / (2 * SIN( Theta ));   Q := 100/D;  Q := SQR( Q );
     IF U-Ulast > Tolg THEN   WRITELN( Lstr );
     WRITE(Lstr, U:10:2, D:10:3, '     *Observed ',Q:8:1);
     IF Nmatch[no] = 0 THEN    WRITELN(Lstr, Spaces(5),'x.x.x.x.x')
        ELSE  WRITELN( Lstr );
     Ulast := U;      No := No + 1;                        GoTo 110;
 Return:  END;            {   List  Matching  lines  }

                 {  Refine lattice parameters         }
 
  PROCEDURE  Solve_For_Shifts;
    VAR
      Deter   :  REAL;
      I, J    :  INTEGER;

  PROCEDURE  Apply_Shifts;
    VAR
      N, M, I, J   :  INTEGER;
    BEGIN
      IF T > (Orp1+1) THEN    Rms := Rms / (T - Orp1);
      FOR I := 1 to Orp1 DO FOR J := 1 to Orp1 DO
                 FYn[i] := Fyn[i] + Ref[i,j] * War[j];
      FOR N := 1 to Orde DO Fyn[n] :=   1E+4 * Fyn[n];
      FOR N := 1 to 3    DO Cel[n] := Cel[n] + Fyn[n];
      IF Orde > 3 THEN      Cel[5] := Cel[5] + Fyn[4];
      IF Orde > 4 THEN      Cel[4] := Cel[4] + Fyn[5];
      IF Orde > 5 THEN      Cel[6] := Cel[6] + Fyn[6];
      Zershf := Fyn[orp1];  ZerTot := ZerTot + ZerShf;
      FOR  M := 1 to MaxObsLns DO   ShfTheta[m] := ShfTheta[m] - ZerShf;
    END;                      {  Apply Shifts  }

    BEGIN
       FOR I := 1 to 8 DO FOR J := 1 to 8 DO     ArrayD[i,j] := Ref[i,j];
       Arrange_Array_Single( Orp1, 8 );         MatInvert( Orp1, Deter );
       Arrange_Array_Double( Orp1, 8 );
       FOR I := 1 to 8 DO FOR J := 1 to 8 DO     Ref[i,j] := ArrayD[i,j];
       IF ABS( Deter ) > 1E-30 THEN    Apply_Shifts;
   END;     {   Solve  for   Shifts     }

 PROCEDURE  Initialize_Refine_Loop;
   BEGIN
       FILLCHAR( Ref, SIZEOF( Ref ), 0 );
       FILLCHAR( Fyn, SIZEOF( Fyn ), 0 );
       FILLCHAR( War, SIZEOF( War ), 0 );
       T := 0;     RMS := 0;    SumEps := 1E-3;     Count := 0;
       IF Druk  THEN WRITELN(Lstr, 'TwoTheta     D       Q       ',
                'H    K    L');
    END;                   {   Initialize  Refine   Loop   }

  PROCEDURE  Complete_Refinement( P : INTEGER );
    VAR
      N     :  INTEGER;
      W, D  :  REAL;

  PROCEDURE  Calculate_Refine_Variables( P, N : INTEGER );
    VAR
      I, J, K, L, H  :  INTEGER;
      Dubl           :  BOOLEAN;
      Q, D, U, W     :  REAL;

  PROCEDURE  Print_Line;
    BEGIN
        D := 100 / SQRT( Q);   U := TwoRadian * Thetc;
        WRITE(Lstr, U:8:2, D:8:3, Q:9:1, H:5, K:5, L:5);
        IF Dubl  THEN WRITELN( Lstr, '    HKL-Combination Used')
                    ELSE WRITELN( Lstr );
    END;                {  Print  Line  }

  PROCEDURE  Retrieve_Values;
    BEGIN
       Q := Qstor[p,n];       L := Lstor[p,n];    Dubl := False;
       H := Hstor[p,n];       K := Kstor[p,n];
       Thetc := Theta_Calculate( Q );
       IF L > 900 THEN BEGIN    L := L - 1000;   Dubl := True;   END;
    END;                  {   Retrieve  Values   }

  PROCEDURE   Vax_Array_Fill;
    BEGIN
       Vax[1] := SQR( H ) * C;  Vax[2] := SQR( K ) * C;
       Vax[3] := SQR( L ) * C;  Vax[4] := H * L * C;
       Vax[5] := K * L * C;     Vax[6] := H * K * C;   Vax[orp1] := 1;
    END;                {    Vax   Array   Fill   }

    BEGIN                           {  Calculate  Refine  Variables   }
       Retrieve_Values;
       IF Druk  THEN     Print_Line;
       Dif := ShfTheta[p] - Thetc;         Count := Count + 1;
       IF SumW <> 0 THEN   W := We[n] / SumW  ELSE   W := 100;
       Rms := Rms + W * SQR( Dif );        T := T + W;
       IF ThetC <> 0 THEN   C := WaveSqr / SIN( 2 * ThetC )
                     ELSE   C := 1;       Vax_Array_Fill;
       FOR J := 1 to Orp1 DO BEGIN     Hulp[j] := W * Vax[j];
           FOR I := 1 to Orp1 DO
                 Ref[i,j] := Ref[i,j] + Vax[i] * Hulp[j];
           War[j] := War[j] + Dif * Hulp[j];  END;
    END;                    {   Calculate  Refine  Variables   }

    BEGIN                                    {   Complete the Refinement  }
      Eps := 1000;      SumW := 0;
      FOR N := 1 to Nrr DO BEGIN
          D := ABS( Qobs - Qstor[p,n] );
          IF D < Eps THEN Eps := D;         IF D > 9 THEN D := 9;
          W := EXP( - SQR( D ) );  We[n] := W;  SumW := SumW + W;  END;
      IF P <= J20 THEN  SumEps := SumEps + Eps;
      FOR N := 1 to Nrr DO  Calculate_Refine_Variables( p,n );
    END;                     {    Complete  The   Refinement  }

   LABEL  1180, Return;
   BEGIN                                       {    R e f i n e   }
      Initialize_Refine_Variables;
      FOR NrefCy := 1 to IrefCy DO BEGIN
         Location := 'Generate Q values.';        Progress_Line;
         Generate_Q;
         IF List2  THEN BEGIN
             List_Matching_Lines;   GoTo Return;  END;
         Initialize_Refine_Loop;
         FOR P := 1 to Nmax DO BEGIN         Nrr := Nmatch[p];
            W := SIN( ShfTheta[p] ) / HalfLambda;
            Qobs := SQR( W );
            IF Druk  THEN BEGIN   D := 100 / W;
               U := TwoRadian * ShfTheta[p];
               WRITE(Lstr, U:8:2, D:8:3, Qobs:9:1,'    *OBSERVED *');  END;
            IF Nrr = 0 THEN BEGIN
                 IF Druk THEN  WRITELN(Lstr, '     Not Indexed');
                                                          GoTo 1180;  END;
            IF Druk THEN   WRITELN( Lstr );
            Complete_Refinement( P );
 1180:  END;                         Solve_for_Shifts;
      END;                  {  Refinement Cycle  }
      Bravais_Lattice;
 Return:    END;            {  Refine   }

  PROCEDURE  Job_Title_Print;
    VAR
       Isp   :  INTEGER;
    BEGIN
       Isp := 81-LENGTH(  JobTitle );
       WRITELN( Lstr );       WRITELN( Lstr );
       WRITELN(Lstr, JobTitle, Spaces(isp),'  >>  Solution  Nr ',Nsol:3);
       WRITELN( Lstr );
    END;                     {  Print Job Title Name line  }

  PROCEDURE  Print_Result;
    VAR
       I  :  INTEGER;
    BEGIN
        WRITELN( Lstr);
        WRITELN(Lstr, 'Reciprocal Lattice Constants After Refinement');
        WRITELN( Lstr );
        FOR I := 1 to 6 DO WRITE(Lstr, cel[i]:9:2 );
        WRITELN( Lstr );   WRITELN( Lstr );
        WRITELN(Lstr, 'Lines Indexed ',Cel[7]:4:0,'     Figure of Merit ',
           Cel[8]:5:1,'      Lattice Type ',Cel[9]:4:0,'      History ',
                 Cel[10]:5:0);        WRITELN( Lstr );
     END;            {  Print Result of Second  Entry   }

  PROCEDURE  Second_Entry;
    VAR
       I  :  INTEGER;
    BEGIN
       Location := 'Second Refine.';                  Progress_Line;
       Refine(1);    IF enl THEN  Druk := True;       Job_Title_Print;
       FOR I := 1 to 10 DO WRITE(Lstr, Cel[i]:9:2);   WRITELN(Lstr);
       WRITELN( Lstr );              Refine(1);       Print_Result;
    END;               {   Second Entry for Refinement   }

  PROCEDURE  Print_Final_Refine;
    BEGIN
      Location := '....Final Refine.';     Progress_Line;
      WRITELN(Lstr,'FINAL pass for  - - - - >   ',JobTitle);   WRITELN( Lstr );
      WRITELN(Lstr,'    The best solution is used to try and index all lines.');
      Refine( 1 );
      WRITELN(Lstr, ' Figure of Merit ', Cel[8]:7:1,
                   '.   Unindexed lines in first twenty indexed is ',NotInd:3);
      WRITELN(Lstr, ' Lines calculated up to 20th indexed line is ', N20:3);
      WRITELN(Lstr,'   Total lines Calculated is ',Ncalc:4);   WRITELN( Lstr );
      IF Ncalc <= 4*Nmax THEN  List2 := True;    Druk := True;    Refine(1);
    END;         {   Print  Final  Refine  }

    LABEL  Return;
    BEGIN                          {   Refine  Parameters  }
       Druk := False;     Prin := False;     List2 := False;
       IF Cel[7] <= 1 THEN BEGIN
          Location := 'Initial Refine.';           Progress_Line;
          Refine(3);                                      GoTo Return;  END;
       IF Soln = False THEN BEGIN   Second_Entry;
          IF Nsol > List THEN                             GoTo Return;
          IF N20 <= 100  THEN BEGIN   List2 := True;    Prin := False;
            Refine(1);      Print_Result;                 GoTo Return;   END;
       END;           Print_Final_Refine;
    Return:  END;                    {   Refine Parameters  }

              {  normalizing of lattice            }
 
  PROCEDURE  Possible_Lattices_Print;
    VAR
       N, I, J  :  INTEGER;
    BEGIN
     WRITE(Lstr, Spaces(12) );
     WRITELN(Lstr,'Possible Lattices based on the above combination of zones');
     WRITELN(Lstr,'Serial    A       B       C       D       E       F',
                              '       Hits  History');
       N4 := LinCombo + 4;  N := 0;
       FOR I := 1 to N4 DO IF Roster[i,7] >= 1 THEN BEGIN
           N := N + 1;
           FOR J := 1 to 8 DO Roster[n,j] := Roster[i,j];
           WRITE(Lstr, N:5);
           FOR J := 1 to 6 DO WRITE(Lstr, Roster[n,j]:8:1 );
           FOR J := 7 to 8 DO WRITE(Lstr, Roster[n,j]:8:0 );
           WRITELN(lstr);       END;
       N4 := N;
    END;               {   Possilbe  Lattices  Print   }

  PROCEDURE  Complete_Lattice_Search;
    VAR
       U        :  REAL;
       I, J, K  :  INTEGER;
    LABEL  Return, 1990;
    BEGIN
       Location := 'Intersect Zones.';  Progress_Line;
       Threed;                Possible_Lattices_Print;
       Location := 'Normalize Zones';   Progress_Line;
       ScreenLine := 'Index of Q values for Lattice Search.';
       ScreenWrite( 1 );    WRITELN;
       IF N4 = 0 THEN BEGIN
          ScreenLine := '> > >    No lattices have been Found.    < < <';
          ScreenWrite( 1 );   ScreenWrite( 11 );
          Mj := -999;  END  ELSE
       FOR J := 1 to N4 DO BEGIN         WRITE(J:4);
          FOR I := 1 to 8 DO  Cel[i] := Roster[j,i];          Normalize3;
          FOR I := 1 to 8 DO  Roster[j,i] := Cel[i];
          IF J <> 1 THEN BEGIN          Mj := j - 1;
             FOR K := 1 to Mj DO BEGIN       U := 0;
                FOR I := 1 to 6 DO  U := U + ABS( Roster[k,i] - Roster[j,i] );
                IF U > Fmineq THEN                                  GoTo 1990;
                IF Roster[k,7] < Roster[j,7] THEN BEGIN
                           Roster[k,7] := 0;  Roster[k,9] := Roster[j,8];  END
                ELSE BEGIN Roster[j,7] := 0;
                           Roster[j,9] := Roster[k,8];  END;  END;   END;
 1990: END;     WRITELN;
    END;                  {   Complete the Lattice Search   }

  PROCEDURE  Normalize_and_Refine_Lattice;
    VAR
       Old       :  ARRAY[1..6] of REAL;
        U        :  REAL;
        I, J, K  :  INTEGER;
    LABEL  2040, 2080, 2100;
    BEGIN
       FOR J := 1 to N4 DO BEGIN
          IF Roster[j,1] < Fmineq THEN                         GoTo 2100;
          IF Roster[j,7] < Fmineq THEN                         GoTo 2100;
          FOR I := 1 to 10 DO   Cel[i] := Roster[j,i]; Cel[10] := Cel[8];
 2040:    Cel[7] := 0;     ZerTot := 0;     Refine_Parameters;
          IF Cel[7] < Println     THEN                         GoTo 2080;
          IF Cel[8] < Printmerit  THEN                         GoTo 2080;
          FOR I := 1 to 6 DO Old[i] := Cel[i];      Normalize3;
          IF Cel[1] < Fminim THEN                              GoTo 2080;
          U := 0;
          FOR I := 1 to 6 DO  U := U + ABS( Cel[i] - Old[i] );
          IF U > 0.1 THEN                                      GoTo 2040;
 2080:    FOR K := 1 to 10 DO  Roster[j,k] := Cel[k];  Roster[j,11] := ZerTot;
 2100: END;
    END;                {    Normalize and Refine Lattices   }

  PROCEDURE  Discard_Lattices;
    VAR
       Message        :  STRING[60];
       ip, J, K       :  INTEGER;

  PROCEDURE  Print_Discard;
    BEGIN
       WRITE(Lstr, J:5);
       FOR Ip := 1 to 10 DO  WRITE(Lstr, Roster[j,ip]:8:1 );
       WRITELN( Lstr,'  --  ', Message );
    END;    {  Print Discarded Lattices   }

  PROCEDURE  ReducedLatticeHead;
    BEGIN
        WRITELN(Lstr, ' Lattices reduced and put into standard form');
        WRITELN(Lstr, Spaces(57),'Lines  Figure    Bravais');
        WRITELN(Lstr,'Serial   A       B       C       D       E       F',
              '       Indxd  of Merit  Lattice  History'); WRITELN(Lstr);
    END;            {   Print Reduced Lattice Header Lines   }

    LABEL  2200;
    BEGIN                                  {  Discard Lattices  }
       Nr := 0;
       IF Enl THEN  ReducedLatticeHead;
       FOR J := 1 to N4 DO BEGIN
           IF Roster[j,1] <= Fmineq THEN BEGIN
               IF Enl THEN BEGIN  Message:='Discarded by Normalize3';
               Print_Discard;  END;                            GOTO 2200;  END;
           IF Roster[j,7] <= Fmineq  THEN BEGIN
               IF Enl THEN BEGIN  Message:='Same as another Lattice';
               Print_Discard; END;                             GoTo 2200;  END;
           IF Roster[j,7] <= Println THEN BEGIN
              IF Enl THEN BEGIN  Message :='Not Enough Lines Indexed.';
              Print_Discard;  END;                             GoTo 2200;  END;
           IF Roster[j,8] <= Printmerit THEN BEGIN
              IF Enl THEN BEGIN  Message :='Figure of Merit too Low.';
              Print_Discard;   END;                            GoTo 2200;  END;
           Nr := Nr + 1;
           FOR K := 1 to 12 DO Roster[nr,k] := Roster[j,k];
           IF Enl THEN BEGIN   Message := '< < <  Lattice  Accepted.  ! ! ! !';
              Print_Discard;        END;
 2200: END;    WRITELN( Lstr );
    END;    {   Discard  Lattices   }

              {  controller for zone refinements   }
 
  PROCEDURE  Roster_Write( L :  INTEGER );
    VAR
       N  :  INTEGER;
    BEGIN
       Theta := TwoRadian * Roster[l,11]; Kkk := ROUND( Roster[l,9] );
       FOR N := 1 to 8 DO WRITE(Lstr, Roster[l,n]:8:1 );
       WRITELN(Lstr,'    ',BravType[kkk]:2,' ',Roster[l,10]:8:0, Theta:10:4);
    END;               {   Roster  Write   }

  PROCEDURE  Roster_Write_Head;
    BEGIN
       WRITELN(Lstr, Spaces(50), 'Lines   Figure  Bravais');
       WRITELN(Lstr,'  Q(A)    Q(B)    Q(C)     Q(D)    Q(E)    Q(F)   Indexed',
                   '  Merit  Lattice History ZeroShift');
    END;                {  Roster  Write  Header   }

  PROCEDURE  Sort_Lines;
    VAR
       Jj, M, N, L, K, H, I  :  INTEGER;
       AA, B, A              :  REAL;
    BEGIN
       WRITELN( Lstr );
       WRITELN(Lstr, 'Trial Lattices After Indexing and LS Refinement');
       Roster_Write_Head;
       Jj := 100;
       FOR L := 1 to Nr DO BEGIN        K := 0;  Aa := 0;
          FOR H := l to Nr DO BEGIN
             M := ROUND( Roster[h,7] );  B := Roster[h,8];
             IF (M-K) > 0 THEN BEGIN
                K := M;  I := H;  Aa := B;  END
              ELSE IF (M-K) = 0 THEN   IF B > Aa THEN BEGIN
                  Aa := B;  I := H;  END;  END;
          IF K < Jj THEN BEGIN     WRITELN( Lstr );  Jj := K;   END;
          FOR N := 1 to 12 DO BEGIN
             A := Roster[i,n];     Roster[i,n] := Roster[l,n];
             Roster[l,n] := A;
             IF N <= 10 THEN Cel[n] := A;  END;
          UnitCell( 0 );  Roster_Write( L );       END;
          WRITELN( Lstr );
    END;                   {   Sort Lines    }

  PROCEDURE  Crystal_System;
     VAR
       Dif          :  REAL;
       M, N, I, K   :  INTEGER;
    BEGIN
       N := 0;   Min := - 1;  Self := True;
       FOR I := 1 to 3 DO IF Nsyst[i] <> 0 THEN N := I;
       IF N <> 0 THEN Self := False;
       IF Self = False THEN BEGIN      N := 0;   M := 0;
          FOR K := 1 to 3 DO BEGIN
             IF Nsyst[k] =  1 THEN N := N + 1;
             IF Nsyst[k] = -1 THEN M := M + 1;   END;
          IF N > 0 THEN Min := 1;
          IF M > 0 THEN Min := 0;
          FOR M := 1 to MaxObsLns DO BEGIN
             Dif := Obs[m] - Obsmt[m];  Obsmt[m] := Obsmt[m] - Dif;
             Obspt[m] := Obspt[m] + Dif;   END;      END;
    END;                    {   Crystal  System   }

  PROCEDURE  Solutions;
    VAR
       Norde, I, J, K     :  INTEGER;
       A                  :  REAL;
    LABEL  3140, 3170, 3200, Return;
    BEGIN
       Nsol := 0;
       FOR I := 1 to NR DO BEGIN
           IF Self = False THEN BEGIN        Norde := 1;
              FOR J := 4 to 6 DO IF ABS( Roster[i,j] ) > Fmineq THEN
                  Norde := Norde + 1;
              IF Norde = 4  THEN Norde := 3;
              IF Nsyst[norde] < Min THEN                        GoTo 3200; END;
           Nsol := Nsol + 1;
           IF I = 1 THEN                                           GoTo 3140;
           A := 0;   K := I - 1;
           FOR J := 1 to 6 DO  A := A + ABS( Roster[k,j] - Roster[i,j] );
           IF A > Fmineq THEN                                      GoTo 3140;
           WRITELN(Lstr, 'Solution', I:2, '  is the same as ',K:2); GoTo 3170;
 3140:     FOR J := 1 to 10 DO  Cel[j] := Roster[i,j];
           ZerTot := Roster[i,11];   Refine_Parameters;
           FOR J := 1 to 10 DO Roster[i,j] := Cel[j];
           Roster[i,11] := ZerTot;
 3170:     IF Nsol = MaxSoln THEN GoTo Return;
 3200: END;
 Return:  END;                      {   Solutions   }

  PROCEDURE  Sort_it_Again;
    VAR
       I, J, L, N  :  INTEGER;
       Aa, A       :  REAL;
    BEGIN
       Nsol := MaxSoln;
       IF Nr < MaxSoln THEN Nsol := Nr;
       FOR I := 1 to Nsol DO  Roster[i,12] := Roster[i,8] / (21-Roster[i,7]);
       FOR I := 1 to Nsol DO BEGIN          Aa := 0;
          FOR J := I to Nsol DO  IF Roster[j,12] > Aa THEN BEGIN
                Aa := Roster[j,12];  L := J;  END;
          FOR N := 1 to 12 DO BEGIN          A := Roster[i,n];
              Roster[i,n] := Roster[l,n];  Roster[l,n] := A;   END;
    END;  END;                 {   Sort it Again}

  PROCEDURE  Lattice_Header;
    BEGIN
       WRITELN(Lstr );
       WRITELN(Lstr, Spaces(20),'The Direct Constants of these Lattices');
       WRITELN(Lstr,'     A        B        C       Alfa     Beta     Gamma',
                         '      Volume');       WRITELN( Lstr );
    END;                {   Lattice Constants Header Lines  }

  PROCEDURE  Bravais_Lattices;
    BEGIN
       BravType[1] := 'A';  BravType[2] := 'B';  BravType[3] := 'C';
       BravType[4] := 'I';  BravType[5] := 'F';  BravType[6] := 'P';
    END;              {   Bravais  Lattice  Types   }

  PROCEDURE  Most_Probable_Solution;
    VAR
       I, J   :  INTEGER;
    BEGIN
       Sort_Lines;
       IF Nr < Maxsoln THEN MaxSoln := Nr;
       WRITELN(Lstr, 'The ',MaxSoln:2,' best solutions follow.');
       Crystal_System;        Solutions;         Sort_it_Again;
       WRITELN(Lstr, 'The ',Nsol:2,' most probable solutions follow.');
       WRITELN(Lstr);           Roster_Write_Head;
       FOR I := 1 to Nsol DO    Roster_Write( i );
       Lattice_Header;
       Location := '*.*.*. Unit Cell .*';    Progress_Line;
       FOR I := 1 to Nsol DO BEGIN
          FOR J := 1 to 10 DO Cel[j] := Roster[i,j];    UnitCell( 1 );  END;
       WRITELN( Lstr );
       FOR I := 1 to Nsol DO BEGIN
          FOR J := 1 to 10 DO  Cel[j] := Roster[i,j];  Symmetry_Test;   END;
    END;                  {   Most  Probable  Solutions   }

  PROCEDURE  Refinement_End;
    VAR
       I  :  INTEGER;
    BEGIN
       Nrun := 2;      Soln := True;
       FOR I := 1 to 10 DO Cel[i] := Roster[1,i];
       ZerTot := Roster[1,11];    Refine_Parameters;
       Roster[1,11] := ZerTot;    Lattice_Header;      UnitCell( 1 );
    END;                   {   Refinement  End  }

    LABEL  Return;
    BEGIN
       Bravais_Lattices;
       IF Soln = False THEN BEGIN
           Complete_Lattice_Search;
           IF Mj = -999 THEN                              GoTo Return; END
        ELSE  N4 := Nr;
       Location := 'Normalize Lattice';       Progress_Line;
       ScreenLine := 'Calculations for Lattice with A, B, C, D, E, F :';
       Screen_Display_Header;  WRITELN;  ScreenWrite( 1 );  WRITELN;
       Soln := False;     Normalize_and_Refine_Lattice;     Discard_Lattices;
       IF Nr = 0 THEN BEGIN
          ScreenLine := 'Not one acceptable Lattice has been found.';
          ScreenWrite( 1 );  ScreenWrite( 11 );  WRITELN;
                                                          GoTo Return;  END;
      Location := 'Best Solutions.';       Progress_Line;
      Most_Probable_Solution;
      IF Roster[1,7] < 17.9 THEN                             GoTo Return;
      IF Roster[1,8] < 10.0 THEN                             GoTo Return;
      Refinement_End;
 Return:  END;                        {  Anito2  }

   PROCEDURE  Clear_the_Lstr_File;
      BEGIN
        WRITELN(Lstr);    WRITELN(Lstr, ChrWdt[Prnt,3], JobTitle );
        WRITELN(Lstr, ChrWdt[Prnt,3], '    * * * *     completed at ',
                     Date, '   . . .   ', time );
        CLOSE( Lstr );
      END;              {   Clear the Lstr File  }

   PROCEDURE  Clear_The_Screen;
      VAR
         I, Line  :  INTEGER;
      BEGIN
        ClrScr;  ScreenLine := PgmTitle;  ScreenWrite( 1 );
        WRITELN( #7 );  WRITELN( #7 );
        ScreenLine := JobTitle;  ScreenWrite( 1 );   WRITELN( #7 );
        WRITELN( #7 );   ScreenLine := '* * * *  Job Completed   * * * *';
        ScreenWrite( 1 );  WRITELN( #7 );  WRITELN( #7 );
        ScreenLine := Date + '       .. .. ..       ' + Time;
        ScreenWrite( 1 );  WRITELN( #7, #7 );  WRITELN( #7 );
        ScreenLine := 'Best  Solution  Real Lattice  Parameters';
        ScreenWrite(1);  WRITELN;
        FOR I := 1 to 3 DO BEGIN
          Line := 2 * i + 14;
          GoToXY(20,line);  WRITE( 'axis[', i:2, ' ] = ',Cel[i]:10:3 );
          GoToXY(50,line);  WRITE( 'angle[',i:2, ' ] =', Cel[i+3]:10:2 ); END;
     END;                  {  Clear the Screen  }

     BEGIN                     {  Main Programme  }
        Initiate_Visser_Algorithm;
        Read_In_The_Control_and_Diffraction_Data;
        IF Nr = 0 THEN BEGIN
           Soln := False;    Nrun := Nrun + 1;
           Crystallographic_Zone_Finding;   END;
        Crystallographic_Lattice_Finding;
        IF Nrun = 1 THEN BEGIN
            Nr := 0;    Soln := False;    Nrun := Nrun + 1;
            Crystallographic_Zone_Finding;
            Crystallographic_Lattice_Finding;         END;
        Clear_the_Lstr_File;             Clear_the_Screen;
     END.



--------------------------------------------------------------------------
END OF FILE