**==aa0001.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*********************************
C**       Program: LCHK         **
C*********************************
C
C
C      LCHK:  Analyse various characteristics of Laue pictures
C
C AUTHOR : Original PERQ version - Mike Elder  1986
C          VAX/VMS version - Ian Clifton   Mar 1988
C          S E R C Daresbury  WA4 4AD UK.
C          Increased portability etc. - John Campbell 1989, 1991
C UPDATES :
C
      INCLUDE 'implicit_none.finc'
C
      COMMON /AMAT/A,AMATX
      COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN
      COMMON /MATRIX/MSYS,MROT,MBEAM,ILGP,ICENT
      COMMON /SPOTS/NR
      COMMON /SYS/ISYS,KSYS
      COMMON /XYZ/IHKL
      REAL CF,DMIN,SPNDLE,WVM,WVS,WVX
      INTEGER ICENT,ILGP,ISYS,MBEAM,MROT,MSYS,NR,NUMA,NUMC
      REAL A(3,3),AMATX(3,3),CCELL(6),CELL(6),PHIS(3)
      INTEGER IHKL(20000),KSYS(3)
C
      INCLUDE 'chanel.finc'
C
C====== FUNCTION TYPES:
C
      INTEGER LUNSTI, LUNSTO
C
C====== LOCALS:
C
      LOGICAL QUIT
      REAL DSTMAX
      INTEGER I,ICHANG,IENT,IFLAG,J
      CHARACTER*1 CHAR
      INTEGER IFAIL
C
C-------------------------------------------------------------------------------
C
C
      ITERM = LUNSTI(0)
      IPRINT = LUNSTO(0)
      IUNIT7 = 3
      CALL TRWSET(ITERM,IPRINT,0,0,0)
      CALL QGINIT(0,IPRINT,0,ITERM,IPRINT)
C
C====== Log file (uses IUNIT6, and units 1&2)
C
      IFAIL = 0
      CALL CCPDPN(-1,'LCHKLOG','NEW','F',0,IFAIL)
      CALL CCPDPN(-2,'LCHKJUNK','NEW','F',0,IFAIL)
C
      IENT = 0
 20   CALL QSCLRA
      CALL QSSCRL(20,24)
      CALL QSREVV
      CALL QSSETC(1,1)
      CALL TSTRNG(
     +      'LCHK - Last update 17/10/91 Authors: M.Elder, I.J. Clifton'
     +      )
      WRITE (1,*) 
     +    ' LCHK - Last update 17/10/91   Authors: M.Elder, I.J.Clifton'
      CALL QSNORV
      CALL SCREEN(IENT,ICHANG)
      CALL QSCLRA
      CALL QSSCRL(1,24)
C
C====== Calculate the master matrix, a, from the cell and phi angles.
C
C
C====== Apply the scale factor to the reciprocal lattice lengths first.
C
      CELL(1) = CELL(1)*WVS
      CELL(2) = CELL(2)*WVS
      CELL(3) = CELL(3)*WVS
      J = 0
      DO 30 I = 1 , 3
         IF (CELL(I).NE.0.) THEN
            J = J + 1
            CCELL(J) = CELL(I)
         END IF
 30   CONTINUE
      CALL ACALC
      CELL(1) = CELL(1)/WVS
      CELL(2) = CELL(2)/WVS
      CELL(3) = CELL(3)/WVS
C
C====== ICHANG denotes a change in a parameter which necessitates
C====== recalculating the reciprocal lattice coords.
C
      IF (ICHANG.NE.0) THEN
C
         DSTMAX = 1./DMIN
C
C====== Put the indices of possible reflections into /xyz/
C
         WRITE (IPRINT,6010)
         CALL HKLGEN(DSTMAX,IFLAG)
         WRITE (IPRINT,6030) NR
         WRITE (IPRINT,6020) NR
      END IF
C
C====== Execute the analysis commands
C
      IF (IENT.EQ.2) THEN
C
C====== Analyse the effect of stepping phi-y and phi-z
C
         CALL SEARCH
      ELSE IF (IENT.EQ.3) THEN
C
C====== Analyse the wavelength variations for symmetry equivs
C
         CALL WAVVAR
      ELSE IF (IENT.EQ.4) THEN
C
C====== Analyse the effect of changing crystal-film distance
C
         CALL VARYCF
      ELSE IF (IENT.EQ.5) THEN
C
C====== Analyse the fraction of reflections recorded against resolution
C
         CALL RESOLN
      ELSE IF (IENT.EQ.6) THEN
C
C====== Analyse the effect of rotating about an offset spindle
C
         CALL MUROT
      ELSE IF (IENT.EQ.7) THEN
C
C====== Analyse the effect of rotating about µ and U
C
         CALL TWOROT
      ELSE
C
C====== Analyse the effect of rotating about the spindle in this orientation
C
         CALL ROTSPI
      END IF
C
      CALL TPROMP('Type <RET> to continue...')
      CALL TREAD(CHAR,QUIT)
      IF (CHAR.EQ.'Q' .OR. CHAR.EQ.'q') STOP
      GO TO 20
