CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C                RT version 2 - (c) Sean Parkin 99.09.09               C
C                                                                      C
C======================================================================C
C                                                                      C
C  A version of the R-tensor program, RT, that works throughout in     C
C  cartesian coordinates.                                              C
C                                                                      C
C  Version for use with SHELXL-2013, may be called from a script.      C
C  or from a command line without further interaction, as follows:     C
C                                                                      C
C  'rt filename' (without extension)                                   C
C                                                                      C
C======================================================================C
C                                                                      C
C  RT is a program to calculate a 3x3 tensor that is used to expand    C
C  commonly used scalar global statistics ('isotropic' or 'spherical'  C
C  approximation) to a form that incorporates 3D spatial information   C
C  ('anisotropic' or 'ellipsoidal' approximation). Input consists of   C
C  a SHELX97 *.res file from which the unit cell and wavelength data   C
C  are read and a SHELX97 *.fcf file created using the LIST 4 option,  C
C  from which the data and indices are extracted. Output consists of   C
C  a file named *.rtn which has one column in addition to the *.fcf    C
C  file that list the 'anisotropic' statistic (e.g. Rrt(hkl))          C
C  calculated using the RT tensor for each reflection.  Other stuff,   C
C  such as the mean anisotropic statistic (which should be the same    C
C  as its global counterpart), the highest and lowest values of (say)  C
C  Rrt(hkl) and values of (e.g.) Rrt(hkl) for each of the cell axes.   C
C  The eigenvalues and eigenvectors of the R-tensor are also found     C
C  and these are given along with their effective R-values.            C
C  This information is also summarized in a separate file *.RT_out.    C
C                                                                      C
C======================================================================C
C                                                                      C
C  To compile the program try 'f77 -o rt rt2.f' or similar.            C
C                                                                      C
C  To run the program you need a *.fcf file. The 'cif' junk need not   C
C  be edited out. You will also need a *.res file - or just the first  C
C  few lines of it i.e. the lines that specify the cell.               C
C                                                                      C
C======================================================================C
C                                                                      C
C  The two routines eigval and eigvec whose calls are commented out    C
C  below are derived from Appendix B of 'Mathematical Techniques in    C
C  Crystallography and Materials Science' by Ted Prince.  They were    C
C  abandoned in favour of JACOBI (adapted from Numerical Recipes in    C
C  Fortran 2nd Ed. by Press, Teukolsky, Vetterling and Flannery) for   C
C  finding the eigenvalues and eigenvectors because they only appear   C
C  to work properly when the angle phi in eigval is smaller than some  C
C  value (unclear to me just what it is, but seems to be about 45      C
C  degrees). If anybody knows why this is, please let me know.  It     C
C  does not seem to be a precision problem.                            C
C                                                                      C
C======================================================================C
C                                                                      C
C  This program along with the algorithm and ideas expressed within    C
C  it can be used freely by anyone on the sole condition that they     C
C  make adequate reference to it.  The most suitable reference is:     C
C                                                                      C
C  Parkin, S. (2000) "Expansion of Scalar Validation Criteria to 3-D:  C
C  The R-Tensor" Acta Cryst. A 56, 157-162.                            C
C                                                                      C
C======================================================================C

      PROGRAM RT

      CHARACTER*80 RESFIL, FCFFIL, OUTFIL, RTNFIL
      INTEGER HKL, NUMDAT, NHIGH, NLOW
      REAL CELLDS, CELLRS, DATMAT, ORT, LAMBDA, FSFMAT, RCRYST, AVFO,
     &   NORMAL, RHSVEC, TENSOR, TVALUE, THLMMT, EIGVAL, EIGVEC,
     &   HKLORT, ORTHKL, EIGHKL
      DIMENSION HKL(100000,3), DATMAT(100000,3), ORT(100000,4),
     &   FSFMAT(100000,3), TVALUE(100000), CELLDS(6), CELLRS(6),
     &   HKLORT(6), ORTHKL(6), NORMAL(6,6), RHSVEC(6), TENSOR(6),
     &   THLMMT(5), EIGVAL(3), EIGVEC(3,3), EIGHKL(3,3)

c S.Parkin 130303: For auto-run operation, do not call BANNER
c      CALL BANNER

      CALL GTNAME (RESFIL, FCFFIL, OUTFIL, RTNFIL)

      CALL GTDATA (FCFFIL, HKL, DATMAT, NUMDAT)

      CALL GTCELL (RESFIL, CELLDS, CELLRS, LAMBDA)

      CALL CONMAT (HKLORT, ORTHKL, CELLDS, CELLRS)

      CALL TRNSFM (HKLORT, HKL, ORT, NUMDAT, LAMBDA)

      CALL FSQTOF (DATMAT, FSFMAT, NUMDAT, RCRYST, AVFO)

      CALL BLDNML (ORT, FSFMAT, NUMDAT, AVFO, RCRYST, NORMAL, RHSVEC)

      CALL CHLSKY (NORMAL, RHSVEC, TENSOR)

      CALL NPDTST (TENSOR)

      CALL CALCTV (TENSOR, NUMDAT, ORT, RCRYST, THLMMT, TVALUE, NHIGH,
     &   NLOW)

C  The following two subroutines have been superseded by JACOBI
C      CALL EIGVAL (TENSOR, EVAL)
C      CALL EIGVEC (TENSOR, EVAL, EVEC)

      CALL JACOBI (TENSOR, EIGVEC, EIGVAL)

      CALL GTINDX (EIGVEC, ORTHKL, EIGHKL)

      CALL APPLYT (NUMDAT, HKL, DATMAT, RCRYST, TVALUE, RTNFIL)

      CALL OUTPUT (NUMDAT, OUTFIL, RTNFIL, RCRYST, AVFO, HKL, NHIGH,
     &   NLOW, THLMMT, EIGVAL, EIGVEC, EIGHKL, TENSOR)

