Title                   :TCBED.ABS
Keywords                :CBED
Computer                :VAX, Macintosh, Machine independent
Operating System        :Any
Programming Language    :Fortran 77
Hardware Requirements   :Line printer
Author(s)               :J.M. Zuo and J.C.H. Spence
Correspondence Address  :Physics Dept, Arizona State Univ., Tempe, AZ 85287, USA
Abstract:
         A Transmission Convergent Beam Electron diffraction program which
simulates the two-dimensional distribution of intensity in CBED disks. The
program includes three-dimensional (HOLZ) diffraction for non-centrosymmetric
crystals with absorption and inclined boundaries. It uses Bloch wave (matrix
diagonalization for general complex matrix) method. A complete description,
source code and test data are given in J. Elec. Micros. Tech., May, 1989
(Zuo et al). An application to GaAs can be found in Phys. Rev. Lett., Vol 61,
P353 (1988). The program produces a numerical output file and a two-digits
intensity map on line printer.
----------------------------------------------------------------------------- 
 
Title                   :TCBED.DOC
Keywords                :CBED
Computer                :VAX, Macintosh, Machine independent
Operating System        :Any
Programming Language    :Fortran 77
Hardware Requirements   :Line printer
Author(s)               :J.M. Zuo and J.C.H. Spence
Correspondence Address  :Physics Dept, Arizona State Univ., Tempe, AZ 85287, USA

A full description of this program including test data is contained in the may,
1989 issue of Journal of Electron Microscopy Technique (Zuo et al)

C
C WRITTEN BY J.M. ZUO, OCT 1987, SEE J.M. ZUO OR J.C.H. SPENCE
C (PHYSICS DEPT, ARIZONA STATE UNIVERSITY, TEMPE, AZ85287, USA,
C TEL (602)965-6486) TO DISCUSS USE OF THIS PROGRAM AND REPORT ANY ERRORS.
C
C THE PURPOSE OF THIS PROGRAM IS TO CALCULATE THE TWO-DIMENSIONAL INTENSITY
C DISTRIBUTION IN TRANSMISSION ELECTRON DIFFRACTION PATTERN FOR A RANGE OF
C INCIDENT BEAM DIRECTIONS. IT THUS CALCULATES CONVERGENT BEAM DISC
C INTENSITY DISTRIBUTIONS FOR ELECTRONS IN 20 Kev - 5 Mev RANGE.
C
C THREE DIMENSIONAL DIFFRACTION (HIGH ORDER LAUE ZONE EFFECTS) ARE
C INCLUDED BY THE RENORMALIZED EIGENVECTOR METHOD. THIS IS DESCRIBED
C IN FOLLOWING PAPERS:
C A.L. Lewis,R.E. Villagrana and A.J.F Metherall,Acta Cryst.
C A34,138 (1978) and J.M. Zuo, J.C.H. Spence, M. O'keffe, Phys. Rev. Letts,
C 61, p353 (1988). More background information can be found in
C Humphreys, Rep. Prog. Phys. 42, p 1825 (1979),
C Spence, "Experimental HREM" Oxford Univ. Press (2nd ed) 1988.
C Zuo and Spence, J. Micros. Techniqe (1988), in press.
C
C THE CRYSTAL IS ASSUMED TO BE PARALLEL SIDED SLAB, WHOSE INCLINATION
C TO THE BEAM MUST BE SPECIFIED. UP TO 100 BEAMS MAY BE INCLUDED. THE
C PROGRAM SHOULD WORK FOR ANY CRYSTAL SYSTEM, BUT HAS ONLY BEEN TESTED
C FOR CUBIC AND HEXAGONAL CRYSTALS.
C
C ABSORPTION IS INCLUDED USING THE GENERAL EIGENVALUES AND
C EIGENVECTORS OF A COMPLEX GENERAL MATRIX. (NOT BY PERTURBATION)
C
C THE SIGN CONVENTION USED IS CONSISTENT WITH AN INCIDENT WAVE
C exp(-2*3.14iK*r), K IS DOWNWARDS TO THE CRYSTAL (WITH ITS SURFACE
C NORMAL UPWARDS). SEE SPENCE (ABOVE) AND Ultramic. 12, p75 (1983)
C FOR DISCUSSION OF SIGN CONVENTION ETC.
C THIS PROGRAM USE THE CRYSTALLOGRAPHIC NOT QUANTUM MECHANICAL SIGN
C CONVENTION, AND IS THEREFORE CONSISTENT WITH THE STRUCTURE FACTORS IN
C THE INTERNATIONAL TABLES FOR X-RAY CRYSTALLOGRAPHY.
C
C THE PROGRAM USES THE EISPACK ROUTINES TO FIND THE EIGENVECTORS AND
C EIGENVALUES OF A COMPLEX GENERAL MATRIX. SEE B.T. SMITH et. al.
C Lecture Notes in Computer Science, Vol. 6. SPRINGER-VERLAG 1976
C SUBROUTINE "INVESE" IS ALSO REQUIRED.
TEST DATA FOR TCBED PROGRAM 
 
