Title:   		LEAST SQUARES UNIT CELL REFINEMENT with INDEXING on the
         		PERSONAL COMPUTER          LSUCRIPC
KeyWords:   		Powder Diffraction, Unit Cell Refinement
Computer:               PC, XY, AT
Operating System:       MS(PC)-DOS
Programming Language:   TURBO Pascal
Hardware:   		192K memory, 1 diskette drive, math coproc.(recommended)
Author:                  Roy Garvey
Correspondence Address:  Department of Chemistry
                         North Dakota State University
                         Fargo,  ND  58105-5516
Abstract:

L S U C R I P C
                A  computer  program implementing the procedure of
           Appleman  and Evans (1973) is now available for the
           IBM-PC.   Output of refinement data has been condensed
           for economy of presentation.  For a discussion of
           methodology, see the USGS-GD-73-003 report number 20 by
           Appleman and Evans (1973).

---------------------------------------------------------------------------
Title:   		LEAST SQUARES UNIT CELL REFINEMENT with INDEXING on the
         		PERSONAL COMPUTER          LSUCRIPC
KeyWords:   		Powder Diffraction, Unit Cell Refinement
Computer:               PC, XY, AT
Operating System:       MS(PC)-DOS
Programming Language:   TURBO Pascal
Hardware:   		192K memory, 1 diskette drive, math coproc.(recommended)
Author:                  Roy Garvey
Correspondence Address:  Department of Chemistry
                         North Dakota State University
                         Fargo,  ND  58105-5516
Documentation:
---------------------------------------------------------------------------

                      INSTRUCTIONS  for  using  LSUCRIPC

              D A T A      I N P U T       Definitions and Specifications

         Parameters are  input  with  spaces separating numeric fields as
    follows:

    FIRST LINE 
         Title describing the data being processed.

    SECOND LINE 
         Unit Cell Parameters as     A  B  C   Alpha  Beta   Gamma
         Control Parameters  Itht    = 0 for defraction data in 2-theta
                                     = 1 for defraction data in 1-theta

                             thtmx   Maximum  diffraction angle (in Itht
                                      units) to be considered in indexing.
                                      Set at 20 if 0 is coded. 

                             Ncyc   Number cycles during which non-indexed
                                    diffraction with angles exceeding ThtMx
                                    are not accepted in least squares.
                                    A zero value is reset to 2 by program.

                             TolMn  Minimum tolerance allowed (in units of
                                     Itht).  If 0 entered, value reset to
                                     0.02o by program.

                             TolMx  Maximum tolerance allowed (in units of
                                     Itht).  If no greater than TolMn, reset
                                     to five times TolMn.  If given as 0, a
                                     value of 0.1o set by program. 

    THIRD LINE Code for Crystallographic System
                   CUBIC  for cubic crystal systems.
                   TETRA  for tetragonal systems.
                   ORTHO  for orthorhombic systems.
                   MONOC  for monoclinic systems.
                   RHOMB  for rhombohedral or trigonal systems.
                   TRICL  for triclinic systems.
                   HEXAG  for hexagonal systems.
                   HEXAR  for hexagonal-rhombohedral systems.

         TheMx    Value comensurate with Itht to be used in determining the
                   minimum interplanar spacing for generated diffraction. 
                   This value may be superceded by given diffraction having 
                   greater Bragg angles.  If 0 and there are no observed 
                   diffraction given, TheMx is set to 90o theta. 

         Pwl     Primary wave length used in computing minimum distance for
                   generater diffraction.  If given as 0,  value of 1.5405 A
                   is assumed.

          The next 27 entries indicate conditions for non-extinction for
    certain diffraction classes.   To be effective,  the integer n entered
    must be 2,3 4 or  6 in the appropriate location;  x,  y,  and z are any
    integers.   If the crystal  type  is HEXAR,  entry 6 will be considered
    to be 3 whether  it  is entered or not by the user.   Each of the
    following 27 items must be  repre sented  by  a  numeric value as
    described above or by 0 as a  place  holder. Each entry must be
    separated form those preceding and following by at  least one space.
                                Class of           Condition for
                 Entry        Diffraction          non-extinction
                   1              hkl                h + k = x.n
                   2              hkl                h + l = x.n
                   3              hkl                k + l = x.n
                   4              hkl       h+k = x.n, h+l = x.n, k+l = x.n
                   5              hkl            h + k + l = x.n
                   6              hkl           -h + k + l = x.n
                   7              hhl                    h = x.n
                   8              hhl                    l = x.n
                   9              hhl                h + l = x.n
                  10              hhl               2h + l = x.n
                  11              0kl                    k = x.n
                  12              0kl                    l = x.n
                  13              0kl                k + l = x.n
                  14              h0l                    h = x.n
                  15              h0l                    l = x.n
                  16              h0l                h + l = x.n
                  17              hk0                    h = x.n
                  18              hk0                    k = x.n
                  19              hk0                h + k = x.n
                  20              hh0                    h = x.n
                  21              h00                    h = x.n
                  22              0k0                    k = x.n
                  23              00l                    l = x.n
                  24              hll                    h = x.n
                  25              hll                    l = x.n
                  26              hll                h + l = x.n
                  27              hll               h + 2l = x.n

              SEmult   Multiple of the standard error of an observation of
                        unit weight to be used in setting tolerances.

    FOURTH LINES
         Diffraction  data.   For  diffraction  that are to be  indexed  by
    the program  by  comparing  corresponding  distances and angles  with 
    those  of theoretical diffraction, h, k, and l must be given as zero. 
    If one or more of  the  diffraction  indices  is  non-zero,   the 
    indices  read  from  the diffraction  entry  will  be maintained 
    throughout  the  run.   Diffraction entries  must be arranged in order
    of increasing Bragg angle (or  decreasing d-spacing).   If  no observed
    diffraction data are given,  the  output  will consist of theoretical
    lines only.  As before, numeric entries are separated by  spaces.   0 
    is entered as a place holder where no data is  input.   Each diffraction
    is represented by the following entries:

              h  k  l   D    Wlc   Obs     Wt

    where h,  k,  l are the corresponding Miller indices.   D is the
    interplanar spacing as computed from the Bragg relationship.  D need not
    be given if OBS is  given.   Wlc  is  the  wavelength (in A)  associated
    with  the  current diffraction.   If given as 0, it is assumed to be the
    same as the wavelength for  the  previous  diffraction  line.   (If  a 
    wavelength  for  the  first diffraction line of a data set is not given,
    the value given or implied  by Pwl  is  used).   Obs is the Bragg angle
    concomitant with Itht and given  as degrees  and fraction thereof.   The
    angle need not be given if D is  given. If  both  theta and D are given,
    D will be recomputed as  the  observed  d- spacing from theta and the
    given or implied wavelength.  WT is the weight to be associated with the
    value of theta for the current diffraction.  If WT is given as 0,  it is
    set to 1.0 by the program.   All entries are set to 0 for the last data
    line.

         A  sample  data  set is included illustrating the use of 0 as  a 
    place holder.

    Cell Refinement of A - ZrPO4    SNOWMASS  WORKSHOP   test data file.
    9.06  5.295 15.414 90 101.75 90 0 0 0 0 0
    MONOC  73  0.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0
     0  0  2 0.0 1.5405981 11.7087 0.0
    -1  1  1  0   0        19.8162  0
     2  0  0  0   0        20.0178  0
     0  1  2  0   0        20.4939  0
     2  0  1  0   0        22.072   0
    -1  1  3  0   0        24.9746  0
     2  0  2  0   0        25.2411  0
    -2  1  2  0   0        27.0218  0
     1  1  3  0   0        27.7465  0
     0  1  4  0   0        29.0032  0
    -3  0  1  0   0        29.5801  0
    -1  1  5  0   0        33.8194  0
     0  2  0  0   0        33.8903  0
    -3  1  1  0   0        34.1718  0
    -3  1  2  0   0        34.655   0
     0  2  2  0   0        35.9646  0
    -3  1  3  0   0        36.1483  0
     3  1  1  0   0        36.301   0
    -2  0  6  0   0        37.2433  0
     1  1  5  0   0        37.3402  0
     0  2  3  0   0        38.4033  0
    -3  1  4  0   0        38.572   0
     3  1  2  0   0        38.7756  0
    -2  2  2  0   0        40.2341  0
     4  0  0  0   0        40.6483  0
    -2  1  6  0   0        41.0746  0
     0  2  4  0   0        41.6506  0
     3  1  3  0   0        42.0089  0
     2  2  2  0   0        42.6968  0
    -4  0  4  0   0        42.8435  0
    -4  1  1  0   0        43.5333  0
    -1  1  7  0   0        44.4877  0
     4  0  2  0   0        44.8181  0
    -2  0  8  0   0        48.2941  0
     1  1  7  0   0        48.4369  0
    -1  3  1  0   0        53.0391  0
    -4  2  2  0   0        53.3706  0
    -5  1  3  0   0        54.0075  0
    -1  2  7  0   0        54.1143  0
    -1  3  3  0   0        55.4481  0
    -4  2  4  0   0        55.7147  0
     1  2  7  0   0        57.5746  0
    -3  2  7  0   0        58.9394  0
     2  1  8  0   0        59.4316  0
    -2  2  8  0   0        60.2872  0
     3  1  7  0   0        60.4909  0
    -5  1  7  0   0        63.1216  0
     1  2  9  0   0        68.7137  0
    -4  1 10  0   0        70.0417  0
     0  4  0  0   0        71.2907  0
    -6  2  2  0   0        71.9717  0
     0  0  0  0   0         0.0     0

         Specific  parameters are identified in the following lines of input
    data. Note that spaces are used as "field Delimiters".

    TITLE LINE
    Cell Refinement of A - ZrPO4    SNOWMASS  WORKSHOP   test data file.

    UNIT CELL PARAMETERS                                    Tolerance
      A     B     C     aLpha   beta   gamma   Itht thtmx  ncyc   mn   mx
    9.06  5.295 15.414   90     101.75  90      0     0     0      0    0
        T
        h
        e        Conditions for non-extinction for certain classes
        m                        1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2
Sysext  x  Pwl 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 Semult