c      WRITE (*,*)
c      WRITE (*,*) 'Program terminated normally. '
      STOP
      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE BANNER

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  Print some identifying information to the screen.                   C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      WRITE (*,*)
      WRITE (*,'(1X,A78)') '+------------------------------------------'
     &   //'----------------------------------+'
      WRITE (*,'(A2,76X,A1)') '|', '|'
      WRITE (*,'(1X,A78)') '|              Program RT version 2 - (c) S'
     &   //'ean Parkin 99.09.09               |'
      WRITE (*,'(A2,76X,A1)') '|', '|'
      WRITE (*,'(1X,A78)') '|        Anisotropic expansion of global sc'
     &   //'alar validation statistics        |'
      WRITE (*,'(A2,76X,A1)') '|', '|'
      WRITE (*,'(1X,A78)') '+------------------------------------------'
     &   //'----------------------------------+'
      WRITE (*,*)

      RETURN
      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE GTNAME (RESFIL, FCFFIL, OUTFIL, RTNFIL)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  Asks user for filename.                                             C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CHARACTER*80 NAME, RESFIL, FCFFIL, OUTFIL, RTNFIL

      CALL GETARG(1, NAME)
      
      if (name.ne.'') then
         RESFIL = NAME(1:INDEX(NAME,' ')-1)//'.res'
         FCFFIL = NAME(1:INDEX(NAME,' ')-1)//'.fcf'
         OUTFIL = NAME(1:INDEX(NAME,' ')-1)//'.rto'
         RTNFIL = NAME(1:INDEX(NAME,' ')-1)//'.rtn'
      else
         WRITE (*, '(1X,A$)') 'Input filename without extension: '
         READ (*, '(A80)') NAME
         RESFIL = NAME(1:INDEX(NAME,' ')-1)//'.res'
         FCFFIL = NAME(1:INDEX(NAME,' ')-1)//'.fcf'
         OUTFIL = NAME(1:INDEX(NAME,' ')-1)//'.rto'
         RTNFIL = NAME(1:INDEX(NAME,' ')-1)//'.rtn'
      endif

      RETURN
      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE GTDATA(FCFFIL, HKL, DATMAT, NUMDAT)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  Reads the *.fcf file and packs arrays with the data.                C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CHARACTER*80 FCFFIL, LINE
      INTEGER HKL, NUMDAT, N, FILEND
      REAL DATMAT
      DIMENSION HKL(100000,3), DATMAT(100000,3)

      N = 1

c  SHELXL-2013 has different format for data lines in a LIST 4 *.fcf
c  100 FORMAT (3I4, 2F12.2, F10.2)
  100 FORMAT (3I4, 2F14.2, F13.2)

      OPEN (11, FILE = FCFFIL, STATUS = 'OLD', ERR = 9999)

 1001 READ (11,'(A80)') LINE
      IF (LINE(55:55).NE.'o') THEN
         GOTO 1001
      ELSE
         READ (LINE, 100) HKL(N,1), HKL(N,2), HKL(N,3),
     &      DATMAT(N,1), DATMAT(N,2), DATMAT(N,3)
      ENDIF

 1002 CONTINUE
      FILEND = ABS(HKL(N,1)) + ABS(HKL(N,2)) + ABS(HKL(N,3))
      IF (FILEND.EQ.0) GOTO 1003
      N = N + 1
      READ (11,100) HKL(N,1), HKL(N,2), HKL(N,3),
     &   DATMAT(N,1), DATMAT(N,2), DATMAT(N,3)
      GOTO 1002

 1003 CONTINUE

      CLOSE (11)
      N = N - 1
      NUMDAT = N

      RETURN

9999  WRITE(*,*)'There is no such file as ', FCFFIL
c      PAUSE
      STOP

      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE GTCELL (RESFIL, CELLDS, CELLRS, LAMBDA)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  Calculates useful direct and reciprocal cell numbers.               C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CHARACTER*80 RESFIL, LINE
      character*4 word
      REAL CELLDS, CELLRS, PI, LAMBDA, A, B, C, AL, BE, GA, AS, BS, CS,
     &   ALS, BES, GAS, CAL, CBE, CGA, SAL, SBE, SGA, CALS, CBES, CGAS,
     &   VOLUME

      DIMENSION CELLDS(6), CELLRS(6)

      PI = 3.14159265354

  100 FORMAT (4X, A76)

      OPEN (12, FILE = RESFIL, STATUS = 'OLD', ERR = 9999)

c S.Parkin 141007: Modification for shelxl-2014, which writes blank 
c lines in the RES file. 
c Note: This fix no longer works for SHELXL2016/6, but see below ...
c1     read (12, '(A4)') word
c      if (word.ne."TITL") then
c         goto 1
c      endif
c      read (12, 100) line

c S.Parkin 161109: Modification for SHELXL2016/6, which again confuses the
c reading in of lambda and the cell dimensions.
1     read (12, '(A80)') line
      word = line(1:4)
      if (word.eq."CELL") then
         close (12)
      else
        goto 1
      endif
      
      OPEN (12, STATUS = 'SCRATCH')
      LINE = LINE(5:80)
      WRITE (12, 100) LINE
      REWIND (12)
      READ (12, *) LAMBDA, CELLDS(1), CELLDS(2), CELLDS(3), CELLDS(4),
     &   CELLDS(5), CELLDS(6)
      CLOSE (12)

      A = CELLDS(1)
      B = CELLDS(2)
      C = CELLDS(3)
      AL = CELLDS(4) * PI / 180.0
      BE = CELLDS(5) * PI / 180.0
      GA = CELLDS(6) * PI / 180.0
      CAL = COS(AL)
      CBE = COS(BE)
      CGA = COS(GA)
      SAL = SIN(AL)
      SBE = SIN(BE)
      SGA = SIN(GA)
      CALS = ((CBE * CGA) - CAL) / (SBE * SGA)
      CBES = ((CAL * CGA) - CBE) / (SAL * SGA)
      CGAS = ((CAL * CBE) - CGA) / (SAL * SBE)
      VOLUME = A * B * C * SQRT(1.0 - (CAL ** 2.0) - (CBE ** 2.0) -
     &         (CGA ** 2.0) + (2.0 * CAL * CBE * CGA))
      AS = B * C * SAL / VOLUME
      BS = A * C * SBE / VOLUME
      CS = A * B * SGA / VOLUME
      ALS = ACOS(CALS)
      BES = ACOS(CBES)
      GAS = ACOS(CGAS)

      CELLRS(1) = AS
      CELLRS(2) = BS
      CELLRS(3) = CS
      CELLRS(4) = ALS * 180.0 / PI
      CELLRS(5) = BES * 180.0 / PI
      CELLRS(6) = GAS * 180.0 / PI

      RETURN