CdS, Bragg cond. at the center of (41-2) and (41-4) disc
2 4
12
0 0 0
4 1 -2
4 1 -4
4 1 0
4 1 2
4 1 -6
4 1 -8
0 0 -2
0 0 -4
0 0 2
0 0 -6
0 0 -8
-1 4 0 -1 4 0
4.1348 4.1348 6.7490 90.0 90.0 120.0
120.0 1 865.0
-1.78552 -0.44638 3.0 0.00 0.00 0.06 0.013894 0.0034734 0.0 50 50
Cd
19.2214 0.5946 17.6444 6.9089 4.461 24.7008 1.6029 87.4825 5.0694
S
6.9053 1.4679 5.2034 22.2151 1.4379 0.2536 1.5863 56.1720 0.8669
1
0 0 -2 0.056818 0.509081
1 0.66667 0.33333 0.00000 1.4
1 0.33333 0.66667 0.50000 1.4
2 0.66667 0.33333 0.37700 1.13
2 0.33333 0.66667 0.87700 1.13
0.57E-01 -0.045 0.13 0.03
1 3
1
0 0 -2 -0.08


 
 OUTPUT LISTING OF TCBED PROGRAM, VERSION 1
 
 CdS, Bragg cond. at the center of (41-2
 
 The unit cell of the crystal is (in Anstroms  and degrees):
 
          a =   4.13480
          b =   4.13480
          c =   6.74900
      alpha =   90.0000
       beta =   90.0000
      gamma =  120.0000
 
 Calculate for  1 thicknesses. Those thickness are:   865.00
 
 The closest zone axis is ( -1  4  0)
 The surface normal of the specimen is ( -1  4  0)
 
 The electron diffraction intensity is calculated for  50 50 different incident
beam directions.
 The initial incident direction is ( -1.78552 -0.44638  3.00000)
 This is the dimensionless component of the incident wavevector in the nearest Z
OLZ.
 In each calculation this is increased by ( 0.00000 0.00000 0.06000) and ( 0.013
89 0.00347 0.00000)
 
 The parameter for X-ray structure factor calculation are:
 
 Cd   48. 19.221399  0.594600 17.644400  6.908900  4.461000 24.700800  1.602900
87.482498  5.069400
 S    16.  6.905300  1.467900  5.203400 22.215099  1.437900  0.253600  1.586300
56.172001  0.866900
 
 Number of structure factors to be refined is  1,which are:
 (  0  0 -2)
 
 The atomic positions in the unit cell and Debye-Waller factors are:
 Cd      0.666670    0.333330    0.000000    1.400000
 Cd      0.333330    0.666670    0.500000    1.400000
 S       0.666670    0.333330    0.377000    1.130000
 S       0.333330    0.666670    0.877000    1.130000
 
 
 U(000)= 0.570E-01, U'(000)/U(000) = -0.045, B=  0.130, C=  0.030
 
 The Ug'/Ug for following reflections are:
      0    0   -2 -0.080000
 
 The volume of the unit cell is   99.92613
 a', b' and c' are :   0.27926   0.27926   0.14817
 
 The incident electron has energy 120.00 kev.
 
 The wavelength is 0.03349229 Angstroms
 
 The indicies and the Ug Values (including absorption)  used in this calculation
 are: 12
 
      Beams           Real         Imaginary      Amplitude        Phase
 (   4   1  -2)    0.68736E-02    0.19481E-02    0.71443E-02    15.823883
 (   4   1  -4)    0.30387E-02   -0.25976E-03    0.30498E-02    -4.885984
 (   4   1   0)    0.98622E-02   -0.11562E-02    0.99298E-02    -6.686487
 (   4   1   2)    0.62244E-02   -0.35068E-02    0.71443E-02   -29.396605
 (   4   1  -6)    0.35266E-02   -0.22105E-02    0.41621E-02   -32.078915
 (   4   1  -8)    0.38135E-02   -0.64182E-03    0.38672E-02    -9.553431
 (   0   0  -2)    0.51828E-01    0.23723E-01    0.57000E-01    24.594271
 (   0   0  -4)    0.14255E-01   -0.27152E-03    0.14258E-01    -1.091180
 (   0   0   2)    0.47398E-01   -0.31661E-01    0.57000E-01   -33.742115
 (   0   0  -6)    0.14306E-01   -0.78553E-02    0.16321E-01   -28.770739
 (   0   0  -8)    0.11884E-01   -0.16854E-02    0.12003E-01    -8.071494
1
 REFLECTION (    4    1   -4  )
 
 THICKNESS     865.00 Anstroms
 
 THE INTENSITY IS SCALED FROM 0 TO 99, MAXIMUM IS  0.127E+00
 
     0  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  8  9  0  1  2  3  4  5  6  7  8  9
 
 0   0  0  0  0  0  0  0  0  0  0  1  2  2  2  2  2  1  0  0  0  0  0  0  1  1
0  0  0  0  2  7 14 24 34 43 45 38 26 16  9
 1   0  0  0  0  0  0  0  0  0  0  1  2  2  1  1  1  0  0  0  0  1  1  2  2  2
1  1  2  4  9 17 27 35 41 43 39 28 15  6  2
 2   0  0  0  0  0  0  0  0  0  0  1  2  2  1  0  0  0  0  0  1  2  3  4  3  2
2  3  7 12 20 30 37 41 39 35 27 15  5  1  2
 3   0  0  0  0  0  0  1  1  0  0  1  2  2  0  0  1  1  0  0  1  3  4  4  2  2
3  8 14 23 32 40 42 38 30 21 13  5  0  2  6
 4   0  0  0  0  1  1  2  2  1  0  0  1  1  0  0  1  1  1  0  1  3  3  2  1  0
4 12 22 32 39 42 38 28 16  8  3  0  2  7 12
 5   0  0  0  0  1  2  3  4  2  0  0  0  0  0  0  1  1  1  0  1  2  3  2  0  0
5 14 26 35 39 36 26 14  5  1  1  3  8 14 16
 6   0  0  0  1  2  3  5  5  4  1  0  0  0  0  0  0  1  0  0  1  4  6  7  4  3
5 13 23 30 30 23 12  3  0  2  6 10 15 18 16
 7   0  0  0  1  2  3  5  6  5  2  0  0  0  0  0  0  0  0  0  4  9 15 17 14  8
6 10 16 19 16  9  2  0  4 10 14 17 20 19 13
 8   1  1  0  1  2  3  4  6  6  3  1  0  0  0  0  0  1  1  4 10 19 28 30 24 14
6  4  7  8  5  1  1  7 15 21 22 22 21 17  9
 9   1  1  1  1  1  1  3  5  6  4  2  0  0  0  2  3  5  7 11 19 30 40 41 32 17
5  1  1  2  3  5 11 19 26 28 26 22 18 12  5
 0   2  2  1  0  0  0  2  4  4  4  2  0  0  2  6 11 14 16 20 27 37 44 43 31 15
3  0  2  6 11 18 25 31 33 30 24 18 13  8  2
 1   2  2  1  0  0  0  1  2  3  3  2  1  2  6 13 20 25 26 27 29 35 38 34 22  8
1  3 10 19 28 34 38 37 33 26 19 13  8  4  1
 2   1  2  2  1  0  1  1  2  2  2  3  3  5 11 21 29 34 33 28 25 24 23 19 10  1
1  9 22 36 44 47 43 35 27 18 12  7  4  2  1
 3   0  1  2  2  2  3  3  3  3  3  4  6 10 17 26 34 37 32 23 14 10  9  8  4  0
4 16 33 48 54 50 39 27 18 11  6  3  2  1  1
 4   0  1  1  2  3  5  5  6  5  5  7 10 14 20 27 32 31 24 13  3  2  6 10 10  9 1
1 23 39 51 52 43 29 17  9  5  3  1  1  1  1
 5   0  0  1  2  4  6  7  7  7  8 10 13 16 20 22 23 19 12  3  0  6 19 29 30 25 2
2 27 37 44 41 29 17  8  4  2  1  0  1  2  1
 6   0  0  0  2  3  5  7  8  8  9 12 14 15 15 14 11  7  2  0  7 24 45 58 56 43 3
2 28 30 31 25 15  7  3  1  1  0  0  1  1  1
 7   0  0  0  1  3  4  5  6  7  9 11 12 12  9  6  3  1  0  6 23 49 74 85 77 57 3
7 25 20 17 11  5  2  1  1  0  0  0  1  1  1
 8   0  0  0  1  2  3  3  3  4  6  8  9  7  4  3  3  4  9 22 43 70 92 99 86 62 3
7 20 11  6  3  1  0  1  0  0  0  0  1  1  0
 9   0  0  0  1  2  2  2  1  1  3  4  5  5  5  7 11 18 26 39 57 78 93 93 79 55 3
1 14  4  1  0  1  1  1  0  0  0  0  1  1  0
 0   0  0  1  1  2  2  1  1  0  1  2  4  7 11 17 26 35 43 51 60 70 75 73 60 41 2
1  7  1  0  1  2  2  1  0  0  0  0  0  0  0
 1   1  1  1  2  2  2  2  2  1  2  4  8 13 20 30 40 49 53 54 52 50 49 46 38 25 1
2  3  0  1  3  3  2  0  0  0  0  0  0  0  0
 2   2  2  1  1  1  1  2  2  3  4  8 14 22 30 39 47 53 53 46 36 27 23 22 19 12
5  0  0  3  4  3  1  0  0  0  0  0  0  0  0
 3   2  2  1  0  0  0  1  2  3  7 13 21 30 37 42 46 47 42 32 19 10  7  8  7  5
1  0  1  3  4  2  0  0  0  0  0  0  0  0  0
 4   1  1  1  0  0  0  0  0  2  8 16 25 33 37 39 38 35 28 17  7  2  1  3  3  2
0  0  2  3  3  1  0  0  0  0  0  0  0  0  0
 5   1  1  2  2  2  1  1  0  1  6 14 23 29 32 31 28 23 16  8  2  0  1  3  4  3
2  2  2  2  2  1  0  0  0  0  0  0  0  0  0
 6   2  3  4  5  6  7  6  4  3  5 10 16 22 24 24 22 19 13  6  1  1  3  6  7  6
5  4  3  2  1  0  0  0  0  0  0  0  0  0  0
 7   7  6  7  9 12 15 16 15 11  6  5  8 13 17 21 23 22 18 11  4  1  3  7 11 11
8  5  3  2  1  0  0  1  1  0  0  0  0  0  0
 8  15 11 10 11 15 22 28 29 22 12  3  1  6 13 21 28 32 28 19  7  0  2  9 15 16 1
2  7  4  2  1  1  1  2  2  1  1  0  0  0  0
 9  24 16 11 10 15 26 38 43 36 21  6  0  3 12 24 35 41 39 26 10  0  2 10 18 19 1
5  9  5  2  1  1  2  3  3  2  1  0  0  0  0
 0  31 20 10  6 11 25 42 50 45 29 12  3  4 13 26 39 46 44 31 13  3  4 12 20 21 1
6  9  5  3  2  2  3  4  4  3  1  0  0  0  0
 1  32 20  8  2  5 19 37 49 46 33 17  7  6 13 24 36 44 42 31 16  6  7 13 18 18 1
4  9  5  4  3  3  4  4  4  2  1  0  0  0  0
 2  29 19  8  0  1 12 27 39 39 31 19 11  9 12 19 28 34 34 26 16  9  9 12 14 13 1
0  7  5  4  3  3  3  3  3  2  1  0  0  0  0
 3  25 19 11  4  2  7 16 24 26 23 17 11  9  9 12 17 21 22 19 13  9  8  8  8  7
5  4  3  3  2  2  2  2  2  1  1  0  0  0  0
 4  22 21 18 13  8  7  8 11 13 13 11  9  7  6  6  8 10 11 11  9  6  5  4  3  2
1  1  1  1  1  1  1  0  1  1  1  0  0  0  0
 5  25 27 27 23 16  9  5  3  4  4  5  4  4  3  2  2  2  3  4  4  3  1  0  0  0
0  0  0  0  0  0  0  0  0  0  1  1  0  0  0
 6  32 35 35 31 22 12  4  1  0  0  1  1  2  2  1  0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0  1  1  1  0  0
 7  41 41 39 32 22 12  5  2  0  0  0  0  0  1  1  1  0  0  0  0  0  0  1  1  2
1  1  1  1  1  1  1  1  1  1  1  1  1  0  0
 8  46 42 35 27 18 10  5  3  2  1  0  0  1  1  2  2  2  1  1  1  2  2  3  3  3
3  3  2  2  3  3  3  2  2  1  1  1  0  0  0
 9  44 35 27 18 11  5  3  2  2  2  1  1  1  1  2  2  2  3  3  4  4  4  4  4  3
3  3  4  4  4  4  3  3  2  1  1  0  0  0  0
 
 
------------------------------------------------------------------------------- 
 
Title                   :TCBED.SRC
Keywords                :CBED
Computer                :VAX, Macintosh, Machine independent
Operating System        :Any
Programming Language    :Fortran 77
Hardware Requirements   :Line printer
Author(s)               :J.M. Zuo and J.C.H. Spence
Correspondence Address  :Physics Dept, Arizona State Univ., Tempe, AZ 85287, USA
Code:
 
C**************************************************************************
C                                                                         *
C    THREE-DIMENSIONAL BLOCH-WAVE TRANSMISSION ELECTRON DIFFRACTION       *
C                             PROGRAM                                     *
C                         WITH ABSORPTION                                 *
C         FOR NON-CENTROSYMMETRIC (OR CENTROSYMMETRIC) CRYSTALS           *
C                                                                         *
C**************************************************************************
C
C WRITTEN BY J.M. ZUO, OCT 1987, SEE J.M. ZUO OR J.C.H. SPENCE
C (PHYSICS DEPT, ARIZONA STATE UNIVERSITY, TEMPE, AZ85287, USA,
C TEL (602)965-6486) TO DISCUSS USE OF THIS PROGRAM AND REPORT ANY ERRORS.
C
C THE PURPOSE OF THIS PROGRAM IS TO CALCULATE THE TWO-DIMENSIONAL INTENSITY
C DISTRIBUTION IN TRANSMISSION ELECTRON DIFFRACTION PATTERN FOR A RANGE OF
C INCIDENT BEAM DIRECTIONS. IT THUS CALCULATES CONVERGENT BEAM DISC
C INTENSITY DISTRIBUTIONS FOR ELECTRONS IN 20 Kev - 5 Mev RANGE.
C
C THREE DIMENSIONAL DIFFRACTION (HIGH ORDER LAUE ZONE EFFECTS) ARE
C INCLUDED BY THE RENORMALIZED EIGENVECTOR METHOD. THIS IS DESCRIBED
C IN FOLLOWING PAPERS:
C A.L. Lewis,R.E. Villagrana and A.J.F Metherall,Acta Cryst.
C A34,138 (1978) and J.M. Zuo, J.C.H. Spence, M. O'keffe, Phys. Rev. Letts,
C 61, p353 (1988). More background information can be found in
C Humphreys, Rep. Prog. Phys. 42, p 1825 (1979),
C Spence, "Experimental HREM" Oxford Univ. Press (2nd ed) 1988.
C Zuo and Spence, J. Micros. Techniqe (1988), in press.
C
C THE CRYSTAL IS ASSUMED TO BE PARALLEL SIDED SLAB, WHOSE INCLINATION
C TO THE BEAM MUST BE SPECIFIED. UP TO 100 BEAMS MAY BE INCLUDED. THE
C PROGRAM SHOULD WORK FOR ANY CRYSTAL SYSTEM, BUT HAS ONLY BEEN TESTED
C FOR CUBIC AND HEXAGONAL CRYSTALS.
C
C ABSORPTION IS INCLUDED USING THE GENERAL EIGENVALUES AND
C EIGENVECTORS OF A COMPLEX GENERAL MATRIX. (NOT BY PERTURBATION)
C
C THE SIGN CONVENTION USED IS CONSISTENT WITH AN INCIDENT WAVE
C exp(-2*3.14iK*r), K IS DOWNWARDS TO THE CRYSTAL (WITH ITS SURFACE
C NORMAL UPWARDS). SEE SPENCE (ABOVE) AND Ultramic. 12, p75 (1983)
C FOR DISCUSSION OF SIGN CONVENTION ETC.
C THIS PROGRAM USE THE CRYSTALLOGRAPHIC NOT QUANTUM MECHANICAL SIGN
C CONVENTION, AND IS THEREFORE CONSISTENT WITH THE STRUCTURE FACTORS IN
C THE INTERNATIONAL TABLES FOR X-RAY CRYSTALLOGRAPHY.
C
C THE PROGRAM USES THE EISPACK ROUTINES TO FIND THE EIGENVECTORS AND
C EIGENVALUES OF A COMPLEX GENERAL MATRIX. SEE B.T. SMITH et. al.
C Lecture Notes in Computer Science, Vol. 6. SPRINGER-VERLAG 1976
C SUBROUTINE "INVESE" IS ALSO REQUIRED.
C
      PROGRAM TCBED
C
C SOME VARIABLES USED IN THIS PROGRAM:
C
C AR,AI : THE MATRIX OF COEFFICIENTS OF BETHE EQUATION, THE DIMENSION
C IS SPECIFIED BY DIMENSION STATEMENTS.
C AX,BX,CX : THE PARAMETERS FOR X-RAY ATOMIC SCATTERING FACTOR (FROM
C INTERNATIONAL TABLES FOR X-RAY CRYSTALLOGRAPHY, Vol 4, p99 (Kynoch
C Press, Birmingham 1974), TABLE 2.2B. GUASSIAN FIT)
C JX,KX,LX: MILLER INDICIES (H,K,L) OF THE REFLECTIONS TO BE INCLUDED
C VR,VI: THE EIGENVALUES
C UGHR,UGHI:(ANGSTROM-2) THE Ug MATRIX, WHICH IS THE SAME FOR ALL TILTS.
C TH: THE THICKNESS OF THE CRYSTAL SLAB IN ANGSTROMS
C RJ,RK,RL: ATOMIC POSITION IN A UNIT CELL (DIMENSIONLESS)
C LAB1: ATOMIC SYMBOL
C ITYPE: SERIAL ORDER OF ATOMIC SPECIES.
C ZT: ATOMIC NUMBER
C DM: DEBYE-WALLER FACTOR FOR EACH ATOM IN THE UNIT CELL.
C Q: THE ABSORPTION PARAMETERS, Ug'/Ug. IF THESE ARE NOT ENTERED,
C THE DEFAULT VALUE OF Q, Q=Ug'/Ug = B*|G|-C*|G|*|G|, WILL BE USED. B,C
C MUST BE ENTERED (See Voss et. al. Z. Naturforsch. A35, p973, 1980)
C HQ,KQ,LQ: (H,K,L) OF REFLECTIONS WHOSE Q IS TO BE ENTERED
C CERTAIN SPECIFIED STRUCTURE FACTORS CAN BE ENTERED TO REPLACE THE VALUES
C TAKEN FROM THE INTERNATIONAL TABLES IF REQUIRED
C KCR,JCR,LCR: THE REFLECTIONS WHOSE STRUCTURE FACTOR IS TO BE REFINED
C SFCRE,SFCIM: THE REFINED STRUCTURE FACTOR Ug (ANGSTROM-2) VALUES
C NOUT: NUMBER OF BEAMS TO BE OUTPUTED
C NTH: NUMBER OF DIFFERENT THICKNESS TO BE CALCULATED AND OUTPUTED
C NXTILT,NYTILT: NUMBER OF POINTS IN X AND Y DIRECTION OF CBED DISK.
C THJ0,THK0,THL0: A INCIDENT BEAM DIRECTION, DEFINED BY A VECTOR FROM
C THE CENTER OF THE LAUE CIRCLE (IN NEAREST ZONE AXIS) TO THE ORIGIN,
C IN (DIMENSIONLESS) UNITS OF RECIPROCAL UNIT CELL VECTORS. THIS INCIDENT
C BEAM DIRECTION MAY BE AT POSITION 1,2,3, OR 4 ON THE CBED DISK, DEPENDING
C ON THE SIGN OF NXTILT AND NYTILT:
C     |---------|   1: NXTILT>0, NYTILT>0
C     |         |   2: NXTILT<0, NYTILT>0
C     2    1    |   3: NXTILT>0, NYTILT<0
C     |         |   4: NXTILT<0, NYTILT<0
C     4----3----|
C XTJ,XTK,XTL: X-DIRECTION INCREMENT OF INCIDENT BEAM (UNITS AS FOR THJ0)
C YTJ,YTK,YTL: Y-DIRECTION INCREMENT OF INCIDENT BEAM (UNITS AS FOR THJ0)
C CCR,CCI EIGENVECTOR MATRIX OF AR,AI. CCR(g,i),CCI(g,i) g NTH BEAM, ITH
C BLOCH WAVE.
C CC0: THE FIRST COLUMN OF THE INVERSE MATRIX OF CCR,CCI
C
C THE OUTPUT:
C
C THIS PROGRAM PRODUCES TWO OUTPUTS, A LISTING, AND THE CBED INTENSITY DATA
C WHICH IS IN THE ARRAY "AIG". THE LISTING FILE PROVIDES MOST INPUT DATA
C AND A 40X40 (MAXIMUM) SCALED ARRAY SHOWING THE TWO-DIMENSIONAL CBED
C INTENSITY IN EACH OF THE DISKS REQUESTED, AS NUMBERS.
C THE SIZE OF AIG IS 65004, THIS CAN STORE 65000 INTENSITY POINTS. THE
C OTHER 4 ARE USED TO STORE NOUT,NTH,NXTILT, AND NYTILT.
C THE OUTPUT DATA FILE IS UNFORMATTED. TO READ THE DATA FILE, USE:
C       READ(UNIT) ITOT,(AIG(I),I=1,ITOT)
C       NYTILT=AIG(ITOT)
C       NXTILT=AIG(ITOT-1)
C       NTH=AIG(ITOT-2)
C       NOUT=AIG(ITOT-3)
C TO GET THE CBED INTENSITY FOR THE I1(TH) REFLECTION, I2(TH) THICKNESS, USE:
C       MXTILT=IABS(NXTILT)
C       MYTILT=IABS(NYTILT)
C       DO 10 IY=1,MYTILT
C       DO 10 IX=1,MXTILT
C       IPST=(IY-1)*MXTILT*NTH*NOUT+(IX-1)*NTH*NOUT+(I1-1)*NOUT+I2
C  10   AAIG(IX,IY)=AIG(IPST)
C THE SIZE OF NOUT,NTH,NXTILT,NYTILT ARE ADJUSTABLE, BUT
C NTH*NOUT*NXTILT*NYTILT MUST BE LESS THAN OR EQUAL TO 65000.
C FOR A LINE SCAN OUTPUT ACROSS A CBED DISK, LET NYTILT=1 (THIS GIVES THE
C INTENSITY WHICH WOULD BE RECORDED BY SCANNING THE CBED PATTERN OVER A POINT
C DETECTOR).
C
C THE INPUT DATA:
C FREE FORMAT IS USED FOR NUMERICAL DATA IN THE INPUT DATA FILE. A SPACE
C MUST SEPERATE EACH NUMBER INPUT.
C
      CHARACTER TEXT*80
      CHARACTER*4 LAB1(20)
      DIMENSION AR(100,100),AI(100,100),VR(100),VI(100),DM(100),
     * Q(50),KCR(50),JCR(50),LCR(50),SFCRE(50),SFCIM(50),
     * UGHR(100,100),UGHI(100,100),CCR(100,100),CCI(100,100),
     * RJ(100),RK(100),RL(100),AX(4,20),BX(4,20),ZT(20),CX(20),
     * ITYPE(100),JX(100),KX(100),LX(100),TH(50),IG(40,40),
     * IOUT(50),AIG(62504),CALFA(3),SALFA(3),ALFA(3),AAIG(100,100)
      REAL KV
      EQUIVALENCE (AR(1,1),AAIG(1,1))
      INTEGER HQ(50),KQ(50),LQ(50)
      COMPLEX C1,C2,C3,X,Y,CC0(100),CVT(100,100),ONE,ZERO
      DATA PI/3.141592654/
C
      NR=10
      NW=12
      ZERO=CMPLX(0.0,0.0)
      ONE=CMPLX(1.0,0.0)
C
C READ NUMBER OF DIFFERENT ATOMS AND TOTAL NUMBER OF ATOMS
C
      READ(NR,10) TEXT
10    FORMAT(80A)
      WRITE (NW,12)
12    FORMAT(' OUTPUT LISTING OF TCBED PROGRAM, VERSION 1',/)
      WRITE (NW,10) TEXT
      READ (NR,*) NTYP,NPOS
C
C READ THE INDICIES OF BEAMS TO BE INCLUDED IN THE DYNAMICAL DIFFRACTION
C CALCULATION.
C
      READ(NR,*) NBEAMS
      DO 20 IBEAM=1,NBEAMS
20    READ(NR,*) JX(IBEAM),KX(IBEAM),LX(IBEAM)
C
C IH,IK,IL, ZONE AXIS DIRECTION (IN REAL SPACE)
C JSH,JSK,ISL THE CRYSTAL SURFACE NORMAL DIRECTION (IN REAL SPACE).
C
      READ(NR,*) IH,IK,IL,JSH,JSK,JSL
C
C A,B,C,ALFA(I) I=1,3 THE REAL SPACE UNIT CELL PARAMETERS
C ALFA(1)=ALPHA, ALFA(2)=BETA, ALFA(3)=GAMMA
C
      READ (NR,*) A,B,C,(ALFA(I),I=1,3)
      WRITE (NW,45) A,B,C,(ALFA(I),I=1,3)
 45   FORMAT (/,1X,'The unit cell of the crystal is (in Anstroms
     * and degrees):',//,10X,'a =',F10.5,
     * /,10X,'b =',F10.5,/,10X,'c =',F10.5,/,6X,'alpha =',F10.4,
     * /,7X,'beta =',F10.4,/,6X,'gamma =',F10.4)