MONOC  73  0.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0

     DIFFRACTION DATA
      h  k  l  D    Wlc      Angle   Weight

      0  0  2 0.0 1.5405981 11.7087   0.0
     -1  1  1  0       0    19.8162    0
      *  *  *  *       *       *       *

     TERMINATOR of H K L DATA INPUT
      0  0  0  0       0       0       0

ERROR  MESSAGES

    Run Time Error Messages that are Fatal with Program Aborted:

      01    Floating point overflow.
      02    Division by zero attempted.
      03    SQRT argument error.
      04    LN argument error.
      10    String length error.
      11    Invalid string index.
      90    Index out of range.
      91    Scalar or subrange out of range.
      92    Out of integer range.
      F0    Overlay file not found.
      FF    Heap / stack collision.

    I / O  Error Messages when I/O checking is active:

      01    File does not exist.
      02    File not opened for input.
      03    File not opened for output.
      04    File not opened.
      10    Error in numeric format.
      20    Operation not allowed on a logical devise.
      21    Not allowed in direct mode.
      22    Assign to STANDARD files not allowed.
      90    Record length mismatch.
      91    Seek beyond end-of-file.
      99    Unexpected end-of-file.
      F0    Disk write error.
      F1    Directory is full.
      F2    File sixe overflow.
      F3    Too many open files.
      FF    File disappeared.


  Special Version of  LSUCRIPC  includes the following new features:

             Ability to pass  < InFile.dat >  and/or  < OutFile.lst >
                file names on the command line.
                   < InFile.dat >  is the file of data and control codes
                used in earlier versions of the program (and described
                in the  LSUCRIPC.doc  file).
                   < OutFile.lst >  is an optional text file to which all
                printed information can be directed.  No longer is it
                necessary to wait for your printer to catch up to the
                compution.
                  If the < OutFile.lst > is to be used, you must also
                include the < InFile.dat > name on the command line as
                the first entry.  If no < OutFile.lst > name is given,
                the program defaults to the system  LST: (PRN: or LPT:)
                device.
                   To read the data file  TEST.dat from drive B: and
                write the printed listing to a file  TEST.out on a
                RAM disk C:  you would type
                           LSUCRIPC  B:TEST.dat  C:TEST.out

            Printer initialization.
                   In the past we have been limiting users to those
                printers using control codes compatible with either
                those of  OKIDATA Microline  or those of the IBM-
                EPSON  series of printers.  For this "Conference"
                version we have included a third option.  By providing
                a  PRINTER.$pr  "Printer Profile" with the decimal
                integer equivalent for the character codes of your
                printer needed to initialize
                     COMPRESSED       EXPANDED          NORMAL
                printing modes, the program will accomodate any
                of a variety of additional printers.  These codes
                are written to the corresponding < OutFile.lst >
                so these files may be printed as desired at a later
                time and retain the text formating of the original.
                   The  PRINTER.$pr  must have three lines of FIVE
                integer numbers (separated by spaces) in a format
                as illustrated below.

           27  48 15  0  0           Compressed Print         Epson Codes
           14   0  0  0  0           Expanded Print
           18   0  0  0  0           Compressed off

-------------------------------------------------------------------------
End of File
Title:   		LEAST SQUARES UNIT CELL REFINEMENT with INDEXING on the
         		PERSONAL COMPUTER          LSUCRIPC
KeyWords:   		Powder Diffraction, Unit Cell Refinement
Computer:               PC, XY, AT
Operating System:       MS(PC)-DOS
Programming Language:   TURBO Pascal
Hardware:   		192K memory, 1 diskette drive, math coproc.(recommended)
Author:                  Roy Garvey
Correspondence Address:  Department of Chemistry
                         North Dakota State University
                         Fargo,  ND  58105-5516
Source Code:
---------------------------------------------------------------------------

PROGRAM Least_Square_Unit_Cell_Refinement_and_Indexing;
   { Based upon algorithm  of Appleman and Evans (1973), USGS-GD-73-003 }
   { "COMMON" Global variables and file definitions }
TYPE
    AbcInt = ARRAY[1..3] of INTEGER;
    AbcXyz = ARRAY[1..3] of REAL;
  String14 = STRING[14];
   DateStr = string[10];
   TimeStr = STRING[5];
   regpack = record
              ax,bx,cx,dx,bp,si,ds,es,flags: integer;
            END;