9999  WRITE(*,*) 'THERE IS NO SUCH FILE AS ', RESFIL
      STOP

      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE CONMAT (HKLORT, ORTHKL, CELLDS, CELLRS)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  Calculates conversion matrices to take the hkl indices based on the C
C  reciprocal cell and give direction cosines in an orthonormal system C
C  and vice versa.                                                     C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      REAL HKLORT, ORTHKL, CELLDS, CELLRS, DTMNNT, PI, C, AS, BS, CS,
     &   CAL, CBES, SBES, CGAS, SGAS
      DIMENSION HKLORT(6), ORTHKL(6), CELLDS(6), CELLRS(6)

      PI = 3.14159265354

      C = CELLDS(3)
      AS = CELLRS(1)
      BS = CELLRS(2)
      CS = CELLRS(3)
      CAL = COS(CELLDS(4) * PI / 180.0)
      CBES = COS(CELLRS(5) * PI / 180.0)
      SBES = SIN(CELLRS(5) * PI / 180.0)
      CGAS = COS(CELLRS(6) * PI / 180.0)
      SGAS = SIN(CELLRS(6) * PI / 180.0)

      HKLORT(1) = AS
      HKLORT(2) = BS * SGAS
      HKLORT(3) = 1.0 / C
      HKLORT(4) = BS * CGAS
      HKLORT(5) = CS * CBES
      HKLORT(6) = -(CS * SBES * CAL)

      DTMNNT = HKLORT(1) * HKLORT(2) * HKLORT(3)

      ORTHKL(1) = 1.0 / HKLORT(1)
      ORTHKL(2) = 1.0 / HKLORT(2)
      ORTHKL(3) = 1.0 / HKLORT(3)
      ORTHKL(4) = - (HKLORT(4) * HKLORT(3)) / DTMNNT
      ORTHKL(5) = ((HKLORT(4) * HKLORT(6)) - (HKLORT(5) * HKLORT(2))) /
     &   DTMNNT
      ORTHKL(6) = - (HKLORT(1) * HKLORT(6)) / DTMNNT

      RETURN
      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE TRNSFM (HKLORT, HKL, ORT, NUMDAT, LAMBDA)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  Transforms the natural coordinate system of the reciprocal lattice  C
C  into an orthonormal coordinate system.  Note that ORT(COUNT,4) is   C
C  not used in RT.                                                     C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER HKL, NUMDAT, COUNT
      REAL HKLORT, ORT, HO, KO, LO, LENGTH, LAMBDA
      DIMENSION HKLORT(6), HKL(100000,3), ORT(100000,4)

      DO 1 COUNT = 1, NUMDAT
         HO = (HKLORT(1) * HKL(COUNT,1)) +
     &        (HKLORT(4) * HKL(COUNT,2)) +
     &        (HKLORT(5) * HKL(COUNT,3))
         KO = (HKLORT(2) * HKL(COUNT,2)) +
     &        (HKLORT(6) * HKL(COUNT,3))
         LO = (HKLORT(3) * HKL(COUNT,3))
         LENGTH = SQRT((HO ** 2.0) + (KO ** 2.0) + (LO ** 2.0))
         ORT(COUNT,1) = HO / LENGTH
         ORT(COUNT,2) = KO / LENGTH
         ORT(COUNT,3) = LO / LENGTH
         ORT(COUNT,4) = (LENGTH * LAMBDA) / 2.0
    1 CONTINUE

      RETURN
      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE FSQTOF (DATMAT, FSFMAT, NUMDAT, RCRYST, AVFO)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  Converts the data given as F^2 in the *.fcf file into F.  Note the  C
C  approximation used for negative F^2.  It gives F as the square root C
C  of 0.7 * sigma(F^2).  In practice this doesn't seem to make much    C
C  difference, but if this is unreasonable, please let me know.        C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER NUMDAT, COUNT
      REAL DATMAT, FSFMAT, RCRYST, SUMDF, SUMF, AVFO
      DIMENSION DATMAT(100000,3), FSFMAT(100000,3)

      SUMDF = 0.0
      SUMF = 0.0

      DO 1 COUNT = 1, NUMDAT
         FSFMAT(COUNT,1) = SQRT(DATMAT(COUNT,1))
         IF (DATMAT(COUNT,2).GE.0.0) THEN
            FSFMAT(COUNT,2) = SQRT(DATMAT(COUNT,2))
         ELSE
            FSFMAT(COUNT,2) = SQRT(0.7 * DATMAT(COUNT,3))
         ENDIF
         FSFMAT(COUNT,3) = DATMAT(COUNT,3) / (2 * FSFMAT(COUNT,2))
         SUMDF = SUMDF + (ABS(FSFMAT(COUNT,2) - FSFMAT(COUNT,1)))
         SUMF = SUMF + FSFMAT(COUNT,2)
    1 CONTINUE

      RCRYST = SUMDF / SUMF
      AVFO = SUMF / NUMDAT

      RETURN
      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE BLDNML (ORT, FSFMAT, NUMDAT, AVFO, RCRYST, NORMAL,
     &   RHSVEC)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  Build the normal equations matrix as a sum of contributions from    C