C
C KV ELECTRON ACCELERATING VOLTAGE IN KILOVOLTS
C
      READ (NR,*) KV,NTH,(TH(I),I=1,NTH)
      WRITE (NW,7) NTH,(TH(I),I=1,NTH)
 7    FORMAT (/,1X,'Calculate for',I3,' thicknesses. Those'
     * ' thickness are: ',5(10f8.2,/))
      WRITE (NW,48) IH,IK,IL,JSH,JSK,JSL
 48   FORMAT (/,1X,'The closest zone axis is (',3I3,')',/,1x,
     * 'The surface normal of the specimen is (',3I3,')')
      READ (NR,*) THJ0,THK0,THL0,XTJ,XTK,XTL,YTJ,YTK,
     * YTL,NXTILT,NYTILT
      WRITE (NW,52) NXTILT,NYTILT,THJ0,THK0,THL0,XTJ,
     * XTK,XTL,YTJ,YTK,YTL
 52   FORMAT (/,1X,'The electron diffraction intensity is calculated
     * for ',2I3,' different incident beam directions.',/,' The initial
     * incident direction is (',3F9.5,')',/,' This is the dimensionless
     * component of the incident wavevector in the nearest ZOLZ.',/,
     * ' In each calculation this is increased by (',3F8.5,
     * ') and (',3f8.5,')')