C
C====== Format statements
C
 6010 FORMAT (
     +     ' Calculate and store the symmetry unique set of reflections'
     +     )
 6020 FORMAT (I6,' Reflections stored')
 6030 FORMAT (///I11)
      END
**==acalc.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       ACALC           **
C***************************
C
C
      SUBROUTINE ACALC
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: Calculate the master matrix, a, from the orientation matrix
C          premultiplied by the phi matrices. Code from IDXREF.
C
      COMMON /AMAT/A,AA
      COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN
      REAL CF,DMIN,SPNDLE,WVM,WVS,WVX
      INTEGER NUMA,NUMC
      REAL A(3,3),AA(3,3),CCELL(6),CELL(6),PHIS(3)
C
C====== LOCALS:
C
      REAL CPX,CPY,CPZ,RAPI,SPX,SPY,SPZ
      REAL AC1(3,3),AC2(3,3),AC3(3,3),AC4(3,3),AC5(3,3),AC6(3,3),C(3,3),
     +     CX(3,3),CY(3,3),CZ(3,3),PC(3,3)
C
C-------------------------------------------------------------------------------
C
C====== Pass the packed cell parameters in ccell to the mats routine
C====== of IDXREF. The AC arrays hold the (unnecessary) derivative matrices.
C
      CALL CLEAR(AA)
      CALL MATS(AA,AC1,AC2,AC3,AC4,AC5,AC6,CCELL)
C
C====== Define all rotation matrices
C
      RAPI = 3.1415926536/180.
      CPX = COS(PHIS(1)*RAPI)
      SPX = SIN(PHIS(1)*RAPI)
      CPY = COS(PHIS(2)*RAPI)
      SPY = SIN(PHIS(2)*RAPI)
      CPZ = COS(PHIS(3)*RAPI)
      SPZ = SIN(PHIS(3)*RAPI)
C
      CALL CLEAR(CX)
      CX(1,1) = 1.0
      CX(2,2) = CPX
      CX(2,3) = -SPX
      CX(3,2) = SPX
      CX(3,3) = CPX
C
      CALL CLEAR(CY)
      CY(1,1) = CPY
      CY(1,3) = SPY
      CY(2,2) = 1.0
      CY(3,1) = -SPY
      CY(3,3) = CPY
C
      CALL CLEAR(CZ)
      CZ(1,1) = CPZ
      CZ(1,2) = -SPZ
      CZ(2,1) = SPZ
      CZ(2,2) = CPZ
      CZ(3,3) = 1.0
C
C====== Matrix a is the transformation from indices to orthogonal
C====== coordinates
C
      CALL MATM33(CY,CX,C)
      CALL MATM33(CZ,C,PC)
      CALL MATM33(PC,AA,A)
C
      RETURN
      END
**==celltest.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       CELLTEST        **
C***************************
C
C
      SUBROUTINE CELLTEST(MSYS,ILGP,ICENT,CELL,LABEL,IOK)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: Check to see if the reciprocal cell parameters are consistent with
C          the lattice type. If inconsistent (eg length of zero in ortho) then
C          give warning message and return iok=1 so user can adjust them. If
C          ok then leave iok and return. To encourage users to understand the
C          principle, set redundant parameters to zero (eg b* and c* in cubic).
C          Note that program assumes b unique for monoclinic, c unique for hex
C          and tetragonal.
C
      INTEGER ICENT,ILGP,IOK,MSYS
      REAL CELL(6)
      CHARACTER LABEL(7,5)*6
C
      INCLUDE 'chanel.finc'
C
C====== LOCALS:
C
      INTEGER I
C
C-------------------------------------------------------------------------------
C
      IF (MSYS.EQ.2) THEN
C
C====== Monoclinic
C
         DO 20 I = 1 , 3
            IF (CELL(I).EQ.0.0) THEN
               WRITE (IPRINT,6010)
               GO TO 150
 
            END IF
 20      CONTINUE
         IF (CELL(5).EQ.0.0) THEN
            WRITE (IPRINT,6020)
            GO TO 150
 
         END IF
         DO 30 I = 4 , 6 , 2
            IF (CELL(I).NE.0.0) THEN
               IF (CELL(I).NE.90.) WRITE (IPRINT,6030)
               CELL(I) = 0.0
            END IF
 30      CONTINUE
         IF (ILGP.NE.1) THEN
            WRITE (IPRINT,6040)
            GO TO 150
 
         END IF
      ELSE IF (MSYS.EQ.3) THEN
C
C====== Orthorhombic
C
         DO 40 I = 1 , 3
            IF (CELL(I).EQ.0.0) THEN
               WRITE (IPRINT,6050)
               GO TO 150
 
            END IF
 40      CONTINUE
         DO 50 I = 4 , 6
            IF (CELL(I).NE.0.0) THEN
               IF (CELL(I).NE.90.) WRITE (IPRINT,6060)
               CELL(I) = 0.0
            END IF
 50      CONTINUE
         IF (ILGP.NE.2) THEN
            WRITE (IPRINT,6070)
            GO TO 150
 
         END IF
      ELSE IF (MSYS.EQ.4) THEN
C
C======== Tetragonal
C
         DO 60 I = 1 , 3 , 2
            IF (CELL(I).EQ.0.0) THEN
               WRITE (IPRINT,6080)
               GO TO 150
 
            END IF
 60      CONTINUE
         IF (CELL(2).NE.0.0) THEN
            IF (CELL(2).NE.CELL(1)) WRITE (IPRINT,6090)
            CELL(2) = 0.0
         END IF
         DO 70 I = 4 , 6
            IF (CELL(I).NE.0.0) THEN
               IF (CELL(I).NE.90.) WRITE (IPRINT,6100)
               CELL(I) = 0.0
            END IF
 70      CONTINUE
         IF (ILGP.NE.3 .AND. ILGP.NE.4) THEN
            WRITE (IPRINT,6110)
            GO TO 150
 
         END IF
      ELSE IF (MSYS.EQ.5) THEN
C
C======== Hexagonal
C
         DO 80 I = 1 , 3 , 2
            IF (CELL(I).EQ.0.0) THEN
               WRITE (IPRINT,6120)
               GO TO 150
 
            END IF
 80      CONTINUE
         IF (CELL(2).NE.0.0) THEN
            IF (CELL(2).NE.CELL(1)) WRITE (IPRINT,6130)
            CELL(2) = 0.0
         END IF
         DO 90 I = 4 , 5
            IF (CELL(I).NE.0.0) THEN
               IF (CELL(I).NE.90.) WRITE (IPRINT,6140)
               CELL(I) = 0.0
            END IF
 90      CONTINUE
         IF (CELL(6).NE.0.) THEN
            IF (CELL(6).NE.60.) WRITE (IPRINT,6130)
            CELL(6) = 0.
         END IF
 
         IF (ILGP.LT.5 .OR. ILGP.GT.9) THEN
            WRITE (IPRINT,6150)
            GO TO 150
 
         END IF
      ELSE IF (MSYS.EQ.6) THEN
C
C====== Rhombohedral
C
         IF (CELL(1).EQ.0.0) THEN
            WRITE (IPRINT,6160)
            GO TO 150
 
         END IF
         DO 100 I = 2 , 3
            IF (CELL(I).NE.0.0) THEN
               IF (CELL(I).NE.CELL(1)) WRITE (IPRINT,6170)
               CELL(I) = 0.0
            END IF
 
 100     CONTINUE
         IF (CELL(4).EQ.0.0) THEN
            WRITE (IPRINT,6180)
            GO TO 150
 
         END IF
         DO 110 I = 5 , 6
            IF (CELL(I).NE.0.0) THEN
               IF (CELL(I).NE.CELL(4)) WRITE (IPRINT,6190)
               CELL(I) = 0.0
            END IF
 110     CONTINUE
         IF (ILGP.LT.7 .OR. ILGP.GT.9) THEN
            WRITE (IPRINT,6200)
            GO TO 150
 
         END IF
      ELSE IF (MSYS.EQ.7) THEN
C
C====== Cubic
C
         IF (CELL(1).EQ.0.0) THEN
            WRITE (IPRINT,6210)
            GO TO 150
 
         END IF
         DO 120 I = 2 , 3
            IF (CELL(I).NE.0.0) THEN
               IF (CELL(I).NE.CELL(1)) WRITE (IPRINT,6220)
               CELL(I) = 0.0
            END IF
 
 120     CONTINUE
         DO 130 I = 4 , 6
            IF (CELL(I).NE.0.0) THEN
               IF (CELL(I).NE.90.) WRITE (IPRINT,6230)
               CELL(I) = 0.0
            END IF
 130     CONTINUE
         IF (ILGP.NE.10 .AND. ILGP.NE.11) THEN
            WRITE (IPRINT,6240)
            GO TO 150
 
         END IF
      ELSE
C
C====== Triclinic
C
         DO 140 I = 1 , 6
            IF (CELL(I).EQ.0.0) THEN
               WRITE (IPRINT,6250)
               GO TO 150
 
            END IF
 140     CONTINUE
         IF (ILGP.NE.0) THEN
            WRITE (IPRINT,6260)
            GO TO 150
 
         END IF
      END IF
C
      IOK = 0
      RETURN
C
 150  IOK = 1
      RETURN
C
C====== Format statements
C
 6010 FORMAT (' Zero cell length in monoclinic case')
 6020 FORMAT (' Zero beta-star in monoclinic case')
 6030 FORMAT (' Monoclinic case must have unique axis b')
 6040 FORMAT (' Laue group should be 2/m for monoclinic')
 6050 FORMAT (' Zero cell length in orthorhombic case')
 6060 FORMAT (' No angles needed for orthorhombic case')
 6070 FORMAT (' Laue group should be mmm for orthorhombic')
 6080 FORMAT (' Zero cell length in tetragonal case')
 6090 FORMAT (' C must be unique axis for tetragonal case')
 6100 FORMAT (' No angles needed for tetragonal case')
 6110 FORMAT (' Laue group should be 4/m or 4/mmm for tetragonal')
 6120 FORMAT (' Zero cell length in hexagonal case')
 6130 FORMAT (' C must be unique axis for hexagonal case')
 6140 FORMAT (' No angles needed for hexagonal case')
 6150 FORMAT (' Laue group should be 6/m or 6/mmm for hexagonal'/
     +        ' 3-bar, 3m1 or 31m for trigonal on hexagonal axes')
 6160 FORMAT (' Zero cell length in rhombohedral case')
 6170 FORMAT (' Only one cell length needed for rhombohedral case')
 6180 FORMAT (' Zero cell angle in rhombohedral case')
 6190 FORMAT (' Only one angle needed for rhombohedral case')
 6200 FORMAT (
     +' Laue group should be 3-bar, 3m1 or 31m for trigonal'' on rhomboh
     +edral axes')
 6210 FORMAT (' Zero cell length in cubic case')
 6220 FORMAT (' Only one cell length needed for cubic case')
 6230 FORMAT (' No angles needed for cubic case')
 6240 FORMAT (' Laue group should be m3 or m3m for cubic')
 6250 FORMAT (' Zero cell parameter in triclinic case')
 6260 FORMAT (' Laue group should be 1-bar for triclinic')
      END
**==clear.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       CLEAR           **
C***************************
C
C
      SUBROUTINE CLEAR(A)
C
      INCLUDE 'implicit_none.finc'
C
      REAL A(9)
      INTEGER I
C
C-------------------------------------------------------------------------------
C
      DO 20 I = 1 , 9
         A(I) = 0.0
 20   CONTINUE
      RETURN
      END
**==gensym.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       GENSYM          **
C***************************
C
C
      SUBROUTINE GENSYM(NSYM,NREF,IH,IK,IL)
C
      INCLUDE 'implicit_none.finc'
C
      COMMON /MATRIX/MSYS,MROT,MBEAM,ILGP,ICENT
      INTEGER IH,IK,IL,NSYM
      INTEGER NREF(3,48)
      INTEGER ICENT,ILGP,MBEAM,MROT,MSYS
C
C====== LOCALS:
C
      INTEGER I,J
C
C-------------------------------------------------------------------------------
C
      NREF(1,1) = IH
      NREF(2,1) = IK
      NREF(3,1) = IL
      NSYM = 1
      IF (ILGP.NE.0) THEN
C
C====== Deal with symmetries higher than triclinic
C
         IF (ILGP.EQ.3) GO TO 30
         IF (ILGP.EQ.4) GO TO 30
         IF (ILGP.EQ.5) GO TO 60
         IF (ILGP.EQ.6) GO TO 60
         IF (ILGP.EQ.7) GO TO 90
         IF (ILGP.EQ.8) GO TO 90
         IF (ILGP.EQ.9) GO TO 90
         IF (ILGP.EQ.10) GO TO 120
         IF (ILGP.EQ.11) GO TO 120
C
C====== 2/m
C
         IF (IK.NE.0) THEN
            IF (IL.NE.0 .OR. IH.NE.0) THEN
               NREF(1,2) = IH
               NREF(2,2) = -IK
               NREF(3,2) = IL
               NSYM = 2
            END IF
         END IF
         IF (ILGP.NE.1) THEN
C
C====== mmm
C
            IF (IH.NE.0 .AND. IL.NE.0) THEN
               DO 20 I = 1 , NSYM
                  NREF(1,I+NSYM) = NREF(1,I)
                  NREF(2,I+NSYM) = NREF(2,I)
                  NREF(3,I+NSYM) = -NREF(3,I)
 20            CONTINUE
               NSYM = 2*NSYM
            END IF
         END IF
      END IF
      GO TO 160
C
C====== 4/m
C
 30   IF (IH.NE.0 .OR. IK.NE.0) THEN
         NREF(1,2) = IK
         NREF(2,2) = -IH
         NREF(3,2) = IL
         NSYM = 2
         IF (IL.NE.0) THEN
            DO 40 I = 1 , NSYM
               NREF(1,I+NSYM) = NREF(1,I)
               NREF(2,I+NSYM) = NREF(2,I)
               NREF(3,I+NSYM) = -NREF(3,I)
 40         CONTINUE
            NSYM = NSYM*2
         END IF
      END IF
      IF (ILGP.NE.3) THEN
C
C====== 4/mmm
C
         IF (IK.NE.0 .AND. IH.NE.IK) THEN
            DO 50 I = 1 , NSYM
               NREF(1,I+NSYM) = NREF(2,I)
               NREF(2,I+NSYM) = NREF(1,I)
               NREF(3,I+NSYM) = NREF(3,I)
 50         CONTINUE
            NSYM = NSYM*2
         END IF
      END IF
      GO TO 160
C
C====== 6/m
C
 60   IF (IH.NE.0 .OR. IK.NE.0) THEN
         NREF(1,2) = -(IH+IK)
         NREF(2,2) = IH
         NREF(3,2) = IL
         NREF(1,3) = IK
         NREF(2,3) = -(IH+IK)
         NREF(3,3) = IL
         NSYM = 3
         IF (IL.NE.0) THEN
            DO 70 I = 1 , NSYM
               NREF(1,I+NSYM) = NREF(1,I)
               NREF(2,I+NSYM) = NREF(2,I)
               NREF(3,I+NSYM) = -NREF(3,I)
 70         CONTINUE
            NSYM = NSYM*2
         END IF
      END IF
      IF (ILGP.NE.5) THEN
C
C====== 6/mmm
C
         IF (IK.NE.0 .AND. IH.NE.IK) THEN
            DO 80 I = 1 , NSYM
               NREF(1,I+NSYM) = NREF(2,I)
               NREF(2,I+NSYM) = NREF(1,I)
               NREF(3,I+NSYM) = NREF(3,I)
 80         CONTINUE
            NSYM = NSYM*2
         END IF
      END IF
      GO TO 160
C
C====== 3bar (hexagonal axes)
C
 90   IF (MSYS.NE.5) THEN
C
C====== 3bar (rhombohedral axes)
C
         IF (MSYS.EQ.6) THEN
            IF (IH.NE.IK .OR. IH.NE.IL) THEN
               NREF(1,2) = NREF(2,1)
               NREF(2,2) = NREF(3,1)
               NREF(3,2) = NREF(1,1)
               NREF(1,3) = NREF(3,1)
               NREF(2,3) = NREF(1,1)
               NREF(3,3) = NREF(2,1)
               NSYM = 3
               IF (ILGP.NE.7) THEN
C
C====== 3m1 and 31m (rhombohedral axes)
C
                  IF (IH.NE.IK) THEN
                     DO 100 I = 1 , NSYM
                        NREF(1,I+NSYM) = NREF(2,I)
                        NREF(2,I+NSYM) = NREF(1,I)
                        NREF(3,I+NSYM) = NREF(3,I)
 100                 CONTINUE
                     NSYM = NSYM*2
                  END IF
               END IF
            END IF
         END IF
         GO TO 160
      ELSE
         IF (IH.NE.0 .OR. IK.NE.0) THEN
            NREF(1,2) = -(IH+IK)
            NREF(2,2) = IH
            NREF(3,2) = IL
            NREF(1,3) = IK
            NREF(2,3) = -(IH+IK)
            NREF(3,3) = IL
            NSYM = 3
         END IF
         IF (ILGP.EQ.7) GO TO 160
C
C====== 3m1 and 31m
C
         IF (IL.EQ.0) THEN
            IF (IH.EQ.IK) GO TO 160
            IF (IK.EQ.0) GO TO 160
 
         ELSE IF (ILGP.EQ.8) THEN
            IF (IK.EQ.0) GO TO 160
 
         ELSE IF (ILGP.EQ.9) THEN
            IF (IH.EQ.IK) GO TO 160
         END IF
      END IF
 
      DO 110 I = 1 , NSYM
         NREF(1,I+NSYM) = NREF(2,I)
         NREF(2,I+NSYM) = NREF(1,I)
         IF (ILGP.EQ.8) THEN
            NREF(3,I+NSYM) = -NREF(3,I)
         ELSE
            NREF(3,I+NSYM) = NREF(3,I)
         END IF
 110  CONTINUE
      NSYM = NSYM*2
      GO TO 160
C
C====== m3
C
 120  IF (IK.NE.0 .OR. IL.NE.0) THEN
         NREF(1,2) = -IH
         NREF(2,2) = IK
         NREF(3,2) = IL
         NSYM = 2
         IF (IK.NE.0) THEN
            DO 130 I = 1 , NSYM
               NREF(1,I+NSYM) = NREF(1,I)
               NREF(2,I+NSYM) = -NREF(2,I)
               NREF(3,I+NSYM) = NREF(3,I)
 130        CONTINUE
            NSYM = NSYM*2
         END IF
      END IF
      IF (IH.NE.IK .OR. IH.NE.IL) THEN
         DO 140 I = 1 , NSYM
            NREF(1,I+NSYM) = NREF(2,I)
            NREF(2,I+NSYM) = NREF(3,I)
            NREF(3,I+NSYM) = NREF(1,I)
            NREF(1,I+2*NSYM) = NREF(3,I)
            NREF(2,I+2*NSYM) = NREF(1,I)
            NREF(3,I+2*NSYM) = NREF(2,I)
 140     CONTINUE
         NSYM = NSYM*3
      END IF
      IF (ILGP.NE.10) THEN
C
C====== m3m
C
         IF (IH.NE.IK) THEN
            DO 150 I = 1 , NSYM
               NREF(1,I+NSYM) = NREF(2,I)
               NREF(2,I+NSYM) = NREF(1,I)
               NREF(3,I+NSYM) = NREF(3,I)
 150        CONTINUE
            NSYM = NSYM*2
         END IF
      END IF
C
 160  DO 180 I = 1 , NSYM
         DO 170 J = 1 , 3
            NREF(J,I+NSYM) = -NREF(J,I)
 170     CONTINUE
 180  CONTINUE
      NSYM = 2*NSYM
C
      RETURN
      END
**==hklgen.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       HKLGEN          **
C***************************
C
C
      SUBROUTINE HKLGEN(DSTMAX,IFLAG)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: Find all reflections within dstar-max in the hemisphere +h, all k,l.
C       In essence the routine flow is to
C       (1) increase h from 0 to max
C       (2) increase k from 0 to max, then
C           decrease k from -1 to min
C       (3) increase l from 0 to max, then
C           decrease l from -1 to min
C  with h varying least, l most rapidly. The max and min values are
C  determined by a dstar which exceeds the limit. Complications are
C  caused by non-orthogonal cells: for given h,k l=0 may exceed dstar-
C  max, but some other l may not. The logic is that if at least one l
C  for given h,k is found within dstar-max, then the l loops end when
C  dstar exceeds dstar-max. If such an l is not found then the loops end
C  in both directions when the calculated dstar is increasing. Similarly
C  for the control of k loops, where the test is for at least one k,l
C  inside for given h.
C
      COMMON /AMAT/A,DUMMY
      COMMON /SPOTS/NR
      COMMON /SYS/ISYS,KSYS
      COMMON /XYZ/IHKL
C
      INCLUDE 'chanel.finc'
C
      REAL DSTMAX
      INTEGER IFLAG
      INTEGER ISYS,NR
      REAL A(3,3),DUMMY(3,3)
      INTEGER IHKL(20000),KSYS(3)
C
C====== LOCALS:
C
      REAL DST2,DST2KM,DST2LM,DSTEST,XG,YG,ZG
      INTEGER H,K,KSTART,L,LIMIT,LSTART
      LOGICAL ABSNT,KIN,KPOS,LIN,LPOS,UNIQ
      CHARACTER IYN*1
      INTEGER IPACK
C
C-------------------------------------------------------------------------------
C
C====== Initialisations
C
      LIMIT = 20000
      DSTEST = DSTMAX*DSTMAX
      NR = 1
      KSTART = 0
      LSTART = 0
C
C====== Start h loop
C
      H = 0
C
C====== Start k loop
C
 20   K = KSTART
      KPOS = .TRUE.
      KIN = .FALSE.
      DST2KM = 1000.
C
C====== Check for this h,k systematically absent, whatever the l value
C
 30   IF (ISYS.NE.0) THEN
         CALL SYSAHK(H,K,ABSNT)
         IF (ABSNT) THEN
            IF (KPOS) THEN
               K = K + 1
 
            ELSE
               K = K - 1
 
            END IF
            GO TO 30
 
         END IF
 
      END IF
C
C====== Start l loop - skip 000
C
      L = LSTART
      IF (H.EQ.0 .AND. K.EQ.0) THEN
         L = 1
      END IF
 
      LPOS = .TRUE.
      LIN = .FALSE.
      DST2LM = 1000.
C
C====== Check for this h,k,l systematically absent
C
 40   IF (ISYS.NE.0) THEN
         CALL SYSABS(H,K,L,ABSNT)
         IF (ABSNT) THEN
            IF (LPOS) THEN
               L = L + 1
 
            ELSE
               L = L - 1
 
            END IF
            GO TO 40
 
         END IF
 
      END IF
C
C====== Calculate dstar squared.
C
      XG = A(1,1)*H + A(1,2)*K + A(1,3)*L
      YG = A(2,1)*H + A(2,2)*K + A(2,3)*L
      ZG = A(3,1)*H + A(3,2)*K + A(3,3)*L
      DST2 = XG**2 + YG**2 + ZG**2
C
C====== Test for within resolution
C
      IF (DST2.LE.DSTEST) THEN
         LIN = .TRUE.
         KIN = .TRUE.
C
C====== Check for this h,k,l in the unique set according to Laue symmetry
C
         CALL UNIQUE(H,K,L,UNIQ)
         IF (UNIQ) THEN
            IHKL(NR) = IPACK(H,K,L)
            NR = NR + 1
            IF (NR.GT.LIMIT) THEN
C
C====== Ran out of array space
C
               WRITE (IPRINT,6010) LIMIT , H , K , L
               WRITE (IUNIT6,6010) LIMIT , H , K , L
               NR = LIMIT
               RETURN
            END IF
         END IF
C
C====== Increment l and continue
C
         IF (LPOS) THEN
            L = L + 1
 
         ELSE
            L = L - 1
 
         END IF
         GO TO 40
C
C====== hkl outside resolution ... test if any hits for this l so far
C
      ELSE IF (LIN) THEN
C
C====== Yes ... reverse l direction
C
         IF (LPOS) THEN
            LPOS = .FALSE.
            L = LSTART - 1
            GO TO 40
C
C====== or increment k if already reversed l
C
         ELSE IF (KPOS) THEN
            K = K + 1
 
         ELSE
            K = K - 1
 
         END IF
         GO TO 30
C
C====== No  ... check if getting closer ... save and continue if so
C
      ELSE IF (DST2.LT.DST2LM) THEN
         DST2LM = DST2
         IF (LPOS) THEN
            L = L + 1
 
         ELSE
            L = L - 1
 
         END IF
         GO TO 40
C
C====== Getting further away ... reverse l and continue
C
      ELSE IF (LPOS) THEN
         LPOS = .FALSE.
         L = LSTART - 1
         GO TO 40
C
C====== No hits for this l so far ... test if any for this k so far
C
      ELSE IF (KIN) THEN
C
C====== Yes ... reverse k direction
C
         IF (KPOS) THEN
            KPOS = .FALSE.
            K = KSTART - 1
            GO TO 30
C
C====== or increment h if already reversed k
C
         ELSE
            H = H + 1
 
         END IF
         GO TO 20
C
C====== No  ... check if getting closer ... save and continue if so
C
      ELSE IF (DST2.LT.DST2KM) THEN
         DST2KM = DST2
         IF (KPOS) THEN
            K = K + 1
 
         ELSE
            K = K - 1
 
         END IF
         GO TO 30
C
C====== Getting further away ... reverse k and continue
C
      ELSE IF (KPOS) THEN
         KPOS = .FALSE.
         K = KSTART - 1
         GO TO 30
C
C====== No hits for this h value so stop, returning iflag=0 to mark end
C
      ELSE IF (H.EQ.0) THEN
         H = H + 1
         GO TO 20
 
      ELSE
         IFLAG = 0
         NR = NR - 1
         RETURN
 
      END IF
C
C====== Format statements
C
 6010 FORMAT (' limit of',I6,' reflections reached at h,k,l=',
     +        3I4/' will have to truncate calculation here')
      END
**==inputv.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       INPUTV          **
C***************************
C
C
      SUBROUTINE INPUTV(PROMPT,VAR,VALUES,NVALS)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose : Extract a value from the VALUES array returned by MENU, or
C           if empty prompt user for the value.
C
      CHARACTER*(*) PROMPT
      REAL VAR, VALUES(*)
      INTEGER NVALS
C
      INCLUDE 'chanel.finc'
C
C====== FUNCTION TYPES:
C
      REAL FLPVAL
C
C====== LOCALS:
C
      LOGICAL QUIT
      CHARACTER*30 TMPSTR
      LOGICAL ERR
C
C-------------------------------------------------------------------------------
C
      IF (NVALS.GT.0) THEN
         VAR = VALUES(1)
      ELSE
         CALL TPROMP(PROMPT)
         CALL TREAD(TMPSTR,QUIT)
         VAR = FLPVAL(TMPSTR,ERR)
      END IF
      END
**==inputw.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       INPUTW          **
C***************************
C
C
      SUBROUTINE INPUTW(PROMPT,WORD,WORDS,NWORDS)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose : Extract a word from the WORDS array returned by MENU, or
C           if empty prompt user for the word.
C
      CHARACTER*(*) PROMPT
      CHARACTER*4 WORD
      CHARACTER*4 WORDS(*)
      INTEGER NWORDS
C
      INCLUDE 'chanel.finc'
C
C====== LOCALS:
C
      LOGICAL QUIT
C
C-------------------------------------------------------------------------------
C
      IF (NWORDS.GT.0) THEN
         WORD = WORDS(1)
      ELSE
         CALL TPROMP(PROMPT)
         CALL TREAD(WORD,QUIT)
      END IF
      CALL CCPUPC(WORD)
      RETURN
      END
**==ipack.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       IPACK           **
C***************************
C
C
      INTEGER FUNCTION IPACK(H,K,L)
C
      INTEGER H,K,L
C
C-------------------------------------------------------------------------------
C
      IPACK = 40000*(H+100) + 200*(K+100) + (L+100)
      RETURN
      END
**==item.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       ITEM            **
C***************************
C
C
      SUBROUTINE ITEM(IOPT,ARRAY,LEN,VIDEO)
 
      INCLUDE 'implicit_none.finc'
C
C Purpose : Displays character string on VDU at given line and column
C
      CHARACTER * (*)  ARRAY
      INTEGER IOPT,LEN
      LOGICAL VIDEO
C
C====== LOCALS:
C
      INTEGER LINE,ICOL
C
C-------------------------------------------------------------------------------
C
      ICOL = (IOPT+16)/17
      LINE = IOPT - 17*(ICOL-1)
      CALL QSSETC(1,24)
      ICOL = ICOL*41 - 40
      LINE = LINE + 1
      CALL QSSETC(ICOL,LINE)
      IF (VIDEO) CALL QSREVV
      CALL TSTRNG(ARRAY(1:LEN))
      CALL QSNORV
      RETURN
      END
**==laue.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       LAUE            **
C***************************
C
C
      SUBROUTINE LAUE(ION)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: Calculate but do not display the Laue photograph.
C          The array ion returned by this routine shows which of the nr
C          reflections in ihkl would appear on the film plot (ion(i)=n>0)
C          or not (ion(i)=0).
C          (n) is the total number of symmetry equivalents for this reflection.
C
      COMMON /AMAT/A,DUMMY
      COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN
      COMMON /RADC/RADI,RADI2
      COMMON /SPOTS/NR
      COMMON /XYZ/IHKL
      INTEGER ION(*)
      REAL CF,DMIN,RADI,RADI2,SPNDLE,WVM,WVS,WVX
      INTEGER NR,NUMA,NUMC
      REAL A(3,3),CCELL(6),CELL(6),DUMMY(3,3),PHIS(3)
      INTEGER IHKL(20000)
C
C====== LOCALS:
C
      REAL C1,DD1,DD2,RF,RMAX,RMAX2,RMIN,RMIN2,S1,X,XF,XG,Y,YF,YG,Z,ZG
      INTEGER I,IH,IK,IL,N,NSYM
      INTEGER NREF(3,48)
C
C-------------------------------------------------------------------------------
C
 
      C1 = SPNDLE*3.1415926536/180.0
      S1 = SIN(C1)
      C1 = COS(C1)
      RMIN = 1./WVX
      RMAX = 1./WVM
      RMIN2 = RMIN*RMIN
      RMAX2 = RMAX*RMAX
C
C====== Loop through all the reflections
C
      DO 30 N = 1 , NR
         ION(N) = 0
C
C====== Unpack reflection indices
C
         CALL UNPACK(IHKL(N),IH,IK,IL)
C
C====== Generate symmetry equivalent reflections
C
         CALL GENSYM(NSYM,NREF,IH,IK,IL)
C
C====== Loop through these symmetry equivalents
C
         DO 20 I = 1 , NSYM
            IH = NREF(1,I)
            IK = NREF(2,I)
            IL = NREF(3,I)
            XG = A(1,1)*IH + A(1,2)*IK + A(1,3)*IL
            YG = A(2,1)*IH + A(2,2)*IK + A(2,3)*IL
            ZG = A(3,1)*IH + A(3,2)*IK + A(3,3)*IL
            X = C1*XG - S1*YG
            IF (X.LT.0.00001) THEN
               Y = S1*XG + C1*YG
               Z = ZG
C
C====== Dists squared of rl point from centre min and max spheres
C
               DD1 = ((RMIN+X)**2+Y*Y+Z*Z)
               IF (DD1.GE.RMIN2) THEN
                  DD2 = ((RMAX+X)**2+Y*Y+Z*Z)
                  IF (DD2.LE.RMAX2) THEN
                     RF = (X*X+Y*Y+Z*Z)/(2*X)
C
C====== rl point can reflect, centre refl sphere at RF,0,0
C
                     IF (CF*(-RF+X).GT.0.) THEN
                        XF = CF*Y/(-RF+X)
                        YF = CF*Z/(-RF+X)
C
C====== Check that point is inside the circle requested by the user
C
                        IF (ABS(XF).LE.RADI) THEN
                           IF (ABS(YF).LE.RADI) THEN
                              IF ((XF*XF+YF*YF).LE.RADI2) THEN
                                 ION(N) = ION(N) + 1
                              END IF
                           END IF
                        END IF
                     END IF
                  END IF
               END IF
            END IF
 20      CONTINUE
 30   CONTINUE
      END
**==lauer.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       LAUER           **
C***************************
C
C
      SUBROUTINE LAUER(ION,WRES)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: Calculate but do not display the Laue photograph.
C          The array ion returned by this routine shows which of the nr
C          reflections in ihkl would appear on the film plot (n) or not (0).
C          (n) is the total number of symmetry equivalents for this reflection.
C          The array wres records the resolution for each.
C
      COMMON /AMAT/A,DUMMY
      COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN
      COMMON /RADC/RADI,RADI2
      COMMON /SPOTS/NR
      COMMON /XYZ/IHKL
      REAL WRES(*)
      INTEGER ION(*)
      REAL CF,DMIN,RADI,RADI2,SPNDLE,WVM,WVS,WVX
      INTEGER NR,NUMA,NUMC
      REAL A(3,3),CCELL(6),CELL(6),DUMMY(3,3),PHIS(3)
      INTEGER IHKL(20000)
C
C====== LOCALS:
C
      REAL C1,DD1,DD2,RF,RMAX,RMAX2,RMIN,RMIN2,S1,X,XF,XG,Y,YF,YG,Z,ZG
      INTEGER I,IH,IK,IL,N,NSYM
      INTEGER NREF(3,48)
C
C-------------------------------------------------------------------------------
C
 
      C1 = SPNDLE*3.1415926536/180.0
      S1 = SIN(C1)
      C1 = COS(C1)
      RMIN = 1./WVX
      RMAX = 1./WVM
      RMIN2 = RMIN*RMIN
      RMAX2 = RMAX*RMAX
C
C====== Loop through all the reflections
C
      DO 30 N = 1 , NR
         ION(N) = 0
C
C====== Unpack reflection indices
C
         CALL UNPACK(IHKL(N),IH,IK,IL)
C
C====== Generate symmetry equivalent reflections
C
         CALL GENSYM(NSYM,NREF,IH,IK,IL)
C
C====== Loop through these symmetry equivalents
C
         DO 20 I = 1 , NSYM
            IH = NREF(1,I)
            IK = NREF(2,I)
            IL = NREF(3,I)
            XG = A(1,1)*IH + A(1,2)*IK + A(1,3)*IL
            YG = A(2,1)*IH + A(2,2)*IK + A(2,3)*IL
            ZG = A(3,1)*IH + A(3,2)*IK + A(3,3)*IL
            X = C1*XG - S1*YG
            IF (X.LT.0.00001) THEN
               Y = S1*XG + C1*YG
               Z = ZG
C
C====== Dists squared of rl point from centre min and max spheres
C
               DD1 = ((RMIN+X)**2+Y*Y+Z*Z)
               IF (DD1.GE.RMIN2) THEN
                  DD2 = ((RMAX+X)**2+Y*Y+Z*Z)
                  IF (DD2.LE.RMAX2) THEN
                     RF = (X*X+Y*Y+Z*Z)/(2*X)
C
C====== rl point can reflect, centre refl sphere at RF,0,0
C
                     IF (CF*(-RF+X).GT.0.) THEN
                        XF = CF*Y/(-RF+X)
                        YF = CF*Z/(-RF+X)
C
C====== Check that point is inside the circle requested by the user
C
                        IF (ABS(XF).LE.RADI) THEN
                           IF (ABS(YF).LE.RADI) THEN
                              IF ((XF*XF+YF*YF).LE.RADI2) THEN
                                 ION(N) = ION(N) + 1
                                 WRES(N) = 1./SQRT(XG**2+YG**2+ZG**2)
                                 GO TO 30
                              END IF
                           END IF
                        END IF
                     END IF
                  END IF
               END IF
            END IF
 
 20      CONTINUE
 30   CONTINUE
      END
**==lauew.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       LAUEW           **
C***************************
C
C
      SUBROUTINE LAUEW(ION,WION)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: Calculate but do not display the Laue photograph.
C          The array ion returned by this routine shows which of the nr
C          reflections in ihkl would appear on the film plot (n) or not (0).
C          (n) is the total number of symmetry equivalents for this reflection.
C          The array wion records the variation in wavelength for each.
C
      COMMON /AMAT/A,DUMMY
      COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN
      COMMON /RADC/RADI,RADI2
      COMMON /SPOTS/NR
      COMMON /XYZ/IHKL
      REAL WION(*)
      INTEGER ION(*)
      REAL CF,DMIN,RADI,RADI2,SPNDLE,WVM,WVS,WVX
      INTEGER NR,NUMA,NUMC
      REAL A(3,3),CCELL(6),CELL(6),DUMMY(3,3),PHIS(3)
      INTEGER IHKL(20000)
C
C====== LOCALS:
C
      REAL C1,DD1,DD2,DIFF,RF,RMAX,RMAX2,RMIN,RMIN2,S1,SUM,X,XF,XG,Y,YF,
     +     YG,Z,ZG
      INTEGER I,IH,IK,IL,J,M,N,NSYM
      REAL WAVE(48)
      INTEGER NREF(3,48)
C
C-------------------------------------------------------------------------------
C
 
      C1 = SPNDLE*3.1415926536/180.0
      S1 = SIN(C1)
      C1 = COS(C1)
      RMIN = 1./WVX
      RMAX = 1./WVM
      RMIN2 = RMIN*RMIN
      RMAX2 = RMAX*RMAX
C
C====== Loop through all the reflections
C
      DO 50 N = 1 , NR
         ION(N) = 0
C
C====== Unpack reflection indices
C
         CALL UNPACK(IHKL(N),IH,IK,IL)
C
C====== Generate symmetry equivalent reflections
C
         CALL GENSYM(NSYM,NREF,IH,IK,IL)
C
C====== Loop through these symmetry equivalents
C
         M = 0
         DO 20 I = 1 , NSYM
            IH = NREF(1,I)
            IK = NREF(2,I)
            IL = NREF(3,I)
            XG = A(1,1)*IH + A(1,2)*IK + A(1,3)*IL
            YG = A(2,1)*IH + A(2,2)*IK + A(2,3)*IL
            ZG = A(3,1)*IH + A(3,2)*IK + A(3,3)*IL
            X = C1*XG - S1*YG
            IF (X.LT.0.00001) THEN
               Y = S1*XG + C1*YG
               Z = ZG
C
C====== Dists squared of rl point from centre min and max spheres
C
               DD1 = ((RMIN+X)**2+Y*Y+Z*Z)
               IF (DD1.GE.RMIN2) THEN
                  DD2 = ((RMAX+X)**2+Y*Y+Z*Z)
                  IF (DD2.LE.RMAX2) THEN
                     RF = (X*X+Y*Y+Z*Z)/(2*X)
C
C====== rl point can reflect, centre refl sphere at rf,0,0
C
                     IF (CF*(-RF+X).GT.0.) THEN
                        XF = CF*Y/(-RF+X)
                        YF = CF*Z/(-RF+X)
C
C====== Check that point is inside the circle requested by the user
C
                        IF (ABS(XF).LE.RADI) THEN
                           IF (ABS(YF).LE.RADI) THEN
                              IF ((XF*XF+YF*YF).LE.RADI2) THEN
                                 ION(N) = ION(N) + 1
                                 M = M + 1
                                 WAVE(M) = -1./RF
                              END IF
                           END IF
                        END IF
                     END IF
                  END IF
               END IF
            END IF
 20      CONTINUE
C
C====== Calculate wave length variation
C
         IF (M.GT.1) THEN
            SUM = 0.
            DO 30 J = 1 , M
               SUM = SUM + WAVE(J)
 30         CONTINUE
            SUM = SUM/M
            DIFF = 0.
            DO 40 J = 1 , M
               DIFF = ABS(WAVE(J)-SUM) + DIFF
 40         CONTINUE
            WION(N) = DIFF/M
         END IF
 50   CONTINUE
      END
**==lgrp.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       LGRP            **
C***************************
C
C
      SUBROUTINE LGRP(MSYS,ILGP)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: Check the crystal system in msys and ask the user for the Laue
C          group (returned in ilgp) if necessary.
C
      INTEGER ILGP,MSYS
C
      INCLUDE 'chanel.finc'
C
C====== LOCALS:
C
      LOGICAL QUIT
      INTEGER ITEMP
C
C-------------------------------------------------------------------------------
C
 
      IF (MSYS.EQ.2) THEN
C
C====== Monoclinic
C
         ILGP = 1
         RETURN
      ELSE IF (MSYS.EQ.3) THEN
C
C====== Orthorhombic
C
         ILGP = 2
         RETURN
      ELSE IF (MSYS.EQ.4) THEN
C
C====== Tetragonal - may be 4/m or 4/mmm
C
 20      CALL TWRITE(
     +              'Select Laue group for symmetry equivalent analysis'
     +              )
         ITEMP = 1
         CALL TASKI('1 = 4/m, 2 = 4/mmm',ITEMP,QUIT)
         IF (ITEMP.NE.1 .AND. ITEMP.NE.2) GO TO 20
         ILGP = 2 + ITEMP
         RETURN
      ELSE IF (MSYS.EQ.5) THEN
C
C====== Hexagonal - may be 6/m or 6/mmm
C====== or 3bar, 3m1 or 31m  for trigonal indexed on hexagonal lattice
C
 30      CALL TWRITE(
     +              'Select Laue group for symmetry equivalent analysis'
     +              )
         ITEMP = 1
         CALL TASKI('1 = 6/m, 2 = 6/mmm, 3 = 3bar, 4 = 3m1, 5 = 31m',
     +              ITEMP,QUIT)
         IF (ITEMP.LT.1 .AND. ITEMP.GT.5) GO TO 30
         ILGP = 4 + ITEMP
         RETURN
      ELSE IF (MSYS.EQ.6) THEN
C
C====== Rhombohedral lattice is trigonal - may be 3bar, 3m1 or 31m
C
 40      CALL TWRITE(
     +              'Select Laue group for symmetry equivalent analysis'
     +              )
         ITEMP = 1
         CALL TASKI('1 = 3bar, 2 = 3m1, 3 = 31m',ITEMP,QUIT)
         IF (ITEMP.LT.1 .OR. ITEMP.GT.3) GO TO 40
         ILGP = 6 + ITEMP
         RETURN
      ELSE IF (MSYS.EQ.7) THEN
C
C====== Cubic
C
 50      CALL TWRITE(
     +              'Select Laue group for symmetry equivalent analysis'
     +              )
         ITEMP = 1
         CALL TASKI('1 = m3, 2 = m3m',ITEMP,QUIT)
         IF (ITEMP.NE.1 .AND. ITEMP.NE.2) GO TO 50
         ILGP = 9 + ITEMP
         RETURN
      ELSE
C
C====== Triclinic
C
         ILGP = 0
         RETURN
      END IF
      END
**==matm33.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       MATM33          **
C***************************
C
C
      SUBROUTINE MATM33(A,B,C)
C
      INCLUDE 'implicit_none.finc'
C
      REAL A(3,3),B(3,3),C(3,3)
C
C====== LOCALS:
C
      INTEGER I,J,K
C
C-------------------------------------------------------------------------------
C
      DO 40 I = 1 , 3
         DO 30 J = 1 , 3
            C(I,J) = 0.0
            DO 20 K = 1 , 3
               C(I,J) = C(I,J) + A(I,K)*B(K,J)
 20         CONTINUE
 30      CONTINUE
 40   CONTINUE
      RETURN
      END
**==mats.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       MATS            **
C***************************
C
C
      SUBROUTINE MATS(AM,AA,AB,AC,AD,AE,AF,CELL)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: Generalised matrix and derivative builder
C          the reciprocal cell parameters come in cell(1-6) but only the
C          unique ones ... cell(1) is astar, cell(2) alphastar for rhombohedral
C          for example.
C          am returns the A matrix ... row 1 is components of reciprocal cell
C          along X, row 2 along Y, row 3 along Z.
C          aa, ab etc are the derivative matrices for A wrt the unique cell
C          parameters in order ... thus aa is dA/dastar, ab dA/alphastar for
C          rhombohedral.
C
      COMMON /MATRIX/MSYS,MROT,MBEAM,IJUNK1,IJUNK2
      REAL AA(3,3),AB(3,3),AC(3,3),AD(3,3),AE(3,3),AF(3,3),AM(3,3),
     +     CELL(6)
      INTEGER IJUNK1,IJUNK2,MBEAM,MROT,MSYS
C
C====== LOCALS:
C
      REAL A,C,COSPSI,COSRHO,FUN,S,SI,SJ,SK
      INTEGER I,J,K,N1,N2,N3
      REAL CC(3),DA(3,3,3),DB(3,3,3),RL(3),SS(3)
C
C-------------------------------------------------------------------------------
C
C
C====== I refers to mrot, the crystal rotation axis, Z
C
      I = IABS(MROT)
      SI = 1.0
      IF (MROT.LT.0) SI = -1.0
C
C====== J refers to mbeam, the xray beam axis, X
C
      J = IABS(MBEAM)
      SJ = 1.0
      IF (MBEAM.LT.0) SJ = -1.0
C
C====== K is therefore the other axis, Y
C
      K = 6 - I - J
C
C====== The sign of k relies on i,j,k and X,Y,Z being right-handed
C
      IF (K.EQ.J+1 .OR. K.EQ.J-2) THEN
         SK = SI*SJ
 
      ELSE
         SK = -SI*SJ
      END IF
C
      IF (MSYS.EQ.2) THEN
C
C====== Monoclinic
C
         RL(1) = CELL(1)
         RL(2) = CELL(2)
         RL(3) = CELL(3)
C
         CC(1) = 0.0
         CC(2) = COS(CELL(4))
         CC(3) = 0.0
C
         SS(1) = 1.0
         SS(2) = SIN(CELL(4))
         SS(3) = 1.0
C
C====== A matrix
C
         AM(1,J) = SJ*RL(J)*SS(K)
         AM(1,K) = SJ*RL(K)*CC(I)
         AM(2,K) = SK*RL(K)*SS(I)*SS(J)
         AM(3,I) = SI*RL(I)
         AM(3,J) = SI*RL(J)*CC(K)
         AM(3,K) = SI*RL(K)*CC(J)
C
C====== Derivatives wrt reciprocal cell lengths
C
         DO 40 N1 = 1 , 3
            DO 30 N2 = 1 , 3
               DO 20 N3 = 1 , 3
                  DA(N1,N2,N3) = 0.
 20            CONTINUE
 30         CONTINUE
 40      CONTINUE
C
         DA(3,I,I) = SI
         DA(1,J,J) = SJ*SS(K)
         DA(3,J,J) = SI*CC(K)
         DA(1,K,K) = SJ*CC(I)
         DA(2,K,K) = SK*SS(I)*SS(J)
         DA(3,K,K) = SI*CC(J)
C
         DO 70 N1 = 1 , 3
            DO 60 N2 = 1 , 3
               DO 50 N3 = 1 , 3
                  IF (DA(N1,N2,N3).NE.0.0) THEN
                     IF (N3.EQ.1) AA(N1,N2) = DA(N1,N2,N3)
                     IF (N3.EQ.2) AB(N1,N2) = DA(N1,N2,N3)
                     IF (N3.EQ.3) AC(N1,N2) = DA(N1,N2,N3)
                  END IF
 50            CONTINUE
 60         CONTINUE
 70      CONTINUE
C
C====== Derivatives wrt reciprocal cell angles (only ss(2) and cc(2))
C
         IF (I.EQ.2) THEN
            AD(1,K) = -SJ*RL(K)*SS(I)
            AD(2,K) = SK*RL(K)*CC(I)
         END IF
 
         IF (J.EQ.2) THEN
            AD(2,K) = SK*RL(K)*CC(J)
            AD(3,K) = -SI*RL(K)*SS(J)
         END IF
 
         IF (K.EQ.2) THEN
            AD(1,J) = SJ*RL(J)*CC(K)
            AD(3,J) = -SI*RL(J)*SS(K)
         END IF
 
         RETURN
      ELSE IF (MSYS.EQ.3) THEN
C
C====== Orthorhombic
C
         RL(1) = CELL(1)
         RL(2) = CELL(2)
         RL(3) = CELL(3)
         AM(1,J) = SJ*RL(J)
         AM(2,K) = SK*RL(K)
         AM(3,I) = SI*RL(I)
C
C====== Derivatives
C
         IF (J.EQ.1) AA(1,1) = SJ
         IF (K.EQ.1) AA(2,1) = SK
         IF (I.EQ.1) AA(3,1) = SI
C
         IF (J.EQ.2) AB(1,2) = SJ
         IF (K.EQ.2) AB(2,2) = SK
         IF (I.EQ.2) AB(3,2) = SI
C
         IF (J.EQ.3) AC(1,3) = SJ
         IF (K.EQ.3) AC(2,3) = SK
         IF (I.EQ.3) AC(3,3) = SI
         RETURN
      ELSE IF (MSYS.EQ.4) THEN
C
C====== Tetragonal
C
         RL(1) = CELL(1)
         RL(2) = CELL(1)
         RL(3) = CELL(2)
         AM(1,J) = SJ*RL(J)
         AM(2,K) = SK*RL(K)
         AM(3,I) = SI*RL(I)
C
C====== Derivatives
C
         IF (J.EQ.1) AA(1,1) = SJ
         IF (K.EQ.1) AA(2,1) = SK
         IF (I.EQ.1) AA(3,1) = SI
C
         IF (J.EQ.2) AA(1,2) = SJ
         IF (K.EQ.2) AA(2,2) = SK
         IF (I.EQ.2) AA(3,2) = SI
C
         IF (J.EQ.3) AB(1,3) = SJ
         IF (K.EQ.3) AB(2,3) = SK
         IF (I.EQ.3) AB(3,3) = SI
         RETURN
      ELSE IF (MSYS.EQ.5) THEN
C
C====== Hexagonal
C
         RL(1) = CELL(1)
         RL(2) = CELL(1)
         RL(3) = CELL(2)
C
         CC(1) = 0.0
         CC(2) = 0.0
         CC(3) = 0.5
C
         SS(1) = 1.0
         SS(2) = 1.0
         SS(3) = SQRT(3.)/2.
C
C====== A matrix
C
         AM(1,J) = SJ*RL(J)*SS(K)
         AM(1,K) = SJ*RL(K)*CC(I)
         AM(2,K) = SK*RL(K)*SS(I)*SS(J)
         AM(3,I) = SI*RL(I)
         AM(3,J) = SI*RL(J)*CC(K)
         AM(3,K) = SI*RL(K)*CC(J)
C
C====== Derivatives wrt reciprocal cell lengths
C
         DO 100 N1 = 1 , 3
            DO 90 N2 = 1 , 3
               DO 80 N3 = 1 , 3
                  DA(N1,N2,N3) = 0.
 80            CONTINUE
 90         CONTINUE
 100     CONTINUE
C
         DA(3,I,I) = SI
         DA(1,J,J) = SJ*SS(K)
         DA(3,J,J) = SI*CC(K)
         DA(1,K,K) = SJ*CC(I)
         DA(2,K,K) = SK*SS(I)*SS(J)
         DA(3,K,K) = SI*CC(J)
C
         DO 130 N1 = 1 , 3
            DO 120 N2 = 1 , 3
               DO 110 N3 = 1 , 3
                  IF (DA(N1,N2,N3).NE.0.0) THEN
                     IF (N3.EQ.1 .OR. N3.EQ.2) AA(N1,N2) = DA(N1,N2,N3)
                     IF (N3.EQ.3) AB(N1,N2) = DA(N1,N2,N3)
                  END IF
 110           CONTINUE
 120        CONTINUE
 130     CONTINUE
         RETURN
      ELSE IF (MSYS.EQ.6) THEN
C
C====== Rhombohedral
C
         A = CELL(1)
         C = COS(CELL(2))
         S = SIN(CELL(2))
         FUN = SQRT(1.-3.*C**2+2.*C**3)
         AM(1,J) = SJ*A*S
         AM(1,K) = SJ*A*C*(1.-C)/S
         AM(2,K) = SK*A*FUN/S
         AM(3,I) = SI*A
         AM(3,J) = SI*A*C
         AM(3,K) = SI*A*C
C
         AA(1,J) = SJ*S
         AA(1,K) = AM(1,K)/A
         AA(2,K) = AM(2,K)/A
         AA(3,I) = SI
         AA(3,J) = SI*C
         AA(3,K) = SI*C
C
         AB(1,J) = SJ*A*C
         AB(1,K) = SJ*A*(C-(1.-C)/S**2)
         AB(2,K) = SK*A*(3.*C*(1.-C)/FUN-FUN*C/S**2)
         AB(3,J) = -SI*A*S
         AB(3,K) = -SI*A*S
         RETURN
      ELSE IF (MSYS.EQ.7) THEN
C
C====== Cubic
C
         A = CELL(1)
         AM(1,J) = SJ*A
         AM(2,K) = SK*A
         AM(3,I) = SI*A
C
         AA(1,J) = SJ
         AA(2,K) = SK
         AA(3,I) = SI
         RETURN
      ELSE
C
C====== Triclinic
C
         RL(1) = CELL(1)
         RL(2) = CELL(2)
         RL(3) = CELL(3)
C
         CC(1) = COS(CELL(4))
         CC(2) = COS(CELL(5))
         CC(3) = COS(CELL(6))
C
         SS(1) = SIN(CELL(4))
         SS(2) = SIN(CELL(5))
         SS(3) = SIN(CELL(6))
C
C====== See Buerger's Xray Crystallography
C====== for the next formulae, p348
C
         COSPSI = (CC(I)-CC(K)*CC(J))/SS(K)
         COSRHO = SQRT(1.-CC(I)**2-CC(J)**2-CC(K)**2+2.*CC(I)*CC(J)
     +            *CC(K))/SS(K)
C
C====== A matrix
C
         AM(1,J) = SJ*RL(J)*SS(K)
         AM(1,K) = SJ*RL(K)*COSPSI
         AM(2,K) = SK*RL(K)*COSRHO
         AM(3,I) = SI*RL(I)
         AM(3,J) = SI*RL(J)*CC(K)
         AM(3,K) = SI*RL(K)*CC(J)
C
C====== Derivatives wrt reciprocal cell lengths
C
         DO 160 N1 = 1 , 3
            DO 150 N2 = 1 , 3
               DO 140 N3 = 1 , 3
                  DA(N1,N2,N3) = 0.
                  DB(N1,N2,N3) = 0.
 140           CONTINUE
 150        CONTINUE
 160     CONTINUE
C
         DA(3,I,I) = SI
         DA(1,J,J) = SJ*SS(K)
         DA(3,J,J) = SI*CC(K)
         DA(1,K,K) = SJ*COSPSI
         DA(2,K,K) = SK*COSRHO
         DA(3,K,K) = SI*CC(J)
C
         DO 190 N1 = 1 , 3
            DO 180 N2 = 1 , 3
               DO 170 N3 = 1 , 3
                  IF (DA(N1,N2,N3).NE.0.0) THEN
                     IF (N3.EQ.1) AA(N1,N2) = DA(N1,N2,N3)
                     IF (N3.EQ.2) AB(N1,N2) = DA(N1,N2,N3)
                     IF (N3.EQ.3) AC(N1,N2) = DA(N1,N2,N3)
                  END IF
 170           CONTINUE
 180        CONTINUE
 190     CONTINUE
C
C====== Derivatives wrt reciprocal cell angles
C
         DB(1,K,I) = -SJ*RL(K)*SS(I)/SS(K)
         DB(2,K,I) = SK*RL(K)*SS(I)*(CC(I)-CC(J)*CC(K))
     +               /(SS(K)**2*COSRHO)
         DB(1,K,J) = SJ*RL(K)*CC(K)*SS(J)/SS(K)
         DB(2,K,J) = SK*RL(K)*SS(J)*(CC(J)-CC(I)*CC(K))
     +               /(SS(K)**2*COSRHO)
         DB(3,K,J) = -SI*RL(K)*SS(J)
         DB(1,J,K) = SJ*RL(J)*CC(K)
         DB(1,K,K) = SJ*RL(K)*(CC(J)-CC(K)*CC(I))/SS(K)**2
         DB(2,K,K) = SK*RL(K)*((CC(K)-CC(I)*CC(J))/(COSRHO*SS(K))-CC(K)
     +               *COSRHO/SS(K))
         DB(3,J,K) = -SI*RL(J)*SS(K)
C
         DO 220 N1 = 1 , 3
            DO 210 N2 = 1 , 3
               DO 200 N3 = 1 , 3
                  IF (DB(N1,N2,N3).NE.0.0) THEN
                     IF (N3.EQ.1) AD(N1,N2) = DB(N1,N2,N3)
                     IF (N3.EQ.2) AE(N1,N2) = DB(N1,N2,N3)
                     IF (N3.EQ.3) AF(N1,N2) = DB(N1,N2,N3)
                  END IF
 200           CONTINUE
 210        CONTINUE
 220     CONTINUE
         RETURN
      END IF
 
      END
**==menu.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       MENU            **
C***************************
C
C
      SUBROUTINE MENU(WORDS,NWORDS,VALUES,NVALS,IOPT)
C
      INCLUDE 'implicit_none.finc'
C
      REAL VALUES(*)
      CHARACTER*4 WORDS(*)
      INTEGER NWORDS,NVALS,IOPT
C
      INCLUDE 'chanel.finc'
C
C====== LOCALS:
C
      LOGICAL ABORT
      CHARACTER*80 STR
      CHARACTER*4 OPTION(34)
      INTEGER KOPT,INOPT,NOPT,I
C
      DATA OPTION/'A   ','B   ','C   ','ALPH','BETA','GAMM','WMIN',
     +     'WMAX','DIST','DMIN','RADI','SPDR','BMDR','PX  ','PY  ',
     +     'PZ  ','SPIN','SYST','SCAL','DATS','DATG','LAUE','LATT',
     +     'ROTA','SEAR','WAVV','VARC','RESO','OFFR','TWOR','STOP',
     +     'HELP','LOG ','    '/
C
      DATA NOPT/34/
C
C-------------------------------------------------------------------------------
C
      INOPT = IOPT
 20   CALL QSSETC(1,24)
      NWORDS = 0
      NVALS = 0
      CALL QSBELL
      WRITE (STR,6010) OPTION(INOPT)
      CALL TPROMP(STR)
      CALL TREAD(STR,ABORT)
      CALL CCPUPC(STR)
      CALL STRTKN(STR,' ,',2,WORDS,NWORDS,VALUES,NVALS)
      IF (NWORDS.GT.0) THEN
         KOPT = 0
         DO 30 I = 1 , NOPT
            IF (WORDS(1).EQ.OPTION(I)) THEN
               KOPT = I
               GO TO 40
            END IF
 30      CONTINUE
      END IF
 40   IF (.NOT.(ABORT)) THEN
         IF (NWORDS.EQ.0) KOPT = INOPT
         IOPT = KOPT
         IF (IOPT.LE.0) THEN
C
C====== Display error message
C
            CALL QSREVV
            WRITE (STR,6020) WORDS(1)
            CALL TWRITE(STR)
            CALL QSNORV
            GO TO 20
         END IF
      END IF
C
      NWORDS = NWORDS - 1
      IF (NWORDS.GE.1) WORDS(1) = WORDS(2)
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('>>> Select option [',A4,']: ')
 6020 FORMAT ('!!! Invalid option: ',A4)
      END
**==murot.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       MUROT           **
C***************************
C
C
      SUBROUTINE MUROT
C
      INCLUDE 'implicit_none.finc'
C
      REAL DEGRAD
      PARAMETER (DEGRAD=3.14159265358979/180.)
      INTEGER MAXONE, MAXALL
      PARAMETER (MAXONE=20000,MAXALL=100000)
C
      COMMON /DATA/CELL, WVM, WVX, WVS, CF, PHIS, SPNDLE, CCELL, NUMA,
     +       NUMC, DMIN
      COMMON /SPOTS/NR
      REAL CF, DMIN, SPNDLE, WVM, WVS, WVX
      INTEGER NR, NUMA, NUMC
      REAL CCELL(6), CELL(6), PHIS(3)
C
      INCLUDE 'chanel.finc'
C
C====== LOCALS:
C
      LOGICAL QUIT
      INTEGER ION(MAXONE)
      INTEGER KSAVE(MAXALL)
      REAL PHISAV(3)
      REAL FRACT, STEP, PHIU, MU, SINPHIU, SINPHIY, COSPHIY
      INTEGER I, ICOUNT, IST, ITOT, J, K, NSTEP
      LOGICAL ISAV
      REAL C(4), SP(20)
C
C-------------------------------------------------------------------------------
C
      WRITE (IPRINT,6010)
      MU = 0.0
      CALL TASKF('Enter spindle offset angle µ (degrees)',MU,1,QUIT)
      MU = MU*DEGRAD
 
      NSTEP = 1
      CALL TASKI('Enter number of steps (counting a spindle)'//
     +           ' setting of 0 as the first)',NSTEP,QUIT)
      STEP = 10.0
      CALL TASKF('Enter step size in degrees',STEP,1,QUIT)
 
      WRITE (IPRINT,6020) MU/DEGRAD
      WRITE (IUNIT6,6020) MU/DEGRAD
      WRITE (IPRINT,6030)
      WRITE (IUNIT6,6030)
C
      DO 20 I = 1 , 3
         PHISAV(I) = PHIS(I)
 20   CONTINUE
      ISAV = .TRUE.
      IST = 0
C
C====== Loop through the spindle steps
C
      DO 80 I = 1 , NSTEP
         PHIU = (I-1)*STEP
         SP(I) = PHIU
         PHIU = PHIU*DEGRAD
         SINPHIU = SIN(PHIU)
         SINPHIY = -SIN(MU)*COS(MU)*(1-COS(PHIU))
         PHIS(2) = ASIN(SINPHIY)
         COSPHIY = COS(PHIS(2))
         PHIS(3) = ASIN(COS(MU)*SINPHIU/COSPHIY)
         PHIS(1) = ASIN(-SIN(MU)*SINPHIU/COSPHIY)
         DO 30 J = 1 , 3
            PHIS(J) = PHIS(J)/DEGRAD + PHISAV(J)
 30      CONTINUE
C
C====== Recalculate the A matrix since phis have changed
C
         CELL(1) = CELL(1)*WVS
         CELL(2) = CELL(2)*WVS
         CELL(3) = CELL(3)*WVS
         J = 0
         DO 40 K = 1 , 3
            IF (CELL(K).NE.0.) THEN
               J = J + 1
               CCELL(J) = CELL(K)
            END IF
 
 40      CONTINUE
 
         CALL ACALC
         CELL(1) = CELL(1)/WVS
         CELL(2) = CELL(2)/WVS
         CELL(3) = CELL(3)/WVS
         CALL LAUE(ION)
C
C====== Find out how many reflections recorded
C
         ICOUNT = 0
         ITOT = 0
         DO 50 J = 1 , 4
            C(J) = 0.
 50      CONTINUE
         DO 60 J = 1 , NR
            K = ION(J)
C
C======= If there is room save the ion array for overlap purposes
C
            IF (ISAV) KSAVE(IST+J) = K
            IF (K.NE.0) THEN
               ITOT = ITOT + K
               IF (K.GT.4) K = 4
               C(K) = C(K) + 1
               ICOUNT = ICOUNT + 1
            END IF
 
 60      CONTINUE
 
         FRACT = FLOAT(ICOUNT)/FLOAT(NR)
         DO 70 J = 1 , 4
            C(J) = C(J)/NR
 70      CONTINUE
         WRITE (IPRINT,6040) PHIU/DEGRAD , FRACT , ICOUNT , (C(J),J=1,4)
     +                       , ITOT
         WRITE (IUNIT6,6040) PHIU/DEGRAD , FRACT , ICOUNT , (C(J),J=1,4)
     +                       , ITOT
         IST = IST + NR
         IF ((IST+NR).GT.MAXALL) ISAV = .FALSE.
 80   CONTINUE
      DO 90 I = 1 , 3
         PHIS(I) = PHISAV(I)
 90   CONTINUE
C
C====== Investigate overlaps?
C
      IF (NSTEP.GT.1) CALL OVERLP(KSAVE,NSTEP,SP)
C
C====== Format statements
C
 6010 FORMAT (' This command analyses the effect of rotating about '/
     +        'the spindle axis at a fixed angle µ.')
 6020 FORMAT (///'  Effect of rotating about spindle at µ = ',F6.1,
     +        ' on the fraction of symmetry'/
     +        '  unique reflections recorded')
 6030 FORMAT (///
     +' Spindle   Fraction recorded  Distribution of symmetry equivalent
     +s'/32X,'  1     2     3     4+'/)
 6040 FORMAT (F6.1,F13.3,' (',I5,')',4X,4F6.3,' (',I5,')')
      END
**==overcn.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       OVERCN          **
C***************************
C
C
      SUBROUTINE OVERCN(ICNT,IREC,KSAVE,NR)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: Count in icnt the reflections in ksave not recorded in irec
C          irec=1 denotes not yet recorded, ksave.ne.0 denotes on film
C
      INTEGER ICNT,NR
      INTEGER KSAVE(NR),IREC(NR)
C
C====== LOCALS:
C
      INTEGER I
C
C-------------------------------------------------------------------------------
C
      ICNT = 0
      DO 20 I = 1 , NR
         IF (IREC(I).EQ.1 .AND. KSAVE(I).NE.0) ICNT = ICNT + 1
 20   CONTINUE
      RETURN
      END
**==overlp.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       OVERLP          **
C***************************
C
C
      SUBROUTINE OVERLP(KSAVE,NSTEP,SP)
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER MAXONE,MAXALL
      PARAMETER (MAXONE=20000,MAXALL=100000)
C
      COMMON /SPOTS/NR
      INTEGER NSTEP
      REAL SP(20)
      INTEGER KSAVE(MAXALL)
      INTEGER NR
C
      INCLUDE 'chanel.finc'
C
C====== LOCALS:
C
      LOGICAL QUIT
      REAL FR,FRLEFT
      INTEGER I,ICNT,IFILM,IMAX,IMAXNO,IOFF,IOFF2,IST,J,K
      REAL OV(20)
      INTEGER KSET(20)
      INTEGER IREC(MAXALL)
C
C-------------------------------------------------------------------------------
C
      I = MAXALL/NR
      IF (I.LT.NSTEP) THEN
         NSTEP = I
         WRITE (IPRINT,6010) NSTEP
         WRITE (IUNIT6,6010) NSTEP
      END IF
 
 20   CALL TWRITE('Analyse the overlap of the symmetry uniques')
      IST = 1
      CALL TASKI('Enter row in above table to use as base (0=end)',IST,
     +           QUIT)
      IF (IST.EQ.0) RETURN
      IF (IST.GT.NSTEP) GO TO 20
      WRITE (IPRINT,6020)
      WRITE (IUNIT6,6020)
C
C====== Zero the array which says which spindle settings used
C====== the value in the array gives the number of the film in the series
C====== ifilm counts through the series
C
      DO 30 I = 1 , NSTEP
         KSET(I) = 0
 30   CONTINUE
      KSET(IST) = 1
      IFILM = 1
C
C====== Set irec array to say all uniques not recorded at start
C
      DO 40 I = 1 , NR
         IREC(I) = 1
 40   CONTINUE
C
C====== Calculate the ones recorded at the first spindle setting
C
      IOFF = (IST-1)*NR + 1
      CALL OVERCN(ICNT,IREC,KSAVE(IOFF),NR)
      CALL OVERST(IREC,KSAVE(IOFF),NR)
      FRLEFT = 1.0
      FR = FLOAT(ICNT)/FLOAT(NR)
      FRLEFT = FRLEFT - FR
      WRITE (IPRINT,6030) SP(IST) , FR , FRLEFT
      WRITE (IUNIT6,6030) SP(IST) , FR , FRLEFT
C
C====== Loop through the remaining spindle settings
C====== imax records max overlap with unrecordeds, imaxno the spindle for it
C
 50   IMAX = -1
      IMAXNO = 0
      DO 60 K = 1 , NSTEP
         IF (KSET(K).EQ.0) THEN
            IOFF = (K-1)*NR + 1
            CALL OVERCN(ICNT,IREC,KSAVE(IOFF),NR)
            IF (ICNT.GT.IMAX) THEN
               IMAX = ICNT
               IMAXNO = K
            END IF
         END IF
 
 60   CONTINUE
C
C====== Check for end - no spindle settings remain
C
      IF (IMAXNO.EQ.0) GO TO 20
C
C====== Output the best overlap spindle setting and repeat
C
      IOFF = NR*(IMAXNO-1) + 1
      CALL OVERST(IREC,KSAVE(IOFF),NR)
      FR = FLOAT(IMAX)/FLOAT(NR)
      FRLEFT = FRLEFT - FR
C
C====== Count the overlaps of this spindle setting with the earlier ones
C
      DO 80 J = 1 , IFILM
         DO 70 K = 1 , NSTEP
            IF (KSET(K).EQ.J) THEN
               IOFF2 = NR*(K-1) + 1
               CALL OVERPR(ICNT,KSAVE(IOFF),KSAVE(IOFF2),NR)
               OV(J) = FLOAT(ICNT)/FLOAT(NR)
               GO TO 80
 
            END IF
 
 70      CONTINUE
 80   CONTINUE
      WRITE (IPRINT,6030) SP(IMAXNO) , FR , FRLEFT , (OV(J),J=1,IFILM)
      WRITE (IUNIT6,6030) SP(IMAXNO) , FR , FRLEFT , (OV(J),J=1,IFILM)
      IFILM = IFILM + 1
      KSET(IMAXNO) = IFILM
C
C====== but only keep going if some unrecordeds remain
C====== or if the last spindle setting yielded some overlaps
C
      IF (IMAX.EQ.0) GO TO 20
      IF (FRLEFT.LE.0.00001) GO TO 20
      GO TO 50
C
C====== Format statements
C
 6010 FORMAT (' *** warning *** only room to save first',I3,' rows')
 6020 FORMAT (//
     +      ' Spndl  Add. fract.  Fract. symm.       Fractional overlap'
     +      /
     +    '        symm. uniq.  uniques not      with previous spindles'
     +    /
     + '         captured     recorded    1      2      3      4      5'
     + /'                        1.000')
 6030 FORMAT (F6.1,4X,F6.3,7X,F6.3,2X,5F7.3/31X,5F7.3/31X,5F7.3/31X,
     +        5F7.3)
      END
**==overpr.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       OVERPR          **
C***************************
C
C
      SUBROUTINE OVERPR(ICNT,KSAVE,KSAVE2,NR)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: Count in icnt the reflections common to ksave and ksave2
C
      INTEGER ICNT,NR
      INTEGER KSAVE(NR),KSAVE2(NR)
C
C====== LOCALS:
C
      INTEGER I
C
C-------------------------------------------------------------------------------
C
      ICNT = 0
      DO 20 I = 1 , NR
         IF (KSAVE(I).NE.0 .AND. KSAVE2(I).NE.0) ICNT = ICNT + 1
 20   CONTINUE
      RETURN
      END
**==overst.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       OVERST          **
C***************************
C
C
      SUBROUTINE OVERST(IREC,KSAVE,NR)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: Set the irec array to show as recorded the reflection in ksave
C
      INTEGER NR
      INTEGER KSAVE(NR),IREC(NR)
C
C====== LOCALS:
C
      INTEGER I
C
C-------------------------------------------------------------------------------
C
 
      DO 20 I = 1 , NR
         IF (IREC(I).EQ.1 .AND. KSAVE(I).NE.0) IREC(I) = 0
 20   CONTINUE
      RETURN
      END
**==resoln.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       RESOLN          **
C***************************
C
C
      SUBROUTINE RESOLN()
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER MAXONE,MAXALL
      PARAMETER (MAXONE=20000,MAXALL=100000)
C
      COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN
      COMMON /SPOTS/NR
      REAL CF,DMIN,SPNDLE,WVM,WVS,WVX
      INTEGER NR,NUMA,NUMC
      REAL CCELL(6),CELL(6),PHIS(3)
C
      INCLUDE 'chanel.finc'
C
C====== LOCALS:
C
      INTEGER ION(MAXALL)
      REAL WRES(MAXALL)
      REAL DMAX,F1,F2,F3,V
      INTEGER ICOUNT,J,K,L
      REAL CON(11)
      INTEGER IBOX(10)
C
      DATA CON/.0,.10,.20,.35,.55,.80,1.2,1.8,2.8,4.5,100./
C
C-------------------------------------------------------------------------------
C
 
      WRITE (IPRINT,6010)
      WRITE (IPRINT,6020)
      WRITE (IUNIT6,6020)
C
      CALL LAUER(ION,WRES)
C
C====== Find out how many reflections recorded
C
      ICOUNT = 0
      DMAX = 0.
      DO 20 J = 1 , NR
         K = ION(J)
         IF (K.NE.0) THEN
            IF (WRES(J).GT.DMAX) DMAX = WRES(J)
            ICOUNT = ICOUNT + 1
         END IF
 20   CONTINUE
      CON(11) = DMAX - DMIN
C
      WRITE (IPRINT,6030) ICOUNT , FLOAT(ICOUNT)/FLOAT(NR)
      WRITE (IUNIT6,6030) ICOUNT , FLOAT(ICOUNT)/FLOAT(NR)
      DO 30 K = 1 , 10
         IBOX(K) = 0
 30   CONTINUE
      DO 60 J = 1 , NR
         IF (ION(J).EQ.0) GO TO 60
         V = WRES(J)
         K = 0
         DO 40 L = 1 , 11
            IF (V.LT.(DMIN+CON(L))) GO TO 50
            K = K + 1
 40      CONTINUE
         K = 10
 50      IBOX(K) = IBOX(K) + 1
 60   CONTINUE
      F2 = 0.0
      F3 = FLOAT(ICOUNT)/FLOAT(NR)
      DO 70 K = 1 , 10
         IF (K.NE.1) F3 = F3 - F1
         F1 = FLOAT(IBOX(K))/FLOAT(NR)
         F2 = F2 + F1
         WRITE (IPRINT,6040) F1 , F2 , F3 , CON(K) + DMIN , CON(K+1)
     +                       + DMIN
         WRITE (IUNIT6,6040) F1 , F2 , F3 , CON(K) + DMIN , CON(K+1)
     +                       + DMIN
 70   CONTINUE
      RETURN
C
C====== Format statements
C
 6010 FORMAT (' This command calculates the resolution distribution',
     +        ' for'/'  the symmetry equivalents in this orientation.')
 6020 FORMAT (///'    Distribution of resolution of reflections'/)
 6030 FORMAT (' Symmetry uniques recorded  =',I6,' (',F4.2,
     +        ')'//' Fract. of unique   Cumulative    Resolution range'/
     +        '   reflections       fraction           '/)
 6040 FORMAT (F12.3,5X,2F7.3,F9.2,' - ',F6.2)
      END
**==rotspi.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       ROTSPI          **
C***************************
C
C
      SUBROUTINE ROTSPI()
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER MAXONE,MAXALL
      PARAMETER (MAXONE=20000,MAXALL=100000)
C
      COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN
      COMMON /SPOTS/NR
      REAL CF,DMIN,SPNDLE,WVM,WVS,WVX
      INTEGER NR,NUMA,NUMC
      REAL CCELL(6),CELL(6),PHIS(3)
C
      INCLUDE 'chanel.finc'
C
C====== LOCALS:
C
      LOGICAL QUIT
      INTEGER ION(MAXONE)
      INTEGER KSAVE(MAXALL)
      REAL FRACT,SPSAV,STEP
      INTEGER I,ICOUNT,IST,ITOT,J,K,NSTEP
      LOGICAL ISAV
      REAL C(4),SP(20)
C
C-------------------------------------------------------------------------------
C
      WRITE (IPRINT,6010)
      NSTEP = 1
      CALL TASKI('Enter number of steps (counting the present spindle'//
     +           ' setting as the first)',NSTEP,QUIT)
      STEP = 10.0
      CALL TASKF('Enter step size in degrees',STEP,1,QUIT)
      WRITE (IPRINT,6020)
      WRITE (IUNIT6,6020)
      WRITE (IPRINT,6030)
      WRITE (IUNIT6,6030)
C
C====== Loop through the spindle steps
C
      SPSAV = SPNDLE
      SPNDLE = SPNDLE - STEP
      ISAV = .TRUE.
      IST = 0
      DO 50 I = 1 , NSTEP
         SPNDLE = SPNDLE + STEP
         SP(I) = SPNDLE
         CALL LAUE(ION)
C
C====== Find out how many reflections recorded
C
         ICOUNT = 0
         ITOT = 0
         DO 20 J = 1 , 4
            C(J) = 0.
 20      CONTINUE
         DO 30 J = 1 , NR
            K = ION(J)
C
C====== If there is room save the ion array for overlap purposes
C
            IF (ISAV) KSAVE(IST+J) = K
            IF (K.NE.0) THEN
               ITOT = ITOT + K
               IF (K.GT.4) K = 4
               C(K) = C(K) + 1
               ICOUNT = ICOUNT + 1
            END IF
 30      CONTINUE
         FRACT = FLOAT(ICOUNT)/FLOAT(NR)
         DO 40 J = 1 , 4
            C(J) = C(J)/NR
 40      CONTINUE
         WRITE (IPRINT,6040) SPNDLE , FRACT , ICOUNT , (C(J),J=1,4) , 
     +                       ITOT
         WRITE (IUNIT6,6040) SPNDLE , FRACT , ICOUNT , (C(J),J=1,4) , 
     +                       ITOT
         IST = IST + NR
         IF ((IST+NR).GT.MAXALL) ISAV = .FALSE.
 50   CONTINUE
      SPNDLE = SPSAV
C
C====== Investigate overlaps?
C
      IF (NSTEP.GT.1) CALL OVERLP(KSAVE,NSTEP,SP)
      RETURN
C
C====== Format statements
C
 6010 FORMAT (' This command analyses the effect of rotating about ',
     +        'the spindle axis.')
 6020 FORMAT (///
     +  '  Effect of rotating about spindle on the fraction of symmetry'
     +  /'  unique reflections recorded')
 6030 FORMAT (///
     +' Spindle   Fraction recorded  Distribution of symmetry equivalent
     +s'/32X,'  1     2     3     4+'/)
 6040 FORMAT (F6.1,F13.3,' (',I5,')',4X,4F6.3,' (',I5,')')
      END
**==screen.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       SCREEN          **
C***************************
C
C
      SUBROUTINE SCREEN(IENT,ICHANG)
C
      INCLUDE 'implicit_none.finc'
C
      COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN
      COMMON /MATRIX/MSYS,MROT,MBEAM,ILGP,ICENT
      COMMON /RADC/RADI,RADI2
      COMMON /SYS/ISYS,KSYS
      INTEGER ICHANG,IENT
      REAL CF,DMIN,RADI,RADI2,RSUB,SPNDLE,SSUB,WVM,WVS,WVX,XSUB,YSUB
      INTEGER ICENT,ILGP,ISYS,MBEAM,MROT,MSYS,NUMA,NUMC
      REAL CCELL(6),CELL(6),PHIS(3)
      INTEGER KSYS(3)
C
      INCLUDE 'chanel.finc'
C
C====== FUNCTION TYPES:
C
      INTEGER LENSTR
C
C====== LOCALS:
C
      LOGICAL QUIT
      REAL RAPI,VALUES(5),VAL
      CHARACTER*4 WORDS(5)
      INTEGER I,IOK,IOPT,J,NVALS,NWORDS
      SAVE IOPT
      LOGICAL LDATEX,LOGG,LSTIEX
      CHARACTER CHARS*4,DATFIL*60,SAVFIL*60
      CHARACTER LABEL(35)*6,NAME(5)*15,PROMPT*50
C
C-------------------------------------------------------------------------------
C
 
      DATA LABEL/'a*    ','b*    ','c*    ','alpha*','beta* ','gamma*',
     +     'wmin  ','wmax  ','dist  ','dmin  ','radi  ','spdr  ',
     +     'bmdr  ','px    ','py    ','pz    ','spin  ','syst  ',
     +     'scal  ','dats  ','datg  ','laue  ','latt  ','rotate',
     +     'search','wavvar','varcf ','resoln','offrot','tworot',
     +     'stop  ','      ','log   ','      ','      '/
C
C====== First entry is ient=0. on subsequent entries, jump to menu.
C
      IF (IENT.GT.0) THEN
         ICHANG = 0
         GO TO 30
 
      END IF
 
      ICHANG = 1
C
C====== Switch logging off initially
C
      LOGG = .FALSE.
      IUNIT6 = 2
C
C====== Set default menu selection to `DATG'
C
      IOPT = 21 - 1
      WORDS(1) = 'DATG'
C
C====== Initialise some values
C
C
C====== Cell parameters to silly values to show which are changed
C
      DO 20 I = 1 , 6
         CELL(I) = 0.0
 20   CONTINUE
C
C====== Matrix to zero. msys=system, mrot=rotation axis, mbeam=beam axis
C
      MSYS = 0
      MROT = 3
      MBEAM = 1
      ISYS = 0
      ICENT = 1
      KSYS(1) = 0
      KSYS(2) = 0
      KSYS(3) = 0
      ILGP = -1
C
C====== Geometry to 0.0
C
      PHIS(1) = 0.0
      PHIS(2) = 0.0
      PHIS(3) = 0.0
      WVM = 0.3
      WVX = 2.6
      CF = 50.
      WVS = 5.0
C
C====== Default spindle setting to 0.0, dmin to 2.5, radius to 59.3 mm.
C
      SPNDLE = 0.0
      DMIN = 2.5
      RADI = 59.3
      RADI2 = RADI*RADI
C
C====== Display menu and get next command
C
 30   CALL UPDMEN(LABEL,LOGG)
      IOPT = IOPT + 1
      CALL MENU(WORDS,NWORDS,VALUES,NVALS,IOPT)
C
C====== Reciprocal cell parameters
C
      IF (IOPT.LE.6) THEN
         IF (IOPT.LE.3) WRITE (PROMPT,6010) LABEL(IOPT)
         IF (IOPT.GT.3) WRITE (PROMPT,6020) LABEL(IOPT)
         CALL INPUTV(PROMPT,CELL(IOPT),VALUES,NVALS)
 
         ICHANG = 1
         GO TO 30
 
      END IF
C
C====== Various limits
C
      IF (IOPT.EQ.7) THEN
         CALL INPUTV('Wavelength minimum (angstroms): ',WVM,VALUES,
     +               NVALS)
         ICHANG = 1
         GO TO 30
 
      END IF
 
      IF (IOPT.EQ.8) THEN
         CALL INPUTV('Wavelength maximum (angstroms): ',WVX,VALUES,
     +               NVALS)
         ICHANG = 1
         GO TO 30
 
      END IF
 
      IF (IOPT.EQ.9) THEN
         CALL INPUTV('Crystal-to-film distance in mm: ',CF,VALUES,NVALS)
         GO TO 30
 
      END IF
 
      IF (IOPT.EQ.10) THEN
         CALL INPUTV('Enter dmin - resolution limit in Angstroms: ',
     +               DMIN,VALUES,NVALS)
         ICHANG = 1
         GO TO 30
 
      END IF
 
      IF (IOPT.EQ.11) THEN
         CALL INPUTV('Enter radius limit for plot - mm: ',RADI,VALUES,
     +               NVALS)
         ICHANG = 1
         RADI2 = RADI*RADI
         GO TO 30
 
      END IF
C
C====== Orientation
C
C
C====== Spindle axis
C
      IF (IOPT.EQ.12) THEN
         CALL TWRITE('Axis 1 = astar, -3 = -cstar etc.')
         CALL INPUTV('Reciprocal axis lying along crystal rotation'//
     +               ' axis: ',VAL,VALUES,NVALS)
         MROT = VAL
         IF (MROT.EQ.0 .OR. MROT.LT.-3 .OR. MROT.GT.3) THEN
            WRITE (IPRINT,6030)
            MROT = 3
 
         END IF
         GO TO 30
 
      END IF
 
C
C====== Beam axis
C
      IF (IOPT.EQ.13) THEN
         CALL TWRITE('Reciprocal axis lying nearest to X-ray beam'//
     +               ' axis in the plane')
         CALL INPUTV('defined by the beam axis and the rotation axis: ',
     +               VAL,VALUES,NVALS)
         MBEAM = VAL
         IF (MBEAM.EQ.0 .OR. MBEAM.LT.-3 .OR. MBEAM.GT.3 .OR. 
     +       IABS(MBEAM).EQ.IABS(MROT)) THEN
            WRITE (IPRINT,6030)
            MBEAM = 1
 
         END IF
         GO TO 30
 
      END IF
C
C====== Phi angles
C
      IF (IOPT.EQ.14) THEN
         CALL INPUTV('Rotation about beam axis, phi(x) in degrees: ',
     +               PHIS(1),VALUES,NVALS)
         GO TO 30
 
      END IF
 
      IF (IOPT.EQ.15) THEN
         CALL INPUTV('Rotation about vertical axis, phi(y) in degrees: '
     +               ,PHIS(2),VALUES,NVALS)
         GO TO 30
 
      END IF
 
      IF (IOPT.EQ.16) THEN
         CALL INPUTV(
     +             'Rotation about horizontal axis, phi(z) in degrees: '
     +             ,PHIS(3),VALUES,NVALS)
         GO TO 30
      END IF
C
C====== Miscellaneous film parameters
C
      IF (IOPT.EQ.17) THEN
         CALL INPUTV('Spindle setting: ',SPNDLE,VALUES,NVALS)
         GO TO 30
 
      END IF
C
C====== Lattice input
C
      IF (IOPT.EQ.18) THEN
         CALL INPUTW('Crystal system: ',CHARS,WORDS,NWORDS)
         IF (CHARS.EQ.'TRIC') THEN
            MSYS = 1
         ELSE IF (CHARS.EQ.'MONO') THEN
            MSYS = 2
         ELSE IF (CHARS.EQ.'ORTH') THEN
            MSYS = 3
         ELSE IF (CHARS.EQ.'TETR') THEN
            MSYS = 4
         ELSE IF (CHARS.EQ.'HEXA') THEN
            MSYS = 5
         ELSE IF (CHARS.EQ.'RHOM') THEN
            MSYS = 6
         ELSE IF (CHARS.EQ.'CUBI') THEN
            MSYS = 7
         ELSE
            WRITE (IPRINT,6040)
            MSYS = 0
            GO TO 30
         END IF
 
         ICHANG = 1
C
C====== May as well get user to set Laue group now
C
         CALL LGRP(MSYS,ILGP)
         GO TO 30
 
      END IF
 
      IF (IOPT.EQ.19) THEN
         CALL INPUTV(
     +              'Scale factor that multiplies recip. cell lengths: '
     +              ,WVS,VALUES,NVALS)
         ICHANG = 1
         GO TO 30
 
      END IF
 
C
C====== Save file for control data
C
      IF (IOPT.EQ.20) THEN
         CALL TPROMP(
     +        'Enter name of file in which to save control data: ')
         CALL TREAD(SAVFIL,QUIT)
C
C====== Look for first blank, append .dat if "." not seen
C
         I = 0
 40      I = I + 1
         IF (SAVFIL(I:I).NE.' ' .AND. SAVFIL(I:I).NE.'.') GO TO 40
         IF (SAVFIL(I:I).EQ.' ') THEN
            SAVFIL(I:I+3) = '.dat'
         END IF
         OPEN (UNIT=IUNIT7,FILE=SAVFIL,STATUS='UNKNOWN',
     +         FORM='FORMATTED')
C
C====== Write the records ... for contents see the comments
C====== above concerning reading them back.
C
         WRITE (IUNIT7,6050) PHIS , CELL
         WRITE (IUNIT7,6060) WVM , WVX , WVS , CF , SPNDLE , DMIN , RADI
         WRITE (IUNIT7,6070) MSYS , MROT , MBEAM , ISYS , KSYS , ILGP , 
     +                       ICENT
         WRITE (IUNIT7,6080)
C
C====== Final mosaic spread record
C
         REWIND (UNIT=IUNIT7)
         CLOSE (UNIT=IUNIT7)
         GO TO 30
      END IF
C
C====== File definition
C
      IF (IOPT.EQ.21) THEN
         CALL TPROMP(
     +        'Enter name of previously saved control data file: ')      
         CALL TREAD(DATFIL,QUIT)
C
C====== Look for first blank, append .dat if "." not seen
C
         I = 0
 50      I = I + 1
         IF (DATFIL(I:I).NE.' ' .AND. DATFIL(I:I).NE.'.') GO TO 50
         IF (DATFIL(I:I).EQ.' ') THEN
            DATFIL(I:I+3) = '.dat'
         END IF
         INQUIRE (FILE=DATFIL,EXIST=LDATEX)
         IF (.NOT.LDATEX) THEN
            WRITE (IPRINT,6090) DATFIL(1:LENSTR(DATFIL))
            GO TO 30
 
         END IF
 
         OPEN (UNIT=IUNIT7,FILE=DATFIL,STATUS='OLD')
C
C====== Record 1: orientation angles, cell parameters
C
         READ (IUNIT7,*) PHIS , CELL
C
C====== Record 2: wavelengths, scale, crystal-to-film distance, dmin, radius.
C
         READ (IUNIT7,*) WVM , WVX , WVS , CF , SPNDLE , DMIN , RADI
         RADI2 = RADI*RADI
C
C====== Record 3: matrix - msys, mrot, mbeam, isys, ksys, ilgp, icent
C
         READ (IUNIT7,*) MSYS , MROT , MBEAM , ISYS , KSYS , ILGP , 
     +                   ICENT
C
         REWIND (UNIT=IUNIT7)
         CLOSE (UNIT=IUNIT7)
         ICHANG = 1
         GO TO 30
 
      END IF
C
C====== Laue group
C
      IF (IOPT.EQ.22) THEN
         IF (MSYS.EQ.0) THEN
            WRITE (IPRINT,6100)
            GO TO 30
 
         END IF
 
         CALL LGRP(MSYS,ILGP)
         ICHANG = 1
         GO TO 30
 
      END IF
 
C
C====== Lattice centering
C
      IF (IOPT.EQ.23) THEN
         CALL INPUTW('Lattice centering: P,A,B,C,I,F or R [P]: ',CHARS,
     +               WORDS,NWORDS)
C
C====== Default lattice centering is p
C
         IF (CHARS.EQ.' ') CHARS = 'P'
C
C====== R only applies to rhombohedral indexed on hexagonal lattice
C
         IF (CHARS.EQ.'R' .AND. MSYS.NE.5) CHARS = 'P'
         IF (CHARS.EQ.'P') THEN
            ISYS = 0
            ICENT = 1
            KSYS(1) = 0
            KSYS(2) = 0
            KSYS(3) = 0
            ICHANG = 1
            GO TO 30
 
         END IF
 
         IF (CHARS.EQ.'A') THEN
            ISYS = 2
            ICENT = 2
            KSYS(2) = 1
            KSYS(3) = 1
            ICHANG = 1
            GO TO 30
 
         END IF
 
         IF (CHARS.EQ.'B') THEN
            ISYS = 2
            ICENT = 3
            KSYS(1) = 1
            KSYS(3) = 1
            ICHANG = 1
            GO TO 30
 
         END IF
 
         IF (CHARS.EQ.'C') THEN
            ISYS = 2
            ICENT = 4
            KSYS(1) = 1
            KSYS(2) = 1
            ICHANG = 1
            GO TO 30
 
         END IF
 
         IF (CHARS.EQ.'I') THEN
            ISYS = 2
            ICENT = 5
            KSYS(1) = 1
            KSYS(2) = 1
            KSYS(3) = 1
            ICHANG = 1
            GO TO 30
 
         END IF
 
         IF (CHARS.EQ.'F') THEN
            ISYS = 4
            ICENT = 6
            ICHANG = 1
            GO TO 30
 
         END IF
 
         IF (CHARS.EQ.'R') THEN
            ISYS = 3
            ICENT = 7
            KSYS(1) = -1
            KSYS(2) = 1
            KSYS(3) = 1
            ICHANG = 1
            GO TO 30
 
         END IF
 
         WRITE (IPRINT,6110)
         ISYS = 0
         ICENT = 1
         KSYS(1) = 0
         KSYS(2) = 0
         KSYS(3) = 0
         GO TO 30
 
      END IF
C
C====== Stop signal.
C
      IF (IOPT.EQ.31) THEN
         CALL QSSCRL(1,24)
         CALL QSCLRA
         STOP
      END IF
C
C====== Help information. (Sorry - removed now)
C
      IF (IOPT.EQ.32) GO TO 30
C
C====== Log file status.
C
      IF (IOPT.EQ.33) THEN
         CALL INPUTW('Output to log file: ',CHARS,WORDS,NWORDS)
         IF (CHARS.EQ.'ON' .OR. CHARS(1:1).EQ.'Y') THEN
            LOGG = .TRUE.
            IUNIT6 = 1
         ELSE IF (CHARS.EQ.'OFF' .OR. CHARS(1:1).EQ.'N') THEN
            LOGG = .FALSE.
            IUNIT6 = 2
         ELSE
            WRITE (IPRINT,6120) CHARS
         END IF
         IOPT = 20
         GO TO 30
      END IF
C
C====== Check things before proceeding
C
      IOK = 0
      IF (CF.EQ.0.0) THEN
         WRITE (IPRINT,6130)
         IOK = 1
      END IF
C
C====== Check cell parameters compatible with lattice type
C
      IF (MSYS.GT.0) CALL CELLTEST(MSYS,ILGP,ICENT,CELL,LABEL,IOK)
C
C====== Copy cell parameters to contiguous array
C
      RAPI = 3.14159265/180.
      NUMC = 0
      NUMA = 0
      J = 0
      DO 60 I = 1 , 6
         IF (CELL(I).NE.0.) THEN
            J = J + 1
            CCELL(J) = CELL(I)
            NUMC = NUMC + 1
            IF (I.GT.3) THEN
               NUMA = NUMA + 1
               CCELL(J) = CCELL(J)*RAPI
            END IF
         END IF
 
 60   CONTINUE
      IF (NUMC.EQ.0) THEN
         WRITE (IPRINT,6140)
         IOK = 1
      END IF
C
C====== Check for matrix
C
      IF (MSYS.EQ.0) THEN
         WRITE (IPRINT,6150)
         IOK = 1
      END IF
 
      IF (ILGP.EQ.-1) THEN
         WRITE (IPRINT,6160)
         IOK = 1
      END IF
 
      IF (IOK.NE.0) GO TO 30
C
C====== Return to main routine to execute the command
C
      IENT = IOPT - 23
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('Recip. lattice length (dimensionless) ',A6,': ')
 6020 FORMAT ('Reciprocal lattice angle (degrees) ',A6,': ')
 6030 FORMAT (' !!! Invalid axis')
 6040 FORMAT (' !!! Invalid crystal lattice')
 6050 FORMAT (3F12.6,3F10.6/3F10.4)
 6060 FORMAT (3F10.6,F10.4,2F10.2,F10.2)
 6070 FORMAT (9I3)
 6080 FORMAT (' 0.2000')
 6090 FORMAT (' !!! File ',A,' does not exist')
 6100 FORMAT (' Please enter crystal lattice type first!')
 6110 FORMAT (' !!! Invalid lattice centering')
 6120 FORMAT (' !!! Illegal choice: ',A4)
 6130 FORMAT (' !!! No value for crystal to film distance')
 6140 FORMAT (' !!! No cell parameters')
 6150 FORMAT (' !!! No lattice information')
 6160 FORMAT (' !!! No laue group')
      END
**==search.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       SEARCH          **
C***************************
C
C
      SUBROUTINE SEARCH()
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER MAXONE,MAXALL
      PARAMETER (MAXONE=20000,MAXALL=100000)
C
      COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN
      COMMON /SPOTS/NR
      REAL CF,DMIN,SPNDLE,WVM,WVS,WVX
      INTEGER NR,NUMA,NUMC
      REAL CCELL(6),CELL(6),PHIS(3)
C
      INCLUDE 'chanel.finc'
C
C====== LOCALS:
C
      LOGICAL QUIT
      INTEGER ION(MAXONE)
      REAL PHIYSV,SPSAV,STEPY,STEPZ
      INTEGER I,ICOUNT,J,K,M,NSTEPY,NSTEPZ
      REAL FR(30)
C
C-------------------------------------------------------------------------------
C
      WRITE (IPRINT,6010)
      WRITE (IUNIT6,6020)
      CALL TWRITE(' ')
      NSTEPY = 1
      CALL TWRITE('Enter number of steps about the vertical axis ')
      CALL TASKI('(counting the present phi-y as the first)',NSTEPY,
     +           QUIT)
      STEPY = 10.0
      CALL TASKF('Enter step size in degrees',STEPY,1,QUIT)
      NSTEPZ = 1
      CALL TWRITE('Enter number of steps about the spindle axis ')
      CALL TASKI('(counting the present phi-z as the first)',NSTEPZ,
     +           QUIT)
      STEPZ = 10.0
      CALL TASKF('Enter step size in degrees',STEPZ,1,QUIT)
C
      SPSAV = SPNDLE
      SPNDLE = 0.
      PHIYSV = PHIS(2)
      WRITE (IPRINT,6030)
      WRITE (IUNIT6,6030)
      WRITE (IPRINT,6040) ((PHIS(3)+(I-1)*STEPZ),I=1,NSTEPZ)
      WRITE (IUNIT6,6040) ((PHIS(3)+(I-1)*STEPZ),I=1,NSTEPZ)
      WRITE (IPRINT,6050)
      WRITE (IUNIT6,6050)
C
C====== Loop through the vertical axis steps
C
      PHIS(2) = PHIS(2) - STEPY
      DO 50 M = 1 , NSTEPY
         PHIS(2) = PHIS(2) + STEPY
C
C====== Recalculate the a matrix since phis have changed
C
         CELL(1) = CELL(1)*WVS
         CELL(2) = CELL(2)*WVS
         CELL(3) = CELL(3)*WVS
         J = 0
         DO 20 I = 1 , 3
            IF (CELL(I).NE.0.) THEN
               J = J + 1
               CCELL(J) = CELL(I)
            END IF
 20      CONTINUE
         CALL ACALC
         CELL(1) = CELL(1)/WVS
         CELL(2) = CELL(2)/WVS
         CELL(3) = CELL(3)/WVS
C
C====== Loop through the spindle steps
C
         SPNDLE = -STEPZ
         DO 40 I = 1 , NSTEPZ
            SPNDLE = SPNDLE + STEPZ
            CALL LAUE(ION)
C
C====== Find out how many reflections recorded
C
            ICOUNT = 0
            DO 30 J = 1 , NR
               K = ION(J)
               IF (K.NE.0) THEN
                  ICOUNT = ICOUNT + 1
               END IF
 30         CONTINUE
            FR(I) = FLOAT(ICOUNT)/FLOAT(NR)
 40      CONTINUE
         WRITE (IPRINT,6060) PHIS(2) , (FR(I),I=1,NSTEPZ)
         WRITE (IUNIT6,6060) PHIS(2) , (FR(I),I=1,NSTEPZ)
 50   CONTINUE
      PHIS(2) = PHIYSV
      SPNDLE = SPSAV
      RETURN
C
C====== Format statements
C
 6010 FORMAT (' This command searches the vertical axis and spindle',
     +        ' axis (phi-y, phi-z)'/
     +        ' to explore the fraction of symmetry',
     +        ' unique reflections recorded.'//
     +        ' Whatever its value now, the',
     +        ' spindle setting will be taken as zero.')
 6020 FORMAT (' Whatever its value now, the',
     +        ' spindle setting will be taken as zero.')
 6030 FORMAT (///'  Fraction of symmetry unique reflections recorded')
 6040 FORMAT (//'       phi-z',9F6.1/12X,9F6.1/12X,9F6.1)
 6050 FORMAT ('  phi-y')
 6060 FORMAT (1X,F6.1,5X,9F6.3/12X,9F6.3/12X,9F6.3)
      END
**==sysabs.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       SYSABS          **
C***************************
C
C
      SUBROUTINE SYSABS(IH,IK,IL,ABSNT)
C
      INCLUDE 'implicit_none.finc'
C
      COMMON /SYS/ISYS,KSYS
      INTEGER IH,IK,IL
      LOGICAL ABSNT
      INTEGER ISYS
      INTEGER KSYS(3)
C
C====== LOCALS:
C
      INTEGER I,J
      INTEGER H(3)
C
C-------------------------------------------------------------------------------
C
      H(1) = IH
      H(2) = IK
      H(3) = IL
      IF (ISYS.EQ.4) THEN
 
         ABSNT = .TRUE.
         DO 20 I = 1 , 2
            IF (MOD(H(I)+H(I+1),2).NE.0) RETURN
 20      CONTINUE
         ABSNT = .FALSE.
         RETURN
      ELSE
C
C====== This version for a,b,c,r,i & f lattice absences
C
         ABSNT = .FALSE.
         J = 0
         DO 30 I = 1 , 3
            J = J + H(I)*KSYS(I)
 30      CONTINUE
         IF (MOD(J,ISYS).NE.0) ABSNT = .TRUE.
         RETURN
      END IF
      END
**==sysahk.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       SYSAHK          **
C***************************
C
C
      SUBROUTINE SYSAHK(IH,IK,ABSNT)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: This variant of sysabs checks to see if a given h,k is absent,
C          whatever the l value
C
      INTEGER IH,IK
      LOGICAL ABSNT
      INTEGER ISYS
      INTEGER KSYS(3)
      COMMON /SYS/ISYS,KSYS
C
C-------------------------------------------------------------------------------
C
      ABSNT = .FALSE.
      IF ((ISYS.EQ.2 .OR. ISYS.EQ.4) .AND. KSYS(3).EQ.0) THEN
         IF (MOD(IH+IK,2).NE.0) ABSNT = .TRUE.
      END IF
 
      RETURN
      END
**==tworot.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       TWOROT          **
C***************************
C
C
      SUBROUTINE TWOROT
C
      INCLUDE 'implicit_none.finc'
C
      REAL DEGRAD
      PARAMETER (DEGRAD=3.14159265358979/180.)
      INTEGER MAXONE, MAXALL
      PARAMETER (MAXONE=20000,MAXALL=100000)
C
      COMMON /DATA/CELL, WVM, WVX, WVS, CF, PHIS, SPNDLE, CCELL, NUMA,
     +       NUMC, DMIN
      COMMON /SPOTS/NR
      REAL CF, DMIN, SPNDLE, WVM, WVS, WVX
      INTEGER NR, NUMA, NUMC
      REAL CCELL(6), CELL(6), PHIS(3)
C
      INCLUDE 'chanel.finc'
C
C====== LOCALS:
C
      LOGICAL QUIT
      INTEGER ION(MAXONE)
      INTEGER KSAVE(MAXALL)
      REAL PHISAV(3), COSPHIY
      REAL FRACT, USTEP, PHIU, MU, MUSTEP
      INTEGER I, IU, IMU, ICOUNT, IST, ITOT, J, K, NU, NMU, NALL, IALL
      LOGICAL ISAV
      REAL C(4), SP(20), ANGSAV(20,2)
C
C-------------------------------------------------------------------------------
C
      WRITE (IPRINT,6010)
      NU = 1
      CALL TASKI('Enter number of phi_U steps (counting a spindle'//
     +           ' setting of 0 as the first)',NU,QUIT)
      USTEP = 10.0
      CALL TASKF('Enter step size in degrees',USTEP,1,QUIT)
      NMU = 1
      CALL TASKI('Enter number of µ steps (counting µ=0 as the first)',
     +           NMU,QUIT)
      MUSTEP = 10.0
      CALL TASKF('Enter step size in degrees',MUSTEP,1,QUIT)
      NALL = NMU*NU
      IF (NALL.GT.20) THEN
         WRITE (IPRINT,6020) NALL
         WRITE (IUNIT6,6020) NALL
      ELSE
         WRITE (IPRINT,6030)
         WRITE (IUNIT6,6030)
         WRITE (IPRINT,6040)
         WRITE (IUNIT6,6040)
C
         DO 20 I = 1 , 3
            PHISAV(I) = PHIS(I)
 20      CONTINUE
         ISAV = .TRUE.
         IST = 0
C
C====== Loop through the spindle steps
C
         DO 90 IU = 1 , NU
            PHIU = (IU-1)*USTEP
            PHIU = PHIU*DEGRAD
            DO 80 IMU = 1 , NMU
               IALL = (IU-1)*NMU + IMU
               SP(IALL) = FLOAT(IALL)
               MU = (IMU-1)*MUSTEP
               MU = MU*DEGRAD
               ANGSAV(IALL,1) = PHIU
               ANGSAV(IALL,2) = MU
               PHIS(2) = ASIN(SIN(MU)*COS(PHIU))
               COSPHIY = COS(PHIS(2))
               PHIS(3) = ASIN(SIN(PHIU)/COSPHIY)
               PHIS(1) = ACOS(COS(MU)/COSPHIY)
               DO 30 J = 1 , 3
                  PHIS(J) = PHIS(J)/DEGRAD + PHISAV(J)
 30            CONTINUE
C
C====== Recalculate the A matrix since phis have changed
C
               CELL(1) = CELL(1)*WVS
               CELL(2) = CELL(2)*WVS
               CELL(3) = CELL(3)*WVS
               J = 0
               DO 40 K = 1 , 3
                  IF (CELL(K).NE.0.) THEN
                     J = J + 1
                     CCELL(J) = CELL(K)
                  END IF
 
 40            CONTINUE
 
               CALL ACALC
               CELL(1) = CELL(1)/WVS
               CELL(2) = CELL(2)/WVS
               CELL(3) = CELL(3)/WVS
               CALL LAUE(ION)
C
C====== Find out how many reflections recorded
C
               ICOUNT = 0
               ITOT = 0
               DO 50 J = 1 , 4
                  C(J) = 0.
 50            CONTINUE
               DO 60 J = 1 , NR
                  K = ION(J)
C
C====== If there is room save the ion array for overlap purposes
C
                  IF (ISAV) KSAVE(IST+J) = K
                  IF (K.NE.0) THEN
                     ITOT = ITOT + K
                     IF (K.GT.4) K = 4
                     C(K) = C(K) + 1
                     ICOUNT = ICOUNT + 1
                  END IF
 
 60            CONTINUE
 
               FRACT = FLOAT(ICOUNT)/FLOAT(NR)
               DO 70 J = 1 , 4
                  C(J) = C(J)/NR
 70            CONTINUE
               WRITE (IPRINT,6050) IALL , FRACT , ICOUNT , (C(J),J=1,4)
     +                             , ITOT
               WRITE (IUNIT6,6050) IALL , FRACT , ICOUNT , (C(J),J=1,4)
     +                             , ITOT
               IST = IST + NR
               IF ((IST+NR).GT.MAXALL) ISAV = .FALSE.
 80         CONTINUE
 90      CONTINUE
         DO 100 I = 1 , 3
            PHIS(I) = PHISAV(I)
 100     CONTINUE
         WRITE (IPRINT,6060)
         WRITE (IUNIT6,6060)
         DO 110 I = 1 , NALL
            WRITE (IPRINT,6070) I , ANGSAV(I,1)/DEGRAD , ANGSAV(I,2)
     +                          /DEGRAD
            WRITE (IUNIT6,6070) I , ANGSAV(I,1)/DEGRAD , ANGSAV(I,2)
     +                          /DEGRAD
 110     CONTINUE
C
C====== investigate overlaps?
C
         IF (NALL.GT.1) CALL OVERLP(KSAVE,NALL,SP)
      END IF
C
C====== Format statements
C
 6010 FORMAT (' This command analyses the effect of changing '/
     +        ' both the spindle axis setting phi_U',
     +        ' and the spindle offset angle µ.')
 6020 FORMAT (' Too many steps altogether (',I3,'), maximum 20')
 6030 FORMAT (///'  Effect of changing phi_U and µ',
     +        ' on the fraction of symmetry'/
     +        '  unique reflections recorded')
 6040 FORMAT (///
     +' Spindle#  Fraction recorded  Distribution of symmetry equivalent
     +s'/32X,'  1     2     3     4+'/)
 6050 FORMAT (I6,F13.3,' (',I5,')',4X,4F6.3,' (',I5,')')
 6060 FORMAT (//' Table of spindle numbers and the corresponding'/
     +        ' phi_U, µ settings'//' Spindle#       phi_U       µ')
 6070 FORMAT (1X,I3,9X,2F8.1)
      END
**==unique.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       UNIQUE          **
C***************************
C
C
      SUBROUTINE UNIQUE(IH,IK,IL,UNIQ)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: This routine returns uniq=true if reflection h,k,l is in the
C          unique set according to the Laue symmetry, otherwise iflag=false.
C
C    To be compatible with the order of reflection generation in hklgen
C    this routine must ensure that it does not reject all l for given k
C    (which causes hklgen to stop incrementing k) if a subsequent k
C    would be in the unique set. Similarly, do not reject all k for given
C    h if a subsequent would be in the unique set.
C
      COMMON /MATRIX/MSYS,MROT,MBEAM,ILGP,ICENT
      INTEGER IH,IK,IL
      LOGICAL UNIQ
      INTEGER ICENT,ILGP,MBEAM,MROT,MSYS
      INTRINSIC IABS
C
C-------------------------------------------------------------------------------
C
 
      UNIQ = .FALSE.
      IF (ILGP.EQ.0) THEN
C
C====== 1bar
C
         IF (IH.LT.0) RETURN
         IF (IH.EQ.0) THEN
            IF (IK.LT.0) RETURN
            IF (IK.EQ.0) THEN
               IF (IL.LT.0) RETURN
            END IF
 
 
         END IF
         GO TO 30
 
      END IF
C
      IF (ILGP.EQ.2) THEN
C
C====== mmm
C
         IF (IH.LT.0) RETURN
         IF (IK.LT.0) RETURN
         IF (IL.LT.0) RETURN
      ELSE IF (ILGP.EQ.3) THEN
C
C====== 4/m
C
         IF (IH.LT.0) RETURN
         IF (IK.LT.0) RETURN
         IF (IH.EQ.0 .AND. IK.NE.0) RETURN
         IF (IL.LT.0) RETURN
      ELSE IF (ILGP.EQ.4) THEN
C
C====== 4/mmm
C
         IF (IH.LT.0) RETURN
         IF (IK.GT.IH) RETURN
         IF (IK.LT.0) RETURN
         IF (IL.LT.0) RETURN
      ELSE IF (ILGP.EQ.5) THEN
C
C====== 6/m
C
         IF (IH.LT.0) RETURN
         IF (IK.LT.0) RETURN
         IF (IL.LT.0) RETURN
         IF (IH.EQ.0 .AND. IK.NE.0) RETURN
      ELSE IF (ILGP.EQ.6) THEN
C
C====== 6/mmm
C
         IF (IH.LT.0) RETURN
         IF (IK.LT.0) RETURN
         IF (IK.GT.IH) RETURN
         IF (IL.LT.0) RETURN
      ELSE IF (ILGP.EQ.7) THEN
C
C====== 3bar ( hexagonal or rhombohedral lattices possible)
C
         IF (MSYS.EQ.5) THEN
            IF (IH.LT.0) RETURN
            IF (IK.LT.0) RETURN
            IF (IH.EQ.0) THEN
               IF (IK.NE.0) RETURN
               IF (IL.LT.0) RETURN
            END IF
 
         ELSE IF (MSYS.EQ.6) THEN
            IF (IH.LE.0) RETURN
            IF (IABS(IL).GT.0) THEN
               IF (IABS(IL).EQ.IH .AND. IK.NE.IH) RETURN
               IF (IABS(IL).LT.IH .AND. 
     +             (IABS(IK).LT.IABS(IL) .OR. IABS(IK).GE.IH)) RETURN
               IF (IABS(IL).GT.IH .AND. 
     +             (IABS(IK).GT.IABS(IL) .OR. IABS(IK).LE.IH)) RETURN
            END IF
 
 
         END IF
      ELSE IF (ILGP.EQ.8) THEN
         GO TO 20
      ELSE IF (ILGP.EQ.9) THEN
         GO TO 20
      ELSE IF (ILGP.EQ.10) THEN
C
C====== m3
C
         IF (IH.LE.0) RETURN
         IF (IL.LT.0) RETURN
         IF (IK.LT.0) RETURN
         IF (IL.GT.0) THEN
            IF (IL.EQ.IH .AND. IK.NE.IH) RETURN
            IF (IL.LT.IH .AND. (IK.LT.IL .OR. IK.GE.IH)) RETURN
            IF (IL.GT.IH .AND. (IK.GT.IL .OR. IK.LE.IH)) RETURN
 
         END IF
      ELSE IF (ILGP.EQ.11) THEN
C
C====== m3m
C
         IF (IH.LE.0) RETURN
         IF (IL.LT.0) RETURN
         IF (IK.LT.0) RETURN
         IF (IK.GT.IH) RETURN
         IF (IK.EQ.0 .AND. IL.NE.0) RETURN
         IF (IK.NE.0 .AND. IK.NE.IH .AND. IL.GE.IK) RETURN
      ELSE
C
C====== 2/m
C
         IF (IH.LT.0) RETURN
         IF (IK.LT.0) RETURN
         IF (IH.EQ.0 .AND. IL.LT.0) RETURN
      END IF
      GO TO 30
C
C====== 3m1 and 31m ( hexagonal or rhombohedral lattices possible)
C
 20   IF (MSYS.EQ.5) THEN
         IF (IH.LT.0) RETURN
         IF (IK.LT.0) RETURN
         IF (IK.GT.IH) RETURN
         IF (IH.EQ.0 .AND. IL.LT.0) RETURN
         IF (IK.EQ.0 .AND. ILGP.EQ.9 .AND. IL.LT.0) RETURN
         IF (IK.EQ.IH .AND. ILGP.EQ.8 .AND. IL.LT.0) RETURN
 
      ELSE IF (MSYS.EQ.6) THEN
         STOP ' Sorry this case not yet implemented'
 
      END IF
C
 30   UNIQ = .TRUE.
      RETURN
      END
**==unpack.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       UNPACK          **
C***************************
C
C
      SUBROUTINE UNPACK(N,H,K,L)
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER H,K,L,N
C
C-------------------------------------------------------------------------------
C
      L = MOD(N,200) - 100
      K = MOD(N/200,200) - 100
      H = N/40000 - 100
      RETURN
      END
**==updmen.f  
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       UPDMEN          **
C***************************
C
C
      SUBROUTINE UPDMEN(LABEL,LOGG)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: Updates option menu with associated values
C
      CHARACTER *6 LABEL(34)
      LOGICAL LOGG
C
C LABEL  : NAMES OF OPTIONS
C
      INCLUDE 'chanel.finc'
C
      COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN
      COMMON /MATRIX/MSYS,MROT,MBEAM,ILGP,ICENT
      COMMON /RADC/RADI,RADI2
      COMMON /SYS/ISYS,KSYS
      REAL CF,DMIN,RADI,RADI2,SPNDLE,WVM,WVS,WVX
      INTEGER ICENT,ILGP,ISYS,MBEAM,MROT,MSYS,NUMA,NUMC
      REAL CCELL(6),CELL(6),PHIS(3)
      INTEGER KSYS(3)
C
C====== LOCALS:
C
      CHARACTER *15 TEMP
      INTEGER I
      CHARACTER MAXIS(7)*7,MSYST(8)*7,MCENT(7)*7,MLAUE(13)*6
C
      DATA MAXIS/'  -cs  ','  -bs  ','  -as  ','       ','  +as  ',
     +     '  +bs  ','  +cs  '/
      DATA MSYST/'       ','tricl. ','monocl.','ortho. ','tetrag.',
     +     'hexag. ','rhomb. ','cubic  '/
      DATA MCENT/'   P   ','   A   ','   B   ','   C   ','   I   ',
     +     '   F   ','   R   '/
      DATA MLAUE/'      ','1-bar ',' 2/m  ',' mmm  ',' 4/m  ','4/mmm ',
     +     ' 6/m  ','6/mmm ','3-bar ',' 3m1  ',' 31m  ','  m3  ',
     +     ' m3m  '/
C
C-------------------------------------------------------------------------------
C
      DO 20 I = 1 , 6
         IF (CELL(I).EQ.0.) THEN
            WRITE (TEMP,6030) LABEL(I)
            CALL ITEM(I,TEMP,15,.FALSE.)
 
         ELSE
            WRITE (TEMP,6040) LABEL(I) , CELL(I)
            CALL ITEM(I,TEMP,15,.FALSE.)
         END IF
C
 20   CONTINUE
C
      WRITE (TEMP,6040) LABEL(7) , WVM
      CALL ITEM(7,TEMP,15,.FALSE.)
C
      WRITE (TEMP,6040) LABEL(8) , WVX
      CALL ITEM(8,TEMP,15,.FALSE.)
C
      WRITE (TEMP,6010) LABEL(9) , CF
      CALL ITEM(9,TEMP,15,.FALSE.)
C
      WRITE (TEMP,6010) LABEL(10) , DMIN
      CALL ITEM(10,TEMP,15,.FALSE.)
C
      WRITE (TEMP,6010) LABEL(11) , RADI
      CALL ITEM(11,TEMP,15,.FALSE.)
      RADI2 = RADI*RADI
C
      WRITE (TEMP,6020) LABEL(12) , MAXIS(MROT+4)
      CALL ITEM(12,TEMP,15,.FALSE.)
C
      WRITE (TEMP,6020) LABEL(13) , MAXIS(MBEAM+4)
      CALL ITEM(13,TEMP,15,.FALSE.)
C
      DO 30 I = 1 , 3
         WRITE (TEMP,6010) LABEL(13+I) , PHIS(I)
         CALL ITEM(13+I,TEMP,15,.FALSE.)
 30   CONTINUE
C
      WRITE (TEMP,6010) LABEL(17) , SPNDLE
      CALL ITEM(17,TEMP,15,.FALSE.)
C
      WRITE (TEMP,6020) LABEL(18) , MSYST(MSYS+1)
      CALL ITEM(18,TEMP,15,.FALSE.)
C
      WRITE (TEMP,6010) LABEL(19) , WVS
      CALL ITEM(19,TEMP,15,.FALSE.)
C
      WRITE (TEMP,6030) LABEL(20)
      CALL ITEM(20,TEMP,15,.TRUE.)
C
      WRITE (TEMP,6030) LABEL(21)
      CALL ITEM(21,TEMP,15,.TRUE.)
C
      WRITE (TEMP,6020) LABEL(22) , MLAUE(ILGP+2)
      CALL ITEM(22,TEMP,15,.FALSE.)
C
      WRITE (TEMP,6020) LABEL(23) , MCENT(ICENT)
      CALL ITEM(23,TEMP,15,.FALSE.)
C
C====== Write out the commands.
C
      DO 40 I = 24 , 32
         WRITE (TEMP,6030) LABEL(I)
         CALL ITEM(I,TEMP,15,.TRUE.)
 40   CONTINUE
C
      IF (LOGG) THEN
         WRITE (TEMP,6020) LABEL(33) , 'ON'
      ELSE
         WRITE (TEMP,6020) LABEL(33) , 'OFF'
      END IF
      CALL ITEM(33,TEMP,15,.TRUE.)
C
C====== Now write out parameters to log file
C
      IF (LOGG) THEN
         DO 50 I = 1 , 6
            IF (CELL(I).EQ.0.) THEN
               WRITE (IUNIT6,6070) LABEL(I)
 
            ELSE
               WRITE (IUNIT6,6080) LABEL(I) , CELL(I)
            END IF
C
 50      CONTINUE
C
         WRITE (IUNIT6,6080) LABEL(7) , WVM
C
         WRITE (IUNIT6,6080) LABEL(8) , WVX
C
         WRITE (IUNIT6,6050) LABEL(9) , CF
C
         WRITE (IUNIT6,6050) LABEL(10) , DMIN
C
         WRITE (IUNIT6,6050) LABEL(11) , RADI
C
         WRITE (IUNIT6,6060) LABEL(12) , MAXIS(MROT+4)
C
         WRITE (IUNIT6,6060) LABEL(13) , MAXIS(MBEAM+4)
C
         DO 60 I = 1 , 3
            WRITE (IUNIT6,6050) LABEL(13+I) , PHIS(I)
 60      CONTINUE
C
         WRITE (IUNIT6,6050) LABEL(17) , SPNDLE
C
         WRITE (IUNIT6,6060) LABEL(18) , MSYST(MSYS+1)
C
         WRITE (IUNIT6,6050) LABEL(19) , WVS
C
         WRITE (IUNIT6,6060) LABEL(22) , MLAUE(ILGP+2)
C
         WRITE (IUNIT6,6060) LABEL(23) , MCENT(ICENT)
      END IF
      RETURN
C
C====== Format statements
C
 6010 FORMAT (A6,F8.2,1X)
 6020 FORMAT (A6,1X,A7,1X)
 6030 FORMAT (A6,9X)
 6040 FORMAT (A6,F8.4,1X)
 6050 FORMAT (1X,A6,F8.2)
 6060 FORMAT (1X,A6,1X,A7)
 6070 FORMAT (1X,A6)
 6080 FORMAT (1X,A6,F8.4)
      END
**==varycf.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       VARYCF          **
C***************************
C
C
      SUBROUTINE VARYCF()
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER MAXONE,MAXALL
      PARAMETER (MAXONE=20000,MAXALL=100000)
C
      COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN
      COMMON /SPOTS/NR
      REAL CF,DMIN,SPNDLE,WVM,WVS,WVX
      INTEGER NR,NUMA,NUMC
      REAL CCELL(6),CELL(6),PHIS(3)
C
      INCLUDE 'chanel.finc'
C
C====== LOCALS:
C
      LOGICAL QUIT
      INTEGER ION(MAXONE)
      REAL CFSAV,FRACT,STEP
      INTEGER I,ICOUNT,ITOT,J,K,NSTEP
C
C-------------------------------------------------------------------------------
C
      WRITE (IPRINT,6010)
      NSTEP = 1
      CALL TASKI('Enter number of steps (counting the present c-to-f'//
     +           ' setting as the first)',NSTEP,QUIT)
      STEP = 5.0
      CALL TASKF('Enter step size in mm',STEP,1,QUIT)
      WRITE (IPRINT,6020)
      WRITE (IUNIT6,6020)
      WRITE (IPRINT,6030)
      WRITE (IUNIT6,6030)
C
C====== Loop through the c-to-f steps
C
      CFSAV = CF
      CF = CF - STEP
      DO 30 I = 1 , NSTEP
         CF = CF + STEP
         CALL LAUE(ION)
C
C====== Find out how many reflections recorded
C
         ICOUNT = 0
         ITOT = 0
         DO 20 J = 1 , NR
            K = ION(J)
            IF (K.NE.0) THEN
               ITOT = ITOT + K
               ICOUNT = ICOUNT + 1
            END IF
 20      CONTINUE
         FRACT = FLOAT(ICOUNT)/FLOAT(NR)
         WRITE (IPRINT,6040) CF , FRACT , ICOUNT , ITOT
         WRITE (IUNIT6,6040) CF , FRACT , ICOUNT , ITOT
 30   CONTINUE
      CF = CFSAV
      RETURN
C
C====== Format statements
C
 6010 FORMAT (' This command analyses the effect of changing the',
     +        ' crystal to film distance.')
 6020 FORMAT (///
     +'  Effect of increasing c-to-f on the fraction of symmetry unique'
     +/'  reflections recorded')
 6030 FORMAT (///' cf(mm)   Fraction of symmetry     Total reflections'/
     +        '            uniques recorded           on film'/)
 6040 FORMAT (F6.1,F13.3,' (',I5,')',12X,I5)
      END
**==wavvar.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       WAVVAR          **
C***************************
C
C
      SUBROUTINE WAVVAR()
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER MAXONE,MAXALL
      PARAMETER (MAXONE=20000,MAXALL=100000)
C
      COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN
      COMMON /SPOTS/NR
      REAL CF,DMIN,SPNDLE,WVM,WVS,WVX
      INTEGER NR,NUMA,NUMC
      REAL CCELL(6),CELL(6),PHIS(3)
C
      INCLUDE 'chanel.finc'
C
C====== LOCALS:
C
      INTEGER ION(MAXALL)
      REAL WION(MAXALL)
      REAL DEL,WMAX
      INTEGER ICOUNT,ITOT,J,K
      INTEGER IBOX(10)
C
C-------------------------------------------------------------------------------
C
      WRITE (IPRINT,6010)
      WRITE (IPRINT,6020)
      WRITE (IUNIT6,6020)
C
      CALL LAUEW(ION,WION)
C
C====== Find out how many reflections recorded
C
      ICOUNT = 0
      WMAX = 0.
      ITOT = 0
      DO 20 J = 1 , NR
         K = ION(J)
         IF (K.GT.1) THEN
            ITOT = ITOT + K
            IF (WION(J).GT.WMAX) WMAX = WION(J)
            ICOUNT = ICOUNT + 1
         END IF
 20   CONTINUE
C
      WRITE (IPRINT,6030) ICOUNT , ITOT
      WRITE (IUNIT6,6030) ICOUNT , ITOT
      DEL = WMAX/10.
      IF (DEL.LT..01) DEL = .01
      DO 30 K = 1 , 10
         IBOX(K) = 0
 30   CONTINUE
      DO 40 J = 1 , NR
         IF (ION(J).GT.1) THEN
            K = WION(J)/DEL + 1
            IF (K.LE.0) K = 1
            IF (K.GT.10) K = 10
            IBOX(K) = IBOX(K) + 1
         END IF
 40   CONTINUE
      DO 50 K = 1 , 10
         WRITE (IPRINT,6040) IBOX(K) , DEL*(K-1) , DEL*K
         WRITE (IUNIT6,6040) IBOX(K) , DEL*(K-1) , DEL*K
 50   CONTINUE
      RETURN
C
C====== Format statements
C
 6010 FORMAT (' This command calculates the wavelength distribution',
     +        ' for'/'  the symmetry equivalents in this orientation.')
 6020 FORMAT (///'    Distribution of wavelength variations'/)
 6030 FORMAT (' Symmetry uniques with duplicates recorded =',
     +        I6/' Total symmetry duplicates recorded        =',
     +        I6//' Number of reflections   Wavelength variation range'/
     +        )
 6040 FORMAT (I10,7X,F19.2,' - ',F4.2)
      END