C  individual reflections.                                             C
C  Note that although the WEIGHT looks like it is not needed, it is    C
C  included so that some other (non-unit) weighting scheme can be      C
C  easily incorporated.                                                C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER NUMDAT, COUNT, X, Y
      REAL ORT, FSFMAT, AVFO, RCRYST, NORMAL, RHSVEC, MATRIX, VECTOR,
     &   HO2, WEIGHT, RHKL, SPDV, WR2S4, WRRS2
      DIMENSION ORT(100000,4), FSFMAT(100000,3), NORMAL(6,6),
     &   RHSVEC(6), MATRIX(6,6), VECTOR(6), HO2(7)

      DO 2 X = 1, 6
         VECTOR(X) = 0.0
         RHSVEC(X) = 0.0
         DO 1 Y = 1, 6
            MATRIX(X, Y) = 0.0
            NORMAL(X, Y) = 0.0
    1    CONTINUE
    2 CONTINUE

      DO 7 COUNT = 1, NUMDAT
         HO2(1) = ORT(COUNT,1) ** 2.0
         HO2(2) = ORT(COUNT,2) ** 2.0
         HO2(3) = ORT(COUNT,3) ** 2.0
         HO2(4) = 2.0 * ORT(COUNT,1) * ORT(COUNT,2)
         HO2(5) = 2.0 * ORT(COUNT,1) * ORT(COUNT,3)
         HO2(6) = 2.0 * ORT(COUNT,2) * ORT(COUNT,3)
         HO2(7) = HO2(1) + HO2(2) + HO2(3)
         WEIGHT = 1.0
C         weight = 1.0 / (fsfmat(count,3) * fsfmat(count,2))
         RHKL = ABS(FSFMAT(COUNT,2) - FSFMAT(COUNT,1)) / AVFO
         WR2S4 = (WEIGHT * RCRYST ** 2.0) / (HO2(7) ** 2.0)
         WRRS2 = (WEIGHT * RCRYST * RHKL) / HO2(7)

         DO 4 X = 1, 6
            VECTOR(X) = WRRS2 * HO2(X)
            SPDV = WR2S4 * HO2(X)
            DO 3 Y = 1, 6
               MATRIX(X,Y) = SPDV * HO2(Y)
    3       CONTINUE
    4    CONTINUE

         DO 6 X = 1, 6
            RHSVEC(X) = RHSVEC(X) + VECTOR(X)
            DO 5 Y = 1, 6
               NORMAL(X,Y) = NORMAL(X,Y) + MATRIX(X,Y)
    5       CONTINUE
    6    CONTINUE

    7 CONTINUE

      RETURN
      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE CHLSKY (NORMAL, RHSVEC, TENSOR)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  Solves the normal equations etc.                                    C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER X, Y, Z
      REAL NORMAL, DIAGNL, RHSVEC, TENSOR, SUM
      DIMENSION NORMAL(6,6), DIAGNL(6), TENSOR(6), RHSVEC(6)

      DO 1 X = 1, 6
         DO 2 Y = X, 6
            SUM = NORMAL(X,Y)
            DO 3 Z = X-1, 1, -1
               SUM = SUM - NORMAL(X,Z) * NORMAL(Y,Z)
    3       CONTINUE
            IF (X.EQ.Y) THEN
               DIAGNL(X) = SQRT(SUM)
            ELSE
               NORMAL(Y,X) = SUM / DIAGNL(X)
            ENDIF
    2    CONTINUE
    1 CONTINUE

      DO 4 X = 1, 6
         SUM = RHSVEC(X)
         DO 5 Z = X-1, 1, -1
            SUM = SUM - NORMAL(X,Z) * TENSOR(Z)
    5    CONTINUE
         TENSOR(X) = SUM / DIAGNL(X)
    4 CONTINUE

      DO 6 X = 6, 1, -1
         SUM = TENSOR(X)
         DO 7 Z = X + 1, 6
            SUM = SUM - NORMAL(Z,X) * TENSOR(Z)
    7    CONTINUE
         TENSOR(X) = SUM / DIAGNL(X)
    6 CONTINUE

      RETURN
      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE NPDTST (TENSOR)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  Tests the R-TENSOR to make sure it is positive definite.  This may  C
C  not be needed, but it is fast and may be useful in rare instances.  C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER X, Y, Z
      REAL TENSOR, MATRIX, DIAGNL, SUM
      DIMENSION TENSOR(6), MATRIX(3,3), DIAGNL(3)

      MATRIX(1,1) = TENSOR(1)
      MATRIX(2,2) = TENSOR(2)
      MATRIX(3,3) = TENSOR(3)
      MATRIX(1,2) = TENSOR(4)
      MATRIX(2,1) = TENSOR(4)
      MATRIX(1,3) = TENSOR(5)
      MATRIX(3,1) = TENSOR(5)
      MATRIX(2,3) = TENSOR(6)
      MATRIX(3,2) = TENSOR(6)

      DO 1 X = 1, 3
         DO 2 Y = X, 3
            SUM = MATRIX(X,Y)
            DO 3 Z = X - 1, 1, -1
               SUM = SUM - MATRIX(X,Z) * MATRIX(Y,Z)
    3       CONTINUE
            IF (X.EQ.Y) THEN
C               IF (SUM.LE.0.0) GOTO 9999
               DIAGNL(X) = SQRT(SUM)
            ELSE
               MATRIX(Y,X) = SUM / DIAGNL(X)
            ENDIF
    2    CONTINUE
    1 CONTINUE

      RETURN