C
C READ IN X-RAY ATOMIC SCATTERING FACTOR FITTING PARAMETERS
C SEE INT. TABLE FOR X-RAY CRYST. KYNOCH PRESS. BIRMINGHAM, U.K
C Vol 4, p99 (1974), TABLE 2.2B, ZT ATOMIC NUMBER
C
      WRITE (NW,53)
53    FORMAT (/,' The parameter for X-ray structure factor calculation
     * are:',/)
      DO 60 I=1,NTYP
      READ (NR,55) LAB1(I)
      READ (NR,*) (AX(IB,I),BX(IB,I),IB=1,4),CX(I)
55    FORMAT(A4)
      ZT(I)=CX(I)
      DO 56 IB=1,4
56    ZT(I)=AX(IB,I)+ZT(I)
      WRITE (NW,54) LAB1(I),ZT(I),(AX(IB,I),BX(IB,I),IB=1,4),CX(I)
54    FORMAT(1X,A4,F4.0,9F10.6)
60    CONTINUE
C
C READ IN THE MODIFIED STRUCTURE FACTORS WHICH WILL BE REFINED.
C THE PHASE IS IN RADIANS
C
      READ (NR,*) NCOREC
      IF (NCOREC.EQ.0) GOTO 78
      DO 70 I=1,NCOREC
      READ (NR,*) JCR(I),KCR(I),LCR(I),SF,PHASE
      SFCRE(I)=SF*COS(PHASE)
      SFCIM(I)=SF*SIN(PHASE)
70    CONTINUE
      WRITE (NW,75) NCOREC,(JCR(I),KCR(I),LCR(I),I=1,NCOREC)
75    FORMAT (/,' Number of structure factors to be refined is',i3,
     * :,',which are:'/,5(' (',3I3,')',:))
C
C READ IN THE TYPES AND POSITIONS OF ATOMS AND DEBYE WALLER FACTORS
C IN A UNIT CELL
C
78    DO 80 I=1,NPOS
80    READ(NR,*) ITYPE(I),RJ(I),RK(I),RL(I),DM(I)
      write (nw,87) (lab1(itype(i)),rj(i),rk(i),rl(i),DM(I),i=1,npos)
87    format(/,' The atomic positions in the unit cell',
     * ' and Debye-Waller factors are:',
     * /,40(1x,a4,4f12.6,/))
C
C READ ABSORPTION COEFFICIENTS AND CONSTANTS FOR PARAMETRIC APPROXIMATION.
C UZERO= U(0,0,0), ABM= U'(0,0,0)/UZERO, BAP=B, CAP=C
C NOTE: ABM IS NEGATIVE, AND BAP AND CAP ARE POSITIVE.
C
      READ (NR,*) UZERO,ABM,BAP,CAP
C
C READ SERIAL NUMBERS OF CBED DISCS TO BE OUTPUT ON FILE AND LISTING.
C
      READ (NR,*) NOUT,(IOUT(I),I=1,NOUT)
      write (nw,96) uzero,abm,BAP,CAP
96    format(/,' U(000)=',e10.3
     * ,', U''(000)/U(000) = ',f6.3,', B= ',F6.3,', C= ',f6.3)
C
C READ IN ABSORPTION coefficients IF SPECIFIED. OTHERWISE USE PARAMETRIC FIT.
C NOTE: Q(I) SHOULD BE NEGATIVE REAL CONSTANTS IN CRYST. SIGN CONVENTION.
C
      READ (NR,*) NAB
      IF(NAB.EQ.0) GOTO 99
      READ (NR,*) (HQ(I),KQ(I),LQ(I),Q(I),I=1,NAB)
      write(nw,98)
98    format(/,' The Ug''/Ug for following reflections are:')
      write(nw,97) (hq(i),kq(i),lq(i),q(i),i=1,nab)
97    FORMAT(2X,3I5,F10.6)
99    DO 100 I=1,3
      ALFA(I)=ALFA(I)*PI/180.0
      SALFA(I)=SIN(ALFA(I))
      CALFA(I)=COS(ALFA(I))
100   CONTINUE
      FACTOR=1.0-CALFA(1)*CALFA(1)-CALFA(2)*CALFA(2)
     * -CALFA(3)*CALFA(3)+2*CALFA(1)*CALFA(2)*CALFA(3)
      FACTOR=SQRT(FACTOR)
      OMEGA=A*B*C*FACTOR
      ASTAR=B*C*SALFA(1)/OMEGA
      BSTAR=A*C*SALFA(2)/OMEGA
      CSTAR=A*B*SALFA(3)/OMEGA
      WRITE (NW,115) OMEGA,ASTAR,BSTAR,CSTAR
115   FORMAT (/,' The volume of the unit cell is ',F10.5,/,
     * ' a'', b'' and c'' are :',3F10.5)
C
C ALAM IS THE RELATIVISTIC ELECTRON WAVELENGTH
C
      ALAM=0.3878314/SQRT(KV*(1.0+0.97846707E-03*KV))
      AL=IH*A
      BL=IK*B
      CL=IL*C
      ZABS=AL*AL+BL*BL+CL*CL+2*AL*BL*CALFA(3)
     * +2*BL*CL*CALFA(1)+2*AL*CL*CALFA(2)
      ZABS=SQRT(ZABS)
      AL1=JSH*A
      BL1=JSK*B
      CL1=JSL*C
      SABS=AL1*AL1+BL1*BL1+CL1*CL1+2.0*AL1*BL1*CALFA(3)
     * +2.0*BL1*CL1*CALFA(1)+2.0*AL1*CL1*CALFA(2)
      SABS=SQRT(SABS)
C
C COSZ : THE COSINE OF ANGLE BETWEEN SURFACE NORMAL AND ZONE AXIS
C
      COSZ=AL*AL1+BL*BL1+CL*CL1+(AL*BL1+BL*AL1)*CALFA(3)
     * +(BL*CL1+CL*BL1)*CALFA(1)+(AL*CL1+CL*AL1)*CALFA(2)
      COSZ=COSZ/(SABS*ZABS)
      CA=BSTAR*CSTAR*(CALFA(2)*CALFA(3)-CALFA(1))/(SALFA(2)*SALFA(3))
      CB=ASTAR*CSTAR*(CALFA(3)*CALFA(1)-CALFA(2))/(SALFA(3)*SALFA(1))
      CC=ASTAR*BSTAR*(CALFA(1)*CALFA(2)-CALFA(3))/(SALFA(1)*SALFA(2))
C
C BIGK IS THE WAVE VECTOR INSIDE THE CRYSTAL
C
      BIGK=1/ALAM
      BIGK=BIGK*BIGK+UZERO
      BIGK=SQRT(BIGK)
      RATIO=(1+1.9569341E-03*KV)/(PI*OMEGA)
      WRITE (NW,46) KV
46    FORMAT (/,1X,'The incident electron has energy ',F6.2,' kev.')
      WRITE (NW,116) ALAM
116   FORMAT (/,' The wavelength is ',F10.8,' Angstroms')
C
C SET UP Ug-g' MATRIX (EXCEPT g=g')
C
      DO 200 IX=1,NBEAMS
      IX1=IX-1
      DO 190 IY=1,IX1
      JXX=JX(IX)-JX(IY)
      KXX=KX(IX)-KX(IY)
      LXX=LX(IX)-LX(IY)
      X1=JXX*ASTAR
      X2=KXX*BSTAR
      X3=LXX*CSTAR