CONST
     Version  : STRING[14] = 'Fargo  87.08.01';
     Ialf     : ARRAY[1..11] of STRING[5] = ('CUBIC','TETRA','ORTHO','MONOC',
                    'RHOMB','TRICL','HEXAG','HEXAR','  2- ','     ','  R  ');
     Iclass   : ARRAY[1..10] of STRING[3] = ('hkl','hhl','0kl','h0l','hk0',
                    'hh0','h00','0k0','00l','hll');
     Icond    : ARRAY[1..10] of STRING[6] = ('  h+k','  h+l','  k+l',' h+k+l',
                    '-h+k+l','  h','  k','  l',' 2h+l',' h+2l');
     Jclass   : ARRAY[1..27] of INTEGER = (1,1,1,1,1,1,2,2,2,2,3,3,3,4,4,4,5,
                    5,5,6,7,8,9,10,10,10,10);
     Jcond    : ARRAY[1..27] of INTEGER = (1,2,3,0,4,5,6,8,2,9,7,8,3,6,8,2,6,
                    7,1,6,6,7,8,6,8,2,10);
     Ihkl     : ARRAY[1..6,1..3] of INTEGER =  (( 1,-1, 0),( 0, 1,-1),
                    ( 1, 0,-1),( 1, 0, 0),( 0, 1, 0),( 0, 0, 1));
     Jhkl     : ARRAY[1..10,1..3] of INTEGER = (( 1, 1, 0),( 1, 0, 1),
                    ( 0, 1, 1),( 1, 1, 1),(-1, 1, 1),( 1, 0, 0),( 0, 1, 0),
                    ( 0, 0, 1),( 2, 0, 1),( 1, 0, 2));
     Ifi      : ARRAY[1..2,1..27] of INTEGER
                  = ((0,0,0,0,0,0,1,1,1,1,4,4,4,5,5,5,6,6,6,1,5,4,4,2,2,2,2),
                     (0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,6,6,5,0,0,0,0));
     Ifj      : ARRAY[1..3,1..27] of INTEGER
                  = ((1,2,3,1,4,5,6,8,2,9,7,8,3,6,8,2,6,7,1,6,6,7,8,6,8,2,10),
                     (0,0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
                     (0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0));
     Test     : ARRAY[1..6,1..7] of INTEGER = ((8,3,1,1,3,1,8),
                     (0,7,2,2,10,2,9),(0,0,3,3,0,3,0),(0,0,0,6,0,4,0),
                     (0,0,0,0,0,5,0), (0,0,0,0,0,6,0));
     HklC     : ARRAY[1..6,1..10] of INTEGER = ((1,2,3,1,2,3,1,1,1,1),
                     (1,2,3,2,3,1,1,1,2,1),(0,0,0,0,0,0,1,2,2,2),
                     (0,0,0,0,0,0,2,2,3,2),(0,0,0,0,0,0,2,3,3,0),
                     (0,0,0,0,0,0,2,3,1,0));
VAR
     Crystal,                    {  Describes the Crystal System            }
     Ialf12, Ja                                          : STRING[5];
     FileName, FileOut                                   : String14;
     N,                          {  Number of diffraction lines observed    }
       Ityp, Jtyp,
     Itht,                       {  Indicator of Theta angle type           }
     Nopt,                       {  Number of Extinction Parameters         }
     Icycle,                     {  Refinement Cycle Counter                }
       Imn,
     Ncyc,                       {  Number of cycles for refinement         }
     dx, cx, M, Irsd, MMn, Nr, Mm,
     Prnt,                       {  Defines printer Type for Initiation     }
     Err,                        {  Global  Error  code                     }
     ScreenX, ScreenY            {  Cursor position for Active Screen       }
                        : INTEGER;
     Vr, Vro,                    {  Reciprocal Cell Volumes                 }
     Vc, Vco,                    {  Direct Cell Volumes                     }
     Q,                          {  1 / d**2                                }
       Fac, Thcmp, D, Consx, Swt, Diff,
     Dobs,                       {  Observed  d  spacing value              }
     Pwl,                        {  Pattern wavelength                      }
     Thtse,                      {  Standard Error for Theta                }
     Thobs,                      {  Observed theta value                    }
       RsdMax, Sum,
     Dmin,                       {  Minimum value of  d  for data           }
     Tolmn,                      {  Tolerances for assigning given line     }
     Tol,                        {    to correspond to a given theoretical  }
     Tolmx,                      {    line calculated in refinement.        }
       Thtmx,
     Themx,
     Semult,
       Vrd, Vcd            : REAL;
     Title                 : STRING[80];
     HKL                   : ABCINT;
     Ao,                         {  Original Unit Cell Distances             }
     Ad,                         {  Direct Cell edge Distances               }
     Angles,                     {  Direct cell Angles in Degrees            }
     Alo,                        {  Direct Cell Angles in minutes            }
     Als,                        {  Saved Values of Direct Cell Angles       }
     Csa,                        {  Cosine of Direct Cell Angles             }
     Sna,                        {  Sines of Direct Cell Angles              }
     X,                          {  Reciprocal Cell edge Distances           }
     Xo,                         {  Saved Values of Reciprocal Cell Edges    }
     Ar,                         {  Angles of the Reciprocal Cell            }
     U,                          {  Cosines of Reciprocal Cell Angles        }
     W,                          {  Sines of the Reciprocal Cell Angles      }
     Uo,                         {  Saved values of Reciprocal Cell Cosines  }
     Dd                    : ABCXYZ;
     LSUCR, LSTR           : TEXT;
     Printer               : BOOLEAN;
     Answr                 : CHAR;
     Ifle, Iop, Jop        : ARRAY[1..27] of INTEGER;
     Ofle, A               : ARRAY[1..27] of REAL;
       B, Cr,
     Sedce,                      {  Standard Errors for Direct Cell          }
     Serce,                      {  Standard Errors for Reciprocal Cell      }
       Z                   : ARRAY[1..6] of REAL;
     Aa, E                 : ARRAY[1..7] of REAL;
     TT, St, Vcme          : ARRAY[1..6,1..6] of REAL;
     Ss                    : ARRAY[1..21] of REAL;
     Dcm                   : ARRAY[1..700] of REAL;
     HKLh                        {  Indices for generated lines              }
                           : ARRAY[1..700,1..3] of INTEGER;
     Obsth,                      {  Observed Experimental theta values       }
     Sw,                         {  Observed Statistical Weight for lines    }
     Wl,                         {  Wavelength of radiation for line         }
     Dsp,                        {  d Spacing for observed line              }
       Dlst                : ARRAY[1..200] of REAL;
     HKLg                        {  Indices for Observed lines               }
                           : ARRAY[1..200,1..3] of INTEGER;
     IoCalcD               : ARRAY[1..200] of INTEGER;
     ChrWdt                : ARRAY[1..3,1..3] of STRING[5];

{ $$$$$       *******  Utility Procedures and Functions  *******      $$$$$$ }

PROCEDURE InitializePrinter;
     VAR
        PrnFile   :  TEXT;

   PROCEDURE  NoPRINTERfile;
      BEGIN
         SOUND( 525 );   DELAY( 750 );   NOSOUND;
         WRITELN('Initialization File  PRINTER.$pr  not Found.');
         ChrWdt[1,1] := Chr(29);     ChrWdt[2,1] := Chr(15);  {compress}
         ChrWdt[1,2] := Chr(31);     ChrWdt[2,2] := Chr(14);  {Expanded}
         WRITELN;
         Write('Select Printer as :   1 = Okidata,  2 = Epson   ');
         READLN(Prnt);
         CASE Prnt OF
               1 : Writeln(LSTR,#24,#27,#56,#31);
               2 : Writeln(LSTR,#24,#27,#48,#14);    END
      END;        {     No  PRINTER  Profile  on Default Drive   }

   PROCEDURE  PRINTERfile;
      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        {   Normal       }
                 READ(PrnFile, Num );
                 ChrWdt[3,i] := ChrWdt[3,i] + CHR( Num );  END;
             READLN( PrnFile );     END;
          CLOSE( PrnFile );   Prnt := 3;
      END;                     {   PRINTER  file found on Default Diskette  }

    BEGIN      {  Initialize Printer  }
       ASSIGN( PrnFile,'PRINTER.$pr' );  {$I-} RESET( PrnFile ); {$I+}
       Err := IOResult;
       IF Err <> 0 THEN   NoPRINTERfile    ELSE    PRINTERfile;
       Writeln(LSTR, ChrWdt[prnt,2], 'Least Squares Unit Cell Refinement');
       Write(LSTR, ChrWdt[prnt,1], ' ':4,'N D S U   version ', Version);
       Writeln(LSTR,' ':6,'after Appleman and Evans (1973).');
       Writeln(LSTR,' ':5,'implementation by',' ':11,'Roy G Garvey');
       Writeln(LSTR,' ':27,'Department  of  CHEMISTRY');
       Writeln(LSTR,' ':23,'NORTH  DAKOTA  STATE  UNIVERSITY');
       Writeln(LSTR,' ':27,'Fargo        58105 - 5516');
       Writeln(LSTR);     WRITELN;
       Write(     'Print Intermediate Reflection Tables ?    Y / N ');
       REPEAT
          Read(Kbd,Answr);  Answr := UpCase( Answr );
       UNTIL Answr in ['Y','N'];
       IF Answr = 'Y' THEN   Printer:= True;  WRITELN;
    END;  {  Initialize 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;
       Hours, Minutes:        string[2];
   BEGIN
       with recpack do    ax := $2c shl 8;
       MsDos(recpack);
       with recpack do  BEGIN
           str(cx MOD 256, Minutes);    str(cx shr 8, Hours);    END;
       Time := Hours + ':' + Minutes;
   END;     {     T i m e    }

PROCEDURE NameOfSystem;
     BEGIN
          ClrScr; GotoXY(20,1);
          Writeln('Least  Squares  Unit  Cell  Refinement');
          GotoXY(29,3);  Writeln('N D S U    version  ', Version);
          GoToXY(65,4);  Writeln('Time = ', Time);
     END;  {  Name of System  }

 PROCEDURE ActiveScreen( ValueToPrint : INTEGER);
    BEGIN
        GotoXY( ScreenX, ScreenY );
        Write( ValueToPrint );
    END;

 FUNCTION Theta( C, D : REAL) : REAL;
    VAR
        Factor  : REAL;
    BEGIN
       Factor := d * d - c * c;
       IF (Factor >= 0) AND (d*d+c*c > 0) THEN
           Theta := 57.29578 * ArcTan( c / sqrt( Factor ))
             ELSE  Theta := 90;
    END;   {  Theta  }

  PROCEDURE  CommandTail;
     BEGIN
          FileName := ParamStr( 1 );    FileOut := ParamStr( 2 );
          IF  FileOut = '' THEN    FileOut := 'LST:';
          WRITELN('  The following information was obtained from the ',
               '"Command Line".');
          WRITELN(' ':8,'LSUCRIPC  <Input File>  <Output File>'); WRITELN;
          Writeln(' InPut = ',FileName,'     OutPut = ',FileOut); WRITELN;
          ASSIGN( LSTR, FileOut );
          IF FileOut = 'LST:' THEN   RESET( LSTR )
                            ELSE   REWRITE( LSTR );
     END;                {    Parse  the  Command  Tail     }

{*****       *****  Data Input Procedures   *****       *****    }

PROCEDURE SymmetryConditions;
   VAR
      L, Icry : INTEGER;

PROCEDURE Four;
   BEGIN
      Ao[2] := Ao[1];
   END;

PROCEDURE Five;
   BEGIN
      FOR L := 1 to 3 DO Angles[l] := 90.0;
   END;  {  Five  }

PROCEDURE Fifteen;
   BEGIN
      Ityp := 2;    N := 2;     Four;
      Five;      Angles[3] := 120.00;
   END;

   BEGIN
      Writeln('S Y M M E T R Y    C O N D I T I O N S');
      Ityp := 0;
      FOR L := 1 to 8 DO IF Crystal = Ialf[l] THEN      Icry := l;
      CASE ICry OF
        1 : BEGIN   Jtyp := 1;   Ityp := 1;  N := 1;
                   Ao[3] := Ao[1];    Four;    Five;    END;
        2 : BEGIN   Jtyp := 5;   Ityp := 5;  N := 2;  Four;   Five;  END;
        3 : BEGIN   Jtyp := 3;   Ityp := 3;  N := 3;  Five;   END;
        4 : BEGIN   Jtyp := 4;   Ityp := 4;  N := 4;
                   Angles[1] := 90.00;   Angles[3] := 90.00;  END;
        5 : BEGIN   Jtyp := 7;   Ityp := 7;  N := 2;
                   FOR L := 2 to 3 DO BEGIN
                       Ao[l] := Ao[1];  Angles[l] := Angles[1]; END; END;
        6 : BEGIN   Ityp := 6;   Jtyp :=6;   N := 6;   END;
        7 : BEGIN   Jtyp := 2;   Fifteen;       END;
        8 : BEGIN   Jtyp := 8;   Fifteen;       END;
      END;  {case  }
    END;             {  Symmetry Conditions Imposed on Direct Cell  }

PROCEDURE InputTheData;
   VAR
       L : INTEGER;

    PROCEDURE  SourceFile;
       BEGIN
          WRITELN( #7 );  Write('Name of Source DATA FILE  >  ', #7 );
          Readln( FileName );
       END;

    PROCEDURE  Does_File_Exist;
       BEGIN
          IF FileName = '' THEN   SourceFile;
          ASSIGN(LSUCR,FileName);    {$I-}   RESET(LSUCR);  {$I+}
          Err := IOResult;
          IF Err <> 0 THEN BEGIN
             WRITELN('  File  ', FileName, 'was not found.');
             SourceFile;  Does_File_Exist; END;
       END;

    PROCEDURE  Parameter_Error;
       BEGIN
          Writeln('Given Parameters in Error');
          FOR L := 1 to 3 DO   Write( Ao[l]:10:5 );          Write('     ');
          FOR L := 1 to 3 DO   Write( Angles[l]:12:4 );    WRITELN;
          Halt;
       END;         {    Parameter  Errors  Detected for Cell   }

   BEGIN
      Does_File_Exist;      Readln(LSUCR, Title);       Write(LSTR, Title);
      Write(LSTR,'    --  >  ',Date,'  . * . * .     ',Time);   Writeln(LSTR);
      FOR L := 1 to 3 DO   Read(LSUCR, Ao[l]);
      FOR L := 1 to 3 DO   Read(LSUCR, Angles[l]);
      Readln(LSUCR, Itht, Thtmx, Ncyc, Tolmn, Tolmx);
      Read(LSUCR, Crystal, Themx, Pwl);
      FOR L := 1 to 27 DO   Read(LSUCR, Iop[l]);
      Readln(LSUCR, Semult);      SymmetryConditions;
      IF (Ityp < 1) OR (Ityp > 7) THEN BEGIN
          Writeln('Crystal System on SYSext given as ',Crystal);
          Halt;      END;
      Sum := 0;
      FOR L := 1 to 3 DO BEGIN
          IF Ao[l] <= 0 THEN   Parameter_Error;
          Alo[l] := Angles[l] * 60.00;
          IF (Alo[l] <= 0) OR (Alo[l] >= 10800) THEN  Parameter_Error;
          Sum := Sum + Alo[l];
          IF Sum > 21599 THEN   Parameter_Error;   END;
   END;                    {  Input the Data  }

 {  @@@@@@      ******   Direct and Reciprocal Cell   ******      @@@@@@  }

PROCEDURE CheckExtinctions(VAR Ih,Ik,Il : INTEGER);
    VAR
        I, J, K, L, Idiv, Ifij, Ifjk, Ians, Iarg : INTEGER;
    LABEL
        12,22,24;
    BEGIN
        Err := 0;
        FOR I := 1 to Nopt DO BEGIN
           L := Jop[i];      Idiv := Iop[l];
           FOR J := 1 to 2 DO BEGIN
              Ifij := Ifi[j,l];
              IF (Ifij < 1) OR (Ifij > 6) THEN                      Goto 12;
              Ians := Ih*Ihkl[ifij,1] + Ik*Ihkl[ifij,2] + Il*Ihkl[Ifij,3];
              IF Ians <> 0 THEN                                     Goto 22;
           END;
       12: FOR K := 1 to 3 DO BEGIN
              IFjk := Ifj[k,l];
              IF (Ifjk < 1) OR (Ifjk > 10) THEN                     GOto 22;
              Iarg := Jhkl[ifjk,1]*Ih + Jhkl[ifjk,2]*Ik + Jhkl[ifjk,3]*Il;
              IF (Iarg MOD Idiv) <> 0 THEN BEGIN
                 Err := 1;                              Goto 24;     END;  END;
    22:  END;
    24:
    END; { Check Extinction }

PROCEDURE OptionSet;
    VAR
        L : INTEGER;
    BEGIN
         Nopt := 0;
         IF Jtyp = 8 THEN Iop[6] := 3;
         FOR L := 1 to 27 DO BEGIN;
              IF Iop[l] IN [2,3,4,6] THEN BEGIN
                  Nopt := Nopt + 1;  Jop[Nopt] := l;    END;     END;
    END;   {  Option Set  }

PROCEDURE ReciprocalCellElements;
     VAR
         I, J, K : INTEGER;
     BEGIN
         FOR I:= 1 to 3 DO BEGIN
             J     := (i MOD 3 ) + 1;  K := (j MOD 3) +1;
             U[i]  := (Csa[j] * Csa[k] - Csa[i]) / (Sna[j] * Sna[k]);
             W[i]  := sqrt(1.0 - U[i] * U[i]);
             X[j]  := 1.0 / (ad[j] * W[i] * Sna[k]);
             Ar[i] := 3437.7468 * (1.5707963 - ArcTan(u[i] / w[i]));  END;
         Vr := x[1] * x[2] * x[3] * Sna[1] * w[2] * w[3];
         Vc := 1 / Vr;
     END;   {  Reciprocal Cell Elements  }

PROCEDURE PrintCellParameters( ABC, Angle : ABCXYZ;  Vol : REAL);
    VAR
        Arg : REAL;
        I   : INTEGER;
    BEGIN
        FOR I := 1 to 3 DO Write(Lstr,ABC[i]:15:7);
        FOR I := 1 to 3 DO BEGIN
           Arg :=  Angle[i]/60;           Write(Lstr,Arg:12:5);        END;
        IF ABC[1] > 1 THEN  Writeln(Lstr,Vol:18:4,' A**3')
             ELSE  WRITELN(Lstr, Vol:18:9,' A**-3');
    END;   {  Print Cell Parameters  }

PROCEDURE DirectCellPrint;
     VAR
          L        : INTEGER;
          VolTest, Angle  :  REAL;
     BEGIN
         Writeln(Lstr);
         Write(Lstr,' ':5,'Cycle   ',Icycle:2,' ':12,'A',' ':14,'B');
         Write(Lstr,' ':14,'C',' ':11,'Alpha',' ':7);
         Writeln(Lstr,'Beta',' ':8,'Gamma',' ':10,'Volume'); Writeln(Lstr);
         Write(Lstr,'    Direct CELL: ');
         VolTest := Ad[1]*Ad[2]*Ad[3]*sna[2];
         WRITELN(' ':20,'Cycle ',Icycle:2,'Direct Cell Volume = ',
                    Voltest:15:5, Angle:10:4 );
         PrintCellParameters( Ad, Als, Vc );         Writeln(Lstr);
         Write(Lstr,'Reciprocal CELL: ');
         PrintCellParameters( X, Ar, Vr );           Writeln(Lstr);
     END;   {  Direct Cell Print  }

  {  ######      ======  h k l  Data Input and Processing  ======      ###### }

PROCEDURE ReadHKLData;
    VAR
        L                  : INTEGER;
        Cnst, Wlc, Obs, Wt : REAL;
    LABEL
        9, 16, 20, 25, 29, 45, 52, 64;
    BEGIN
        M := 0;
        WHILE NOT Eof(LSUCR) DO BEGIN
    9:     Readln(Lsucr, HKL[1], HKL[2], HKL[3], D, Wlc, Obs, Wt);
           IF Wlc > 0 THEN   Pwl := Wlc;
           IF Wt <= 0 THEN    Wt := 1;
           IF Obs > 0 THEN                                            Goto 20;
           IF D   > 0 THEN                                            Goto 25;
   16:     Write('Entry with indices: ');
           FOR L := 1 to 3 DO   Write(HKL[l]:5);
           Writeln('  D = ',D:10:6);
           Writeln('  Wl = ',Pwl:9:6,'  Obs Angle = ',
                   Obs:9:4,'  Wght = ',Wt:10:5,'  Not Processed');
           IF (HKL[1] = 0) AND (HKL[2] = 0) AND (HKL[3] = 0) THEN     Goto 64;
                                                                      Goto 9;
  20:      Thobs := Obs;   IF Itht = 0 THEN    Thobs := 0.5*Thobs;
           IF Thobs >= 90 THEN                                        Goto 16;
           Dobs  := 0.5*Pwl/Sin(0.017453293*Thobs);                   Goto 29;
  25:      Cnst  := 0.5*Pwl;       Thobs := Theta (Cnst,D);
           IF (Thobs <= 0) OR (Thobs >= 90) THEN                      Goto 16;
           Dobs  := D;
  29:      IF Nopt < 1 THEN                                           Goto 45;
           IF (HKL[1] = 0) AND (HKL[2] = 0) AND (HKL[3] = 0) THEN     Goto 45;
           CheckExtinctions(HKL[1],HKL[2],HKL[3]);
           IF Err = 0 THEN                                            Goto 45;
               Write('Reflection with H K L of  ');
               FOR L := 1 to 3 DO   Write(HKL[l]:5);     Writeln;
               Writeln(' not processed - Should not exist',
                        ' according to SYSext');                      Goto  9;
  45:      M := M + 1;
           IF M <= 200 THEN                                           Goto 52;
             Writeln('Over  200  Observations given  -  Only',
                      ' first 200 Used');           M := 200;         Goto 64;
  52:      Obsth[m] := Thobs;      Dsp[m] := Dobs;     Wl[m] := Pwl;
           Sw[m]    := sqrt(Wt);
           FOR L := 1 to 3 DO   HKLg[m,l] := HKL[l];      END;  {  while  }
  64:
   END;   {  Read  h k l  }

PROCEDURE InitializeParameters;
   VAR
       I, J, K, L, II : INTEGER;
       Arad, Themxo   : REAL;
   BEGIN
       Writeln('I N I T I A L I Z E     P A R A M E T E R S');
       Itht := Itht MOD 2;
       IF Itht = 0 THEN BEGIN
          Thtmx := 0.5*Thtmx;       Tolmn := 0.5*Tolmn;
          Tolmx := 0.5*Tolmx;       Themx := 0.5*Themx;       END;
       IF Thtmx <= 0     THEN    Thtmx := 20;
       IF Ncyc  <  1     THEN     Ncyc := 2;
       IF Tolmn <= 0     THEN    Tolmn := 0.02;
       IF Tolmx <= Tolmn THEN    Tolmx := 5*Tolmn;
       IF Pwl   <= 0     THEN      Pwl := 1.5405;
       IF Semult <= 0    THEN   Semult := 2;
       Ialf12 := Ialf[9+Itht];        OptionSet;
       Imn    := Round(((N+3)*N)/2);           {  from SYMMETRY  }
       FOR L:= 1 to 3 DO BEGIN
          Ad[l]  := Ao[l];  Als[l] := Alo[l];     Arad := 0.00029088821*Alo[l];
          Csa[l] := Cos(Arad);       Sna[l] := Sin(Arad);       END;
       ReciprocalCellElements;           DirectCellPrint;
       FOR L := 1 to 3 DO BEGIN   Uo[l] := U[l];     Xo[l] := X[l];   END;
       Vro := Vr;    Vco := Vc;     ReadHKLData;        Mm := M;
       IF M < 1 THEN   Themxo := 0
                ELSE   Themxo := Obsth[m] + Tolmx;
       IF Themxo > Themx THEN    Themx := Themxo;
       IF (Themx <= 0) OR (Themx > 90) THEN Themx := 90;
       Dmin := 0.5*Pwl/Sin(0.017453293*Themx);       Writeln(Lstr);
       Writeln(Lstr,'1-Theta Angles  Thtmx = ',Thtmx:6:1,'  Ncyc = ',
             Ncyc:2,'  Tolmn =',Tolmn:7:4,'  Tolmx =',Tolmx:7:4,'  Themx =',
             Themx:6:1,'  Dmin =',Dmin:10:6,'   ',Crystal,' System');
       Write(Lstr,Nopt:2,' Conditions for non-Extinctions Called for');
       Writeln(Lstr,'     Class Condition(s)');
       IF Nopt > 0 THEN   BEGIN
           FOR L := 1 to nopt DO BEGIN
               J := Jop[l];   I := Jclass[j];    K := Jcond[j];
               WRITE(Lstr, ' ':50);
              IF K IN [1..10] THEN
                 Writeln(Lstr, Iclass[i],Icond[k],' = ',Iop[j]:1,'n')
              ELSE   Writeln(Lstr, 'hkl     h+k= ',Iop[j]:1,'n, h+l= ',
                                   Iop[j]:1,'n, k+l= ',Iop[j]:1,'n');   END;
          END;       Writeln(Lstr);
    END;           {  Initialize Parameters }

{    END  Data Input.    Begin Generation of Calculated Diffraction Data   }

 {  %%%%%%      ------    Generate Q values for h k l     ------      %%%%%% }

PROCEDURE Skweze;
    VAR
       Ih, Jh                            : ARRAY[1..3] of INTEGER;
       I, J, K, L, Jmi, J2, J1, K2, Ieql : INTEGER;
       Tol                               : REAL;

PROCEDURE TestHKL;
    VAR
        I, Intry, J, K, L, Isum : INTEGER;
    LABEL
        5, 20, 32, 34;
    BEGIN
       FOR I := 1 to 6 DO BEGIN
           Intry := Test[i,ityp];
           IF (Intry < 1) OR (Intry > 10) THEN             Goto 32;
           Isum := 0;         J := 1;
     5:    K := Hklc[j,intry];
           IF (K < 1) OR (K > 3) THEN                      Goto 20;
           L := Hklc[j+1,intry];      J    := J + 2;
           Isum := Isum + Ih[k]*Ih[l] - Jh[k]*Jh[l];
           IF J < 6 THEN                                   Goto  5;
    20:    IF Isum <> 0 THEN BEGIN   Ieql := 1;            Goto 34;
       END;       END;
    32:  Ieql   := 0;
    34:
    END;                           {  Test H K L  }

    LABEL
        12, 19, 27, 35, 40, 50;
    BEGIN                              { Skweze }
        Write('S K W E Z E   Procedure   J =');
        ScreenX := WhereX + 2;   ScreenY := WhereY;     Tol := 1E-5;
        I := 1;
  12:   J := I + 1;
        IF J > Nr THEN                                          Goto 50;
        IF (Dcm[i] - Dcm[j]) >= Tol THEN                        Goto 40;
        FOR K := 1 to 3 DO BEGIN
            Ih[k] := HKLh[i,k];   Jh[k] := HKLh[j,k];    END;
  19:   TestHKL;        ActiveScreen(j);
        IF Ieql <> 0 THEN                                       Goto 27;
        J := J + 1;
        IF J > Nr THEN                                          Goto 27;
        IF (Dcm[i] - Dcm[j]) >= TOl THEN                        Goto 27;
        FOR K := 1 to 3 DO Jh[k] := HKLh[j,k];                  Goto 19;
  27:   Jmi := J - I - 1;   J2  := Nr - Jmi;  J1  := i + 1;
        IF J2 < J1 THEN                                         Goto 35;
        FOR K := J1 to J2 DO BEGIN
            K2 := K + Jmi;        Dcm[k] := Dcm[k2];
            FOR L := 1 to 3 DO HKLh[k,l] := HKLh[k2,l];   END;
  35:   Nr := J2;
  40:   I  := I + 1;                                            Goto 12;
  50:   Writeln;
    END;   {  Skweze  }

PROCEDURE Qcomputation(HKLr : AbcInt);
     BEGIN
         Q := (Z[1]*HKLr[1] + Z[4]*HKLr[2] + Z[6]*HKLr[3])*HKLr[1]
             +(Z[4]*HKLr[1] + Z[2]*HKLr[2] + Z[5]*HKLr[3])*HKLr[2]
             +(Z[6]*HKLr[1] + Z[5]*HKLr[2] + Z[3]*HKLr[3])*HKLr[3];
     END;                      {  Q Computation  }

PROCEDURE SortD;
     VAR
          J, K, I : INTEGER;
     BEGIN
         Write('S O R T D     Nr = ');
         ScreenX := WhereX + 5;   ScreenY := WhereY;    I := Nr;
         WHILE I > 1 DO BEGIN
            I := I-1;
            FOR J := 1 to I DO BEGIN
               IF Dcm[j] < Dcm[j+1] THEN BEGIN
                  Fac := Dcm[j];         Dcm[j] := Dcm[j+1];
                  Dcm[j+1] := Fac;
                  FOR K:= 1 to 3 DO BEGIN
                        HKL[k]   := HKLh[j,k]; HKLH[j,k] := HKLh[J+1,k];
                        HKLh[J+1,k] := HKL[k];               END;
                   ActiveScreen(j);  END;   END;   END;   Writeln;
     END;                        {  Sort D spacings  }

PROCEDURE GenerateMillerIndices;
   VAR
       I, J, K, L, Il, Ik, Ik1, Ik2, Ih, Ih1, Ih2 : INTEGER;
       Qmax, Td, Tn, Dscrs, Dscrm, Gt, Ht, Den    : REAL;
   LABEL
      15,31,38,40,42,44,46,50,56,57,58,60;
   BEGIN
      Write('M I L L E R     I N D I C E S   Procedure    Nr =');
      ScreenX := WhereX + 4;       ScreenY := WhereY;
      Qmax    := 1/(Dmin*Dmin);         Nr := 0;
      FOR I := 1 to 3 DO BEGIN
          J := (I MOD 3) + 1;       K := (J MOD 3) + 1;
          Z[i] := SQR( X[i] );
          Z[i+3] := X[i] * X[j] * U[k];    END;
      E[1] := Z[1];     E[2] := Z[2];     E[3] := Z[4];
      Td   := E[1] * E[2] - E[3] * E[3];    Il := 0;
  15: HKL[3] := Il;
      E[4] := (Z[5] + Z[5]) * HKL[3];
      E[5] := (Z[6] + Z[6]) * HKL[3];
      E[6] := Z[3]*HKL[3]*HKL[3] - Qmax;
      Tn := E[3]*E[5] - E[1]*E[4];
      Dscrs := Tn*Tn - (4*E[1]*E[6] - E[5]*E[5])*Td;
      ActiveScreen(Nr);
      IF Dscrs < 0 THEN                                Goto 60;
      Dscrm := Sqrt(Dscrs);
      Ik1 := Trunc((Tn - Dscrm) * 0.5 / Td);
      Ik2 := Trunc((Tn + Dscrm) * 0.5 / Td);
      CASE Jtyp OF
         1:          IF Ik1 < Il THEN Ik1 := Il;
         2,3,4,5,8 : IF Ik1 < 0  THEN Ik1 := 0;
         6,7:
      END;  { case }
      Ik := Ik2;
      IF Ik < Ik1 THEN                                  Goto 58;
  31: HKL[2] := Ik;
      GT := (E[3] + E[3])*HKL[2] + E[5];
      HT := (Z[5]*HKL[3] + Z[2]*HKL[2])*HKL[2] + (Z[5]*HKL[2] +
                   Z[3]*HKL[3])*HKL[3] - Qmax;
      Den := Z[1] + Z[1];
      Dscrs  := Gt*Gt - (Den+Den)*HT;
      IF Dscrs < 0 THEN                                 Goto 57;
      Dscrm := sqrt(Dscrs);
      Ih1 := Trunc(-(Gt+Dscrm)/Den);
      Ih2 := Trunc((Dscrm-GT)/Den);
      CASE Jtyp OF
         1,2,5 :  Goto 42;
         3     :  Goto 40;
         4,6,7 :  Goto 44;
         8     :  Goto 38;
     END;  { case }
  38: IF (Il MOD 3) = 0  THEN                           Goto 42;
  40: IF Ih1 < 0  THEN Ih1 := 0;                        Goto 44;
  42: IF Ih1 < Ik THEN Ih1 := Ik;
  44: Ih := Ih2;
      IF Ih < Ih1 THEN                                  Goto 57;
  46: HKL[1] := Ih;
      IF Nopt < 1   THEN                                Goto 50;
      CheckExtinctions(Ih, Ik, Il);
      IF Err  <> 0   THEN                               Goto 56;
  50: Qcomputation(HKL);
      IF Q <= 0   THEN                                  GOto 56;
      Nr := Nr + 1;           Dcm[Nr] := 1/sqrt(q);
      FOR J := 1 to 3 DO   HKLh[nr,j] := HKL[j];
      IF Nr < 700  THEN                                 Goto 56;
      SortD;      Skweze;
      IF Nr < 700  THEN                                 Goto 56;
      Writeln(Lstr,'More than  700  Reflections ',
                   'exist for Dmin = ',Dmin:10:5);      Goto 60;
  56: Ih := Ih - 1;      IF Ih >= Ih1 THEN              Goto 46;
  57: Ik := Ik - 1;      IF Ik >= Ik1 THEN              Goto 31;
  58: Il := Il + 1;                                     Goto 15;
  60: Writeln;
  END;  { Generate Miller Indices }

{   Begin Procedures of the  MAIN  COMPUTATION  CYCLE   * * * * * * * * * * }

 {  &&&&&&      ===     Generate D Spacings Data     ===      &&&&&&   }

PROCEDURE Coordinates;
    VAR
         J : INTEGER;
    BEGIN
        Writeln('C O O R D I N A T E S    Procedure');
        CASE Ityp OF
        1 : FOR J := 1 to 3 DO BEGIN   Cr[j] := B[1];  Cr[j+3] := 0;   END;
        2 : BEGIN
              FOR J := 1 to 2 DO   Cr[j] := B[1];        Cr[3] := B[2];
              FOR J := 4 to 6 DO   Cr[j] := 0;            END;
        3 : FOR J := 1 to 3 DO BEGIN   Cr[j] := B[j];  Cr[j+3] := 0;   END;
        4 : BEGIN
              FOR J := 1 to 3 DO   Cr[j] := B[j];
              Cr[4] := 0;          Cr[5] := B[4];        Cr[6] := 0;   END;
        5 : BEGIN
              FOR J := 1 to 2 DO   Cr[j] := B[1];        Cr[3] := B[2];
              FOR J := 4 to 6 DO   Cr[j] := 0;            END;
        6 : FOR J := 1 to 6 DO   Cr[j] := B[j];
        7 : FOR J := 1 to 3 DO BEGIN  Cr[j] := B[1];  Cr[j+3] := B[2];  END;
       END;  { CASE }
 END;                             {  Coordinates  }

PROCEDURE PositiveDefiniteInversion;
    VAR
        Ie, Np1, Np2, Inm, I, J, Ii, Ij, L, K : INTEGER;
    BEGIN
        Writeln(' ':13, 'Positive Definite Inversion');
        Ie := 0;      Np1 := N + 1;      Np2 := N + 2;
        Inm := Trunc(((N + 3)*N)/2) - 1;      Ii := 1;
        FOR I:= 1 to N DO BEGIN
            B[i] := A[ii];         J := I + 1;
            WHILE J <= N DO BEGIN
               Sum := 0;      Ij := Ii + J - I;      L := J - 1;
               FOR K := i to l DO BEGIN
                   Sum := Sum - A[ij]*B[k];  Ij := Ij + Np1 - k;   END;
               B[j] := A[ij]*Sum;        J := J + 1;   END;   { While }
           J := N;           Ij := Inm;
           REPEAT
                Sum := B[j];          K := N;
                WHILE K > J DO BEGIN
                   Sum := Sum - A[ij]*B[k];   Ij := Ij - 1;
                   K := K - 1;              END;  { while }
                B[j] := A[ij] * Sum;            J := J - 1;
                IF J >= I THEN Ij := Ij - 2;
           UNTIL J < I;
           FOR J := i to n DO BEGIN    Ie := Ie + 1;    A[ie] := B[j];   END;
           Ii := Ii + Np2 - I;       END;
    END;   {  Positive Definate Inversion  }

PROCEDURE DirectCellParameters;
     VAR
         I, J, K, L    : INTEGER;
         Angle, Dabc   : REAL;
     LABEL
         18, 19;
     BEGIN
        Vro := Vr;
        Vr  := Sqrt(1- U[1]*U[1]- U[2]*U[2]- U[3]*U[3]+
                    2*U[1]*U[2]*U[3]) *X[1]*X[2]*X[3];
        Vco := Vc;     Vc := 1.0 / Vr;
        FOR I := 1 to 3 DO BEGIN
                J := (I MOD 3) + 1;     K := (J MOD 3) + 1;
            Ao[i]  := Ad[i];    Ad[i] := X[j]*X[k]*W[i]/Vr;
            Alo[i] := Als[i];
            IF Ityp > 5 THEN                              Goto 18;
            IF (I = 2) AND (Ityp = 4) THEN                Goto 18;
            IF (I = 3) AND (Ityp = 2) THEN                Goto 18;
            Csa[i] := 0;                                  Goto 19;
    18:     Csa[i] := (U[j]*U[k] - U[i])/(W[j]*W[k]);
    19:     Sna[i] := Sqrt( 1 - Csa[i]*Csa[i] );
            Als[i] := 3437.7468 *
                      (1.5707963 - Arctan( Csa[i]/Sna[i] ));
            Ar[i] :=  3437.7468 *
                      (1.5707963 - Arctan( U[i] / W[i] ));  END;
            DirectCellPrint;            Vcd := Vc - Vco;
            Writeln(Lstr,'*--  Direct Cell     Corrections: ');
            WRITE(Lstr, ' ':17 );
            FOR L := 1 to 3 DO BEGIN
               Dabc := Ad[l] - Ao[l];    Write(Lstr, Dabc:15:7);   END;
            FOR L := 1 to 3 DO BEGIN
               Angle := (Als[l] - Alo[l]) / 60.00;
               Write(Lstr, Angle:12:5);       END;
            Writeln(Lstr, Vcd:18:6);   Writeln(Lstr);
    END;    {  Direct Cell Parameters  }

PROCEDURE SortReflections;
    VAR
        I, J, K : INTEGER;
        Wlc     : REAL;
    BEGIN
       Write('S O R T   R E F L E C T I O N S    I =');
       ScreenX := WhereX + 5;   ScreenY := WhereY;
       I := M;
       IF I > 1 THEN  REPEAT
          I := I - 1;
          FOR J := 1 to I DO IF Dsp[j] < Dsp[j+1] THEN BEGIN
             Dobs  := Dsp[j];     Dsp[j] := Dsp[j+1];     Dsp[j+1] := Dobs;
             Thobs := Obsth[j]; Obsth[j] := Obsth[j+1]; Obsth[j+1] := Thobs;
             Wlc   := Wl[j];       Wl[j] := Wl[j+1];       Wl[j+1] := Wlc;
             Swt   := Sw[j];       Sw[j] := Sw[j+1];       Sw[j+1] := Swt;
             FOR K := 1 to 3 DO BEGIN
                HKL[k] := HKLg[j,k];      HKLg[j,k] := HKLg[j+1,k];
                HKLg[j+1,k] := HKL[k];    END;     ActiveScreen(I);   END;
         UNTIL I = 1;
      Writeln;
   END;  {  Sort Reflections  }

PROCEDURE Dspacing;
     VAR
         I, J, K, Ic : INTEGER;
     BEGIN
         FOR I:= 1 to 3 DO BEGIN
             J := (i MOD 3)+1;             K := (j MOD 3)+1;
             Z[i] := X[i]*X[i];       Z[i+3] := X[i]*X[j]*U[k];       END;
         FOR Ic := 1 to Nr DO BEGIN
             FOR J:= 1 to 3 DO   HKL[j] := HKLh[ic,j];  Qcomputation(HKL);
             IF Q > 0 THEN      Dcm[ic] := 1/sqrt(q)
                      ELSE      Dcm[ic] := 0;            END;
     END;                     {  D Spacing  }

  {  $$$$$$     ***        Theoretical D Spacings       ***      $$$$$$   }

PROCEDURE StandErrors;
    VAR
        I, Iec, Ii, J, K, L, Ij12, Ij  : INTEGER;
        Factr, Angle                   : REAL;
    BEGIN
        FOR I := 1 to 21 DO   Ss[i] := 0;
        Factr := Thtse*Thtse;      PositiveDefiniteInversion;
        Iec   := Trunc(((N + 1)*N)/2);
        FOR I := 1 to Iec DO   A[i] := A[i]*Factr;
        CASE Ityp OF
          1 : BEGIN
               FOR I := 1 to 3 DO   Ss[i] := A[1];
               Ss[7] := A[1];   Ss[8] := A[1];  Ss[12] := A[1];    END;
       5, 2 : BEGIN
               Ss[1] := A[1];   Ss[2] := A[1];    Ss[3] := A[2];
               Ss[7] := A[1];   Ss[8] := A[2];   Ss[12] := A[3];   END;
          3 : BEGIN
               FOR I:= 1 to 3 DO   Ss[i] := A[i];
               Ss[12] := A[6];               END;
          4 : BEGIN
               FOR I := 1 to 3 DO       Ss[i] := A[i];    Ss[5]  := A[4];
               FOR I := 7 to 8 DO       Ss[i] := A[i-2];
               FOR I := 7 to 9 DO   Ss[2*i-4] := A[i];
                    Ss[19] := A[10];              END;
          6 : FOR I := 1 to 21 DO       Ss[i] := A[i];
          7 : BEGIN
               FOR I := 1 to 3 DO       Ss[i] := A[1];
               Ss[7] := A[1];           Ss[8] := A[1];    Ss[12] := A[1];
               FOR I :=  4 to  6 DO     Ss[i] := A[2];
               FOR I :=  9 to 11 DO     Ss[i] := A[2];
               FOR I := 13 to 15 DO     Ss[i] := A[2];
               FOR I := 16 to 21 DO     Ss[i] := A[3];               END;
         END;  {  case  }
       Ii := 1;
       FOR I := 1 to 3 DO BEGIN
          Serce[i] := sqrt(Ss[ii]);     Ii := Ii + 7 - I;     END;
       FOR I := 4 to 6 DO BEGIN
          Serce[i] := 3437.7468 * Sqrt(Ss[ii]) / W[i-3];
          Ii := Ii + 7 - I;       END;
       FOR I := 1 to 3 DO    Sna[i] := Sqrt(1-Csa[i]*Csa[i]);
       FOR I := 1 to 3 DO FOR J := 1 to 3 DO BEGIN
          IF I = J THEN     Tt[i,j] := 0
            ELSE            Tt[i,j] := -Ad[i]/X[i];     END;
       FOR I := 4 to 6 DO FOR J := 1 to 3 DO   Tt[i,j] := 0;
       FOR I := 1 to 3 DO BEGIN
           J := (I MOD 3) + 1;       K := (J MOD 3) + 1;
           Dd[i] := - Ad[j]*Ad[k]*X[j]*X[k];        END;
       FOR I:= 1 to 3 DO FOR J := 4 to 6 DO BEGIN
           IF (I-J+3) = 0 THEN BEGIN
              K := (J MOD 3) + 1;      L := (K MOD 3) + 1;
              Tt[i,j] := Dd[j-3] * Ad[j-3] * Csa[k] * Csa[l];
           END ELSE
              Tt[i,j] := Dd[j-3] * Ad[j-3] * Csa[k] * Csa[l];  END;
       FOR I := 1 to 3 DO BEGIN
           J := (I MOD 3) + 1;    K := (J MOD 3) + 1;
           Dd[i] := 1/(W[j]*W[k]);            END;
       FOR I := 4 to 6 DO FOR J := 4 to 6 DO BEGIN
          IF I = J THEN   Tt[i,j] := Dd[i-3] / Sna[j-3]
          ELSE BEGIN      Ij12 := 12-I-J;
             Tt[i,j] := Dd[i-3]*Csa[Ij12]/Sna[j-3]; END;  END;
       FOR I := 1 to 6 DO FOR J := 1 to 6 DO BEGIN
          Sum := 0;          Ij  := i;
          FOR K := 1 to 6 DO BEGIN
             SUm := Sum + Ss[ij]*Tt[j,k];
             IF I > K THEN   Ij := Ij+6-K
                 ELSE        Ij := Ij + 1;          END;
          St[i,j] := Sum;       END;
       FOR I := 1 to 6 DO  FOR J := 1 to 6 DO BEGIN
              Sum := 0;
              FOR K := 1 to 6 DO   Sum := Sum + Tt[i,k]*St[k,j];
              Vcme[i,j] := Sum;           END;
       FOR I := 1 to 6 DO  BEGIN
           IF Vcme[i,i] < 0 THEN Vcme[i,i] := 0;
           Sedce[i] := Sqrt(Vcme[i,i]);       END;
       FOR I := 4 to 6 DO   Sedce[i] := 3437.7468 * Sedce[i];
       Factr := Ad[1]*Ad[2]*Ad[3];       Factr := Factr*Factr/Vc;
       FOR I := 1 to 3 DO BEGIN
           J := (I MOD 3) +1;   K := (J MOD 3) + 1;  Cr[i] := Vc/Ad[i];
           Cr[i+3] := Factr*Sna[i]*(Csa[i]-Csa[j]*Csa[k]);       END;
       Factr := 0;
       FOR I := 1 to 6 DO FOR J := 1 to 6 DO
          Factr := Factr + Cr[i]*Cr[j]*Vcme[i,j];
       IF Factr > 0 THEN    Factr := Sqrt(Factr)
           ELSE             Factr := 0;
       Writeln(Lstr,'$--  Direct Cell   Standard   ERRorS ');
       WRITE(Lstr, ' ':17 );
       FOR I := 1 to 3 DO    Write(Lstr,Sedce[i]:15:7);
       FOR I := 4 to 6 DO  BEGIN
           Angle := Sedce[i] / 60.0;  Write(Lstr, Angle:12:5);   END;
       Writeln(Lstr,Factr:18:6);       Writeln(Lstr);
       WRITELN(Lstr,'#--  Reciprocal Cell    Standard  ERRorS');
       WRITE(Lstr, ' ':17);
       FOR I := 1 to 3 DO    Write(Lstr, Serce[i]:15:9);
       FOR I := 4 to 6 DO  BEGIN
           Angle := Serce[i] / 60.0;  Write(Lstr, Angle:12:7);  END;
       WRITELN(Lstr);         Writeln(Lstr);
    END;  { StandardError }

PROCEDURE TheoreticalReflections;
    VAR
        D1, D2, Wl1, Wt, Wl2, WlC,
          Wlchelp, Df, Thclc         : REAL;
        I2, Igp, Ihp, L              : INTEGER;
        IHKL, JHKL                   : ARRAY[1..3] of INTEGER;

PROCEDURE PrintLine2;
    BEGIN
        Write(Lstr, Ihp:3);
        FOR L := 1 to 3 DO    Write(Lstr,Ihkl[l]:5);
        Writeln(Lstr,'  ',D:13:6, Wlc:25:6, Thcmp:16:5);
    END;

PROCEDURE PrintLine1;
    BEGIN
      Write(Lstr,'***');
      FOR L := 1 to 3 DO   Write(Lstr,Jhkl[l]:5);
      Writeln(Lstr,Ja:4, Df:11:6, Dobs:13:6, Wlc:12:6, Thclc:16:5, Thobs:15:4,
           Diff:16:5, Wt:14:5);
    END;

    LABEL 20, 22, 32, 40, 60, 70, 75, 99;
    BEGIN                                     {  Theoretical Reflections  }
        Writeln(Lstr,Title,
         '        HKL Listing -   ***  Refers to FIXed,    R  to Rejects    ',
              Crystal:8);  Writeln(Lstr);
        Wl1 := Pwl;        I2  := 0;
        IF Mm > 0 THEN BEGIN
           Wl1 := Wl[1];    I2 := Abs(IoCalcD[1]);     D1 := Dsp[1];   END;
        Wl2 :=Wl1;
        If Nr > 0 THEN   D1 := Dcm[1];
        D2 := D1;        Igp := 1;        Ihp := 1;
        Write(Lstr,'  N    H    K    L       D  Calc        D  Obs      ');
        Write(Lstr,'Lambda   ',Ialf12:5,'Theta Calc ',Ialf12:5,'Theta Obs ');
        Writeln(Lstr,Ialf12:5,'Theta Diff      Weight');
    20: IF Ihp <= Nr THEN                                           Goto 22;
        IF Igp >  Mm THEN                                           Goto 99;
                                                                    Goto 40;
    22: FOR L := 1 to 3 DO   Ihkl[l] := HKLh[ihp,l];
        D   := Dcm[ihp];
        IF Ihp = I2 THEN                                            Goto 40;
        Wlc := Wl1;
        IF Abs(D-D2) < Abs(D1-D) THEN   Wlc := Wl2;
        Wlchelp := 0.5*Wlc;        Thcmp := Theta(Wlchelp,D);
        IF Itht = 0 THEN    Thcmp := Thcmp + Thcmp;
        PrintLine2;
    32: Ihp := Ihp + 1;                                             Goto 20;
    40: FOR L := 1 to 3 DO    Jhkl[l] := HKLg[igp,l];
        Df    := Dlst[igp];     Wt := Sw[igp]*Sw[igp];     Dobs := Dsp[igp];
        Thobs := Obsth[igp];   Wlc := Wl[igp];          Wlchelp := 0.5*Wlc;
        Thclc := Theta(Wlchelp,Df);                        Diff := Thclc-Thobs;
        IF Itht = 0 THEN BEGIN
            Thclc := Thclc + Thclc;      Thobs := Thobs + Thobs;
            Diff  := Diff + Diff;        END;
       Ja := Ialf[10];
       IF IoCalcD[igp] < 0 THEN   JA := Ialf[11];
       IF (Jhkl[1] <> 0) OR (Jhkl[2] <> 0) OR (Jhkl[3] <> 0) THEN   Goto 70;
       Write(Lstr, Ihp:3);
       FOR L := 1 to 3 DO   Write(Lstr, Ihkl[l]:5);
       Writeln(Lstr,Ja:5, Df:11:6, Dobs:13:6, Wlc:12:6, Thclc:16:5, Thobs:15:4,
          Diff:16:5, Wt:14:5);
   60: Igp := Igp + 1;    D1 := D2;     Wl1 := Wl2;
       IF Igp > Mm THEN   I2 := 0
          ELSE BEGIN
             IF Abs(IoCalcD[igp]) = I2 THEN                          Goto 40;
             I2 := Abs(IoCalcD[igp]);       D2 := Dcm[i2];
             Wl2 := Wl[igp];         END;  { else }                  Goto 32;
      Wlchelp := 0.5*Wlc;
  70: Thcmp   := Theta(Wlchelp,D);
      IF Itht = 0 THEN   Thcmp := Thcmp + Thcmp;
      IF Igp  < 2 THEN                                               Goto 75;
      IF Abs(IoCalcD[igp-1]) <> I2 THEN                              Goto 75;
      PrintLine1;                                                    Goto 60;
  75: PrintLine1;                                                    Goto 60;
  99: I2 := 1;
 END;  { Theoretical Reflections }

PROCEDURE SquareRootMethod;
    VAR
       I, J, K, L, Il, I1, I2, Ij, Ii : INTEGER;
       Sum, Ai                        : REAL;
    LABEL
       999;
    BEGIN
       I2 := N + 1;
       FOR I := 1 to N DO BEGIN
          J := N;      K := I;       Sum := 0;      Il := I-1;
          IF Il > 0 THEN FOR L := 1 to Il DO BEGIN
             Fac := A[k];    Sum := Sum - Fac*Fac;
             K := K + J;     J := J - 1;      END;
          Sum  := Sum + A[k];       Err := i;
          IF Sum <= 0 THEN   Goto 999;
          FAc  := 1 / Sqrt(Sum);        A[k] := Fac;     I1 := I + 1;
          FOR J := I1 to I2 DO BEGIN
             Sum := 0;       K := N;      Ij := j;       Ii := I;
             IF Il > 0 THEN FOR L := 1 to Il DO BEGIN
                Ai := A[ii];       Sum := Sum - Ai*A[ij];
                Ii := Ii + K;       Ij := Ij + K;      K := K - 1;     END;
             Sum := Sum + A[ij];     A[ij] := Fac*Sum;         END;    END;
       J := Trunc(((N+3)*N)/2);      I := N;
       WHILE I >= 1 DO BEGIN
           Sum := A[j];     J := J - 1;      K := I + 1;     L := N;
           WHILE L >= K  DO BEGIN
               Sum := Sum - A[j]*B[l];       J := J - 1;     L := L - 1;  END;
               B[i] := A[j]*Sum;             J := J - 1;     I := I - 1;  END;
       Err := 0;
  999:
  END;  { Square Root Method }

  {  ******      ======      Finish Computations     ======      ******   }

PROCEDURE HKLPrint;
    VAR
        Wt, Tc,Ob, Td  : REAL;
        L              : INTEGER;
     BEGIN
         Tc := Thcmp;      Wt := Swt*Swt;      Ob := Thobs;      Td := Diff;
         IF Itht = 0 THEN BEGIN
             Tc := Tc+Tc;        Ob := Ob+Ob;        Td := Td+Td;       END;
         Write(Lstr, M:4);
         FOR L:= 1 to 3 DO   Write(Lstr, HKL[l]:4);
         Writeln(Lstr, JA:4, D:11:6, Dobs:13:6,
                 Pwl:12:6, Tc:16:5, Ob:15:5, Td:16:5, Wt:14:5);
     END;                    {  H K L  Print  }

PROCEDURE DerivativesofQ;
     VAR
        J, K, L :INTEGER;
     BEGIN
        FOR J := 1 to 3 DO BEGIN
            K := (j MOD 3)+1;            L := (k MOD 3) +1;
            B[j] :=Fac*(X[j]*HKL[j]*HKL[j] + HKL[j]*HKL[l]*X[l]*U[k]
                    + HKL[j]*HKL[k]*X[k]*U[l]);
            B[j+3] := Fac*HKL[k]*HKL[l]*X[k]*X[l]        END;
     END;                     {  Derivatives of Q  }

PROCEDURE ObservationEquation;
    VAR
       J               : INTEGER;
       Divisor, Eqcns  : REAL;
    BEGIN
        M  := M + 1;        Divisor := Cos(0.017453293 * Thcmp);
        IF Divisor > 0 THEN BEGIN
             Fac := (57.29578 * D * Consx * Swt) / Divisor;
             Eqcns := Swt*Diff;      DerivativesOfQ;
             E[n+1] := -Eqcns;
             CASE Ityp OF
                1 : E[1] := B[1] + B[2] + B[3];
                2 : BEGIN   E[1] := B[1] + B[2];   E[2] := B[3];  END;
               3 : FOR J := 1 to 3 DO      E[j] := B[j];
               4 : BEGIN
                      FOR J := 1 to 3 DO   E[j] := B[j];
                      E[4] := B[5];                 END;
               5 : BEGIN
                      E[1] := B[1] + B[2];    E[2] := B[3];    END;
               6 : FOR J := 1 to 6 DO    E[j] := B[j];
               7 : FOR J := 1 to 2 DO    E[j] := B[3*j-2] + B[3*j-1] + B[3*j];
            END; { case }
         END ELSE    M := M - 1;
    END;   {  Observation  Equation  }

PROCEDURE NormalEquations;
     VAR
          Np1, J, K, L : INTEGER;
     BEGIN
         L := 0;         Np1 := n + 1;
         FOR J:= 1 to N DO  FOR K := j to np1 DO BEGIN
             L := l + 1;     A[l] := A[l] + E[j]*E[k];   END;
        Sum := Sum + E[np1]*E[np1];
        IF abs(E[np1]) > abs(Rsdmax) THEN BEGIN
             Irsd := M;       Rsdmax := E[np1];    END;
     END;    {  Normal Equations  }

PROCEDURE SelectReflections;
    VAR
        Tolp, Cnst, Db, Thtry, Dmndf                : REAL;
        J, K, L, Ihp, Jhp, Igp, Isum, Jadrs, Kadrs  : INTEGER;
    LABEL
        28,32,40,48,50,56,60,64,66,68,70,90,98,99,100,999;
    BEGIN
        Write('S E L E C T     R E F L E C T I O N S');
        ScreenX := WhereX + 4;     ScreenY := WhereY;
        Icycle  := Icycle +1;         Tolp := Tol;
        IF Itht = 0 THEN              Tolp := Tolp + Tolp;
        IF Printer THEN
          Writeln(Lstr,'Cycle ',Icycle:3,'    ',Title,'    ',Ialf12:3,
                       'Theta Tolerance = ',Tolp:9:5,'     ',Crystal);
        M := 0;       Irsd := 0;       Sum := 0;      Rsdmax := 0;
        IF Printer THEN BEGIN
           Write(Lstr,'  N    H   K   L        D Calc        D Obs');
           Write(Lstr,'       Lambda    ',Ialf12:5,'Theta Calc ',Ialf12:5);
           Writeln(Lstr,'Theta Obs ',Ialf12:5,'Theta Diff       Weight');  END;
        IF MM < 1 THEN    Goto 999;
        FOR L := 1 to Imn DO    A[l] := 0;        Ihp := 1;
        FOR Igp := 1 to MM DO BEGIN
           ActiveScreen(Igp);           Ja  := Ialf[10];
           Swt := Sw[igp];      Thobs := Obsth[igp];       Isum := 0;
           FOR L := 1 to 3 DO BEGIN
              HKL[l] := HKLg[igp,l];
              IF HKL[l] <> 0 THEN      Isum := Isum + 1;   END;
           Pwl := Wl[igp];    Dobs := Dsp[igp];
           Dmndf := Abs(Dcm[Ihp]-Dobs);
           Jhp := Ihp + 1;
    28:    IF Jhp > Nr THEN                                       Goto 32;
           IF Abs(Dcm[jhp]-Dobs) > Dmndf THEN                     Goto 32;
                 Dmndf := Abs(Dcm[jhp]-Dobs);
                 Jhp := Jhp + 1;                                  Goto 28;
    32:    Ihp := Jhp - 1;       IoCalcD[igp] := Ihp;
           IF Isum <> 0 THEN                                      Goto 70;
           D := Dcm[ihp];     Dlst[igp] := D;    Consx := 0.5*Pwl;
           IF D <= Consx THEN                                     Goto 99;
           Thcmp := Theta(Consx,D);
           IF Abs((Thobs-Thcmp)*Swt) > Tol THEN                   Goto 99;
           J := Igp - 1;      Jadrs := 50;
    40:    IF (J < 1) OR (J > MM) THEN                            Goto 48;
           IF Abs(Obsth[j]-Thobs)*Swt < Tolmn THEN                Goto 99;
           Cnst  := 0.5*Wl[j];
           IF D <= Cnst THEN                                      Goto 48;
           Thtry := Theta(Cnst,D);
           IF Abs((Obsth[j]-Thtry)*Sw[j]) <= Tol THEN             Goto 99;
    48:    IF Jadrs = 56 THEN                                     Goto 56;
    50:    J := igp + 1;      Jadrs := 56;                        Goto 40;
    56:    K := ihp - 1;      Kadrs := 66;
    60:    IF (K < 1) OR (K > Nr) THEN                            Goto 64;
           Db := Dcm[k];
           IF Db <= Consx THEN                                    Goto 64;
           Thtry := Theta(Consx,Db);
           IF Abs(Thtry-Thcmp)*Swt < Tolmn THEN                   Goto 99;
           IF Abs(Thtry-Thobs)*Swt <= TOl  THEN                   Goto 99;
    64:    IF Kadrs = 68 THEN                                     Goto 68;
    66:    K := Ihp + 1;      Kadrs := 68;                        Goto 60;
    68:    IF (Thobs > Thtmx) AND (Icycle <= Ncyc) THEN           Goto 100;
           FOR L := 1 to 3 DO   HKL[l] := HKLh[ihp,l];            Goto 90;
    70:    Consx := 0.5*Pwl;        Qcomputation(HKL);
           IF Q <= 0 THEN                                         Goto 99;
           D     := 1/sqrt(Q);     Dlst[igp] := D;
           IF Consx >= D THEN                                     Goto 99;
           Thcmp := Theta(Consx,D);
           IF Abs((Thcmp-Thobs)*Swt) > Tol THEN   JA := Ialf[11];
    90:    Diff  := Thcmp-Thobs;
           IF JA = Ialf[11] THEN                                  Goto 98;
           ObservationEquation;
           IF Printer THEN   HKLPrint;
           NormalEquations;                                       Goto 100;
    98:    HKLPrint;
    99:    IoCalcD[igp] := -IoCalcD[igp];
   100:
      END; { for }
   999: Writeln;
   END;  { Select Reflections }

 { @@@@@@      =====       Main Computation Cycle      =====      @@@@@@  }

PROCEDURE MainComputationCycle;
    VAR
         I, J, K  : INTEGER;
         Crmax    : REAL;

PROCEDURE NoSolution;
    BEGIN
       Writeln(Lstr);
       Writeln(Lstr, '* * *      N O     S O L U T I O N       * * * *');
       Writeln(Lstr);       TheoreticalReflections;       Halt;
    END;   { No Solution }

PROCEDURE ResetTol;
   BEGIN
      IF (Icycle >= 10) OR ((Icycle > Ncyc) AND (Tol <= Tolmn)) THEN
              NoSolution;
      Tol := 0.5*Tol;       Ncyc := 0;
      IF Tol < Tolmn THEN    Tol := Tolmn;
   END;

   PROCEDURE  Singular_Equation;
       BEGIN
           Write('Singular Normal Equation Matrix - Row ',Err:2);
           Writeln(' Linearly Dependent on Previous Rows');
           ResetTol;
       END;

    PROCEDURE  Corrections( J : INTEGER );
       VAR
          K  :  INTEGER;
       BEGIN
          FOR K := 1 to J DO BEGIN
                 X[k] := X[k] - Cr[k];
                 U[k] := U[k] - Cr[k+3];      END;
       END;              {   Corrections   }

    PROCEDURE  Correct_Reciprocal_Cell_Edge( J : INTEGER );
       BEGIN
          Writeln('Corrected value of Reciprocal Cell Edge ',j:2,
                     ' less than or equal to zero');
          Corrections( j );                ResetTol;
       END;              {  Correct Reciprocal Cell Edge  }

    PROCEDURE  Correct_Reciprocal_Cell_Cosine( J : INTEGER );
       BEGIN
           Writeln('Corrected Value of Reciprocal Cell Cosine ',J:2,
                     ' >= 1.0 in Absolute Value');
           Corrections( j );             ResetTol;
       END;         {   Correct Reciprocal Cell Cosine   }

    LABEL   80;
    BEGIN
       Writeln('Your Job has reached    M A I N    Computation Cycle');
 80:   SelectReflections;
       IF Mm  < N THEN   NoSolution;
       Mmn := M - N;
       IF Mmn < 0 THEN BEGIN            ResetTol;         Goto 80;  END
          ELSE BEGIN     J := N + 1;      K := J;
            FOR I := 1 to N DO BEGIN
              Aa[i] := A[k];    J := J - 1;
              K := K + J;   END;  END;
        SquareRootMethod;
        IF Err <> 0 THEN BEGIN  Singular_Equation;        Goto 80;  END;
        Coordinates;
        FOR J := 1 to 3 DO BEGIN
           X[j] := X[j] + Cr[j];       U[j] := U[j] + Cr[j+3];
           W[j] := SQRT( 1 - U[j]*U[j] );
           IF X[j] <= 0.00 THEN BEGIN
               Correct_Reciprocal_Cell_Edge( j );         Goto 80;  END;
           IF Abs(U[j]) - 1.0 >= 0.00 THEN BEGIN
              Correct_Reciprocal_Cell_Cosine( j );        Goto 80;  END; END;
       DirectCellParameters;
       IF Mmn > 0 THEN BEGIN
           FOR I := 1 to N DO   Sum := Sum - Aa[i]*B[i];
           IF Sum < 0 THEN      Sum := 0;
           Thtse := sqrt(Sum/Mmn);  Tol   := Semult*Thtse;
               END ELSE             Tol   := 0.5*Tol;
       IF Tol < Tolmn THEN    Tol := Tolmn;
       IF Tol > Tolmx THEN    Tol := Tolmx;
       IF Icycle <= Ncyc THEN BEGIN    Dspacing;    SortD;  Goto 80; END
        ELSE BEGIN          Crmax := 0;
          FOR I := 1 to N DO IF Abs(B[i]) > Crmax THEN
              Crmax := Abs(B[i]);
          IF (Crmax >= 0.5E-6) AND (Icycle < 10) THEN BEGIN
             Dspacing;    SortD;                            Goto 80; END;
          IF Icycle >= 10 THEN
             Writeln(Lstr,'N O T    C O N V E R G E N T   in 10 Cycles');
          IF Printer = False THEN BEGIN
                Writeln(Lstr);    Printer := True;
                Writeln(Lstr,'F I N A L    V A L U E S    for  ',Title);
                FOR I := 1 to 30 DO Write(Lstr,'*   ');  Writeln(Lstr);
                DirectCellParameters;          END;
          IF Mmn > 0 THEN    StandErrors;
          TheoreticalReflections;              END;
    END;                           {  Main Computation Cycle  }

BEGIN         {    Main Program     }
    TextMode( bw80 );      ICycle := 0;
    Err := 0;      Printer := False;      NameOfSystem;     CommandTail;
    InitializePrinter;        InputTheData;        InitializeParameters;
    GenerateMillerIndices;    SortD;        Skweze;      Tol  := TolMx;
    MainComputationCycle;
    WRITELN(LSTR);  WRITELN(LSTR,'eoj  .. ',Title,'    . .  >  ',Time);
    WRITELN(LSTR,'            LSUCRIPC   v.  ',Version);
    CLOSE( LSTR );    Writeln;
    Writeln(#7,'   Your data has been Processed.  Call Again.',#7);
    SOUND(789);  DELAY(1000);  NOSOUND;   Writeln;
END.

-------------------------------------------------------------------------
End of File