C9999 PAUSE 'R-TENSOR IS NON-POSITIVE DEFINITE '
C     STOP

      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE CALCTV (TENSOR, NUMDAT, ORT, RCRYST, THLMMT, TVALUE,
     &   NHIGH, NLOW)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  Calculates the T values for each datapoint in the *.fcf file.       C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER NUMDAT, COUNT, NHIGH, NLOW
      REAL TENSOR, ORT, RCRYST, THLMMT, TVALUE, HO2, SMT
      DIMENSION TENSOR(6), ORT(100000,4), THLMMT(5), TVALUE(100000),
     &   HO2(7)

      SMT = 0.0
      THLMMT(1) = 1.0
      THLMMT(2) = 1.0

      DO 2 COUNT = 1, NUMDAT
         HO2(1) = ORT(COUNT,1) ** 2.0
         HO2(2) = ORT(COUNT,2) ** 2.0
         HO2(3) = ORT(COUNT,3) ** 2.0
         HO2(4) = 2.0 * ORT(COUNT,1) * ORT(COUNT,2)
         HO2(5) = 2.0 * ORT(COUNT,1) * ORT(COUNT,3)
         HO2(6) = 2.0 * ORT(COUNT,2) * ORT(COUNT,3)
         HO2(7) = HO2(1) + HO2(2) + HO2(3)

         TVALUE(COUNT) = ((TENSOR(1) * HO2(1)) + (TENSOR(2) * HO2(2)) +
     &      (TENSOR(3) * HO2(3)) + (TENSOR(4) * HO2(4)) +
     &      (TENSOR(5) * HO2(5)) + (TENSOR(6) * HO2(6))) / HO2(7)

         IF (TVALUE(COUNT).GT.THLMMT(1)) THEN
            THLMMT(1) = TVALUE(COUNT)
            NHIGH = COUNT
         ELSE IF (TVALUE(COUNT).LT.THLMMT(2)) THEN
            THLMMT(2) = TVALUE(COUNT)
            NLOW = COUNT
         ENDIF
         SMT = SMT + TVALUE(COUNT)
    2 CONTINUE

      THLMMT(1) = TVALUE(NHIGH)
      THLMMT(2) = TVALUE(NLOW)
      THLMMT(3) = SMT / NUMDAT
      THLMMT(4) = TVALUE(NUMDAT / 2)
      THLMMT(5) = (TVALUE(NUMDAT / 4) + (2.0 * TVALUE(NUMDAT / 2)) +
     &   TVALUE(3 * NUMDAT / 4)) / 4.0

      RETURN
      END

C=======================================================================
C
C----------------------------------------------------------------------C
C  The two routines eigval and eigvec that are commented out here are  C
C  derived from Prince 'Mathematical Techniques in Crystallography and C
C  Materials Science'.  They were abandoned in favour of JACOBI for    C
C  finding the eigenvalues and eigenvectors.  They only seem to work   C
C  properly when the angle phi in eigval is smaller than some value.   C
C  If  anybody knows why this is, please let me know.  It does not     C
C  seem to be a precision problem.  The variables eval and evec here   C
C  are named eigval and eigvec in the program itself, which are the    C
C  names I gave to the subroutines, but that is not the problem.       C
C----------------------------------------------------------------------C
C
C      subroutine eigval (TENSOR, eval)
C
C      real a, b, c, TENSOR, eval, q, r, phi, cc, pi
C      dimension TENSOR(6), eval(3)
C
C      pi = 3.14159265354
C
C      a = 0.0 - (TENSOR(1) + TENSOR(2) + TENSOR(3))
C      b = ((TENSOR(1) * TENSOR(2)) - (TENSOR(4) ** 2.0)) +
C     &    ((TENSOR(1) * TENSOR(3)) - (TENSOR(5) ** 2.0)) +
C     &    ((TENSOR(2) * TENSOR(3)) - (TENSOR(6) ** 2.0))
C      c = -((TENSOR(1) * TENSOR(2) * TENSOR(3)) +
C     &     (2.0 * TENSOR(4) * TENSOR(5) * TENSOR(6)) -
C     &     (TENSOR(1) * TENSOR(6) ** 2.0) -
C     &     (TENSOR(2) * TENSOR(5) ** 2.0) -
C     &     (TENSOR(3) * TENSOR(4) ** 2.0))
C
C      q = ((a ** 2.0) / 3.0) - b
C      r = -(((2.0 * a ** 3.0) / 27.0) - ((a * b) / 3.0) + c)
C      phi = acos(sqrt((27.0 * r * r) / (4.0 * q ** 3.0)))
C      write (*,*) 'PHI EQUALS ', phi, phi * 180.0 / pi
C
C      cc = sqrt((4.0 * q) / 3.0)
C
C      eval(1) = cc * cos(phi / 3.0) - (a / 3.0)
C      eval(2) = cc * cos(((2.0 * pi) - phi) / 3.0) - (a / 3.0)
C      eval(3) = cc * cos(((2.0 * pi) + phi) / 3.0) - (a / 3.0)
C
C      return
C      end
C
C=======================================================================
C
C-----------------------------------------------------------------------
C
C      subroutine eigvec (TENSOR, eval, evec)
C
C      integer count
C      real TENSOR, eval, evec, x1, y1, a11, a22, a33, a12, a13, a23,
C     &   delta, d
C      dimension TENSOR(6), eval(3), evec(3, 3), tttttt(3)
C
C      a11 = TENSOR(1)
C      a22 = TENSOR(2)
C      a33 = TENSOR(3)
C      a12 = TENSOR(4)
C      a13 = TENSOR(5)
C      a23 = TENSOR(6)
C
C      do 1 count = 1, 3
C
C         delta = ((a11 - eval(count)) * (a22 - eval(count))) -
C     &      a12 ** 2.0
C
C         x1 = (0.0 - (a13 * (a22 - eval(count))) +
C     &      (a12 * a23)) / delta
C         y1 = (0.0 - (a23 * (a11 - eval(count))) +
C     &      (a12 * a13)) / delta
C
C        d = sqrt((x1 ** 2.0) + (y1 ** 2.0) + 1.0)
C
C         evec(count, 1) = x1 / d
C         evec(count, 2) = y1 / d
C         evec(count, 3) = 1.0 / d
C
C    1 continue
C
C      return
C      end
C
C========================================================================