C
C SSQURE IS (1/2d(h,k,l))**2 = (sin(theta Bragg)/lambda)**2
C
      SSQURE=0.25*(X1*X1+X2*X2+X3*X3+
     * 2*(JXX*KXX*CC+JXX*LXX*CB+KXX*LXX*CA))
      GL=2.0*SQRT(SSQURE)
      JUDGE=ABS(JXX)+ABS(KXX)+ABS(LXX)
      IF (JUDGE.EQ.0) GOTO 190
C
C GET ABSORPTION COEFFICIENTS USING EMPRICAL FORMULA OF Voss et al.
C
      ABG=-BAP*GL+CAP*GL*GL
      ABP=ABG
      DO 110 IA=1,NAB
      IF (JXX.NE.HQ(IA)) GOTO 105
      IF (KXX.NE.KQ(IA)) GOTO 105
      IF (LXX.NE.LQ(IA)) GOTO 105
      ABP=Q(IA)
      GOTO 117
105   IF (JXX.NE.-HQ(IA)) GOTO 110
      IF (KXX.NE.-KQ(IA)) GOTO 110
      IF (LXX.NE.-LQ(IA)) GOTO 110
      ABP=Q(IA)
      GOTO 117
110   CONTINUE
117   SFRE=0.0
      SFIM=0.0
      IF (NCOREC.EQ.0) GOTO 121
      DO 120 IR=1,NCOREC
      IDEN=1
      IF(JXX.EQ.JCR(IR).AND.KXX.EQ.KCR(IR)
     * .AND.LXX.EQ.LCR(IR)) GOTO 160
       JJ=-JXX
       KK=-KXX
       LL=-LXX
       IDEN=-1
      IF(JJ.EQ.JCR(IR).AND.KK.EQ.KCR(IR)
     * .AND.LL.EQ.LCR(IR)) GOTO 160
120   CONTINUE
121   CONTINUE
      DO 140 IP=1,NPOS
      PHASE=2.0*PI*(JXX*RJ(IP)+KXX*RK(IP)+LXX*RL(IP))
      ITP=ITYPE(IP)
      SAVE=CX(ITP)
C
C CALCULATE THE X-RAY ATOMIC SCATTERING FACTOR.
C SAVE IS X-RAY SCATTERING FACTOR IN "NUMBER OF ELECTRONS"
C SUM GAUSSIAN SERIES OF INT. TABLES VOL 4 P71 (1974), EQN (3).
C
      DO 130 IB=1,4
130   SAVE=SAVE+AX(IB,ITP)*EXP(-BX(IB,ITP)*SSQURE)
C
C GET ELECTRON ATOMIC SCATTERING FACTOR HERE USING MOTT-BETHE FORMULA
C
      SAVE=0.023933754*(ZT(ITP)-SAVE)/SSQURE
C
C ADD IN THE DEBYE-WALLER FACTORS HERE
C
      SAVE=SAVE*EXP(-DM(IP)*SSQURE)
C
C STRUCTURE FACTOR OF ELECTRONS Fg
C
      SFRE=SFRE+SAVE*COS(PHASE)
      SFIM=SFIM+SAVE*SIN(PHASE)
140   CONTINUE
C
C GET Ug HERE (ANGSTROMS-2)
C
      SFRE=SFRE*RATIO
      SFIM=SFIM*RATIO
C
C ADD IN THE ABSORPTION PART AND STORE IN THE MATRIX
C
      UGHR(IX,IY)=SFRE-ABP*SFIM
      UGHI(IX,IY)=SFIM+ABP*SFRE
      UGHR(IY,IX)=SFRE+ABP*SFIM
      UGHI(IY,IX)=-SFIM+ABP*SFRE
      GOTO 190
160   D1=SFCRE(IR)
      D2=SFCIM(IR)*IDEN
      UGHR(IX,IY)=D1-ABP*D2
      UGHI(IX,IY)=D2+ABP*D1
      UGHR(IY,IX)=D1+ABP*D2
      UGHI(IY,IX)=-D2+ABP*D1
190   CONTINUE
200   CONTINUE
      WRITE (NW,230) NBEAMS
230   FORMAT (/,' The indicies and the Ug Values (including absorption)
     *  used in this calculation are:',i3,//,6X,'Beams',11X,'Real',
     * 9X,'Imaginary',6X,'Amplitude',8x,'Phase')
      DO 236 I=1,NBEAMS
      AMP=SQRT(UGHR(I,1)**2.0+UGHI(I,1)**2.0)
      IF (AMP.LT.1.0E-5) GOTO 236
      PHASE=180.0*ASIN(UGHI(I,1)/AMP)/PI
      IF (UGHR(I,1).LT.0.0) THEN
        SAVE=180.0-ABS(PHASE)
        PHASE=SIGN(SAVE,PHASE)
      END IF
      WRITE (NW,235) JX(I),KX(I),LX(I),UGHR(I,1),UGHI(I,1),AMP,PHASE
235   FORMAT (' (',3I4,')',3E15.5,F13.6)
236   CONTINUE
C
C ABM IS THE MEAN ABSORPTION COEFFICIENTS
C
      ABM=ABM*UZERO
      NXCENT=1
      NYCENT=1
      IF (NXTILT.GT.0) NXCENT=NXTILT/2
      IF (NYTILT.GT.0) NYCENT=NYTILT/2
      MXTILT=IABS(NXTILT)
      MYTILT=IABS(NYTILT)
      IPST=0
