Title                  :Stereographic Projection Program Version 4.00
Keywords               :Stereographic Projection, CBED, SAD, TEM, AEM, 
HOLZ
Computer               :IBM PC/XT/AT or compatible
Operating System       :PC-DOS
Programming Language   :QUICKBASIC 4.5
Hardware Requirements  :Hewlett Packard HP7470A Plotter required.
Author(s)              :John R. Porter
Correspondence Address :Rockwell International Science Center,
                       :Thousand Oaks, CA 91360.
Abstract:

This QuickBASIC program plots stereographic projections, CBED and
HOLZ line simulations, and performs axis/angle pair calculations.  
Output is to the screen or a Hewlett Packard HP 7470A Plotter.  Stereographic 
projections are drawn to scale for subsequent manipulations with a standard 
18cm. Wulff net and can be plotted for cubic, hexagonal, tetragonal and 
monoclinic crystal systems in any orientation.  Plotted poles can be 
restricted to those allowed by structure factor considerations (for certain 
structures) or restricted to those in certain Laue zones.  The program 
leads the user through menus to determine the structure and orientation  
for the projection, which is then plotted accordingly.  Version 4.00 is 
significantly enhanced compared to earlier versions.  

EDITORS NOTE:  THIS IS A NEW UPDATED VERSION OF STERPROJ- IBM  
(FEB.1990)
-----------------------------------------------------------------------------
Title                  :Stereographic Projection Program Version 4.00
Keywords               :Stereographic Projection, CBED, SAD, TEM, AEM, HOLZ
Computer               :IBM PC/XT/AT or compatible
Operating System       :PC-DOS
Programming Language   :QUICKBASIC version 4.5
Hardware Requirements  :Hewlett Packard HP7470A Plotter required.
Author(s)              :John R. Porter
Correspondence Address :Rockwell International Science Center,
                       :Thousand Oaks, CA 91360.

SP.EXE is a program for an IBM PC which plots most combinations of
stereographic projections, CBED and HOLZ line simulations and axis/angle pair
calculations for cubic, tetragonal, hexagonal and monoclinic crystals.  Version
4.0 of the program has been significantly enhanced compared with the previous
program available through EMMPDL.  In particular, screen output has been
improved, more crystal systems are incorporated, an easier to use menu
interface has been designed and the program runs faster.

The program is written in QUICKBASIC and requires a Hewlett Packard HP 7470A
Plotter connected to an IBM PC or compatible via COM1:.  The communications
parameters should be as described on p. 35 of the HP7470A Plotter
Interconnection Guide (i.e. set dip-switches on plotter for 9600 baud, 7 data
bits, 1 stop bit, no parity - or change program line 4250).  The plotter must
be connected to the computer via a null modem cable.  For screen output, the
IBM color graphics adapter or an enhanced graphics adapter is required.

The ASCII file SP.BAA can be compiled by the Microsoft Quickbasic compiler
Version 4.5.  Before compiling, select the type of graphics adapter in the
beginning of the program (the default is EGA).

To run the compiled program, at the DOS prompt, type SP<RETURN>.

The program opens to a menu and the selections should be answered in turn
(asterisks indicate questions to be answered).  The program saves information
about crystals in data files with the extension .XTL.  Sample .XTL files are at
the end of this documentation file or the program will make them by selecting
menu choice 2.  Therefore, initially select a crystal by choosing 1 or 2.
Lattice parameter changes can be made in choice 3.  Use item 4 to select
whether output is to the screen or the plotter.  In item 5, select from an
Axis/Angle pair computation, a blank circle (for manual stereographic
projection plots), a CBED pattern simulation, a HOLZ line pattern simulation or
a Stereographic projection.  Not all options are available for output to the
screen, however.

Depending on the selection, further choices need to be made.  There are two
ways to input orientation data for stereographic projections. One option is to
enter the pole which is to appear at the center of a projection and a second
pole which is to appear on the horizontal axis.  More useful, perhaps, to
electron microscopists is the second option which allows you to enter two poles
which are to appear on the circumference of a projection - these correspond to
indexed spots on a diffraction pattern.  Having selected the way the data is to
be entered, the actual orientations are entered in choice 6 from the main menu.

Of interest to people working with convergent beam diffraction patterns is the
ability to plot only those poles which appear in the zero, first and second
order Laue zones.  The program provides this as an option.  The structure
factor parameter "p" is explained in "Practical Analytical Electron Microscopy
in Materials Science" by D. B. Williams on p. 128.  Basically, for certain
orientations of some crystals, structure factor considerations mean that entire
Laue zones will be forbidden.  "p" is a multiplicity factor which tells the
program to consider only every second or third zone.  "p" should be 1 or 2 in
cubic crystals and 1 or 3 in hexagonal crystals.

In item 7, the program will ask you for the highest value of h, k, (i) and l to
plot (3 is recommended for a first try for a stereographic projection, 20 for a
diffraction pattern) and for primitive hexagonal systems the increment for
changing l values (i.e. an increment of 2 will plot l's of say -4, -2, 0, 2, 4
etc.).

Finally, the program will give you the option of manually entering additional
poles.  Experiment!  After confirming that your entry parameters are correct,
the program will proceed to plot the projection on the plotter.  The projection
will be 18 cm in diameter but this can be altered by anyone with knowledge of
BASIC.  Wulff nets and Polar nets of this size can be obtained from N. P. Nies,
969 Skyline Drive, Laguna Beach, CA 92651.  Tel. (714) 494-4619.

Remember that if you make a mistake during program entry, you can reselect a
menu item.  The options are extensive, but one example set of answers is listed
here:

The first choice in the program is to select the crystal.  Select "2<RTN>".
Enter "Aluminum<RTN>", "c"<RTN>", "fcc<RTN>" and "4.05<RTN>" in response to the
questions (do not type the "'s).  The program will store the crystal data in a
file called ALUMINUM.XTL and another time the data can be recalled in choice 1.

Choose 4 and select "s<RTN>" for output to screen.

Choose 5 and select "p<RTN>" for a stereographic projection.  Then "2<RTN>" to
select the data input method (see above).  Then "<RTN>" to NOT select 1st. and
2nd. HOLZ poles only.  (Try that another time).

Now the program requires the Miller indices for a pole which is to appear on
the circumference, in the vicinity of the East pole.  h, k and l will be
selected separately.  Then the exact position is specified by entering the
angle the input pole makes with the horizontal axis by rotation about the
center.  This angle should be less than 90 degrees and counterclockwise is
positive.  Next, the indices of a second pole, located anticlockwise from the
first pole, should be entered.  This uniquely defines the geometry of the
projection.  Note that defaulting any answers in this part of the program
corresponds to entering zero.

Therefore, choose 6 and select "2<RTN>2<RTN><RTN>" for h1k1l1, "<RTN>" for the
rotation angle and "1<RTN>-3<RTN>1<RTN>" for h2k2l2.

Choose 7 and select "2<RTN>" for max hkl.  Select "N<RTN>" to NOT select more
poles.

Choose 8 and a stereographic projection will be plotted for the 1-1-4
orientation selected.  The 220 and 1-31 poles will be on the circumference.
Select "N<RTN>" to return to the main menu.

Go back to item 5 and select "h<RTN>" for HOLZ line simulation.  Default the
next two questions with "<RTN><RTN>"  Go to 7 and select a max hkl of 15 and a
HOLZ parameter, p, of 2.  Replot in 8 and get a HOLZ line simulation for the
1-1-4 orientation.  The choices are numerous!

You can color code for higher order Laue zones - this selection will cause the
program only to plot poles in the zero, first and second order Laue zones.
When plotted, the ZOLZ poles will be black open circles, the FOLZ poles will be
red and the SOLZ poles filled black circles.  Defaulting the question (i.e.
just pressing <RETURN> will plot all poles in the specified range).

Finally - will you wish to plot more poles (say the 7 -3 32 pole)?  If so,
answer Y and follow screen instructions.

Other pathways through the menus are similar.  Of interest could be the option
that appears when selecting a projection for alumina (hexagonal) which plots
facet planes.  For this option, the poles corresponding to the planes on which
alumina is known to facet are plotted instead of a specified range of h, k and
l.

A little more specialized is the axis / angle pair pathway.  Here data for two
crystal orientations is entered, corresponding to two grains. Key poles for
grain 1 are plotted as open circles and for grain 2 as closed circles.  Now you
have to select two poles from grain 1 (a and b) which are to rotate to
equivalent poles (a' and b') in grain 2.  All four poles should be on the
as-plotted projection.  Follow the menu instructions to choose poles a, b, a'
and b'.  Be careful to rotate only to equivalent poles.  i.e. for the hexagonal
{11.0} poles, only rotate red to red and black to black (note that {00.1} poles
are always equivalent and there is no corresponding restriction).

The simulated CBED pattern menus provide similar choices.  The large plotted
circles represent the positions of the HOLZ rings and the small circles the
positions of the reflections.  If no reflections are plotted, there are two
possible explanations that don't involve program error.  First, the specified
range of h, k and l may be too small (up to 50 can be needed in some cases - in
which case plotting can take a long time while the computer completes 1,000,000
iterations of a loop).  Second, if a structure factor calculation is being
performed, the parameter p may need to be entered as 2 or 3.

The following notes make reference to figures 1 and 2.  These figures should be
plotted by printing them to the plotter, using the DOS PRINT command.  The
ASCII code for each figure is at the end of this file between asterisks.  Use
your word processor to mark each block and save each block as an ASCII file
named FIGURE.1 and FIGURE.2 respectively.  Now, there is a small problem.  Each
C in the figure codes must be converted to a <CTRL>C and each H converted to a
<CTRL>H.  Again, use your word processor (Sidekick works well).  Resave the
code.  Then, at the DOS prompt, type PRINT FIGURE.1 and PRINT FIGURE.2.
Respond AUX when asked for the printing device (make sure PRINT.COM is present
on your disk).  If this doesn't work for any reason - write to the address
below for copies of the figures.

The procedures used in writing the program are explained with reference to
figures 1 and 2.  To plot a pole, e.g. hkl in figure 1, a number of values need
to be calculated.  First, the direction u3v3w3 at the center of the projection
is established as a reference together with two poles at the circumference,
h1k1l1 and h2k2l2.  u1v1w1 and u2v2w2 are then the zone axes of h1k1l1 and
h2k2l2 with the plane normal to u3v3w3 (defined as h3k3l3).

For each pole, hkl, the angle the line from the center to hkl makes with the
horizontal x-axis, rho, needs to be determined.  This is done by measuring the
angle uvw (zone axis of hkl and h3k3l3) makes with u1v1w1 (zone axis of h1k1l1
with h3k3l3).  However, rho is calculated from an arccosine formula, which
cannot differentiate between rho and (2 pi - rho).  Therefore a pole in the
upper right hand quadrant is defined, h2k2l2, and angles rho 1 and rho 2 are
calculated in similar ways.  Then the true value of rho can be established by
the requirement that rho = rho 1 + rho 2.

Having established rho, the radial distance of the pole from the center needs
to be established.  Figure 2 shows the geometry of the projection.  The radial
distance is the diameter multiplied by the tangent of theta / 2, where theta is
the interplanar angle between hkl and h3k3l3.

During the execution of the program, all potential poles within a user
specified range are selected.  If certain criterion are met, the pole is
plotted according to the geometrical procedures outlined above.

The extra features of the program (simulated CBED patterns, axis / angle pair
determinations, etc.) make use of the same routines.  Remember, however, that a
simulated diffraction pattern will need to be plotted out to higher maximum
values of h, k and l (probably 20 or higher if poles in HOLZ rings are to be
plotted).

The program itself is extensively documented.  For more information, examine
the program listing.

In computereze, this program is classified as "Shareware".  In other words, if
you find the program useful, I would like you or your institution to buy it.
In return for a nominal fee ($50.00), the author will provide updates as they
become available, in both compiled and source code form, or make minor
modifications for specific purposes.  Suggestions for changes (or corrections)
would be welcomed.  Contact:-

         John R. Porter,
         1096 Calle Las Trancas,
         Thousand Oaks, CA 91360

         (805) 373 4702 - Daytime.

Figure.1 code:
*****
in;sp 1;sc -1400,1400,-1000,1000;
pa 0,0;vs 15;ci 1000;vs;
sr 0.5,1;
pa 0,0;pd 1000,0;pu 0,0;pd 707,707;pu 0,0;pd 0,1000;pu;
pa 0,0;pd -707,707;pu;ci10,pa -850,707;lbu v wC;pa -850,692;lb 2 2 2C;
pa 0,0;ci 10;pa -180,0;lb u v wC;pa -180,-15;lb  3 3 3C;
pa 1000,0; ci 10; lb h k lC;pa 1000,-15; lb  1 1 1C;
pa 707,707; ci 10; lb h k lC;pa 707,692; lb  2 2 2C;
pa 0,-40; lb h k lC;pa 0,-55; lb  3 3 3C;
pa 0,1000; ci 10; pa 0,970; lb u v wC;pa 0,955;lb  1 1 1C;
pa 200,0; pd; aa 0,0,45;pu;
pa 0,0; pd -300,-200;pu; ci 10; lb hklC;pu
pa 0,0; pd 554,-832;pu;ci 10;lb uvwC;pu
pa 0,200;pd; aa 0,0 213;pu -350,350;pd;aa 0,0,168;pu
pa 200,100;lb rho 1C;pa -270,50; lbrhoC;
pa 0,300;pd;aa 0,0,45;pu;pa -220,300;lbrho 1C;pa -620,0;lbrho 2C;
pa -1300,900; sr;lb Figure 1C;
pa -1400,-1000;sp 0;
*****
Figure.2 code
*****
in;sp 1;sc -1400,1400,-1000,1000;
pa 0,0;vs 15;ci 500;vs;
sr 0.5,1;
pa -1100,500;pd 1100,500;pu 0,0;pd 353,353;pu 0,-500;
pd 0,500;pu;pa -500,0; pd 500,0; pa 0,-500;pd -1000,500;pu 0,-500;pd 1000,500;
pu 0,-500;pd 414,500;pu;pa 5,-350;lb0H-/2C; pa 20,50;lbOH-C;
pa 0,100; pd;ar 0,-100,-45;pu 0,-250; pd;ar 0,-250,-22;pu;
pa -50,525;lbh k lC;pa -50,510;lb 3 3 3C;
pa 400,525;lbhklC;pa 950,525;lbh k lC;pa 950,510;lb 1 1 1C;
pu 0,0;ci 10; pa -1300,900; sr;lb Figure 2C; pa -1400,-1000;;sp 0;
*****
The following should be saved as ALUMINA.XTL
*****
"Alumina"
"H"
"Alumina"
4.758
12.991
*****
The following should be saved as ZINC.XTL
*****
"Zinc"
"H"
"cph"
2.664
4.946
*****

Title                  :Stereographic Projection Program Version 4.00
Keywords               :Stereographic Projection, CBED, SAD, TEM, AEM, HOLZ
Computer               :IBM PC/XT/AT or compatible
Operating System       :PC-DOS
Programming Language   :QUICKBASIC 4.5
Hardware Requirements  :Hewlett Packard HP7470A Plotter required.
Author(s)              :John R. Porter
Correspondence Address :Rockwell International Science Center,
                       :Thousand Oaks, CA 91360.

SOURCE CODE


REM Stereographic Projection Program for Hewlett Packard Plotter
REM IBM PC/XT/AT or compatible
REM Compile with QUICKBASIC 4.5
REM Hewlett Packard Plotter and EGA recommended for output
REM John R. Porter
REM Rockwell International Science Center, Thousand Oaks, CA 91360
REM Version 4.09
REM
REM Contact address:-
REM
REM             John R. Porter
REM             Rockwell International Science Center
REM             P. O. Box 1085
REM             1049 Camino Dos Rios
REM             Thousand Oaks, CA 91360
REM
REM             (805) 373 4702 - daytime
   DEFINT N: DEFDBL A-M: DEFDBL O-Z
   display$ = "EGA" ' This variable must be either "EGA" or "CGA"
   IF display$ = "CGA" THEN
      disp1 = 2
      disp2 = 9
   ELSE
      disp1 = 9
      disp2 = 8
   END IF
   DIM a(15, 4), B%(3), C(4, 2), D(7), E(9)
   h = 0: k = 0: i = 0: l = 0: nmax = 0: LMAX% = 0: LZ = 0: PI = 4 * ATN(1): NPN = 1: flag = 0: LINCR% = 1: RH = 0: RV = 0: CA = 0: LPA = 0: LPC = 0: NC = 0: SN = .000001
   RD = 0: GMAX = 0: G1MIN = 0: G1MAX = 0: G2MIN = 0: G2MAX = 0
   C$ = CHR$(3): LNF$ = CHR$(10): CR$ = CHR$(13): HOLZ$ = "": HOLZR$ = "": LABEL$ = "": LF$ = "20": plot% = 0
    DIM NB9(61)
    DIM NB8(61)
    DIM NB7(61)
    DIM NB6(61)
    DIM NB5(61)
    DIM NB4(61)
    DIM NB3(61)
    DIM NB2(61)
    DIM NB1(61)
    DIM N0(61)
    DIM N1(61)
    DIM N2(61)
    DIM N3(61)
    DIM N4(61)
    DIM N5(61)
    DIM N6(61)
    DIM N7(61)
    DIM N8(61)
    DIM N9(61)
   CLS
   SCREEN disp1, 0
   IF display$ = "EGA" THEN COLOR 14, 1
    WINDOW (0, 0)-(640, 200)
   LINE (8, disp2)-(13, disp2)
    LOCATE 25, 1: PRINT 9; : GET (8, 1)-(16, 9), NB9
    LOCATE 25, 1: PRINT 8; : GET (8, 1)-(16, 9), NB8
    LOCATE 25, 1: PRINT 7; : GET (8, 1)-(16, 9), NB7
    LOCATE 25, 1: PRINT 6; : GET (8, 1)-(16, 9), NB6
    LOCATE 25, 1: PRINT 5; : GET (8, 1)-(16, 9), NB5
    LOCATE 25, 1: PRINT 4; : GET (8, 1)-(16, 9), NB4
    LOCATE 25, 1: PRINT 3; : GET (8, 1)-(16, 9), NB3
    LOCATE 25, 1: PRINT 2; : GET (8, 1)-(16, 9), NB2
    LOCATE 25, 1: PRINT 1; : GET (8, 1)-(16, 9), NB1
    CLS
    LOCATE 25, 1: PRINT 0; : GET (8, 1)-(16, 9), N0
    LOCATE 25, 1: PRINT 1; : GET (8, 1)-(16, 9), N1
    LOCATE 25, 1: PRINT 2; : GET (8, 1)-(16, 9), N2
    LOCATE 25, 1: PRINT 3; : GET (8, 1)-(16, 9), N3
    LOCATE 25, 1: PRINT 4; : GET (8, 1)-(16, 9), N4
    LOCATE 25, 1: PRINT 5; : GET (8, 1)-(16, 9), N5
    LOCATE 25, 1: PRINT 6; : GET (8, 1)-(16, 9), N6
    LOCATE 25, 1: PRINT 7; : GET (8, 1)-(16, 9), N7
    LOCATE 25, 1: PRINT 8; : GET (8, 1)-(16, 9), N8
    LOCATE 25, 1: PRINT 9; : GET (8, 1)-(16, 9), N9
    SCREEN 0, 0, 0
10  COLOR 3, 1, 1: CLS : PRINT : PRINT "         Stereographic Projection for Hewlett Packard Plotter"
    PRINT : PRINT : PRINT "Choose from:"
    PRINT : IF NMEN1 = 0 AND NMEN2 = 0 AND NMEN3 = 0 THEN PRINT "* ";  ELSE PRINT "  ";
    PRINT "1 - Retrieve crystal data file"; TAB(50); CRYSTAL$
    IF NMEN1 = 0 AND NMEN2 = 0 AND NMEN3 = 0 THEN PRINT "* ";  ELSE PRINT "  ";
    PRINT "2 - Make crystal data file"
    IF NMEN1 = 0 AND NMEN2 = 0 AND NMEN3 = 0 THEN PRINT "* ";  ELSE PRINT "  ";
    PRINT "3 - Modify crystal data"
    IF NMEN4 = 0 THEN PRINT "* ";  ELSE PRINT "  ";
    PRINT "4 - Select output device"; TAB(50); device$
    IF NMEN5 = 0 THEN PRINT "* ";  ELSE PRINT "  ";
    PRINT "5 - Select output type"; TAB(50); TYPE$
    IF NMEN6 = 0 THEN PRINT "* ";  ELSE PRINT "  ";
    PRINT "6 - Input crystal orientation"; TAB(50); H1$; K1$; I1$; L1$; TAB(65); H2$; K2$; I2$; L2$
    IF NMEN7 = 0 THEN PRINT "* ";  ELSE PRINT "  ";
    PRINT "7 - Select output variables"
    PRINT "  8 - Plot"
    PRINT "  0 - Quit"
    PRINT : INPUT "Make selection: ", NSEL
   ON NSEL GOTO 20, 30, 60, 70, 90, 180, 230, 270
   IF NSEL = 0 THEN PRINT : INPUT "Are you sure (Y/N): ", ans$: IF ans$ = "y" OR ans$ = "Y" THEN SYSTEM ELSE 10
   GOTO 10
20  CLS : PRINT : PRINT "         Stereographic Projection for Hewlett Packard Plotter"
   'ON ERROR GOTO 21
   NMEN1 = 1
   PRINT : PRINT
   PRINT
   FILES "*.xtl"
   PRINT : PRINT : INPUT "Enter crystal: ", CRYSTAL$
    OPEN CRYSTAL$ + ".xtl" FOR INPUT AS #2
   INPUT #2, CRYSTAL$: INPUT #2, st$: INPUT #2, SF$: INPUT #2, LPA
   IF st$ = "M" OR st$ = "O" THEN INPUT #2, LPB
   IF st$ <> "C" THEN INPUT #2, LPC
   IF st$ = "M" THEN INPUT #2, BETA
   CLOSE #2
   'ON ERROR GOTO 0
   IF SF$ = "" THEN SF% = 0
   IF SF$ = "fcc" THEN SF% = 1
   IF SF$ = "fct" THEN SF% = 2
   IF SF$ = "ooe" THEN SF% = 3
   IF SF$ = "dc" THEN SF% = 4
   IF SF$ = "bcc" THEN SF% = 5
   IF SF$ = "bct" THEN SF% = 6
   IF SF$ = "cph" THEN SF% = 7
   IF SF$ = "Alumina" THEN SF% = 11
   IF SF$ = "Alpha SiC" THEN SF% = 12 ' SG 186
   IF SF$ = "Beta SiC" THEN SF% = 1 ' SG 219
   IF SF$ = "Rutile" THEN SF% = 136
   IF SF$ = "Anatase" THEN SF% = 141
   IF SF$ = "Alpha Silicon Nitride" THEN SF% = 12 ' SG 159
   IF SF$ = "Beta Silicon Nitride" THEN SF% = 176
   GOTO 10
'21 PRINT : INPUT "Error in file name.  Try again - Press <RTN> to continue. ", a$
'   RESUME 20
30  CLS : PRINT : PRINT "         Stereographic Projection for Hewlett Packard Plotter"
   NMEN2 = 1
   PRINT : PRINT : INPUT "Enter crystal: ", CRYSTAL$
   OPEN CRYSTAL$ + ".xtl" FOR OUTPUT AS #2
40  PRINT : PRINT "Enter crystal structure:": PRINT : PRINT "C(ubic)": PRINT "T(etragonal)": PRINT "H(exagonal)": PRINT "O(rthorhombic)": PRINT "M(onoclinic)": PRINT : INPUT "Make selection: "; st$
   IF st$ = "c" THEN st$ = "C"
   IF st$ = "t" THEN st$ = "T"
   IF st$ = "h" THEN st$ = "H"
   IF st$ = "o" THEN st$ = "O"
   IF st$ = "m" THEN st$ = "M"
   IF st$ = "C" OR st$ = "H" OR st$ = "M" OR st$ = "O" OR st$ = "T" THEN 50 ELSE 40
50  WRITE #2, CRYSTAL$: WRITE #2, st$
   SF$ = "": SF% = 0
   IF st$ = "C" THEN PRINT : INPUT "Enter structure factor (P(rimitive), fcc, bcc or dc(diamond): "; SF$
   IF st$ = "T" THEN PRINT : INPUT "Enter structure factor - P(rimitive), fct, bct or ooe: "; SF$
   IF st$ = "H" THEN PRINT : INPUT "Enter structure factor - P(rimitive), cph or A(lumina): "; SF$
   IF SF$ = "p" OR SF$ = "P" THEN SF$ = "": SF% = 0
   IF SF$ = "FCC" THEN SF$ = "fcc": SF% = 1
   IF SF$ = "BCC" THEN SF$ = "bcc": SF% = 5
   IF SF$ = "DC" THEN SF$ = "dc": SF% = 4
   IF SF$ = "FCT" THEN SF$ = "fct": SF% = 2
   IF SF$ = "BCT" THEN SF$ = "bct": SF% = 6
   IF SF$ = "CPH" THEN SF$ = "cph": SF% = 7
   IF SF$ = "a" OR SF$ = "A" THEN SF$ = "Alumina": SF% = 11
   WRITE #2, SF$
   PRINT : INPUT "Enter lattice parameter, a: ", LPA: WRITE #2, LPA
   IF st$ = "M" OR st$ = "O" THEN INPUT "Enter lattice parameter, b: ", LPB: WRITE #2, LPB
   IF st$ = "T" OR st$ = "H" OR st$ = "M" OR st$ = "O" THEN INPUT "Enter lattice paramete, c: ", LPC: WRITE #2, LPC
   IF st$ = "M" THEN INPUT "Enter angle, beta: ", BETA: BETA = BETA * PI / 180: WRITE #2, BETA
   CLOSE #2
   GOTO 10
60  LP = 0: PRINT "Lattice parameter 'a' is currently "; LPA; : INPUT " Enter new value: ", LP: IF LP > 0 THEN LPA = LP
   LP = 0: IF st$ = "O" OR st$ = "M" THEN PRINT "Lattice parameter 'b' is currently "; LPB; : INPUT " Enter new value: ", LP: IF LP > 0 THEN LPB = LP
   LP = 0: IF st$ <> "C" THEN PRINT "Lattice parameter 'c' is currently "; LPC; : INPUT " Enter new value: ", LP: IF LP > 0 THEN LPC = LP
   LP = 0: IF st$ = "M" THEN PRINT "'beta' currently "; BETA; : INPUT " Enter new value: ", LP: IF LP > 0 THEN BETA = LP
   NMEN3 = 1
   GOTO 10
70  CLS : PRINT : PRINT "         Stereographic Projection for Hewlett Packard Plotter"
   NMEN4 = 1: NMEN6 = 0
80  PRINT : PRINT : INPUT "Output to screen or plotter (S or P): ", SPB$
   IF SPB$ = "s" OR SPB$ = "S" THEN SPB$ = "S": device$ = "Screen"
   IF SPB$ = "p" OR SPB$ = "P" THEN SPB$ = "P": device$ = "Plotter"
   IF SPB$ = "b" OR SPB$ = "B" THEN SPB$ = "B": device$ = "Screen and plotter"
   IF SPB$ <> "S" AND SPB$ <> "P" AND SPB$ <> "B" THEN 80
   IF SPB$ = "S" THEN 10
   GOSUB 1000
   CLOSE #1
   GOTO 10
90  CLS : PRINT : PRINT "         Stereographic Projection for Hewlett Packard Plotter"
   NMEN5 = 1
100 PRINT : PRINT "Select:": PRINT
   IF SPB$ = "P" THEN IF st$ = "C" OR st$ = "T" OR st$ = "H" THEN PRINT "Axis/angle pair plot (A)"
   IF SPB$ = "P" THEN PRINT "Blank circle (B)"
   PRINT "Diffraction Pattern (D)": PRINT "HOLZ line pattern (H) ": PRINT "Stereographic Projection (P)"
   IF SPB$ = "P" THEN IF st$ = "C" OR SF% = 11 THEN PRINT "Unit triangle (T)"
   PRINT : INPUT "Make selection: ", DPSP$
   IF DPSP$ = "d" OR DPSP$ = "D" THEN DPSP$ = "D": TYPE$ = "CBED simulation"
   IF DPSP$ = "h" OR DPSP$ = "H" THEN DPSP$ = "H": TYPE$ = "HOLZ line simulation"
   IF DPSP$ = "p" OR DPSP$ = "P" THEN DPSP$ = "P": TYPE$ = "Stereographic Projection"
   IF DPSP$ = "a" OR DPSP$ = "A" THEN DPSP$ = "A": AAP = 0: NMEN7 = 1: TYPE$ = "Axis/angle pair"
   IF DPSP$ = "t" OR DPSP$ = "T" THEN DPSP$ = "T": NMEN6 = 1: NMEN7 = 1: TYPE$ = "Unit triangle"
   IF DPSP$ = "b" OR DPSP$ = "B" THEN DPSP$ = "B": NMEN6 = 1: NMEN7 = 1: TYPE$ = "Blank Circle"
   IF DPSP$ = "A" OR DPSP$ = "B" OR DPSP$ = "D" OR DPSP$ = "H" OR DPSP$ = "P" OR DPSP$ = "T" THEN 110 ELSE 100
110 IF DPSP$ = "B" THEN 170
   IF DPSP$ = "D" OR DPSP$ = "H" THEN PRINT : INPUT "Operating voltage? (default 120) ", KV: PRINT : INPUT "Convergence half-angle (mrad)? (default 10) ", CONANG: PRINT
   IF DPSP$ = "D" THEN INPUT "Camera length? (default 1000) ", CAMLEN
   IF KV = 0 THEN KV = 120
   K120 = SQR(KV * 1000 + KV * KV) / 12.3' reciprocal wavelength
   IF CONANG = 0 THEN CONANG = 10
   IF DPSP$ = "H" THEN CAMLEN = 100000! / CONANG
   IF CAMLEN = 0 THEN CAMLEN = 1000
   SCALE = 10 ' CONVERT TO PLOTTER UNITS
   LAMDAL = CAMLEN / K120 * SCALE
   CONANG = CONANG / 1000' radians
   DRRL = CONANG * K120 ' Disc radius - recip lattice
   DR = DRRL * LAMDAL
   IF DPSP$ = "H" THEN HOLZ$ = "H"
   IF DPSP$ = "P" OR DPSP$ = "A" OR DPSP$ = "T" THEN DR = 10
   IF DPSP$ = "D" OR DPSP$ = "H" THEN IP$ = "2": IPUH$ = "H": GOTO 160
   IF DPSP$ = "A" OR DPSP$ = "T" THEN 10
120 PRINT : INPUT "Select center pole (1) or diffraction pattern data (2): ", IP$
   IF IP$ = "1" OR IP$ = "2" THEN 130 ELSE 120
130 IF IP$ = "1" AND st$ <> "C" THEN PRINT : INPUT "Choose to input direction indices, uvw, or plane normals, hkl (U/H): ", IPUH$ ELSE 140
   IF IPUH$ = "u" THEN IPUH$ = "U"
   IF IPUH$ = "h" THEN IPUH$ = "H"
   IF IPUH$ <> "U" AND IPUH$ <> "H" THEN 130
140 IF st$ = "H" THEN PRINT : INPUT "Choose to output direction indices, uvw, or plane normals, hkl (U/H): ", opuh$ ELSE 150
   IF opuh$ = "u" THEN opuh$ = "U"
   IF opuh$ = "h" THEN opuh$ = "H"
   IF opuh$ <> "U" AND opuh$ <> "H" THEN 140
   IF opuh$ = "U" THEN 200
150 IF DPSP$ <> "D" AND DPSP$ <> "H" THEN IF IP$ = "2" OR st$ = "C" OR IPUH$ = "U" THEN PRINT : INPUT "Select (H) to color code for HOLZ patterns: ", HOLZ$' routine needs u3,v3,w3
160 IF DPSP$ = "D" THEN HOLZ$ = "H": PRINT : INPUT "Select (H) to display HOLZ rings ", HOLZR$: IF HOLZR$ = "h" THEN HOLZR$ = "H"
   IF HOLZR$ <> "H" AND HOLZR$ <> "" THEN 160
   IF HOLZ$ = "h" THEN HOLZ$ = "H"
   IF HOLZ$ <> "H" AND HOLZ$ <> "" THEN 150
170 GOTO 10
180 IF DPSP$ = "A" THEN AAP = 0: NC = 0
190 IF DPSP$ = "A" THEN AAP = AAP + 1
   IF DPSP$ = "T" THEN PRINT : PRINT "This selection not needed for this output selection": PRINT : INPUT "Press <RTN> to continue", ans$: GOTO 10
200 COLOR 3, 1, 1: CLS : PRINT : PRINT "         Stereographic Projection for Hewlett Packard Plotter"
   NMEN6 = 1
   IF DPSP$ = "D" THEN POLE$ = " reflection" ELSE POLE$ = " pole"
   IF DPSP$ = "A" THEN IP$ = "2": PRINT : PRINT "Enter data for grain "; AAP
210 IF IPUH$ = "U" THEN PRINT : PRINT "Enter direction indices for"; POLE$; " at center of stereogram": PRINT : INPUT ; "u = ", U3: INPUT ; ", v = ", V3: INPUT " , w = ", W3: T3 = -(U3 + V3)
   IF IPUH$ = "U" THEN PRINT : PRINT "Enter direction indices for"; POLE$; " on right hand horizontal axis": PRINT : INPUT ; "u = ", U4: INPUT ; ", v = ", V4: INPUT " , w = ", W4: T4 = -(U4 + V4): GOTO 220
   IF st$ <> "H" AND IP$ = "1" THEN PRINT : PRINT "Enter Miller indices for"; POLE$; " at center of stereogram"' Pole at center of projection - cubic
   IF st$ = "H" AND IP$ = "1" THEN PRINT : PRINT "Enter Miller-Bravais indices for"; POLE$; " at center of stereogram"' Pole at center of projection  - hexagonal
   IF st$ <> "H" AND IP$ = "2" THEN PRINT : PRINT "Enter Miller indices for"; POLE$; " near RH horizontal axis"
   IF st$ = "H" AND IP$ = "2" THEN PRINT : PRINT "Enter Miller-Bravais indices for"; POLE$; " near RH horizontal axis"
   PRINT : INPUT ; "h = ", H1: INPUT ; " , k = ", K1: INPUT " , l = ", L1
   I1 = -(H1 + K1)
   H1$ = STR$(H1): K1$ = STR$(K1): I1$ = STR$(I1): L1$ = STR$(L1)
   IF st$ <> "H" THEN I1$ = ""
   IF IP$ = "2" THEN PRINT : PRINT "Enter angle of"; POLE$; : INPUT " with horizontal axis (+ve is anticlockwise): ", THETA' applies to diffraction pattern data input only
   THETA = THETA * PI / 180
   IF st$ <> "H" AND IP$ = "1" THEN PRINT : PRINT "Enter Miller indices for"; POLE$; " on the right hand horizontal axis"
   IF st$ = "H" AND IP$ = "1" THEN PRINT : PRINT "Enter Miller-Bravais indices for"; POLE$; " on the right hand horizontal axis"
   IF st$ <> "H" AND IP$ = "2" THEN PRINT : PRINT "Enter Miller indices for"; POLE$; " anticlockwise from first "; POLE$
   IF st$ = "H" AND IP$ = "2" THEN PRINT : PRINT "Enter Miller-Bravais indices for"; POLE$; " anticlockwise from first"; POLE$
   PRINT : INPUT ; "h = ", H2: INPUT ; " , k = ", K2: INPUT " , l = ", L2
   I2 = -(H2 + K2)
   H2$ = STR$(H2): K2$ = STR$(K2): I2$ = STR$(I2): L2$ = STR$(L2)
   IF st$ <> "H" THEN I2$ = ""
   IF AAP = 2 THEN 280
220 GOTO 10
230 IF DPSP$ = "A" OR DPSP$ = "T" THEN PRINT : PRINT "This selection not needed for this output type": PRINT : INPUT "Press <RTN> to continue", ans$: GOTO 10
   NMEN7 = 1
   CLS : PRINT : PRINT "         Stereographic Projection for Hewlett Packard Plotter"
240 IF DPSP$ = "P" AND SF% = 11 THEN PRINT : INPUT "Enter <F> for facet plane indices: ", FACET$
   IF FACET$ = "f" THEN FACET$ = "F"
   IF FACET$ <> "" AND FACET$ <> "F" THEN 240
   IF FACET$ = "F" THEN EXTRA$ = "N": HOLZ$ = "": GOTO 10
   IF DPSP$ = "A" THEN 10
   IF st$ <> "H" THEN PRINT : INPUT "Enter maximum value for h, k, l: ", nmax
   IF st$ = "H" THEN PRINT : INPUT "Enter maximum value for h, k, i: ", nmax
   IF st$ = "H" THEN PRINT : INPUT "Enter maximum value for l: ", LMAX%
   IF st$ = "H" THEN IF SF% = 0 AND HOLZ$ <> "H" THEN PRINT : INPUT "Enter increment for l: ", LINCR%
   IF HOLZ$ = "H" AND SF% > 0 THEN PRINT : INPUT "Enter HOLZ structure parameter, p: ", NP
   IF DPSP$ = "H" AND SF% = 0 THEN PRINT : INPUT "Enter HOLZ ring contributing to pattern: ", NP
   IF NP = 0 THEN NP = 1
   IF nmax = 0 THEN nmax = 2
   IF LMAX% = 0 THEN LMAX% = 2
   IF LINCR% = 0 THEN LINCR% = 2
   IF HOLZ$ = "H" OR SF% > 0 THEN LINCR% = 1
   IF st$ <> "H" THEN LMAX% = nmax
   IF DPSP$ = "H" THEN EXTRA$ = "N": GOTO 260
250 PRINT : PRINT "Will you wish to plot more"; POLE$; "s? (Y/N): "; : INPUT "", EXTRA$
   IF EXTRA$ = "y" THEN EXTRA$ = "Y"
   IF EXTRA$ = "n" THEN EXTRA$ = "N"
   IF EXTRA$ = "N" OR EXTRA$ = "Y" THEN 260 ELSE 250
260 GOTO 10
270 IF SPB$ <> "S" THEN GOSUB 1010
   IF SPB$ = "P" AND DPSP$ = "B" THEN GOSUB 10380
   IF SPB$ = "P" AND DPSP$ = "B" THEN GOTO 350
   IF NMEN1 = 0 AND NMEN2 = 0 AND NMEN3 = 0 THEN INPUT "Crystal must be selected - Press <RTN> to continue", ans$: GOTO 10
   IF NMEN4 = 0 THEN INPUT "Output device must be selected - Press <RTN> to continue", ans$: GOTO 10
   IF NMEN5 = 0 THEN INPUT "Output style must be selected - Press <RTN> to continue", ans$: GOTO 10
   IF NMEN6 = 0 THEN INPUT "Orientation must be selected - Press <RTN> to continue", ans$: GOTO 10
   IF NMEN7 = 0 THEN INPUT "Output variables must be set - Press <RTN> to continue", ans$: GOTO 10
   CA = LPC / LPA
   IF CA = 0 THEN CA = 1
   AC2 = 3 / (4 * CA * CA)' 3/4 of inverse square of c/a ratio
   ACSQ = 1 / (CA * CA)' inverse square of c/a ratio
   ASQ = LPA * LPA: BSQ = LPB * LPB: CSQ = LPC * LPC: COSBETA = COS(BETA): SINBETA = SIN(BETA)
   CASQO3 = CA * CA / 3
   IF DPSP$ = "T" THEN 380
   IF IPUH$ = "U" THEN GOSUB 1330
   IF IP$ = "1" AND IPUH$ <> "U" THEN GOSUB 1340
280 IF IP$ = "2" THEN GOSUB 1350 ' returns u3v3w3 ((Miller) integers) and h3k3l3 (not integers). uvw are Miller-Bravais zone axis for Hexagonal.
   GOSUB 8610 ' compute lhsphi
   IF st$ = "H" THEN HPSQ = NP * NP / (LPA * LPA * (U3 * U3 + V3 * V3 - U3 * V3) + LPC * LPC * W3 * W3)' H parameter squared
   IF st$ = "C" THEN HPSQ = NP * NP / (LPA * LPA * (U3 * U3 + V3 * V3 + W3 * W3))
   IF st$ = "T" THEN HPSQ = NP * NP / (LPA * LPA * (U3 * U3 + V3 * V3) + LPC * LPC * W3 * W3)
   IF st$ = "O" THEN HPSQ = NP * NP / (LPA * LPA * U3 * U3 + LPB * LPB * V3 * V3 + LPC * LPC * W3 * W3)
   IF st$ = "M" THEN HPSQ = NP * NP / (LPA * LPA * U3 * U3 + LPB * LPB * V3 * V3 + LPC * LPC * W3 * W3 + 2 * U3 * W3 * LPA * LPC * COSBETA)
   IF DPSP$ = "D" THEN G1SQ = 2 * K120 * SQR(HPSQ) - HPSQ: G2SQ = G1SQ * 2: GMAX = 1.5: G1MIN = G1SQ * .8: G1MAX = G1SQ * 1.2: G2MIN = G2SQ * .8: G2MAX = G2SQ * 1.2
   IF DPSP$ = "H" THEN G1SQ = 2 * K120 * SQR(HPSQ) - HPSQ: G1MIN = (SQR(G1SQ) - DRRL) * (SQR(G1SQ) - DRRL) * 1.001: G1MAX = (SQR(G1SQ) + DRRL) * (SQR(G1SQ) + DRRL) * .999: x = SQR(G1SQ) / (2 * K120): Y = ATN(x / SQR(1 - x * x)): HR = CAMLEN * SCALE  _
* TAN(2 * Y)
   IF AAP = 2 THEN 290
   RH = 0: RV = 0: LABEL$ = ""
   GOSUB 10380 ' Draw circle and tick marks or center spot
   IF IPUH$ = "U" THEN h = 0: k = 0: l = 0
   IF IPUH$ = "H" THEN h = H3: k = K3: l = L3
   IF st$ = "C" THEN h = H3: k = K3: l = L3
290 IF IP$ = "2" THEN h = U: k = V: i = -(h + k): l = W' center pole - plane with indices uvw
   H4 = h: K4 = k: L4 = l: I4 = -(H4 + K4)' used to prevent subsequent plotting
   IF SF% > 0 THEN GOSUB 9680 ' Only plot allowed reflections
   IF plot% > 0 THEN IF DPSP$ <> "A" AND DPSP$ <> "D" AND DPSP$ <> "H" THEN IF IP$ = "1" AND IPUH$ = "H" OR st$ = "C" THEN GOSUB 8970 ' Plots point and prints label at center
   h = H1: k = K1: l = L1: i = -(h + k)
   GOSUB 8660 ' Calculates radial distance RD from center
   IF RD > 1005 THEN PRINT "Second"; POLE$; " must be less than 90 degrees from first"; POLE$; " - re-enter": GOTO 210
   RHO = 0
   IF st$ = "H" THEN GOSUB 1530 ELSE GOSUB 1540 ' Calculates zone axis for pole at center and pole on x-axis
   U1 = U: V1 = V: W1 = W
   GOSUB 1560 ' compute lhsden
   IF IP$ = "1" THEN h = U + H1: k = V + K1: l = W + L1 ' Pole in upper right hand quadrant
   IF IP$ = "2" THEN h = H2: k = K2: : l = L2: i = -(h + k)
   IF st$ = "H" THEN GOSUB 1530 ELSE GOSUB 1540
   GOSUB 1550 ' Calculate rho1 between u1v1w1 and u2v2w2
   RHO1 = RHO
   U2 = U: V2 = V: W2 = W
   GOSUB 8430
   IF DPSP$ = "P" AND IP$ = "2" THEN GOSUB 1390 ' Plot input poles first
   flag = 1: LZ = 0' Allows coloring for HOLZ
   IF SPB$ = "P" THEN COLOR 4, 0, 0: CLS
   IF FACET$ = "F" THEN GOSUB 1130
   IF FACET$ = "F" THEN 300
   IF DPSP$ = "P" THEN GOSUB 1020
   IF DPSP$ = "D" THEN IF W3 <> 0 THEN GOSUB 1040 ELSE IF V3 <> 0 THEN GOSUB 1060 ELSE GOSUB 1080
   IF DPSP$ = "H" THEN IF W3 <> 0 THEN GOSUB 1100 ELSE IF V3 <> 0 THEN GOSUB 1110 ELSE GOSUB 1120
   IF DPSP$ = "A" THEN GOSUB 1140
   IF DPSP$ = "A" AND AAP = 1 THEN 190
   IF DPSP$ = "A" THEN 340
300 IF SPB$ = "P" THEN COLOR 3, 1, 1: CLS : PRINT : PRINT "         Stereographic Projection for Hewlett Packard Plotter"
   IF EXTRA$ = "N" THEN 320
   EXTRA$ = "EXTRA"
   PRINT : PRINT "Input additional "; POLE$; "s - input 0,0,0 to end"
310 PRINT : INPUT ; "h = ", h: INPUT ; " , k = ", k: INPUT " , l = ", l: i = -(h + k)
   IF h = 0 AND k = 0 AND l = 0 THEN 320
   NPM = 0: GOSUB 9360
   flag = 0
   IF DPSP$ <> "D" THEN GOSUB 1470 ELSE LZ = (h * U3 + k * V3 + l * W3) / NP: IF LZ > 2 OR LZ < 0 THEN 310 ELSE GOSUB 1500
   GOTO 310
320 IF LINCR% = 1 THEN 330 ' true for all except hexagonal
   h = 0: k = 0: i = 0
   flag = 1
   IF U3 = 0 AND V3 = 0 THEN 330
   FOR l% = -1 TO 1 STEP 2
   l = l%
   IF HOLZ$ = "H" THEN GOSUB 9340
   IF NPM > 2 THEN NPM = 2: plot% = 1
   GOSUB 1470
   NEXT
330 NPM = 0: GOSUB 9360
   IF st$ = "H" THEN U = (2 * U3 - V3): V = (2 * V3 - U3): W = W3 * 3: GOSUB 9570 ' return Miller-Bravais index for final label
   IF st$ = "H" THEN U3 = U: V3 = V: T3 = -(U3 + V3): W3 = W
   h = CINT(U3 * 100) / 100: k = CINT(V3 * 100) / 100: i = CINT(T3 * 100) / 100: l = CINT(W3 * 100) / 100
   IF DPSP$ = "D" THEN B$ = " Zone Axis" ELSE B$ = " Projection"
   IF SPB$ <> "P" THEN
      IF display$ = "EGA" THEN COLOR 3, 1
      IF st$ = "H" THEN PRINT B$; " "; U3; " "; V3; " "; -(U3 + V3); " "; W3 ELSE PRINT B$; " "; U3; " "; V3; " "; W3
   END IF
   IF SPB$ = "S" THEN INPUT "Press return to continue ", ans$: GOTO 370
   PRINT #1, "PA 800,900;"
   PRINT #1, "SR;" ' default type size
   LABEL$ = "TITLE": LF$ = "30": GOSUB 8970
   PRINT #1, "LB"; B$ + CR$ + LNF$ + LNF$ + C$
   LA$ = STR$(CINT(LPA * 1000) / 1000): LB$ = STR$(CINT(LPB * 1000) / 1000): LC$ = STR$(CINT(LPC * 1000) / 1000): BETA$ = STR$(CINT(BETA / PI * 1800) / 10)
   IF st$ = "M" THEN PRINT #1, "pr 100,0;lba = "; LA$ + CR$ + LNF$; "b = "; LB$ + CR$ + LNF$; "c = "; LC$ + CR$ + LNF$; "beta = "; BETA$ + C$; "pr 3,20;lbo"; C$
   IF DPSP$ = "P" THEN CA = CINT(CA * 1000) / 1000: IF st$ = "H" OR st$ = "T" THEN B$ = STR$(CA): GOSUB 10620
   IF DPSP$ = "D" OR DPSP$ = "H" THEN IF st$ = "H" OR st$ = "T" THEN PRINT #1, "pr 100,0;lba = "; LA$ + CR$ + LNF$; "c = "; LC$ + C$
   IF DPSP$ = "D" OR DPSP$ = "H" THEN IF st$ = "C" THEN PRINT #1, "pr 100,0;lba = "; LA$ + C$
   IF DPSP$ = "D" OR DPSP$ = "H" THEN KV$ = STR$(KV): PRINT #1, "PA 800,-900;lb"; KV$; " kV"; C$
   PRINT #1, "PA -1380,870;"
   PRINT #1, "LB"; CRYSTAL$ + C$
   IF HOLZ$ = "H" AND DPSP$ <> "H" THEN GOSUB 10640 ' label HOLZ colors
340 IF DPSP$ = "A" THEN IF BNS$ = "y" OR BNS$ = "Y" THEN GOSUB 10670 ' label rotation axis
350 GOSUB 10730
360 CLOSE #1
370 SCREEN 0, 0, 0: LF$ = "20": GOTO 10
380 IP$ = "1"' Start of unit triangle routine
   POLE$ = "pole"
   H1 = 0: K1 = 0: I1 = 0: L1 = 1: H2 = 0: K2 = 1: I2 = -1: L2 = 0: CA = 0: EXTRA$ = "Y": IF st$ = "C" THEN L2 = 1
   IF st$ = "C" THEN
385   INPUT "Input data manually or from a *.UTD file (M/U)?", utd$
   IF utd$ <> "m" AND utd$ <> "M" AND utd$ <> "u" AND utd$ <> "U" THEN 385
   ELSE
   utd$ = "M"
   END IF
   IF utd$ = "u" OR utd$ = "U" THEN
      opuh$ = "H": lab = 2
   ELSE
      IF st$ = "H" THEN
390      PRINT : INPUT "Choose to output direction indices, uvw, or plane normals, hkl (U/H): ", opuh$
	 IF opuh$ = "u" THEN opuh$ = "U"
	 IF opuh$ = "h" THEN opuh$ = "H"
	 IF opuh$ <> "U" AND opuh$ <> "H" THEN 390
      END IF
410   PRINT : PRINT "Do you wish to label "; POLE$; "s? (Y/N): "; : INPUT "", ans$
      IF ans$ = "Y" OR ans$ = "y" THEN lab = 1
      IF ans$ = "N" OR ans$ = "n" THEN lab = 2
      IF ans$ <> "Y" AND ans$ <> "y" AND ans$ <> "N" AND ans$ <> "n" THEN GOTO 410
      IF lab = 1 THEN INPUT "Angles or indices (A/I): ", ans$
      IF ans$ = "A" OR ans$ = "a" THEN lab = 3
   END IF
   GOSUB 1340 ' Calculate h3, k3, l3
   GOSUB 8610
   IF st$ = "H" THEN GOSUB 11190 ' Plot and label unit triangle - hexagonal
   IF st$ = "C" THEN GOSUB 11280 ' Cubic
   h = H1: k = K1: l = L1: i = -(h + k)
   RHO = 0
   IF st$ = "H" THEN GOSUB 1530 ELSE GOSUB 1540 ' Calculates zone axis for pole at center and pole on x-axis
   U1 = U: V1 = V: W1 = W
   GOSUB 1560
   h = U + H1: k = V + K1: l = W + L1 ' Pole in upper right hand quadrant
   IF st$ = "H" THEN GOSUB 1530 ELSE GOSUB 1540
   GOSUB 1550 ' Calculate rho1 between u1v1w1 and u2v2w2
   RHO1 = RHO
   U2 = U: V2 = V: W2 = W
   GOSUB 8430
   IF utd$ = "u" OR utd$ = "U" THEN
      CLS : PRINT : PRINT "         Stereographic Projection for Hewlett Packard Plotter"
      PRINT : PRINT
      FILES "*.utd"
      PRINT : PRINT : INPUT "Enter data file name: ", ut$
      PRINT : PRINT
      OPEN ut$ + ".utd" FOR INPUT AS #2
      INPUT #2, nmax
      FOR n = 1 TO nmax
      INPUT #2, h: INPUT #2, k: INPUT #2, l
      IF h * h + k * k + l * l < .99 OR h * h + k * k + l * l > 1.01 THEN PRINT h, k, l: END
      flag = 0: plot% = 1
      h = ABS(h): k = ABS(k): l = ABS(l)
      IF h > k THEN x = h: h = k: k = x
      IF k > l THEN x = k: k = l: l = x
      IF h > k THEN x = h: h = k: k = x
      GOSUB 1470
      NEXT
      CLOSE #2
   ELSE
      IF SPB$ = "P" THEN CLS : PRINT : PRINT "         Stereographic Projection for Hewlett Packard Plotter"
      EXTRA$ = "EXTRA"
      PRINT : PRINT "Input "; POLE$; "s - input 0,0,0 to end"
420   PRINT : INPUT ; "h = ", h: INPUT ; " , k = ", k: INPUT " , l = ", l: i = -(h + k)
      IF h = 0 AND k = 0 AND l = 0 THEN 440
      IF st$ = "H" AND l < 0 THEN h = -h: k = -k: l = -l: i = -(h + k)
430   flag = 0: plot% = 1
      GOSUB 1470
      IF st$ = "H" AND opuh$ = "U" THEN l = LTEMP
      IF st$ = "H" AND plot% = 0 THEN k = h: h = i: i = -(h + k): GOTO 430
      GOTO 420
   END IF
440 PRINT
   IF st$ = "H" THEN GOSUB 11240
   IF st$ = "C" THEN GOSUB 11320
   CLOSE #1
   GOTO 10
REM *************************************************************
REM BEGINNING OF SUBROUTINES
1000 PRINT : PRINT : INPUT "Set up plotter and press return to continue ", B$
1010 CLOSE #1: OPEN "COM2:9600,S,7,1,RS,CS65535,DS,CD" FOR RANDOM AS #1' This line is IBM specific and sets up communication parameters at 9600 baud for output port #1
     GOSUB 10200 ' Initialize plotter
     RETURN
1020 FOR h = -nmax TO nmax
     FOR k = -nmax TO nmax
     IF SPB$ = "P" THEN LOCATE 25, 65, 0: PERCENT% = (h + nmax + (k + nmax) / (2 * nmax + 1)) * 100 / (2 * nmax + 1): PRINT PERCENT%
     HK = h * U3 + k * V3
     FOR l% = -LMAX% TO LMAX% STEP LINCR%
     l = l%
    IF HOLZ$ = "H" THEN LZ = l * W3 + HK: IF LZ > NP * 2 THEN 1030 ELSE IF LZ < 0 THEN 1030 ELSE LZ = LZ / NP
    GOSUB 1460
1030 NEXT
     NEXT
     NEXT
     RETURN
1040 VW = V3 / W3 ' CBED routine - w3 <> 0
     FOR LZ = 0 TO 2
     GOSUB 9340
     Z = LZ * NP
     FOR h = -nmax TO nmax
     ZH = (Z - h * U3) / W3
     FOR k = -nmax TO nmax
     l = ZH - VW * k
    IF l <> CINT(l) THEN 1050
    IF l > LMAX% THEN 1050 ELSE IF l < -LMAX% THEN 1050
    GOSUB 1500 ' CBED routine
1050 NEXT
     NEXT
     NEXT
     RETURN
1060 FOR LZ = 0 TO 2 ' v3 <> 0: w3 = 0
     GOSUB 9340
     Z = LZ * NP
     FOR h = -nmax TO nmax
     k = (Z - h * U3) / V3
    IF k <> CINT(k) THEN 1070
     FOR l = -LMAX% TO LMAX%
    GOSUB 1500
1070 NEXT
     NEXT
     NEXT
     RETURN
1080 FOR LZ = 0 TO 2 ' v3 = 0:w3 = 0
     GOSUB 9340
     Z = LZ * NP
     h = Z / U3
    IF h <> CINT(h) THEN 1090
     FOR k = -nmax TO nmax
     FOR l = -LMAX% TO LMAX%
    GOSUB 1500
1090 NEXT
     NEXT
     NEXT
     RETURN
1100 LZ = 1 ' Start of HOLZ line subroutines
     VW = V3 / W3 ' HOLZ routine - w3 <> 0
     GOSUB 9340
     FOR h = -nmax TO nmax
     ZH = (NP - h * U3) / W3
     FOR k = -nmax TO nmax
     l = ZH - VW * k
    IF l = CINT(l) AND ABS(l) <= LMAX% THEN GOSUB 1510 ' HOLZ routine
    NEXT
     NEXT
     RETURN
1110 LZ = 1 ' w3 = 0
     GOSUB 9340 ' check "p"
     FOR h = -nmax TO nmax
     k = (NP - h * U3) / V3
    IF k = CINT(k) THEN
      FOR l = -LMAX% TO LMAX%
      GOSUB 1510 ' HOLZ routine
      NEXT
    END IF
    NEXT
     RETURN
1120 LZ = 1 ' v3 = 0: w3 = 0
     GOSUB 9340
     h = NP / U3
    IF h = CINT(h) THEN
      FOR k = -nmax TO nmax
      FOR l = -LMAX% TO LMAX%
      GOSUB 1510
    NEXT
     NEXT
    END IF
    RETURN
1130 h = 0: k = 0: l = 6: GOSUB 1460' Basal plane - facet plane poles for alumina
    l = -6: GOSUB 1460
    h = 1: k = 1: l = 0: GOSUB 1460 ' Prism planes
    h = -2: GOSUB 1460
    h = 1: k = -2: GOSUB 1460
    h = -1: k = -1: GOSUB 1460
    h = 2: k = -1: GOSUB 1460
    h = -1: k = 2: GOSUB 1460
    h = -1: k = 0: l = 2: GOSUB 1460 ' Rhombohedral planes
    h = 1: k = -1: GOSUB 1460
    h = 0: k = 1: GOSUB 1460
    h = 1: k = 0: l = -2: GOSUB 1460
    h = -1: k = 1: GOSUB 1460
    h = 0: k = -1: GOSUB 1460
    h = 1: k = 0: l = 4: GOSUB 1460  ' 10-14 type planes
    h = -1: k = 1: GOSUB 1460
    h = 0: k = -1: GOSUB 1460
    h = -1: k = 0: l = -4: GOSUB 1460
    h = 1: k = -1: GOSUB 1460
    h = 0: k = 1: GOSUB 1460
    h = 1: k = 0: l = 1: GOSUB 1460 ' 10-11 type planes
    h = -1: k = 1: GOSUB 1460
    h = 0: k = -1: GOSUB 1460
    k = 1: l = -1: GOSUB 1460
    h = -1: k = 0: GOSUB 1460
    h = 1: k = -1: GOSUB 1460
     RETURN
1140 IF st$ = "H" THEN GOSUB 1160 ELSE GOSUB 1170 ' Axis/angle pair subroutine
1150 IF AAP = 2 THEN
      PRINT : INPUT "Do you wish to locate a rotation axis? (Y/N): ", BNS$
      IF BNS$ = "y" OR BNS$ = "Y" THEN GOSUB 1200 ELSE IF BNS$ <> "n" AND BNS$ <> "N" THEN 1150
    END IF
    RETURN
1160 NPN = 0: h = -1: k = -1: i = 2: l = 0: GOSUB 1470' plot pole
    GOSUB 1180
     NC = NC + 1
    h = -1: k = 2: i = -1: l = 0: GOSUB 1470
    GOSUB 1180
     NC = NC + 1
    h = 2: k = -1: i = -1: l = 0: GOSUB 1470
    GOSUB 1180
     NC = NC + 1
    h = 0: k = 0: i = 0: l = 1: GOSUB 1470
    GOSUB 1180
     NC = NC + 1
    NPN = 1: h = 1: k = -2: i = 1: l = 0: GOSUB 1470
    GOSUB 1180
     NC = NC + 1
    h = 1: k = 1: i = -2: l = 0: GOSUB 1470
    GOSUB 1180
     NC = NC + 1
    h = -2: k = 1: i = 1: l = 0: GOSUB 1470
    GOSUB 1180
     NC = NC + 1
    h = 0: k = 0: i = 0: l = -1: GOSUB 1470
    GOSUB 1180
     NC = NC + 1
     RETURN
1170 NPN = AAP - 1: h = 0: k = 0: l = 1: GOSUB 1470
    GOSUB 1180
     NC = NC + 1
    l = -1: GOSUB 1470
    GOSUB 1180
     NC = NC + 1
    k = 1: l = 0: GOSUB 1470
    GOSUB 1180
     NC = NC + 1
    k = -1: GOSUB 1470
    GOSUB 1180
     NC = NC + 1
    h = 1: k = 0: GOSUB 1470
    GOSUB 1180
     NC = NC + 1
    h = -1: GOSUB 1470
    GOSUB 1180
     NC = NC + 3
     RETURN
1180 a(NC, 0) = RHO + THETA: a(NC, 1) = PHI' conversion to standard cubic hkl
     IF COS(a(NC, 0)) > SN THEN a = (1 - COS(a(NC, 0)) * COS(a(NC, 0))) / (COS(a(NC, 0)) * COS(a(NC, 0))): IF ABS(a) < SN THEN a = 0
     IF COS(a(NC, 1)) > SN THEN B = (1 - COS(a(NC, 1)) * COS(a(NC, 1))) / (COS(a(NC, 1)) * COS(a(NC, 1))): IF ABS(B) < SN THEN B = 0
    IF COS(a(NC, 0)) <= SN THEN a(NC, 2) = SQR(B / (1 + B)): a(NC, 3) = 0: a(NC, 4) = SQR(1 / (1 + B)): GOTO 1190
    IF COS(a(NC, 1)) <= SN THEN a(NC, 2) = SQR(a / (1 + a)): a(NC, 3) = SQR(1 / (1 + a)): a(NC, 4) = 0: GOTO 1190
     a(NC, 2) = SQR((B * a) / ((1 + B) * (1 + a))): a(NC, 3) = SQR(B / ((1 + B) * (1 + a))): a(NC, 4) = SQR(1 / (1 + B))' equivalent h k l in standard cubic projection
1190 IF a(NC, 0) > 2 * PI THEN a(NC, 2) = -a(NC, 2)
     IF a(NC, 0) > 0 AND a(NC, 0) < PI THEN a(NC, 2) = -a(NC, 2)
     IF a(NC, 0) > PI / 2 AND a(NC, 0) < 3 * PI / 2 THEN a(NC, 3) = -a(NC, 3)
     IF a(NC, 1) > PI / 2 THEN a(NC, 4) = -a(NC, 4)
     RETURN
1200 IF st$ <> "H" THEN 1240
1210 CLS : PRINT "Select poles for rotation axis"' hexagonal
     PRINT : PRINT "Pole 1: -1 -1  2  0": PRINT "Pole 2: -1  2 -1  0": PRINT "Pole 3:  2 -1 -1  0": PRINT "Pole 4:  0  0  0  1": PRINT "Pole 5:  1 -2  1  0": PRINT "Pole 6:  1  1 -2  0": PRINT "Pole 7: -2  1  1  0": PRINT "Pole 8:  0  0  0 -1"
     PRINT : INPUT "Select pole a from grain 1 (1-8) ", B%(0): INPUT "Select pole b from grain 1 (1-8) ", B%(1): INPUT "Select pole a' from grain 2 (1-8) ", B%(2): INPUT "Select pole b' from grain 2 (1-8) ", B%(3)
    IF B%(0) < 1 OR B%(0) > 8 THEN 1210
    IF B%(1) < 1 OR B%(1) > 8 THEN 1210
    IF B%(2) < 1 OR B%(2) > 8 THEN 1210
    IF B%(3) < 1 OR B%(3) > 8 THEN 1210
    GOSUB 1270
1220 PRINT : INPUT "Do you want to plot another? (Y/N): ", ans$
    IF ans$ = "y" OR ans$ = "Y" THEN GOTO 1210 ELSE IF ans$ = "n" OR ans$ = "N" THEN 1230 ELSE 1220
1230 GOTO 1260
1240 CLS : PRINT "Select poles for rotation axis"' cubic and tetragonal
    PRINT : PRINT "Pole 1:  0  0  1": PRINT "Pole 2:  0  0 -1": PRINT "Pole 3:  0  1  0": PRINT "Pole 4:  0 -1  0": PRINT "Pole 5:  1  0  0": PRINT "Pole 6: -1  0  0"
     PRINT : INPUT "Select pole a from grain 1 (1-6) ", B%(0): INPUT "Select pole b from grain 1 (1-6) ", B%(1): INPUT "Select pole a' from grain 2 (1-6) ", B%(2): INPUT "Select pole b' from grain 2 (1-6) ", B%(3)
    GOSUB 1270
1250 PRINT : INPUT "Do you want to plot another? (Y/N): ", ans$
    IF ans$ = "y" OR ans$ = "Y" THEN GOTO 1240 ELSE IF ans$ = "n" OR ans$ = "N" THEN 1260 ELSE 1250
1260 RETURN
1270 C(0, 0) = a(B%(0) - 1, 2) - a(B%(2) + 7, 2): C(0, 1) = a(B%(0) - 1, 3) - a(B%(2) + 7, 3): C(0, 2) = a(B%(0) - 1, 4) - a(B%(2) + 7, 4)' (a - a prime) standard cubic hkl
     C(1, 0) = a(B%(1) - 1, 2) - a(B%(3) + 7, 2): C(1, 1) = a(B%(1) - 1, 3) - a(B%(3) + 7, 3): C(1, 2) = a(B%(1) - 1, 4) - a(B%(3) + 7, 4)' (b - b prime) standard cubic hkl
     C(2, 0) = C(0, 1) * C(1, 2) - C(0, 2) * C(1, 1): C(2, 1) = C(0, 2) * C(1, 0) - C(0, 0) * C(1, 2): C(2, 2) = C(0, 0) * C(1, 1) - C(0, 1) * C(1, 0)
     IF C(2, 2) < 0 THEN C(2, 0) = -C(2, 0): C(2, 1) = -C(2, 1): C(2, 2) = -C(2, 2)
    IF st$ = "C" OR st$ = "T" THEN 1280
     D(0) = (a(2, 2) * C(2, 0) + a(2, 3) * C(2, 1) + a(2, 4) * C(2, 2)) / SQR((a(2, 2) * a(2, 2) + a(2, 3) * a(2, 3) + a(2, 4) * a(2, 4)) * (C(2, 0) * C(2, 0) + C(2, 1) * C(2, 1) + C(2, 2) * C(2, 2)))' cos angle hkl^(2-1-10)
     D(1) = (a(1, 2) * C(2, 0) + a(1, 3) * C(2, 1) + a(1, 4) * C(2, 2)) / SQR((a(1, 2) * a(1, 2) + a(1, 3) * a(1, 3) + a(1, 4) * a(1, 4)) * (C(2, 0) * C(2, 0) + C(2, 1) * C(2, 1) + C(2, 2) * C(2, 2)))' cos angle hkl^(-12-10)
     D(2) = (a(3, 2) * C(2, 0) + a(3, 3) * C(2, 1) + a(3, 4) * C(2, 2)) / SQR((a(3, 2) * a(3, 2) + a(3, 3) * a(3, 3) + a(3, 4) * a(3, 4)) * (C(2, 0) * C(2, 0) + C(2, 1) * C(2, 1) + C(2, 2) * C(2, 2)))' cos angle hkl^(0001)
     D(3) = -(D(0) + D(1)) ' i index - hexagonal
     D(4) = D(2) * CA ' l index - hexagonal
    GOTO 1290
1280 D(0) = (a(4, 2) * C(2, 0) + a(4, 3) * C(2, 1) + a(4, 4) * C(2, 2)) / SQR((a(4, 2) * a(4, 2) + a(4, 3) * a(4, 3) + a(4, 4) * a(4, 4)) * (C(2, 0) * C(2, 0) + C(2, 1) * C(2, 1) + C(2, 2) * C(2, 2)))' cos angle hkl^(100)
     D(1) = (a(2, 2) * C(2, 0) + a(2, 3) * C(2, 1) + a(2, 4) * C(2, 2)) / SQR((a(2, 2) * a(2, 2) + a(2, 3) * a(2, 3) + a(2, 4) * a(2, 4)) * (C(2, 0) * C(2, 0) + C(2, 1) * C(2, 1) + C(2, 2) * C(2, 2)))' cos angle hkl^(010)
     D(3) = 1 ' redundant index
     D(4) = (a(0, 2) * C(2, 0) + a(0, 3) * C(2, 1) + a(0, 4) * C(2, 2)) / SQR((a(0, 2) * a(0, 2) + a(0, 3) * a(0, 3) + a(0, 4) * a(0, 4)) * (C(2, 0) * C(2, 0) + C(2, 1) * C(2, 1) + C(2, 2) * C(2, 2)))' cos angle hkl^(001)-cubic
     IF st$ = "T" THEN D(4) = D(4) * CA ' tetragonal - l index
1290 IF ABS(D(0)) < ABS(D(1)) AND ABS(D(0)) < ABS(D(3)) AND ABS(D(0)) < ABS(D(4)) AND ABS(D(0)) > .1 THEN D(5) = ABS(D(0)): FOR J% = 0 TO 4: D(J%) = D(J%) / D(5): NEXT
     IF ABS(D(1)) < ABS(D(0)) AND ABS(D(1)) < ABS(D(3)) AND ABS(D(1)) < ABS(D(4)) AND ABS(D(1)) > .1 THEN D(5) = ABS(D(1)): FOR J% = 0 TO 4: D(J%) = D(J%) / D(5): NEXT
     IF ABS(D(3)) < ABS(D(0)) AND ABS(D(3)) < ABS(D(1)) AND ABS(D(3)) < ABS(D(4)) AND ABS(D(3)) > .1 THEN D(5) = ABS(D(3)): FOR J% = 0 TO 4: D(J%) = D(J%) / D(5): NEXT
     IF ABS(D(4)) < ABS(D(0)) AND ABS(D(4)) < ABS(D(1)) AND ABS(D(4)) < ABS(D(3)) AND ABS(D(4)) > .1 THEN D(5) = ABS(D(4)): FOR J% = 0 TO 4: D(J%) = D(J%) / D(5): NEXT
     FOR J% = 0 TO 4
     D(J%) = CINT(D(J%) * 1000) / 1000
     NEXT
     IF st$ = "H" THEN PRINT : PRINT "h = "; D(0), "k = "; D(1), "i = "; D(3), "l = "; D(4): PRINT
     IF st$ = "C" OR st$ = "T" THEN PRINT : PRINT "h = "; D(0), "k = "; D(1), "l = "; D(4): PRINT
     COSPHI = C(2, 2) / SQR(C(2, 0) * C(2, 0) + C(2, 1) * C(2, 1) + C(2, 2) * C(2, 2))
    IF COSPHI >= 1 AND COSPHI < 1.001 THEN PHI = 0: RD = 0: GOTO 1300
    IF COSPHI <= -1 AND COSPHI > -1.001 THEN PHI = PI: RD = 100000: GOTO 1300
    PHI = PI / 2 - ATN(COSPHI / SQR(1 - COSPHI * COSPHI))
    RD = 1000 * TAN(PHI / 2)
1300 IF C(2, 1) = 0 THEN RHO = 0: GOTO 1310
     COSRHO = C(2, 1) / SQR(C(2, 0) * C(2, 0) + C(2, 1) * C(2, 1))
    IF COSRHO >= 1 AND COSRHO < 1.001 THEN RHO = 0: GOTO 1310
    IF COSRHO <= -1 AND COSRHO > -1.001 THEN RHO = PI: GOTO 1310
    RHO = PI / 2 - ATN(COSRHO / SQR(1 - COSRHO * COSRHO)): IF C(2, 0) > 0 THEN RHO = 2 * PI - RHO
1310 RH = RD * COS(RHO) ' Calculate plotting co-ordinates
     RV = RD * SIN(RHO)
     IF RH < 1 AND RH > -1 THEN RH = 0
     IF RV < 1 AND RV > -1 THEN RV = 0
     GOSUB 10770
     C(3, 0) = a(B%(0) - 1, 3) * C(2, 2) - a(B%(0) - 1, 4) * C(2, 1)' a ^ axis (cross product)
     C(3, 1) = a(B%(0) - 1, 4) * C(2, 0) - a(B%(0) - 1, 2) * C(2, 2)
     C(3, 2) = a(B%(0) - 1, 2) * C(2, 1) - a(B%(0) - 1, 3) * C(2, 0)
     C(4, 0) = a(B%(2) + 7, 3) * C(2, 2) - a(B%(2) + 7, 4) * C(2, 1)' a prime ^ axis (cross product)
     C(4, 1) = a(B%(2) + 7, 4) * C(2, 0) - a(B%(2) + 7, 2) * C(2, 2)
     C(4, 2) = a(B%(2) + 7, 2) * C(2, 1) - a(B%(2) + 7, 3) * C(2, 0)
     COSALPHA = (C(3, 0) * C(4, 0) + C(3, 1) * C(4, 1) + C(3, 2) * C(4, 2)) / SQR((C(3, 0) * C(3, 0) + C(3, 1) * C(3, 1) + C(3, 2) * C(3, 2)) * (C(4, 0) * C(4, 0) + C(4, 1) * C(4, 1) + C(4, 2) * C(4, 2)))' rotation angle
    IF COSALPHA >= 1 AND COSALPHA < 1.001 THEN ALPHA = 0: GOTO 1320
    IF COSALPHA <= -1 AND COSALPHA > -1.001 THEN ALPHA = PI: GOTO 1320
    ALPHA = PI / 2 - ATN(COSALPHA / SQR(1 - COSALPHA * COSALPHA))
1320 ALPHA% = ALPHA * 1800 / PI
     ALPHA$ = STR$(ALPHA% / 10): GOSUB 10790
     PRINT "The rotation angle = "; ALPHA% / 10; " degrees"
     RETURN
1330 IF st$ = "T" THEN H3 = U3: K3 = V3: L3 = W3 / ACSQ: I3 = -(H3 + K3): H1 = U4: K1 = V4: I1 = -(H1 + K1): L1 = W4 / ACSQ   ' Pole at center of projection and on right hand axis
    IF st$ <> "H" THEN RETURN
     H3 = U3: K3 = V3: L3 = W3 / (2 * AC2): I3 = -(H3 + K3): H1 = U4: K1 = V4: I1 = -(H1 + K1): L1 = W4 / (2 * AC2): H2 = H1: K2 = K1: I2 = I1: L2 = L1    ' Pole at center of projection and on right hand axis - uvw input
     U = 2 * U3 + V3: V = 2 * V3 + U3: W = W3: GOSUB 9620
     U3 = U: V3 = V: T3 = -(U3 + V3): W3 = W
    RETURN
1340 H3 = H1: L3 = L1: I3 = I1: K3 = K1: H1 = H2: L1 = L2: I1 = I2: K1 = K2 ' See Figure 1
     IF st$ = "C" THEN U3 = H3: V3 = K3: W3 = L3 ' Beam direction - used in HOLZ subroutines and final label
     IF st$ = "H" THEN U3 = H3: V3 = K3: W3 = L3 * 2 * AC2: T3 = -(U3 + V3) ' Beam direction - used in final label
     IF st$ = "T" THEN U3 = H3: V3 = K3: W3 = L3 * ACSQ: T3 = -(U3 + V3) ' Beam direction - used in final label
     U = U3: V = V3: W = W3
     RETURN
1350 h = H2: k = K2: l = L2: H3 = H1: K3 = K1: L3 = L1 ' temporary assignment
    GOSUB 1540 ' Calculate beam direction, u3v3w3 - Miller indices
1360 GOSUB 9130 ' reduce to smallest integers - returns NN
    IF NN > 0 THEN U = U / NN: V = V / NN: W = W / NN: T = -(U + V): GOTO 1360
     U3 = U: V3 = V: W3 = W: T3 = -(U3 + V3) ' Miller indices used for all crystal systems
    IF st$ = "H" THEN GOSUB 1530 ' Miller-Bravais indices
1370 IF st$ = "H" THEN GOSUB 9130 ' reduce to smallest integers - returns NN
    IF st$ = "H" AND NN > 0 THEN U = U / NN: V = V / NN: W = W / NN: T = -(U + V): GOTO 1370
     IF st$ = "C" THEN GOSUB 9660 ' indices of plane normal to u3v3w3
     IF st$ = "H" THEN GOSUB 9420
     IF st$ = "T" THEN GOSUB 9440
     IF st$ = "M" THEN GOSUB 9460
1380 IF SF% = 0 THEN GOSUB 9210
    IF SF% = 0 AND NN > 0 THEN h = h / NN: k = k / NN: l = l / NN: i = -(h + k) / NN: GOTO 1380
     H3 = h: K3 = k: L3 = l: I3 = -(H3 + K3)
     RETURN
1390 IF SF% = 0 THEN
      h = H4: k = K4: l = L4: i = -(h + k) ' Plot poles near center
      IF st$ = "H" THEN GOSUB 1530 ELSE GOSUB 1540 ' Calculate zone axis
      IF st$ <> "C" THEN IF U <> 0 OR V <> 0 OR W <> 0 THEN GOSUB 1490 ELSE RD = 0: GOSUB 8970 ' Plot point near center
      IF st$ = "C" THEN GOSUB 8970
    END IF
    IF IPUH$ = "U" THEN RETURN
     h = H1: k = K1: l = L1: i = -(h + k) ' Plot input poles
    IF SF% > 0 THEN
      GOSUB 9680
      IF plot% = 0 THEN 1420 ELSE 1410
    END IF
1400 GOSUB 9210 ' reduce hkl's
    IF NN > 0 THEN h = h / NN: k = k / NN: l = l / NN: i = -(h + k): GOTO 1400
1410 GOSUB 1480 ' plot poles
1420 H1 = h: K1 = k: L1 = l
     h = H2: k = K2: l = L2: i = -(h + k) ' Plot input poles
    IF SF% > 0 THEN
      GOSUB 9680
      IF plot% = 0 THEN 1450 ELSE 1440
    END IF
1430 IF DPSP$ <> "D" THEN GOSUB 9210
    IF NN > 0 THEN h = h / NN: k = k / NN: l = l / NN: i = -(h + k): GOTO 1430
1440 GOSUB 1480 ' plot poles
1450 H2 = h: K2 = k: L2 = l
     RETURN
1460 plot% = 1: IF h = 0 AND k = 0 AND l = 0 THEN RETURN ' Main subroutine
     IF st$ <> "H" THEN i = 0 ELSE i = -(h + k): IF i > nmax OR i < -nmax THEN RETURN
    IF FACET$ = "F" THEN 1480
     IF SF% = 0 THEN GOSUB 9210 ' Check indices for divisibility
     IF SF% = 0 AND NN > 0 THEN RETURN
     IF SF% > 0 THEN GOSUB 9680 ' Only plot allowed reflections
     IF plot% = 0 THEN RETURN
1470 IF opuh$ <> "U" THEN ' Second entry point
      IF DPSP$ = "A" THEN 1480
      IF IP$ = "2" THEN
	 IF h = H1 AND k = K1 AND l = L1 THEN RETURN ' already plotted
	 IF h = H2 AND k = K2 AND l = L2 THEN RETURN ' already plotted
      END IF
      IF h = H4 AND k = K4 AND l = L4 THEN RETURN ' already plotted
    END IF
    IF st$ = "H" AND EXTRA$ <> "EXTRA" THEN IF i > nmax OR i < -nmax THEN RETURN
1480 IF opuh$ = "U" THEN LTEMP = l: l = 2 * CA * CA * LTEMP / 3 ' Find equivalent plane indices
   IF st$ = "H" THEN GOSUB 1530 ELSE GOSUB 1540 ' Calculate zone axis
1490 GOSUB 8660 ' Calculates radial distance
    IF DPSP$ = "P" THEN IF RD > 1005 OR PHI = PI THEN RETURN ' plot would be outside projection
    GOSUB 1550 ' Calculates angle RHO between zone axes
     GOSUB 8350 ' Calculate zone axis with pole 3
    IF U = 0 AND V = 0 AND W = 0 THEN 1520
     GOSUB 9090 ' Determine sign of rho
    IF DPSP$ = "A" THEN IF RD > 1005 OR PHI = PI THEN RETURN ' plot would be outside projection
     IF st$ = "H" AND lab > 0 AND RHO > PI * 2 / 3 THEN plot% = 0: RETURN
     IF st$ = "H" AND lab > 0 THEN IF RHO > PI / 3 THEN RHO = 2 * PI / 3 - RHO
    IF flag <> 0 THEN ' Prevents coloring for HOLZ
      IF HOLZ$ = "" THEN IF DPSP$ = "A" THEN GOSUB 10750 ELSE GOSUB 9290 ' Change pen color
      IF HOLZ$ = "H" THEN GOSUB 9340
    END IF
    GOSUB 8770 ' Calculates RH and RV
     IF opuh$ = "U" THEN l = LTEMP ' Convert back to direction indices
     GOSUB 8970 ' Plot and label
     RETURN
1500 plot% = 1: IF st$ <> "H" THEN i = 0 ELSE i = -(h + k): IF i > nmax OR i < -nmax THEN RETURN ' CBED subroutine
     IF SF% > 0 THEN GOSUB 9680 ' Only plot allowed reflections
     IF plot% = 0 THEN RETURN
     IF st$ = "H" THEN GSQ = (4 * (h * h + h * k + k * k) / (3 * LPA * LPA)) + (l * l / (LPC * LPC))' Hexagonal diffraction pattern
     IF st$ = "C" THEN GSQ = (h * h + k * k + l * l) / (LPA * LPA)
     IF st$ = "T" THEN GSQ = (h * h + k * k) / (LPA * LPA) + l * l / (LPC * LPC)
     IF st$ = "O" THEN GSQ = h * h / (LPA * LPA) + k * k / (LPB * LPB) + l * l / (LPC * LPC)
     IF st$ = "M" THEN GSQ = h * h / (LPA * LPA * SINBETA * SINBETA) + k * k / (LPB * LPB) + l * l / (LPC * LPC * SINBETA * SINBETA) - (2 * h * l * COSBETA) / (LPA * LPC * SINBETA * SINBETA)
     IF GSQ = 0 THEN RETURN
     GHSQ = GSQ - (HPSQ * LZ * LZ)
     IF LZ = 0 THEN IF GHSQ > GMAX THEN RETURN
     IF LZ = 1 THEN IF GHSQ < G1MIN OR GHSQ > G1MAX THEN RETURN
     IF LZ = 2 THEN IF GHSQ < G2MIN OR GHSQ > G2MAX THEN RETURN
     x = SQR(GSQ) / (2 * K120): Y = ATN(x / SQR(1 - x * x)): RD = CAMLEN * SCALE * TAN(2 * Y)
    IF st$ = "H" THEN GOSUB 1530 ELSE GOSUB 1540 ' Calculate zone axis
    GOSUB 1550 ' Calculates angle RHO between zone axes
     GOSUB 8350 ' Calculate zone axis with pole 3
     GOSUB 9090 ' Determine sign of rho
     GOSUB 8770 ' Calculates RH and RV
     GOSUB 8970 ' Plot and label
     RETURN
1510 plot% = 1:  IF st$ <> "H" THEN i = 0 ELSE i = -(h + k): IF i > nmax OR i < -nmax THEN RETURN' HOLZ subroutine
     IF SF% > 0 THEN GOSUB 9680 ' Only plot allowed reflections
     IF plot% = 0 THEN RETURN
     IF st$ = "H" THEN GSQ = (4 * (h * h + h * k + k * k) / (3 * LPA * LPA)) + (l * l / (LPC * LPC))' Hexagonal diffraction pattern
     IF st$ = "C" THEN GSQ = (h * h + k * k + l * l) / (LPA * LPA)
     IF st$ = "T" THEN GSQ = (h * h + k * k) / (LPA * LPA) + l * l / (LPC * LPC)
     IF st$ = "O" THEN GSQ = h * h / (LPA * LPA) + k * k / (LPB * LPB) + l * l / (LPC * LPC)
     IF st$ = "M" THEN GSQ = h * h / (LPA * LPA * SINBETA * SINBETA) + k * k / (LPB * LPB) + l * l / (LPC * LPC * SINBETA * SINBETA) - (2 * h * l * COSBETA) / (LPA * LPC * SINBETA * SINBETA)
     IF GSQ = 0 THEN RETURN
     IF GSQ < G1MIN OR GSQ > G1MAX THEN RETURN
    x = SQR(GSQ) / (2 * K120): Y = ATN(x / SQR(1 - x * x)): RD = CAMLEN * SCALE * TAN(2 * Y)
    IF st$ = "H" THEN GOSUB 1530 ELSE GOSUB 1540 ' Calculate zone axis
    GOSUB 1550 ' Calculates angle RHO between zone axes
     GOSUB 8350 ' Calculate zone axis with pole 3
     GOSUB 9090 ' Determine sign of rho
     GOSUB 8770 ' RH and RV
     GOSUB 8820 ' HOLZ line positions
    LOCATE 1, 1: PRINT h; k; l
     GOSUB 10570 ' Plot lines
     RETURN
1520 GOSUB 10750
     RHO = 0: PHI = 0: RD = 0: RH = 0: RV = 0
     GOSUB 10750
     GOSUB 8970
     RETURN
1530 U = l * (2 * K3 + H3) - L3 * (2 * k + h)' Calculates zone axis - hexagonal - Miller Bravais Notation
     V = L3 * (2 * h + k) - l * (2 * H3 + K3)
     W = 3 * (H3 * k - h * K3)
     T = -(U + V)
     RETURN
1540 U = K3 * l - k * L3 ' Calculates zone axis - all systems in Miller notation
     V = L3 * h - l * H3
     W = H3 * k - h * K3
     RETURN
1550 IF st$ = "C" THEN GOSUB 8320
    IF st$ = "T" THEN GOSUB 1580
    IF st$ = "H" THEN GOSUB 1570
    IF st$ = "M" THEN GOSUB 1590
    IF COSRHO >= 1 AND COSRHO < 1.001 THEN
      RHO = 0
    ELSEIF COSRHO <= -1 AND COSRHO > -1.001 THEN
      RHO = PI
    ELSE
      RHO = PI / 2 - ATN(COSRHO / SQR(1 - COSRHO * COSRHO))
    END IF
    RETURN
1560 IF st$ = "H" THEN LHSDEN = U1 * U1 + V1 * V1 + U1 * V1 + CASQO3 * W1 * W1
     IF st$ = "T" THEN LHSDEN = ACSQ * (U1 * U1 + V1 * V1) + W1 * W1
     IF st$ = "M" THEN LHSDEN = ASQ * U1 * U1 + BSQ * V1 * V1 + CSQ * W1 * W1 + 2 * LPA * LPC * U1 * W1 * COSBETA
     IF st$ = "C" THEN LHSDEN = U1 * U1 + V1 * V1 + W1 * W1
     IF st$ = "O" THEN LHSDEN = ASQ * U1 * U1 + BSQ * V1 * V1 + CSQ * W1 * W1
     RETURN
1570 IF U * U + V * V + U * V + CASQO3 * W * W = 0 THEN COSRHO = 0: RETURN' hexagonal
     COSRHO = (U1 * U + V1 * V + (U1 * V + V1 * U) / 2 + CASQO3 * W1 * W) / SQR(LHSDEN * (U * U + V * V + U * V + CASQO3 * W * W))
     RETURN
1580 IF U * U + V * V + W * W = 0 THEN COSRHO = 0: RETURN
     COSRHO = (ACSQ * (U1 * U + V1 * V) + W1 * W) / SQR(LHSDEN * (ACSQ * (U * U + V * V) + W * W))
     RETURN
1590 DENOMSQ = LHSDEN * (ASQ * U * U + BSQ * V * V + CSQ * W * W + 2 * LPA * LPC * U * W * COSBETA)
     IF DENOMSQ = 0 THEN COSRHO = 0: RETURN ' monoclinic
     COSRHO = (ASQ * U1 * U + BSQ * V1 * V + CSQ * W1 * W + LPA * LPC * (W1 * U + U1 * W) * COS(BETA)) / SQR(DENOMSQ)
     RETURN
8320 IF U * U + V * V + W * W = 0 THEN COSRHO = 0: RETURN' cubic
     COSRHO = (U1 * U + V1 * V + W1 * W) / SQR(LHSDEN * (U * U + V * V + W * W))
     RETURN
8350 IF st$ = "C" THEN GOSUB 8580
     IF st$ = "T" THEN GOSUB 8510
     IF st$ = "H" THEN GOSUB 8480
     IF st$ = "M" THEN GOSUB 8540
     IF COSRHO2 >= 1 AND COSRHO2 < 1.001 THEN RHO2 = 0: GOTO 8420
     IF COSRHO2 <= -1 AND COSRHO2 > -1.001 THEN RHO2 = PI: GOTO 8420
     RHO2 = PI / 2 - ATN(COSRHO2 / SQR(1 - COSRHO2 * COSRHO2))
8420 RETURN
8430 IF st$ = "H" THEN LHSDEN2 = U2 * U2 + V2 * V2 + U2 * V2 + CASQO3 * W2 * W2
     IF st$ = "T" THEN LHSDEN2 = ACSQ * (U2 * U2 + V2 * V2) + W2 * W2
     IF st$ = "M" THEN LHSDEN2 = ASQ * U2 * U2 + BSQ * V2 * V2 + CSQ * W2 * W2 + 2 * LPA * LPC * U2 * W2 * COSBETA
     IF st$ = "C" THEN LHSDEN2 = U2 * U2 + V2 * V2 + W2 * W2
     IF st$ = "O" THEN LHSDEN2 = ASQ * U2 * U2 + BSQ * V2 * V2 + CSQ * W2 * W2
     RETURN
8480 IF U * U + V * V + U * V + CASQO3 * W * W = 0 THEN COSRHO2 = 0: RETURN' hexagonal
     COSRHO2 = (U2 * U + V2 * V + (U2 * V + V2 * U) / 2 + CASQO3 * W2 * W) / SQR(LHSDEN2 * (U * U + V * V + U * V + CASQO3 * W * W))
     RETURN
8510 IF U * U + V * V + W * W = 0 THEN COSRHO2 = 0: RETURN
     COSRHO2 = (ACSQ * (U2 * U + V2 * V) + W2 * W) / SQR(LHSDEN2 * (ACSQ * (U * U + V * V) + W * W))
     RETURN
8540 DENOMSQ = LHSDEN2 * (ASQ * U * U + BSQ * V * V + CSQ * W * W + 2 * LPA * LPC * U * W * COSBETA)
     IF DENOMSQ = 0 THEN COSRHO2 = 0: RETURN ' monoclinic
     COSRHO2 = (ASQ * U2 * U + BSQ * V2 * V + CSQ * W2 * W + LPA * LPC * (W2 * U + U2 * W) * COSBETA) / SQR(DENOMSQ)
     RETURN
8580 IF U * U + V * V + W * W = 0 THEN COSRHO2 = 0: RETURN' cubic
     COSRHO2 = (U2 * U + V2 * V + W2 * W) / SQR(LHSDEN2 * (U * U + V * V + W * W))
     RETURN
8610 IF st$ = "H" THEN LHSPHI = H3 * H3 + K3 * K3 + H3 * K3 + AC2 * L3 * L3
     IF st$ = "C" THEN LHSPHI = H3 * H3 + K3 * K3 + L3 * L3
     IF st$ = "T" THEN LHSPHI = H3 * H3 + K3 * K3 + L3 * L3 * ACSQ
     IF st$ = "M" THEN LHSPHI = H3 * H3 / ASQ + K3 * K3 * SINBETA * SINBETA / BSQ + L3 * L3 / CSQ - 2 * H3 * L3 / (LPA * LPC) * COSBETA
     RETURN
8660 IF st$ = "H" THEN COSPHI = (H3 * h + K3 * k + (H3 * k + K3 * h) / 2 + AC2 * L3 * l) / SQR(LHSPHI * (h * h + k * k + h * k + AC2 * l * l))' hexagonal
     IF st$ = "C" THEN COSPHI = (H3 * h + K3 * k + L3 * l) / SQR(LHSPHI * (h * h + k * k + l * l))' cubic
     IF st$ = "T" THEN COSPHI = ((H3 * h) + (K3 * k) + (L3 * l * ACSQ)) / SQR(LHSPHI * (h * h + k * k + l * l * ACSQ))' tetragonal
     IF st$ <> "M" THEN 8720
     DENOMSQ = LHSPHI * (h * h / ASQ + k * k * SINBETA * SINBETA / BSQ + l * l / CSQ - 2 * h * l * COSBETA / (LPA * LPC))
     COSPHI = (H3 * h / ASQ + K3 * k * SINBETA * SINBETA / BSQ + L3 * l / CSQ - (L3 * h + l * H3) * COSBETA / (LPA * LPC)) / SQR(DENOMSQ)
8720 IF COSPHI >= 1 AND COSPHI < 1.001 THEN PHI = 0: RD = 0: GOTO 8760
     IF COSPHI <= -1 AND COSPHI > -1.001 THEN PHI = PI: RD = 100000: GOTO 8760
     PHI = PI / 2 - ATN(COSPHI / SQR(1 - COSPHI * COSPHI))
     RD = 1000 * TAN(PHI / 2)
8760 RETURN
8770 RHOTHETA = RHO + THETA
     RH = RD * COS(RHOTHETA) ' Calculate plotting co-ordinates
     RV = RD * SIN(RHOTHETA)
     IF RH < 1 AND RH > -1 THEN RH = 0
     IF RV < 1 AND RV > -1 THEN RV = 0
     RETURN
8820 RHOTHETA = RHO + THETA: IF RHOTHETA < SN THEN 8920
     E(0) = (HR - RD) / SIN(RHOTHETA)' y intercept of HOLZ line
     E(1) = TAN(PI / 2 + RHOTHETA)'slope
     E(2) = E(1) * E(1) + 1: E(3) = 2 * E(0) * E(1): E(4) = E(0) * E(0) - (DR * DR)' a b & c of quadratic
     E(5) = E(3) * E(3) - 4 * E(2) * E(4)
     E(6) = (-E(3) + SQR(E(5))) / (2 * E(2))'x1
     E(7) = E(0) + E(1) * E(6)' y1
     E(8) = (-E(3) - SQR(E(5))) / (2 * E(2))'x2
     E(9) = E(0) + E(1) * E(8)' y2
     RETURN
8920 E(6) = HR - RD
     E(8) = E(6)
     E(7) = SQR(DR * DR - E(6) * E(6))
     E(9) = -E(7)
     RETURN
8970 h$ = STR$(h) ' Prepares plotting label
     HS$ = LEFT$(h$, 1): HI$ = RIGHT$(h$, LEN(h$) - 1)
     k$ = STR$(k)
     KS$ = LEFT$(k$, 1): KI$ = RIGHT$(k$, LEN(k$) - 1)
     i$ = STR$(i)
     IS$ = LEFT$(i$, 1): II$ = RIGHT$(i$, LEN(i$) - 1)
     IF st$ <> "H" THEN IS$ = "": II$ = ""
     l$ = STR$(l)
     LS$ = LEFT$(l$, 1): LI$ = RIGHT$(l$, LEN(l$) - 1)
     IF LEN(l$) = 3 THEN LS$ = LS$ + LS$
     IF LABEL$ = "TITLE" THEN GOSUB 10920 ELSE GOSUB 10810 ' Plotter instruction
     RETURN
9090 IF RHO < RHO1 THEN RHO2 = -RHO2 ' provides sign for angle obtained from arccosine
     IF RHO - RHO1 - RHO2 > .01 THEN RHO = 2 * PI - RHO
     IF RHO - RHO1 - RHO2 < -.01 THEN RHO = 2 * PI - RHO
     RETURN
9130 NN = 0 ' reduces uvw's to smallest integer values
     FOR n = 5 TO 2 STEP -1
     IF U / n <> INT(U / n) THEN 9190
     IF V / n <> INT(V / n) THEN 9190
     IF W / n <> INT(W / n) THEN 9190
     NN = n
9190 NEXT
     RETURN
9210 NN = 0 ' eliminates the plotting of reducable hkl's
     FOR n = 5 TO 2 STEP -1
     IF h / n <> INT(h / n) THEN 9270
     IF k / n <> INT(k / n) THEN 9270
     IF l / n <> INT(l / n) THEN 9270
     NN = n
9270 NEXT
     RETURN
9290 NPM = h MOD 2 ' Change pen color for even and odd h values
     IF NPM = NPN THEN GOTO 9330
     NPN = NPM ' Pen number
     GOSUB 10750
9330 RETURN
9340 IF LZ <> CINT(LZ) THEN COLOR 13: LOCATE 1, 1: PRINT CHR$(7); "An inappropriate value of p is being used": PRINT "Color coding would be incorrect": PRINT : PRINT : GOTO 360
     NPM = LZ
9360 IF NPM = 2 THEN NPM = 0 ' Change pen color for HOLZ
     NPM = NPM + 1: IF NPM = 2 THEN NPM = 0
     IF NPM = NPN THEN GOTO 9410
     NPN = NPM ' Pen number: 0 - red, 1 - black
     GOSUB 10750
9410 RETURN
9420 h = U: k = V: l = 2 * CASQO3 * W: i = -(h + k) ' indices of plane normal to uvw - hexagonal
     RETURN
9440 h = U: k = V: l = W * (CA * CA)' indices of plane normal to uvw - tetragonal
     RETURN
9460 h = U * ASQ + W * LPC * LPA * COSBETA: k = V * BSQ: l = U * LPC * LPA * COSBETA + W * CSQ ' monoclinic
     RETURN
9570 IF U = CINT(U) AND V = CINT(V) AND W = CINT(W) THEN 9620
     UA = ABS(U): VA = ABS(V): WA = ABS(W)
     IF UA > VA AND UA > WA THEN U = CINT(9 * U / UA): V = CINT(9 * V / UA): W = CINT(9 * W / UA)
     IF VA > UA AND VA > WA THEN U = CINT(9 * U / VA): V = CINT(9 * V / VA): W = CINT(9 * W / VA)
     IF WA > UA AND WA > VA THEN U = CINT(9 * U / WA): V = CINT(9 * V / WA): W = CINT(9 * W / WA)
9620 GOSUB 9130 ' Reduce to smallest integers
     IF NN > 0 THEN U = U / NN: V = V / NN: W = W / NN: GOTO 9620
     T = -(U + V)
     RETURN
9660 h = U: k = V: l = W ' indices of plane normal to uvw - cubic
     RETURN
9680 GOSUB 9931
     IF plot% = 0 THEN RETURN
     IF DPSP$ = "D" OR DPSP$ = "H" THEN RETURN
     GOSUB 9210 ' Check for divisibility
     IF NN = 0 THEN RETURN
     h = h / NN: k = k / NN: l = l / NN
     GOSUB 9931
     h = h * NN: k = k * NN: l = l * NN
     IF plot% > 0 THEN plot% = 0: RETURN
     GOSUB 9931
     RETURN
9931 IF SF% < 5 THEN GOSUB 9940 ' Only plot allowed reflections
     IF SF% = 5 OR SF% = 6 THEN GOSUB 9990
     IF SF% = 7 THEN GOSUB 10010 ' cph
     IF SF% = 11 THEN GOSUB 10030 ' Alumina
     IF SF% = 12 OR SF% = 176 THEN GOSUB 10090 ' SG 159,176,186
     IF SF% = 36 THEN GOSUB 10121 ' Silicon oxynitride
     IF SF% = 136 THEN GOSUB 10130 ' Rutile
     IF SF% = 141 THEN GOSUB 10181 ' Anatase
     RETURN
9940 MH = ABS(h MOD 2): MK = ABS(k MOD 2): ML = ABS(l MOD 2) ' fcc structure factor
     IF MH = MK AND MH = ML THEN plot% = 1 ELSE plot% = 0
     IF SF% = 4 THEN IF h MOD 2 = 0 AND k MOD 2 = 0 AND l MOD 2 = 0 AND (h + k + l) MOD 4 <> 0 THEN plot% = 0
     IF SF% = 3 THEN IF MH = 1 AND MK = 1 AND ML = 0 THEN plot% = 1 ' tetragonal zirconia structure factor
     RETURN
9990 IF (h + k + l) MOD 2 = 0 THEN plot% = 1 ELSE plot% = 0 ' bcc and bct
      RETURN
10010 IF (h + 2 * k) MOD 3 = 0 AND l MOD 2 <> 0 THEN plot% = 0 ELSE plot% = 1' cph
      RETURN
10030 IF (-h + k + l) MOD 3 = 0 THEN plot% = 1 ELSE plot% = 0' Alumina structure factor
      i = -(h + k)
      IF h = 0 AND l MOD 2 <> 0 THEN plot% = 0
      IF k = 0 AND l MOD 2 <> 0 THEN plot% = 0
      IF i = 0 AND l MOD 2 <> 0 THEN plot% = 0
      RETURN
10090 plot% = 1: IF h = 0 AND k = 0 AND l MOD 2 <> 0 THEN plot% = 0' SG 159,176,186
      i = -(h + k)
      IF SF% <> 176 THEN IF h = k OR h = i OR k = i THEN IF l MOD 2 <> 0 THEN plot% = 0
      RETURN
10121 IF (h + k) MOD 2 <> 0 THEN plot% = 0 ELSE plot% = 1
      IF k = 0 AND l MOD 2 <> 0 THEN plot% = 0
      RETURN
10130 plot% = 1: IF h = 0 THEN IF (k + l) MOD 2 <> 0 THEN plot% = 0' SG 136
      IF k = 0 THEN IF (h + l) MOD 2 <> 0 THEN plot% = 0
      IF h = 0 AND k = 0 THEN IF l MOD 2 <> 0 THEN plot% = 0 ELSE plot% = 1
      IF h = 0 AND l = 0 THEN IF k MOD 2 <> 0 THEN plot% = 0 ELSE plot% = 1
      IF k = 0 AND l = 0 THEN IF h MOD 2 <> 0 THEN plot% = 0 ELSE plot% = 1
      RETURN
10181 plot% = 1: IF (h + k + l) MOD 2 <> 0 THEN plot% = 0'SG 141
      IF l = 0 THEN IF h MOD 2 <> 0 THEN plot% = 0
      IF l = 0 THEN IF k MOD 2 <> 0 THEN plot% = 0
      IF h = k THEN IF (h + k + l) MOD 4 <> 0 THEN plot% = 0
      RETURN
      REM Plotter instruction subroutines
10200 PRINT #1, "IN;" ' Initialize HP Plotter
      COLOR 12
      PRINT #1, "OI;"
      WHILE EOF(1) = -1: WEND
      INPUT #1, id$
      IF id$ <> "7470A" AND id$ <> "7475A" THEN PRINT : PRINT "You are not using an HP plotter": PRINT : INPUT "Press <RTN> to continue ", B$: CLOSE #1: GOTO 10
      device$ = id$ + " Plotter"
      PRINT #1, "oe;"
      WHILE EOF(1) = -1: WEND
      INPUT #1, a$
      IF a$ <> "0" THEN PRINT : PRINT "HP plotter error": IF a$ = "8" THEN PRINT : PRINT "Lower pinch wheel"
     IF a$ <> "0" THEN PRINT : INPUT "Press <RTN> to continue ", B$: CLOSE #1: GOTO 10
       PRINT #1, CHR$(27) + ".E"
      WHILE EOF(1) = -1: WEND
      INPUT #1, a$
     IF a$ <> "0" THEN PRINT : PRINT "Communication error": PRINT : INPUT "Press <RTN> to continue ", B$: GOTO 10
      IF id$ = "7475A" THEN PRINT #1, "SC -1380,1380,-995,995;" ELSE PRINT #1, "SC -1400,1400,-1000,1000;"
      COLOR 3
      RETURN
10380 IF SPB$ = "S" THEN 10500 ELSE PRINT #1, "SP 1;PA 0,0;" ' Draw circumference
      IF DPSP$ = "D" THEN PRINT #1, "CI 10,30; CI 13,30;"
      IF DPSP$ = "H" THEN PRINT #1, "CI "; DR; ";": GOTO 10490
      PRINT #1, "VS15;" ' Pen velocity
      IF DPSP$ = "D" THEN PRINT #1, "VS; SR 0.4,0.8;": IF HOLZR$ = "" THEN RETURN ELSE PRINT #1, "CI "; CINT(LAMDAL * SQR(G1SQ)); ";CI "; CINT(LAMDAL * SQR(G2SQ)); ";": GOTO 10490
      PRINT #1, "CI 1000;"
      PRINT #1, "VS;" ' Pen velocity
      PRINT #1, "pa -50,0;pd;pa 50,0;pu;pa 0,-50;pd;pa 0,50;pu;"
      PRINT #1, "pa -1000,0;pd;pa -950,0;pu;pa 950, 0;pd;pa 1000, 0;pu;"
      PRINT #1, "pa 0,-1000;pd;pa 0, -950;pu;pa 0, 950;pd;pa 0, 1000;pu;"
      PRINT #1, "SR 0.5,1;" ' Set type size
10490 IF SPB$ = "P" THEN RETURN
10500 SCREEN disp1, 0: WINDOW (-1450, -1070)-(1450, 1140)
     IF display$ = "EGA" THEN COLOR 14, 1
   CIRCLE (0, 0), DR
      IF DPSP$ = "H" THEN RETURN
      IF DPSP$ = "P" THEN CIRCLE (0, 0), 1000: LINE (-1000, 0)-(-970, 0): LINE (970, 0)-(1000, 0): LINE (0, -1000)-(0, -950): LINE (0, 950)-(0, 1000): RETURN
      IF HOLZR$ = "" THEN RETURN
      CIRCLE (0, 0), LAMDAL * SQR(G1SQ): CIRCLE (0, 0), LAMDAL * SQR(G2SQ)
      RETURN
10570 IF SPB$ = "S" THEN 10600
      PRINT #1, "pa "; CINT(E(6)); ","; CINT(E(7)); "; pd; pa "; CINT(E(8)); ","; CINT(E(9)); "; pu;"
      IF SPB$ = "P" THEN RETURN
10600 LINE (E(6), E(7))-(E(8), E(9))
      RETURN
10620 PRINT #1, "LB  c/a = "; B$ + C$
      RETURN
10640 PRINT #1, "PA -1380,-600;LBp ="; NP; C$' label for HOLZ symbols
      PRINT #1, "PA -1380,-700;CI 10,30;LB ZOLZ"; C$; "PA -1380,-900;PD;PU;CI 5,45;CI 7,30;CI 10,30; LB SOLZ"; C$; "PA -1380,-800; SP2; CI 8,30; LB FOLZ"; C$
      RETURN
10670 PRINT #1, "PA -1380,-600;LBh = "; D(3); C$' label for rotation axis
      PRINT #1, "PA -1380,-700;LBk = "; D(4); C$
      IF st$ = "H" THEN PRINT #1, "PA -1380,-800;LBi = "; D(5); C$; "PA -1380,-900;LBl = "; D(6); C$
      IF st$ = "C" OR st$ = "T" THEN PRINT #1, "PA -1380,-800;LBl = "; D(6); C$
      CLS
      RETURN
10730 PRINT #1, "PA -1400,-1000;SP 0;" ' End of plotting
      RETURN
10750 IF SPB$ = "S" THEN RETURN ELSE IF NPN = 0 THEN PRINT #1, "SP2;" ELSE PRINT #1, "SP1;"
      RETURN
10770 PRINT #1, "PA"; RH; ","; RV; ";pd;pu;ci 20,15;"
      RETURN
10790 PRINT #1, "pr 20,20;LB"; ALPHA$ + C$; "pr 10,30;ci 7,15;"
      RETURN
10810 IF SPB$ = "S" THEN 10970
      IF LZ = 2 OR AAP = 2 THEN PRINT #1, "PA"; RH; ","; RV; ";PD;PU;CI 4,45;CI 7,30"
      IF lab = 0 THEN
	 IF LZ = 1 THEN
	    PRINT #1, "PA"; RH; ","; RV; ";CI 8,30;PR 20,30;"' LAB value determines size of plotted circles
	 ELSE
	    PRINT #1, "PA"; RH; ","; RV; ";CI 10,30;PR 20,30;"
	 END IF
      END IF
      IF DPSP$ = "T" AND st$ = "C" THEN PRINT #1, "PA"; RH; ","; RV; ";CI 2,45;PR 5,7;": LF$ = "5": GOTO 10880
      IF DPSP$ = "D" THEN IF h > 9 OR h < -9 OR k > 9 OR k < -9 OR i > 9 OR i < -9 OR l > 9 OR l < -9 THEN 10960
      IF lab > 0 THEN PRINT #1, "PA"; RH; ","; RV; "; CI 7,30;"
10880 IF lab = 2 THEN 10960
      IF lab = 3 THEN INPUT "Enter rotation angle "; ANG: PH = 5: IF ANG < 0 THEN ANG = -ANG: PH = -60' label with rotation angles for axis/angle pairs
      IF lab = 3 THEN ANG$ = STR$(ANG): PRINT #1, "pr "; PH; ",0;lb"; ANG$ + C$: RETURN
      IF lab > 0 THEN PRINT #1, "PR 20,30;"
10920 PRINT #1, "LB" + HS$ + KS$ + IS$ + LS$ + CR$ + C$
      PRINT #1, "PR 0,-"; LF$; ";"
      PRINT #1, "LB" + HI$ + KI$ + II$ + LI$ + C$
      IF LEN(l$) = 3 THEN PRINT #1, "PR 0,-7; LB" + CHR$(8) + CHR$(8) + CHR$(95) + CHR$(95) + C$
10960 IF SPB$ = "P" THEN RETURN
10970 IF ABS(RH) > 1301 OR ABS(RV) > 1001 THEN RETURN
   CIRCLE (RH, RV), DR, LZ + 11
      IF st$ <> "H" THEN i = 0
      IF h > 9 OR h < -9 OR k > 9 OR k < -9 OR i > 9 OR i < -9 OR l > 9 OR l < -9 THEN RETURN
      NH = 20
      ON h + 10 GOSUB 11081, 11082, 11083, 11084, 11085, 11090, 11100, 11110, 11120, 11130, 11140, 11150, 11160, 11170, 11171, 11172, 11173, 11174, 11175
      NH = 60
      ON k + 10 GOSUB 11081, 11082, 11083, 11084, 11085, 11090, 11100, 11110, 11120, 11130, 11140, 11150, 11160, 11170, 11171, 11172, 11173, 11174, 11175
      NH = 100
      IF st$ = "H" THEN ON i + 10 GOSUB 11081, 11082, 11083, 11084, 11085, 11090, 11100, 11110, 11120, 11130, 11140, 11150, 11160, 11170, 11171, 11172, 11173, 11174, 11175
      IF st$ = "H" THEN NH = 140
      ON l + 10 GOSUB 11081, 11082, 11083, 11084, 11085, 11090, 11100, 11110, 11120, 11130, 11140, 11150, 11160, 11170, 11171, 11172, 11173, 11174, 11175
      RETURN
11081 PUT (RH + NH, RV + 20), NB9, OR: RETURN
11082 PUT (RH + NH, RV + 20), NB8, OR: RETURN
11083 PUT (RH + NH, RV + 20), NB7, OR: RETURN
11084 PUT (RH + NH, RV + 20), NB6, OR: RETURN
11085 PUT (RH + NH, RV + 20), NB5, OR: RETURN
11090 PUT (RH + NH, RV + 20), NB4, OR: RETURN
11100 PUT (RH + NH, RV + 20), NB3, OR: RETURN
11110 PUT (RH + NH, RV + 20), NB2, OR: RETURN
11120 PUT (RH + NH, RV + 20), NB1, OR: RETURN
11130 PUT (RH + NH, RV + 20), N0, OR: RETURN
11140 PUT (RH + NH, RV + 20), N1, OR: RETURN
11150 PUT (RH + NH, RV + 20), N2, OR: RETURN
11160 PUT (RH + NH, RV + 20), N3, OR: RETURN
11170 PUT (RH + NH, RV + 20), N4, OR: RETURN
11171 PUT (RH + NH, RV + 20), N5, OR: RETURN
11172 PUT (RH + NH, RV + 20), N6, OR: RETURN
11173 PUT (RH + NH, RV + 20), N7, OR: RETURN
11174 PUT (RH + NH, RV + 20), N8, OR: RETURN
11175 PUT (RH + NH, RV + 20), N9, OR: RETURN
      RETURN
11190 PRINT #1, "IN; SC -600,1365,-250,1180; SP 1; PA 0,0: PD; PA 1000,0; AA 0,0,60; PA 0,0; PU;" ' Alumina unit triangle
      PRINT #1, "CI 7,30; PA -90,-30;SR 0.5,1; LB0001" + C$ + "PA 1000,0; CI 7,30; PA 1020,-10;LB  -" + C$ + "PA 1020,-30;LB0110" + C$
      PRINT #1, "PA 866,500; CI 7,30; PA 880,540; LB- -" + C$ + "PA 880,520; LB1210" + C$
      PRINT #1, "PA 500,866; CI 7,30; PA 520,900; LB-" + C$ + "PA 520,880; LB1100" + C$
      RETURN
11240 PRINT #1, "SR; PA 900,1000; LBUnit triangle" + C$
      IF SF% = 11 THEN PRINT #1, "PA 900,900; LBAlumina" + C$
      PRINT #1, "PA -600,-250; SP 0;"
      RETURN
11280 PRINT #1, "IN; SC -100,600,-400,100; SP 1; PA 0,0: PD; PA 414,0; AR -1414,0,-15,1; PA 0,0; PU;" ' Cubic unit triangle
      PRINT #1, "PA -60,10;SR 1.5,3; LB001" + C$ + "PA 430,10; LB011" + C$
      PRINT #1, "PA 385,-395; LB111" + C$
      RETURN
11320 PRINT #1, "SR; PA 450,80; LBUnit triangle" + C$
      PRINT #1, "SR; PA 450,60; LBCubic" + C$
      PRINT #1, "PA -100,-400; SP 0;"
      RETURN