C-----------------------------------------------------------------------

      SUBROUTINE JACOBI (TENSOR, EIGVEC, EIGVAL)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This subroutine performs the Jacobi rotation and sorts the          C
C  eigenvalues and eigenvectors (adapted from Numerical Recipes in     C
C  Fortran 2nd Ed. by Press, Teukolsky, Vetterling and Flannery).      C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER NUMROT, I, J, K, P, Q
      REAL TENSOR, MATRIX, EIGVEC, EIGVAL, BLD, ACC, SINPHI,
     &   COSPHI, ARP, ARQ, SUM, T, TAU, THETA, THRESH, SWAP
      DIMENSION TENSOR(6), MATRIX(3,3), EIGVEC(3,3), EIGVAL(3), BLD(3),
     &   ACC(3)

C  Fill matrix with the appropriate values from TENSOR, set NUMROT
C  to zero, initialize EIGVEC to the identity matrix, initialize
C  EIGVAL and BLD to the diagonal of MATRIX, and fill ACC with zeros.

      MATRIX(1,1) = TENSOR(1)
      MATRIX(2,2) = TENSOR(2)
      MATRIX(3,3) = TENSOR(3)
      MATRIX(1,2) = TENSOR(4)
      MATRIX(2,1) = TENSOR(4)
      MATRIX(1,3) = TENSOR(5)
      MATRIX(3,1) = TENSOR(5)
      MATRIX(2,3) = TENSOR(6)
      MATRIX(3,2) = TENSOR(6)

      NUMROT = 0
      DO 1 P = 1, 3
         DO 2 Q = 1, 3
            EIGVEC(P,Q) = 0.0
2        CONTINUE
         EIGVEC(P,P) = 1.0
         EIGVAL(P) = MATRIX(P,P)
         BLD(P) = MATRIX(P,P)
         ACC(P) = 0.0
1     CONTINUE

C  Jacobi proper starts here.

      DO 3 I = 1, 50
         SUM = 0.0
         DO 4 P = 1, 2
            DO 5 Q = P+1, 3
               SUM = SUM + ABS(MATRIX(P,Q))
5           CONTINUE
4        CONTINUE
         IF (SUM.EQ.0.0) GO TO 1000
         IF (I.LT.4) THEN
            THRESH = 0.02 * SUM
         ELSE
            THRESH = 0.0
         ENDIF
         DO 6 P = 1, 3
            DO 7 Q = P+1, 3
               ARP = 100.0 * ABS(MATRIX(P,Q))
               IF ((I.GT.4).AND.(ABS(EIGVAL(P))+ARP.EQ.ABS(EIGVAL(P)))
     &            .AND.(ABS(EIGVAL(Q))+ARP.EQ.ABS(EIGVAL(Q)))) THEN
                  MATRIX(P,Q) = 0.0
               ELSE IF (ABS(MATRIX(P,Q)).GT.THRESH) THEN
                  ARQ = EIGVAL(Q) - EIGVAL(P)
                  IF (ABS(ARQ)+ARP.EQ.ABS(ARQ)) THEN
                     T = MATRIX(P,Q) / ARQ
                  ELSE
                     THETA = 0.5 * ARQ / MATRIX(P,Q)
                     T = 1.0 / (ABS(THETA) + SQRT(1.0 + THETA ** 2.0))
                     IF (THETA.LT.0.0) T = 0.0 - T
                  ENDIF
                  COSPHI = 1.0 / SQRT(1.0 + T ** 2.0)
                  SINPHI = T * COSPHI
                  TAU = SINPHI / (1.0 + COSPHI)
                  ARQ = T * MATRIX(P,Q)
                  ACC(P) = ACC(P) - ARQ
                  ACC(Q) = ACC(Q) + ARQ
                  EIGVAL(P) = EIGVAL(P) - ARQ
                  EIGVAL(Q) = EIGVAL(Q) + ARQ
                  MATRIX(P,Q) = 0.0
                  DO 8 J = 1, P-1
                     ARP = MATRIX(J,P)
                     ARQ = MATRIX(J,Q)
                     MATRIX(J,P) = ARP - SINPHI * (ARQ + ARP * TAU)
                     MATRIX(J,Q) = ARQ + SINPHI * (ARP - ARQ * TAU)
8                 CONTINUE
                  DO 9 J = P+1, Q-1
                     ARP = MATRIX(P,J)
                     ARQ = MATRIX(J,Q)
                     MATRIX(P,J) = ARP - SINPHI * (ARQ + ARP * TAU)
                     MATRIX(J,Q) = ARQ + SINPHI * (ARP - ARQ * TAU)
9                 CONTINUE
                  DO 10 J = Q+1, 3
                     ARP = MATRIX(P,J)
                     ARQ = MATRIX(Q,J)
                     MATRIX(P,J) = ARP - SINPHI * (ARQ + ARP * TAU)
                     MATRIX(Q,J) = ARQ + SINPHI * (ARP - ARQ * TAU)
10                CONTINUE
                  DO 11 J = 1, 3
                     ARP = EIGVEC(J,P)
                     ARQ = EIGVEC(J,Q)
                     EIGVEC(J,P) = ARP - SINPHI * (ARQ + ARP * TAU)
                     EIGVEC(J,Q) = ARQ + SINPHI * (ARP - ARQ * TAU)
11                CONTINUE
                  NUMROT = NUMROT + 1
               ENDIF
7           CONTINUE
6        CONTINUE
         DO 12 P = 1, 3
            BLD(P) = BLD(P) + ACC(P)
            EIGVAL(P) = BLD(P)
            ACC(P) = 0.0
12       CONTINUE
3     CONTINUE

1000  CONTINUE