C
C DO LOOP FOR EACH INCIDENT BEAM DIRECTION
C
      DO 410 JTILT=1,MYTILT
      DO 400 ITILT=1,MXTILT
      THJ=THJ0+(ITILT-NXCENT)*XTJ+(JTILT-NYCENT)*YTJ
      THK=THK0+(ITILT-NXCENT)*XTK+(JTILT-NYCENT)*YTK
      THL=THL0+(ITILT-NXCENT)*XTL+(JTILT-NYCENT)*YTL
      XJ=THJ*ASTAR
      XK=THK*BSTAR
      XL=THL*CSTAR
      XYK=XJ*XJ+XK*XK+XL*XL+2*(THJ*THK*CC+THJ*THL*CB+THK*THL*CA)
      BIGKZ=-SQRT(BIGK*BIGK-XYK)
      BIGKN=BIGKZ*COSZ+(XJ*JSH*A+XK*JSK*B+XL*JSL*C)/SABS
      DO 300 IX=1,NBEAMS
      GN=(JX(IX)*JSH+KX(IX)*JSK+LX(IX)*JSL)/SABS
      II=IX-1
      DO 240 IY=1,II
      GNSTAR=(JX(IY)*JSH+KX(IY)*JSK+LX(IY)*JSL)/SABS
      DIVIDE=SQRT(1+(GN/BIGKN))*SQRT(1+(GNSTAR/BIGKN))
      AR(IX,IY)=UGHR(IX,IY)/DIVIDE
      AI(IX,IY)=UGHI(IX,IY)/DIVIDE
240   CONTINUE
      II=IX+1
      DO 245 IY=II,NBEAMS
      GNSTAR=(JX(IY)*JSH+KX(IY)*JSK+LX(IY)*JSL)/SABS
      DIVIDE=SQRT(1.0+GN/BIGKN)*SQRT(1.0+GNSTAR/BIGKN)
      AR(IX,IY)=UGHR(IX,IY)/DIVIDE
      AI(IX,IY)=UGHI(IX,IY)/DIVIDE
245   CONTINUE
C
C FIND EXCITATION ERRORS
C
      XJ=2*THJ+JX(IX)
      XK=2*THK+KX(IX)
      XL=2*THL+LX(IX)
      SAVE=XJ*JX(IX)*ASTAR*ASTAR+XK*KX(IX)*BSTAR*BSTAR
     * +XL*LX(IX)*CSTAR*CSTAR+(XJ*KX(IX)+XK*JX(IX))*CC
     * +(XJ*LX(IX)+XL*JX(IX))*CB+(XK*LX(IX)+XL*KX(IX))*CA
      GZZ=(JX(IX)*IH+KX(IX)*IK+LX(IX)*IL)/ZABS
C
C                2       2
C  SAVE=2*Kz*Sg=K - (K+g) =-2K*g-g*g
C  K is downwards so that its z-component is negative.
C
      SAVE=-2.0*BIGKZ*GZZ-SAVE
      SAVE=SAVE/(1+(GN/BIGKN))
      AR(IX,IX)=SAVE
      AI(IX,IX)=ABM
300   CONTINUE
C
C EIGENVALUES AND EIGENVECTORS
C
      CALL EISPACK(100,NBEAMS,AR,AI,VR,VI,CCR,CCI,IERR)
      IF (IERR.EQ.0) GOTO 350
      WRITE (NW,355)
355   FORMAT ('  EISPACK FAILS')
      STOP 'FAIL'
350   CONTINUE
C
C CALCULATE THE FIRST COLUMN OF C-INVERSE
C
      DO 500 IJ=1,NBEAMS
      DO 500 JI=1,NBEAMS
      GN=(JX(IJ)*JSH+KX(IJ)*JSK+LX(IJ)*JSL)/SABS
      CCR(IJ,JI)=CCR(IJ,JI)/SQRT(1.0+GN/BIGKN)
      CCI(IJ,JI)=CCI(IJ,JI)/SQRT(1.0+GN/BIGKN)
500   CVT(IJ,JI)=CMPLX(CCR(IJ,JI),CCI(IJ,JI))
      CALL INVERSE(100,NBEAMS,CVT,CC0,IERR)
      IF (IERR.EQ.0) THEN
        WRITE (NW,357)
357     FORMAT ('  INVERSE FAILS')
        STOP 'FAIL'
      END IF
C
C CALCULATE THE INTENSITY
C
      DO 380 ITH=1,NTH
      DO 380 IRL=1,NOUT
      IX=IOUT(IRL)
      RPT=0.0
      UPT=0.0
      DO 370 II=1,NBEAMS
      C1=CMPLX(CCR(IX,II),CCI(IX,II))
      PHASE=VR(II)*TH(ITH)*PI/BIGKN
      CF=COS(PHASE)
      SF=SIN(PHASE)
      C2=CMPLX(CF,SF)
      C3=CC0(II)
      DECAY=PI*VI(II)*TH(ITH)/BIGKN
      DECAY=EXP(-DECAY)
      C1=C1*C2*C3*DECAY
      RPT=RPT+REAL(C1)
      UPT=UPT+AIMAG(C1)
370   CONTINUE
      IPST=IPST+1
      AIG(IPST)=RPT*RPT+UPT*UPT
380   CONTINUE
400   CONTINUE
410   CONTINUE
      ITOT=IPST+4
C
C OUTPUT
C
      AIG(ITOT-3)=NOUT
      AIG(ITOT-2)=NTH
      AIG(ITOT-1)=NXTILT
      AIG(ITOT)=NYTILT
640   WRITE (16) ITOT,(AIG(I),I=1,ITOT)
C
C MAKE 40X40 LISTING
C
      MXT=MXTILT
      MYT=MYTILT
      IF (MXTILT.GT.40) MXT=40
      IF (MYTILT.GT.40) MYT=40
      WRITE(NW,1365)
      DO 1600 I5=1,NOUT
      DO 1500 I1=1,NTH
      IX=IOUT(I5)
      WRITE(NW,1250) JX(IX),KX(IX),LX(IX)
1250  FORMAT(' REFLECTION (',3I3,' )',/)
      WRITE(NW,1300) TH(I1)
1300  FORMAT(' THICKNESS ',F6.2,' Anstroms',/)
      DO 1310 I3=1,MYT
      DO 1310 I2=1,MXT
      IPST=(I3-1)*MXTILT*NTH*NOUT+(I2-1)*NTH*NOUT+(I1-1)*NOUT+I5
1310  AAIG(I2,I3)=AIG(IPST)
      AMX=0.0
      DO 1320 I3=1,MYT
      DO 1320 I2=1,MXT
1320  IF(AAIG(I2,I3).GT.AMX) AMX=AAIG(I2,I3)
      WRITE(NW,1330) AMX
1330  FORMAT(' THE INTENSITY IS SCALED FROM 0 TO 99, MAXIMUM IS '
     * ,E10.3,/)
      DO 1340 I3=1,MYT
      DO 1340 I2=1,MXT
1340  IG(I2,I3)=IFIX(99*AAIG(I2,I3)/AMX)
      WRITE(NW,1345) ((I2,I2=0,9),I3=1,4)
1345  FORMAT(3X,40I3,/)
      DO 1380 I3=1,MYT
      I4=MYT-I3+1
      I6=I3-(I3-1)/10*10-1
1380  WRITE(NW,1360) I6,(IG(I2,I4),I2=1,MXT)
1360  FORMAT(I2,X,40I3)
C
C GO TO A NEW PAGE HERE
C
      IF (I5.EQ.NOUT.AND.I1.EQ.NTH) GOTO 1500
      WRITE (NW,1365)
1365  FORMAT('1')
1500  CONTINUE
1600  CONTINUE
      STOP
      END
c
c
c
        subroutine inverse(nm,n,a,b,ierr)
c
c   to calculate the first row of c-inverse by Gaussian elimination
c   and back substitution
c     by j.m. zuo oct. 1987
c
        complex a(nm,n),t,zero,b(n)
        integer ipvt(100),kp1,jpvt,ierr,nm1
        ierr=100
        zero=cmplx(0.0,0.0)
        if (n.eq.1) goto 35
        amx=0.0
c
c   pivoting for the first row
c
        do 5 j=1,n
          if (cabs(a(1,j)).gt.amx) then
                amx=cabs(a(1,j))
                jpvt=j
          else
          end if
5       continue
        do 10 i=1,n
          t=a(i,jpvt)
          a(i,jpvt)=a(i,n)
          a(i,n)=t
10      continue
c
c   interchange firt row and last row
c
        do 12 i=1,n
          t=a(1,i)
          a(1,i)=a(n,i)
          a(n,i)=t
12      continue
c
c   Gaussian elimination with partial pivoting
c
        nm1=n-1
        do 35 k=1,nm1
          kp1=k+1
          m=k
          if (kp1.eq.nm1.or.k.eq.nm1) goto 15
          do 15 i=kp1,nm1
                if(cabs(a(i,k)).gt.cabs(a(m,k))) m=i
15        continue
          ipvt(k)=m
          t=a(m,k)
          a(m,k)=a(k,k)
          a(k,k)=t
          if (cabs(t).le.1.0e-07) goto 35
          do 20 i=kp1,n
                a(i,k)=-a(i,k)/t
20        continue
          do 30 j=kp1,n
                t=a(m,j)
                a(m,j)=a(k,j)
                a(k,j)=t
                if (t.eq.zero) goto 30
            do 25 i=kp1,n
                a(i,j)=a(i,j)+a(i,k)*t
25          continue
30        continue
35      continue
        if (cabs(a(n,n)).le.1.0e-07) ierr=0
c
c   back substitution to get x
c
        b(n)=1.0/a(n,n)
        do 50 kb=2,n
          k=n-kb+1
          kp1=k+1
          b(k)=zero
          do 40 j=kp1,n
40          b(k)=b(k)-a(k,j)*b(j)
          b(k)=b(k)/a(k,k)
50      continue
        t=b(jpvt)
        b(jpvt)=b(n)
        b(n)=t
        return
        end
c
c
c
c
        subroutine  eispack(nm,n,ar,ai,wr,wi,vr,vi,ierr)
c
c   eigenvalues and eigenvectors of a complex general matrix
c   calling EISPACK subroutine (B. T. Smith, Lecture notes in computer
c   science, Vol 6, Springer-Verlag 1976)
c   implemented by J.M. Zuo sept. 1987
c
        integer nm,n,is1,is2,iv1(100)
        real ar(nm,n),ai(nm,n),wr(n),wi(n),vr(nm,n),vi(nm,n),fv1(100)
        real fv2(100),fv3(100)
c
        call cbal(nm,n,ar,ai,is1,is2,fv1)
        call corth(nm,n,is1,is2,ar,ai,fv2,fv3)
        call comqr2(nm,n,is1,is2,fv2,fv3,ar,ai,wr,wi,vr,vi,ierr)
        if (ierr.ne.0) return
        call cbabk2(nm,n,is1,is2,fv1,n,vr,vi)
c
c   Normalize the eigenvectors
c
        do 160 i=1,n
          sum=0.0
          do 100 j=1,n
                term=vr(j,i)**2+vi(j,i)**2
                sum=sum+term
100       continue
          sum=sqrt(sum)
          do 120 j=1,n
                vr(j,i)=vr(j,i)/sum
                vi(j,i)=vi(j,i)/sum
120       continue
160     continue
        return
        end
 
 
        subroutine cbal(nm,n,ar,ai,low,igh,scale)
c
        integer i,j,k,l,m,n,jj,nm,igh,low,iexc
        real ar(nm,n),ai(nm,n),scale(n)
        real c,f,g,r,s,b2,radix
        real abs
        logical noconv
        radix=2.0
        b2=radix*radix
        k=1
        l=n
        go to 100
20      scale(m)=j
        if (j.eq.m) goto 50
        do 30 i=1,l
          f=ar(i,j)
          ar(i,j)=ar(i,m)
          ar(i,m)=f
          f=ai(i,j)
          ai(i,j)=ai(i,m)
          ai(i,m)=f
30      continue
        do 40 i=k,n
          f=ar(j,i)
          ar(j,i)=ar(m,i)
          ar(m,i)=f
          f=ai(j,i)
          ai(j,i)=ai(m,i)
          ai(m,i)=f
40      continue
50      goto (80,130) iexc
80      if(l.eq.1) goto 280
        l=l-1
100     do 120 jj=1,l
          j=l+1-jj
          do 110 i=1,l
            if(i.eq.j) goto 110
            if(ar(j,i).ne.0.0.or.ai(j,i).ne.0.0) goto 120
110       continue
          m=l
          iexc=1
          go to 20
120     continue
        goto 140
130     k=k+1
140     do 170 j=k,l
          do 150 i=k,l
            if(i.eq.j) goto 150
            if(ar(i,j).ne.0.0.or.ai(i,j).ne.0.0) goto 170
150     continue
        m=k
        iexc=2
        goto 20
170     continue
        do 180 i=k,l
180       scale(i)=1.0
190     noconv=.false.
        do 270 i=k,l
          c=0.0
          r=0.0
          do 200 j=k,l
            if(j.eq.i) goto 200
            c=c+abs(ar(j,i))+abs(ai(j,i))
            r=r+abs(ar(i,j))+abs(ai(i,j))
200     continue
        g=r/radix
        f=1.0
        s=c+r
210     if(c.ge.g) goto 220
        f=f*radix
        c=c*b2
        goto 210
220     g=r*radix
230     if(c.lt.g) goto 240
        f=f/radix
        c=c/b2
        goto 230
240     if((c+r)/f.ge.0.95*s) goto 270
        g=1.0/f
        scale(i)=scale(i)*f
        noconv=.true.
        do 250 j=k,n
          ar(i,j)=ar(i,j)*g
          ai(i,j)=ai(i,j)*g
250     continue
        do 260 j=1,l
          ar(j,i)=ar(j,i)*f
          ai(j,i)=ai(j,i)*f
260     continue
270     continue
        if(noconv) goto 190
280     low=k
        igh=l
        return
        end
 
c
c    corth reduces a complex general matrix to upper hessenberg
c    form using unitary similarity transformations.
c
        subroutine corth(nm,n,low,igh,ar,ai,ortr,orti)
c
        integer i,j,m,n,ii,jj,la,mp,igh,nm,kp1,low
        real ar(nm,n),ai(nm,n),ortr(igh),orti(igh)
        real f,g,h,fi,fr,scale,sqrt,cabs,abs
        complex cmplx
c
        la=igh-1
        kp1=low+1
        if(la.lt.kp1) goto 200
        do 180 m=kp1,la
          h=0.0
          ortr(m)=0.0
          orti(m)=0.0
          scale=0.0
          do 90 i=m,igh
90              scale=scale+abs(ar(i,m-1))+abs(ai(i,m-1))
          if (scale.eq.0.0) goto 180
          mp=m+igh
          do 100 ii=m,igh
                i=mp-ii
                ortr(i)=ar(i,m-1)/scale
                orti(i)=ai(i,m-1)/scale
                h=h+ortr(i)*ortr(i)+orti(i)*orti(i)
100       continue
          g=sqrt(h)
          f=cabs(cmplx(ortr(m),orti(m)))
          if (f.eq.0.0) goto 103
          h=h+f*g
          g=g/f
          ortr(m)=(1.0+g)*ortr(m)
          orti(m)=(1.0+g)*orti(m)
          goto 105
103       ortr(m)=g
          ar(m,m-1)=scale
105       do 130 j=m,n
                fr=0.0
                fi=0.0
            do 110 ii=m,igh
                i=mp-ii
                fr=fr+ortr(i)*ar(i,j)+orti(i)*ai(i,j)
                fi=fi+ortr(i)*ai(i,j)-orti(i)*ar(i,j)
110         continue
                fr=fr/h
                fi=fi/h
            do 120 i=m,igh
                ar(i,j)=ar(i,j)-fr*ortr(i)+fi*orti(i)
                ai(i,j)=ai(i,j)-fr*orti(i)-fi*ortr(i)
120         continue
130       continue
          do 160 i=1,igh
                fr=0.0
                fi=0.0
            do 140 jj=m,igh
                j=mp-jj
                fr=fr+ortr(j)*ar(i,j)-orti(j)*ai(i,j)
                fi=fi+ortr(j)*ai(i,j)+orti(j)*ar(i,j)
140         continue
                fr=fr/h
                fi=fi/h
            do 150 j=m,igh
                ar(i,j)=ar(i,j)-fr*ortr(j)-fi*orti(j)
                ai(i,j)=ai(i,j)+fr*orti(j)-fi*ortr(j)
150         continue
160       continue
          ortr(m)=scale*ortr(m)
          orti(m)=scale*orti(m)
          ar(m,m-1)=-g*ar(m,m-1)
          ai(m,m-1)=-g*ai(m,m-1)
180     continue
200     return
        end
 
c
c   comqr2 use QR method to calculate eigenvalues and eigenvectors
c   from eispack.
c
      subroutine comqr2(nm,n,low,igh,ortr,orti,hr,hi,wr,wi,zr,zi,ierr)
        integer i,j,k,l,m,n,en,ii,jj,ll,nm,nn,igh,ip1,
     *    its,low,lp1,enm1,iend,ierr,min0
        real hr(nm,n),hi(nm,n),wr(n),wi(n),zr(nm,n),zi(nm,n),
     *    ortr(igh),orti(igh)
        real si,sr,ti,tr,xr,xi,yi,yr,zzi,zzr,norm,machep,abs
        complex z3,csqrt
c
c   machep, relative accuracy of single precision real number operation
c   on VAX 750
c
        machep=1.19e-07
        ierr=0
        do 100 i=1,n
          do 100 j=1,n
            zr(i,j)=0.0
            zi(i,j)=0.0
            if (i.eq.j) zr(i,j)=1.0
100     continue
        iend=igh-low-1
        if(iend) 180,150,105
105     do 140 ii=1,iend
          i=igh-ii
          if(ortr(i).eq.0.0.and.orti(i).eq.0.0) goto 140
          if(hr(i,i-1).eq.0.0.and.hi(i,i-1).eq.0.0) goto 140
          norm=hr(i,i-1)*ortr(i)+hi(i,i-1)*orti(i)
          ip1=i+1
          do 110 k=ip1,igh
                ortr(k)=hr(k,i-1)
                orti(k)=hi(k,i-1)
110       continue
          do 130 j=i,igh
                sr=0.0
                si=0.0
            do 115 k=i,igh
                sr=sr+ortr(k)*zr(k,j)+orti(k)*zi(k,j)
                si=si+ortr(k)*zi(k,j)-orti(k)*zr(k,j)
115         continue
                sr=sr/norm
                si=si/norm
            do 120 k=i,igh
                zr(k,j)=zr(k,j)+sr*ortr(k)-si*orti(k)
                zi(k,j)=zi(k,j)+sr*orti(k)+si*ortr(k)