C  Jacobi is done, now sort the eigenvalues and eigenvectors into
C  descending order.

      DO 13 I = 1, 2
         K = I
         SWAP = EIGVAL(I)
         DO 14 J = I + 1, 3
            IF (EIGVAL(J).GE.SWAP) THEN
               K = J
               SWAP = EIGVAL(J)
            ENDIF
14       CONTINUE
         IF (K.NE.I) THEN
            EIGVAL(K) = EIGVAL(I)
            EIGVAL(I) = SWAP
            DO 15 J = 1, 3
               SWAP = EIGVEC(J,I)
               EIGVEC(J,I) = EIGVEC(J,K)
               EIGVEC(J,K) = SWAP
15          CONTINUE
         ENDIF
13    CONTINUE

      RETURN
      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE GTINDX (EIGVEC, ORTHKL, EIGHKL)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  Converts orthogonal coordinates (direction cosines) of the eigen-   C
C  values and eigenvectors back into HKL (fractional) indices.         C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER COUNT
      REAL EIGVEC, ORTHKL, EIGHKL
      DIMENSION EIGVEC(3,3), ORTHKL(6), EIGHKL(3,3)

      DO 1 COUNT = 1, 3
         EIGHKL(1,COUNT) = (ORTHKL(1) * EIGVEC(1,COUNT)) +
     &                     (ORTHKL(4) * EIGVEC(2,COUNT)) +
     &                     (ORTHKL(5) * EIGVEC(3,COUNT))
         EIGHKL(2,COUNT) = (ORTHKL(2) * EIGVEC(2,COUNT)) +
     &                     (ORTHKL(4) * EIGVEC(3,COUNT))
         EIGHKL(3,COUNT) = (ORTHKL(3) * EIGVEC(3,COUNT))
    1 CONTINUE

      RETURN
      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE APPLYT (NUMDAT, HKL, DATMAT, RCRYST, TVALUE, RTNFIL)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  Writes a file *.rtn that is the same as *.fcf apart from it has an  C
C  an extra column which lists the effective R-values obtained by the  C
C  R-TENSOR for each datapoint present in the dataset.                 C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CHARACTER*80 RTNFIL
      INTEGER NUMDAT, HKL, COUNT
      REAL DATMAT, TVALUE, RCRYST
      DIMENSION HKL(100000,3), DATMAT(100000,3), TVALUE(100000)

  100 FORMAT (3I4, 2F14.2, F13.2, 2X, F8.5)

C Modification by H.Hope for setting BBEdit filetype on Mac
C Removed by S.Parkin 02/08/30 for g77 OSX version
C      OPEN (12, FILE = RTNFIL, STATUS = 'UNKNOWN' ,CREATOR='R*ch')
      OPEN (12, FILE = RTNFIL, STATUS = 'UNKNOWN')

      DO 1 COUNT  = 1, NUMDAT
         WRITE (12, 100) HKL(COUNT,1), HKL(COUNT,2), HKL(COUNT,3),
     &      DATMAT(COUNT,1), DATMAT(COUNT,2), DATMAT(COUNT,3),
     &      TVALUE(COUNT) * RCRYST
    1 CONTINUE

      CLOSE (12)

      RETURN
      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE OUTPUT (NUMDAT, OUTFIL, RTNFIL, RCRYST, AVFO, HKL,
     &   NHIGH, NLOW, THLMMT, EIGVAL, EIGVEC, EIGHKL, TENSOR)

      CHARACTER*80 NAME, OUTFIL, RTNFIL, LINE
      INTEGER NUMDAT, HKL, NHIGH, NLOW, COUNT
      REAL RCRYST, AVFO, THLMMT, EIGVAL, EIGVEC, EIGHKL, TENSOR
      DIMENSION HKL(100000,3), THLMMT(5), EIGVAL(3), EIGVEC(3,3),
     &   EIGHKL(3,3), TENSOR(6)

      NAME = RTNFIL(1:INDEX(RTNFIL,'.')-1)
      TRACE = TENSOR(1) + TENSOR(2) + TENSOR(3)
      geotra = (TENSOR(1) * TENSOR(2) * TENSOR(3)) ** (1./3.)
      DTMNNT = (TENSOR(1) * TENSOR(2) * TENSOR(3)) -
     &   (TENSOR(1) * TENSOR(6) ** 2.0) +
     &   (TENSOR(4) * TENSOR(5) * TENSOR(6)) -
     &   (TENSOR(3) * TENSOR(4) ** 2.0) +
     &   (TENSOR(4) * TENSOR(5) * TENSOR(6)) +
     &   (TENSOR(2) * TENSOR(5) ** 2.0)

C Mofification by H.Hope for setting filetype to BBEdit
C Removed by S.Parkin 02/08/30 for OSX version
C      OPEN (13, FILE = OUTFIL, STATUS = 'UNKNOWN',CREATOR='R*ch')
      OPEN (13, FILE = OUTFIL, STATUS = 'UNKNOWN')

      WRITE (13,'(1X,A78)') '+========================================='
     &   //'===================================+'
      WRITE (13,'(A2,76X,A1)') '|', '|'
      WRITE (13,'(1X,A78)') '|        Output from program RT version 2 '
     &   //'- (c) Sean Parkin 99.09.09         |'
      WRITE (13,'(1X,A78)') '|                [Parkin, S. (2000) Acta C'
     &   //'ryst. A56, 157-162]                |'
      WRITE (13,'(A2,76X,A1)') '|', '|'
      WRITE (13,'(1X,A78)') '+========================================='
     &   //'===================================+'
      WRITE (13,'(A2,76X,A1)') '|', '|'
      WRITE (13,'(1X,A39,A10,A29)') '|   Results of the R-TENSOR applie'
     &   //'d to ',NAME,'                            |'
      WRITE (13,'(A2,76X,A1)') '|', '|'
      WRITE (13,'(1X,A30,I8,A40)') '|      Number of reflections = ',
     &   NUMDAT,'                                       |'
      WRITE (13,'(1X,A30,F8.4,A40)') '|      Conventional R1       = ',
     &   RCRYST,'                                       |'
      WRITE (13,'(1X,A30,F8.2,A40)') '|      Average Fobs          = ',
     &   AVFO,'                                       |'
      WRITE (13,'(A2,76X,A1)') '|', '|'
      WRITE (13,'(1X,A78)') '+-----------------------------------------'
     &   //'-----------------------------------+'
      WRITE (13,'(A2,76X,A1)') '|', '|'
      WRITE (13,'(1X,A78)') '|   The R-TENSOR is                       '
     &   //'                                   |'
      WRITE (13,'(A2,76X,A1)') '|', '|'
      WRITE (13,'(A2,16X,3(F10.5,2X),24X,A1)') '|', TENSOR(1),
     &   TENSOR(4), TENSOR(5), '|'
      WRITE (13,'(A2,16X,3(F10.5,2X),24X,A1)') '|', TENSOR(4),
     &   TENSOR(2), TENSOR(6), '|'
      WRITE (13,'(A2,16X,3(F10.5,2X),24X,A1)') '|', TENSOR(5),
     &   TENSOR(6), TENSOR(3), '|'
      WRITE (13,'(A2,76X,A1)') '|', '|'
      WRITE (13,'(1X,A78)') '+-----------------------------------------'
     &   //'-----------------------------------+'
      WRITE (13,'(A2,76X,A1)') '|', '|'
      WRITE (13,'(1X,A78)') '|   Eigenvalues of the R-TENSOR give extre'
     &   //'ma of possible R-values,           |'
      WRITE (13,'(1X,A78)') '|   eigenvectors are expressed as fraction'
     &   //'al (h,k,l) indices                 |'
      WRITE (13,'(A2,76X,A1)') '|', '|'
      WRITE (13,'(1X,A25,F8.4,A8,3F8.3,A13)') '|      max. R-value    ='
     &   //' ',EIGVAL(1) * RCRYST, '  along ', EIGHKL(1,1),
     &   EIGHKL(2,1), EIGHKL(3,1), '            |'
      WRITE (13,'(1X,A25,F8.4,A8,3F8.3,A13)') '|      mid. R-value    ='
     &   //' ',EIGVAL(2) * RCRYST, '  along ', EIGHKL(1,2),
     &   EIGHKL(2,2), EIGHKL(3,2), '            |'
      WRITE (13,'(1X,A25,F8.4,A8,3F8.3,A13)') '|      min. R-value    ='
     &   //' ',EIGVAL(3) * RCRYST, '  along ', EIGHKL(1,3),
     &   EIGHKL(2,3), EIGHKL(3,3), '            |'
      WRITE (13,'(1X,A25,F8.4,A45)') '|      ellipticity     = ',
     &   EIGVAL(1) / EIGVAL(3), '                                      '
     &   //'      |'
      WRITE (13,'(A2,76X,A1)') '|', '|'
      WRITE (13,'(1X,A78)') '|   For directions represented in the data'
     &   //'set, the R-tensor gives            |'
      WRITE (13,'(A2,76X,A1)') '|', '|'
      WRITE (13,'(1X,A35,F8.4,A8,3I4,A15)') '|      highest effective R'
     &   //'-value = ', THLMMT(1) * RCRYST, '  along ', HKL(NHIGH,1),
     &   HKL(NHIGH,2), HKL(NHIGH,3), '              |'
      WRITE (13,'(1X,A35,F8.4,A8,3I4,A15)') '|      lowest effective R-'
     &   //'value  = ', THLMMT(2) * RCRYST, '  along ', HKL(NLOW,1),
     &   HKL(NLOW,2), HKL(NLOW,3), '              |'
      WRITE (13,'(1X,A35,F8.4,A35)') '|      mean effective R-value    '
     &   //'= ', THLMMT(3) * RCRYST, '                                |'
      WRITE (13,'(1X,A35,F8.4,A35)') '|      median effective R-value  '
     &   //'= ', THLMMT(4) * RCRYST, '                                |'
      WRITE (13,'(1X,A35,F8.4,A35)') '|      trimean effective R-value '
     &   //'= ', THLMMT(5) * RCRYST, '                                |'
      WRITE (13,'(1X,A35,F8.4,A35)') '|      ratio R(high)/R(low)      '
     &   //'= ', THLMMT(1) / THLMMT(2), '                              '
     &   //'  |'
      WRITE (13,'(A2,76X,A1)') '|', '|'
      WRITE (13,'(1X,A78)') '|   Axial R-values, which need not be repr'
     &   //'esented in the dataset are         |'
      WRITE (13,'(A2,76X,A1)') '|', '|'
      WRITE (13,'(1X,A45,F8.4,A25)') '|      R-value along a*, (h00) re'
     &   //'flections = ', TENSOR(1) * RCRYST, '                       '
     &   //' |'
      WRITE (13,'(1X,A45,F8.4,A25)') '|      R-value along b*, (0k0) re'
     &   //'flections = ', TENSOR(2) * RCRYST, '                       '
     &   //' |'
      WRITE (13,'(1X,A45,F8.4,A25)') '|      R-value along c*, (00l) re'
     &   //'flections = ', TENSOR(3) * RCRYST, '                       '
     &   //' |'
      WRITE (13,'(A2,76X,A1)') '|', '|'
      WRITE (13,'(1X,A78)') '+========================================='
     &   //'===================================+'

      REWIND (13)

c S.Parkin 130303: For auto-run operation, suppress output to screen by
c commenting out the next five lines.
      WRITE (*,*)
      DO 1 COUNT = 1, 46
         READ (13,'(A80)') LINE
         WRITE (*,*) LINE
    1 CONTINUE

      CLOSE (13)

      RETURN
      END

C=======================================================================