120         continue
130       continue
140     continue
150     l=low+1
        do 170 i=l,igh
          ll=min0(i+1,igh)
          if(hi(i,i-1).eq.0.0) goto 170
          norm=cabs(cmplx(hr(i,i-1),hi(i,i-1)))
          yr=hr(i,i-1)/norm
          yi=hi(i,i-1)/norm
          hr(i,i-1)=norm
          hi(i,i-1)=0.0
          do 155 j=i,n
                si=yr*hi(i,j)-yi*hr(i,j)
                hr(i,j)=yr*hr(i,j)+yi*hi(i,j)
                hi(i,j)=si
155       continue
          do 160 j=1,ll
                si=yr*hi(j,i)+yi*hr(j,i)
                hr(j,i)=yr*hr(j,i)-yi*hi(j,i)
                hi(j,i)=si
160       continue
          do 165 j=low,igh
                si=yr*zi(j,i)+yi*zr(j,i)
                zr(j,i)=yr*zr(j,i)-yi*zi(j,i)
                zi(j,i)=si
165       continue
170     continue
180     do 200 i=1,n
          if (i.ge.low.and.i.le.igh) goto 200
          wr(i)=hr(i,i)
          wi(i)=hi(i,i)
200     continue
        en=igh
        tr=0.0
        ti=0.0
220     if (en.lt.low) goto 680
        its=0
        enm1=en-1
240     do 260 ll=low,en
          l=en+low-ll
          if (l.eq.low) goto 300
        if (abs(hr(l,l-1)).le.machep*(abs(hr(l-1,l-1))
     *  +abs(hi(l-1,l-1))+abs(hr(l,l))+abs(hi(l,l)))) goto 300
260     continue
300     if (l.eq.en) goto 660
        if (its.eq.30) goto 1000
        if (its.eq.10.or.its.eq.20) goto 320
        sr=hr(en,en)
        si=hi(en,en)
        xr=hr(enm1,en)*hr(en,enm1)
        xi=hi(enm1,en)*hr(en,enm1)
        if (xr.eq.0.0.and.xi.eq.0.0) goto 340
        yr=(hr(enm1,enm1)-sr)/2.0
        yi=(hi(enm1,enm1)-si)/2.0
        z3=csqrt(cmplx(yr**2-yi**2+xr,2.0*yr*yi+xi))
        zzr=real(z3)
        zzi=aimag(z3)
        if (yr*zzr+yi*zzi.ge.0.0) goto 310
        zzr=-zzr
        zzi=-zzi
310     z3=cmplx(xr,xi)/cmplx(yr+zzr,yi+zzi)
        sr=sr-real(z3)
        si=si-aimag(z3)
        goto 340
320     sr=abs(hr(en,enm1))+abs(hr(enm1,en-2))
        si=0.0
340     do 360 i=low,en
          hr(i,i)=hr(i,i)-sr
          hi(i,i)=hi(i,i)-si
360     continue
        tr=tr+sr
        ti=ti+si
        its=its+1
        lp1=l+1
        do 500 i=lp1,en
          sr=hr(i,i-1)
          hr(i,i-1)=0.0
          norm=sqrt(hr(i-1,i-1)*hr(i-1,i-1)+hi(i-1,i-1)
     *      *hi(i-1,i-1)+sr*sr)
          xr=hr(i-1,i-1)/norm
          wr(i-1)=xr
          xi=hi(i-1,i-1)/norm
          wi(i-1)=xi
          hr(i-1,i-1)=norm
          hi(i-1,i-1)=0.0
          hi(i,i-1)=sr/norm
c
          do 490 j=i,n
                yr=hr(i-1,j)
                yi=hi(i-1,j)
                zzr=hr(i,j)
                zzi=hi(i,j)
                hr(i-1,j)=xr*yr+xi*yi+hi(i,i-1)*zzr
                hi(i-1,j)=xr*yi-xi*yr+hi(i,i-1)*zzi
                hr(i,j)=xr*zzr-xi*zzi-hi(i,i-1)*yr
                hi(i,j)=xr*zzi+xi*zzr-hi(i,i-1)*yi
490       continue
500     continue
        si=hi(en,en)
        if(si.eq.0.0) goto 540
        norm=cabs(cmplx(hr(en,en),si))
        sr=hr(en,en)/norm
        si=si/norm
        hr(en,en)=norm
        hi(en,en)=0.0
        if(en.eq.n) goto 540
        ip1=en+1
        do 520 j=ip1,n
          yr=hr(en,j)
          yi=hi(en,j)
          hr(en,j)=sr*yr+si*yi
          hi(en,j)=sr*yi-si*yr
520     continue
540     do 600 j=lp1,en
          xr=wr(j-1)
          xi=wi(j-1)
          do 580 i=1,j
                yr=hr(i,j-1)
                yi=0.0
                zzr=hr(i,j)
                zzi=hi(i,j)
                if(i.eq.j) goto 560
                yi=hi(i,j-1)
                hi(i,j-1)=xr*yi+xi*yr+hi(j,j-1)*zzi
560             hr(i,j-1)=xr*yr-xi*yi+hi(j,j-1)*zzr
                hr(i,j)=xr*zzr+xi*zzi-hi(j,j-1)*yr
                hi(i,j)=xr*zzi-xi*zzr-hi(j,j-1)*yi
580       continue
          do 590 i=low,igh
                yr=zr(i,j-1)
                yi=zi(i,j-1)
                zzr=zr(i,j)
                zzi=zi(i,j)
                zr(i,j-1)=xr*yr-xi*yi+hi(j,j-1)*zzr
                zi(i,j-1)=xr*yi+xi*yr+hi(j,j-1)*zzi
                zr(i,j)=xr*zzr+xi*zzi-hi(j,j-1)*yr
                zi(i,j)=xr*zzi-xi*zzr-hi(j,j-1)*yi
590       continue
600     continue
        if(si.eq.0.0) goto 240
        do 630 i=1,en
          yr=hr(i,en)
          yi=hi(i,en)
          hr(i,en)=sr*yr-si*yi
          hi(i,en)=sr*yi+si*yr
630     continue
        do 640 i=low,igh
          yr=zr(i,en)
          yi=zi(i,en)
          zr(i,en)=sr*yr-si*yi
          zi(i,en)=sr*yi+si*yr
640     continue
        goto 240
660     hr(en,en)=hr(en,en)+tr
        wr(en)=hr(en,en)
        hi(en,en)=hi(en,en)+ti
        wi(en)=hi(en,en)
        en=enm1
        goto 220
680     norm=0.0
        do 720 i=1,n
          do 720 j=i,n
                norm=norm+abs(hr(i,j))+abs(hi(i,j))
720     continue
        if(n.eq.1.or.norm.eq.0.0) goto 1001
        do 800 nn=2,n
          en=n+2-nn
          xr=wr(en)
          xi=wi(en)
          enm1=en-1
          do 780 ii=1,enm1
                i=en-ii
                zzr=hr(i,en)
                zzi=hi(i,en)
                if(i.eq.enm1) goto 760
                ip1=i+1
            do 740 j=ip1,enm1
                zzr=zzr+hr(i,j)*hr(j,en)-hi(i,j)*hi(j,en)
                zzi=zzi+hr(i,j)*hi(j,en)+hi(i,j)*hr(j,en)
740         continue
760             yr=xr-wr(i)
                yi=xi-wi(i)
                if (yr.eq.0.0.and.yi.eq.0.0) yr=machep*norm
                z3=cmplx(zzr,zzi)/cmplx(yr,yi)
                hr(i,en)=real(z3)
                hi(i,en)=aimag(z3)
780       continue
800     continue
        enm1=n-1
        do 840 i=1,enm1
          if(i.ge.low.and.i.le.igh) goto 840
          ip1=i+1
          do 820 j=ip1,n
                zr(i,j)=hr(i,j)
                zi(i,j)=hi(i,j)
820       continue
840     continue
        do 880 jj=low,enm1
          j=n+low-jj
          m=min0(j-1,igh)
          do 880 i=low,igh
                zzr=zr(i,j)
                zzi=zi(i,j)
            do 860 k=low,m
                zzr=zzr+zr(i,k)*hr(k,j)-zi(i,k)*hi(k,j)
                zzi=zzi+zr(i,k)*hi(k,j)+zi(i,k)*hr(k,j)
860         continue
                zr(i,j)=zzr
                zi(i,j)=zzi
880     continue
        goto 1001
1000    ierr=en
1001    return
        end
 
 
        subroutine cbabk2(nm,n,low,igh,scale,m,zr,zi)
c
        real scale(n),zr(nm,m),zi(nm,m),s
        if (igh.eq.low) goto 120
        do 110 i=low,igh
          s=scale(i)
          do 100 j=1,m
                zr(i,j)=zr(i,j)*s
                zi(i,j)=zi(i,j)*s
100       continue
110     continue
120     do 140 ii=1,n
          i=ii
          if (i.ge.low.and.i.le.igh) goto 140
          if (i.lt.low) i=low-ii
          k=scale(i)
          if (k.eq.i) goto 140
          do 130 j=1,m
                s=zr(i,j)
                zr(i,j)=zr(k,j)
                zr(k,j)=s
                s=zi(i,j)
                zi(i,j)=zi(k,j)
                zi(k,j)=s
130       continue
140     continue
        return
        end