      PROGRAM THKL
C
C     TRANSFORM hkl AND/OR TEST SYSTEMATIC ABSENCES.
C
C     COMPILE WITH /CONTINUATIONS=50.
C
      CHARACTER ATIME*8,ADATE*9,TITLE*80,FILE*80
      DIMENSION CELL(6),T(3,3),G(3,3),GINV(3,3)
      DIMENSION JKL(3,13)
      DATA JKL / 1,0,0, 0,2,0, 0,0,3,
     & 1,2,0, -1,2,0, 1,0,3, -1,0,3, 0,2,3, 0,-2,3,
     & 1,2,3, -1,2,3, 1,-2,3, 1,2,-3/
      COMMON /STATS/ CUTOFF,NA(34),NP(34),XA(34),XP(34)
      DATA CELL,T,CUTOFF,SMAX,NA,NP,XA,XP,NTOT,NSIG /155*0/
      data ii,a1,a2,a3,a4,xtime /6*0/
      CALL TIME (ATIME)
      CALL DATE (ADATE)
      IO1=10
      IO2=20
      ILP=60
      OPEN (UNIT=IO1,FILE='thkl.dat',STATUS='OLD')
      READ (IO1,'(A)') TITLE
      READ (IO1,'(A)') FILE
      READ (IO1,*,END=5) ((T(I,J),J=1,3),I=1,3)
      READ (IO1,*,END=5) CELL
      READ (IO1,*,END=5) CUTOFF,SMAX
    5 CLOSE (UNIT=IO1)
      OPEN (UNIT=IO1,FILE=FILE,STATUS='OLD',FORM='UNFORMATTED')
      OPEN (UNIT=ILP,FILE='thkl.lp',STATUS='NEW')
      DET=T(1,1)*T(2,2)*T(3,3)+T(1,2)*T(2,3)*T(3,1)+T(1,3)*T(2,1)*T(3,2)
     & -T(3,1)*T(2,2)*T(1,3)-T(3,2)*T(2,3)*T(1,1)-T(3,3)*T(2,1)*T(1,2)
      IF (DET.NE.0) THEN
        WRITE (ILP,6000) ATIME,ADATE,TITLE
        WRITE (ILP,6100) FILE
        WRITE (ILP,6001) ((T(I,J),J=1,3),I=1,3),CELL
        CALL NEWABC (CELL,T)
        WRITE (ILP,6002) CELL
        DO I=1,13
          J=JKL(1,I)
          K=JKL(2,I)
          L=JKL(3,I)
          IH=J
          IK=K
          IL=L
          CALL NEWHKL (IH,IK,IL,T)
          WRITE (ILP,6003) J,K,L,IH,IK,IL
        END DO
        IF (DET.LT.0) THEN
          WRITE (ILP,6004)
          TYPE 6004
        END IF
        OPEN (UNIT=IO2,FILE='data.newhkl',STATUS='NEW')
        TYPE 6600
        WRITE (ILP,6000) ATIME,ADATE,TITLE
        WRITE (ILP,6006)
      END IF
 6000 FORMAT ('1PROGRAM THKL.  ',A8,', ',A9,'.  ',A)
 6100 FORMAT ('0REFLECTION DATA FILE:  ',A)
 6001 FORMAT ('0UNIT CELL AND INDEX TRANSFORMATION:'/
     &'0( NEW A )   ( T11 T12 T13 ) ( OLD A )'/
     &' ( NEW B ) = ( T21 T22 T23 ) ( OLD B )'/
     &' ( NEW C )   ( T31 T32 T33 ) ( OLD C )'/
     &'0( NEW H )   ( T11 T12 T13 ) ( OLD H )'/
     &' ( NEW K ) = ( T21 T22 T23 ) ( OLD K )'/
     &' ( NEW L )   ( T31 T32 T33 ) ( OLD L )'/
     &'0(T11 T12 T13)   (',3F7.3,' )'/
     &' (T21 T22 T23) = (',3F7.3,' )'/
     &' (T31 T32 T33)   (',3F7.3,' )'/
     &'0OLD CELL:         A         B         C     ALPHA      BETA',
     &'     GAMMA'/' --- ---- ',3F10.4,3F10.3)
 6002 FORMAT (
     &'0NEW CELL:         A         B         C     ALPHA      BETA',
     &'     GAMMA'/' --- ---- ',3F10.4,3F10.3/
     &'0OLD  H  K  L   NEW  H  K  L'/
     &' ---  -  -  -   ---  -  -  -')
 6003 FORMAT (' ',3X,3I3,6X,3I3)
 6004 FORMAT ('0ATTENTION!  THE GIVEN T-MATRIX CORRESPONDS TO A LEFT-HAN
     &DED TRANSFORMATION.'/' ---------')
 6006 FORMAT ('0OLD  H  K  L   NEW  H  K  L         Y  SIGMA(Y)     XTIM
     &E'/' ---  -  -  -   ---  -  -  -         -  --------     -----')
 6005 FORMAT (' ',3X,3I3,6X,3I3,2F10.2,F10.3)
 6600 FORMAT ('0DATA WITH TRANSFORMED hkl WILL BE WRITTEN TO A NEW REFLE
     &CTION FILE NAMED'/' ''data.newhkl''.')
      IF (CUTOFF.LE.0) CUTOFF=3
      IF (SMAX.GT.0) CALL METRIC (CELL,G,GINV)
      IEND=0
    1 CALL READ1 (IO1,IEND,II,IH,IK,IL,A1,A2,A3,A4,Y,SIGY,XTIME)
      IF (IEND.NE.0) GO TO 9
      IF (SIGY.LE.0) GO TO 1
      NTOT=NTOT+1
      IF (Y/SIGY.GE.CUTOFF) NSIG=NSIG+1
      IF (DET.NE.0) THEN
        IPRINT=0
        CALL DPRINT (II,IH,IK,IL,IPRINT)
        IF (IPRINT.NE.0) THEN
          J=IH
          K=IK
          L=IL
        END IF
        CALL NEWHKL (IH,IK,IL,T)
        CALL WRITE1 (IO2,II,IH,IK,IL,A1,A2,A3,A4,Y,SIGY,XTIME)
        IF (IPRINT.NE.0) WRITE (ILP,6005) J,K,L,IH,IK,IL,Y,SIGY,XTIME
      END IF
      IF (SMAX.GT.0.AND.SINTHL(IH,IK,IL,GINV).GT.SMAX) GO TO 1
      CALL HKLTEST (IH,IK,IL,Y,SIGY)
      GO TO 1
    9 CONTINUE
      WRITE (ILP,6000) ATIME,ADATE,TITLE
      WRITE (ILP,6100) FILE
      WRITE (ILP,6200) CUTOFF,NTOT,NSIG
      IF (SMAX.GT.0) WRITE (ILP,6300) SMAX
      IF (DET.NE.0) WRITE (ILP,6400)
      WRITE (ILP,6009) (NA(I)+NP(I),
     & FLOAT(NA(I))/MAX(1,NA(I)+NP(I)),XA(I)/MAX(1,NA(I)),
     & FLOAT(NP(I))/MAX(1,NA(I)+NP(I)),XP(I)/MAX(1,NP(I)),I=1,34)
 6200 FORMAT ('0POSSIBLE SYSTEMATIC ABSENCES:'/
     &'0ABSENT IF FSQ .LT. CUTOFF*SIGMA(FSQ), WHERE CUTOFF = ',F5.2,'.'/
     &'0N0 = ',I7,' TOTAL MEASUREMENTS'/
     &' N1 = ',I7,' MEASUREMENTS WITH FSQ .GE. CUTOFF*SIGMA(FSQ)')
 6300 FORMAT ('0hkl TESTED IF SIN(THETA)/LAMBDA < ',F5.3,'.')
 6400 FORMAT ('0TESTS REFER TO NEW hkl.')
 6009 FORMAT ('0F = N(ABSENT OR PRESENT)/N(MEASURED)'/'0Z = <FSQ/SIGMA(F
     &SQ)>'/
     &'0CLASS                        MEASURED      ABSENT        PRESENT
     &'/
     &'                                            F     Z       F     Z
     &'/
     &' -----                        --------     ---   ---     ---   --
     &-'/
     &'0hkl GENERAL CONDITIONS:'/
     &'0hkl,h+k=2n+1               ',I10,2(2X,F6.2,F6.1)/
     &' hkl,k+l=2n+1               ',I10,2(2X,F6.2,F6.1)/
     &' hkl,l+h=2n+1               ',I10,2(2X,F6.2,F6.1)/
     &' hkl,h+k,k+l(,l+h)=2n+1     ',I10,2(2X,F6.2,F6.1)/
     &' hkl,h+k+l=2n+1             ',I10,2(2X,F6.2,F6.1)/
     &' hkl,-h+k+l=3n+1 or 2       ',I10,2(2X,F6.2,F6.1)/
     &' hkl,h-k+l=3n+1 or 2        ',I10,2(2X,F6.2,F6.1)/
     &' hkl,h-k=3n+1 or 2          ',I10,2(2X,F6.2,F6.1)/
     &'0hk0 ZONAL CONDITIONS:'/
     &'0hk0,h=2n+1                 ',I10,2(2X,F6.2,F6.1)/
     &' hk0,k=2n+1                 ',I10,2(2X,F6.2,F6.1)/
     &' hk0,h+k=2n+1               ',I10,2(2X,F6.2,F6.1)/
     &' hk0,h+k=4n+1, 2, or 3      ',I10,2(2X,F6.2,F6.1)/
     &'00kl ZONAL CONDITIONS:'/
     &'00kl,k=2n+1                 ',I10,2(2X,F6.2,F6.1)/
     &' 0kl,l=2n+1                 ',I10,2(2X,F6.2,F6.1)/
     &' 0kl,k+l=2n+1               ',I10,2(2X,F6.2,F6.1)/
     &' 0kl,k+l=4n+1, 2, or 3      ',I10,2(2X,F6.2,F6.1)/
     &'0h0l ZONAL CONDITIONS:'/
     &'0h0l,l=2n+1                 ',I10,2(2X,F6.2,F6.1)/
     &' h0l,h=2n+1                 ',I10,2(2X,F6.2,F6.1)/
     &' h0l,l+h=2n+1               ',I10,2(2X,F6.2,F6.1)/
     &' h0l,l+h=4n+1, 2, or 3      ',I10,2(2X,F6.2,F6.1)/
     &'0hhl, hkh, AND hkk ZONAL CONDITIONS:'/
     &'0hhl,l=2n+1                 ',I10,2(2X,F6.2,F6.1)/
     &' hhl,2h+l=4n+1, 2, or 3     ',I10,2(2X,F6.2,F6.1)/
     &'0hkh,k=2n+1                 ',I10,2(2X,F6.2,F6.1)/
     &' hkh,2h+k=4n+1, 2, or 3     ',I10,2(2X,F6.2,F6.1)/
     &'0hkk,h=2n+1                 ',I10,2(2X,F6.2,F6.1)/
     &' hkk,2k+h=4n+1, 2, or 3     ',I10,2(2X,F6.2,F6.1)/
     &'0h00, 0k0, AND 00l AXIAL CONDITIONS:'/
     &'0h00,h=2n+1                 ',I10,2(2X,F6.2,F6.1)/
     &' h00,h=4n+1, 2, or 3        ',I10,2(2X,F6.2,F6.1)/
     &'00k0,k=2n+1                 ',I10,2(2X,F6.2,F6.1)/
     &' 0k0,k=4n+1, 2, or 3        ',I10,2(2X,F6.2,F6.1)/
     &'000l,l=2n+1                 ',I10,2(2X,F6.2,F6.1)/
     &' 00l,l=4n+1, 2, or 3        ',I10,2(2X,F6.2,F6.1)/
     &' 00l,l=3n+1 or 2            ',I10,2(2X,F6.2,F6.1)/
     &' 00l,l=6n+1, 2, 3, 4, or 5  ',I10,2(2X,F6.2,F6.1))
      STOP 'PROGRAM THKL FINIS'
      END
C-----------------------------------------------------------------------
      SUBROUTINE DPRINT (II,IH,IK,IL,IPRINT)
C
C     IPRINT = 1 FOR H00, 0K0, 00L, HH0, 0KK, H0H, AND HHH.
C     IPRINT = 0 OTHERWISE.
C
      IF (II.LT.0) GO TO 10
      J=ABS(IH)
      K=ABS(IK)
      L=ABS(IL)
      IF (J.EQ.0.AND.K.EQ.0) GO TO 11
      IF (J.EQ.0.AND.L.EQ.0) GO TO 11
      IF (K.EQ.0.AND.L.EQ.0) GO TO 11
      IF (J.EQ.K.AND.L.EQ.0) GO TO 11
      IF (J.EQ.L.AND.K.EQ.0) GO TO 11
      IF (K.EQ.L.AND.J.EQ.0) GO TO 11
      IF (J.EQ.K.AND.K.EQ.L) GO TO 11
   10 IPRINT=0
      RETURN
   11 IPRINT=1
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE HKLTEST (H,K,L,Y,SIGY)
C
C     CONDITIONS LIMITING POSSIBLE REFLECTIONS
C
C     hkl,h+k=2n
C     hkl,k+l=2n
C     hkl,l+h=2n
C     hkl,h+k,k+l(,l+h)=2n
C     hkl,h+k+l=2n
C     hkl,-h+k+l=3n
C     hkl,h-k+l=3n
C     hkl,h-k=3n
C
C     hk0,h=2n
C     hk0,k=2n
C     hk0,h+k=2n
C     hk0,h+k=4n
C
C     0kl,k=2n
C     0kl,l=2n
C     0kl,k+l=2n
C     0kl,k+l=4n
C
C     h0l,l=2n
C     h0l,h=2n
C     h0l,l+h=2n
C     h0l,l+h=4n
C
C     hh(-2h)l,l=2n
C     h(-h)0l,l=2n
C     hhl,l=2n (rhombohedral axes)
C     hhl,l=2n
C
C     hkh,k=2n
C     hkk,h=2n
C     hhl,2h+l=4n
C     hkh,2h+k=4n
C     hkk,2k+h=4n
C
C     h00,h=2n
C     h00,h=4n
C
C     0k0,k=2n
C     0k0,k=4n
C
C     00l,l=2n
C     00l,l=4n
C     000l,l=2n
C     000l,l=3n
C     000l,l=6n
C
      INTEGER H
      COMMON /STATS/ CUTOFF,NA(34),NP(34),XA(34),XP(34)
      IF (SIGY.LE.0) RETURN
      T=Y/SIGY
C
C     GENERAL hkl
C
      IF (MOD(H+K,2).NE.0)    CALL TEST (1,T)
      IF (MOD(K+L,2).NE.0)    CALL TEST (2,T)
      IF (MOD(L+H,2).NE.0)    CALL TEST (3,T)
      IF (MOD(H+K,2).NE.0.OR.MOD(K+L,2).NE.0) CALL TEST (4,T)
      IF (MOD(H+K+L,2).NE.0)  CALL TEST (5,T)
      IF (MOD(-H+K+L,3).NE.0) CALL TEST (6,T)
      IF (MOD(H-K+L,3).NE.0)  CALL TEST (7,T)
      IF (MOD(H-K,3).NE.0)    CALL TEST (8,T)
C
C     ZONAL hk0, 0kl, AND h0l
C
      IF (L.EQ.0) THEN
        IF (MOD(H,2).NE.0)   CALL TEST ( 9,T)
        IF (MOD(K,2).NE.0)   CALL TEST (10,T)
        IF (MOD(H+K,2).NE.0) CALL TEST (11,T)
        IF (MOD(H+K,4).NE.0) CALL TEST (12,T)
      ELSE IF (H.EQ.0) THEN
        IF (MOD(K,2).NE.0)   CALL TEST (13,T)
        IF (MOD(L,2).NE.0)   CALL TEST (14,T)
        IF (MOD(K+L,2).NE.0) CALL TEST (15,T)
        IF (MOD(K+L,4).NE.0) CALL TEST (16,T)
      ELSE IF (K.EQ.0) THEN
        IF (MOD(L,2).NE.0)   CALL TEST (17,T)
        IF (MOD(H,2).NE.0)   CALL TEST (18,T)
        IF (MOD(L+H,2).NE.0) CALL TEST (19,T)
        IF (MOD(L+H,4).NE.0) CALL TEST (20,T)
      END IF
C
C     ZONAL hhl, hkh, AND hkk
C
      IF (ABS(K).EQ.ABS(H)) THEN
        IF (MOD(L,2).NE.0)     CALL TEST (21,T)
        IF (MOD(2*H+L,4).NE.0) CALL TEST (22,T)
      ELSE IF (ABS(L).EQ.ABS(H)) THEN
        IF (MOD(K,2).NE.0)     CALL TEST (23,T)
        IF (MOD(2*H+K,4).NE.0) CALL TEST (24,T)
      ELSE IF (ABS(L).EQ.ABS(K)) THEN
        IF (MOD(H,2).NE.0)     CALL TEST (25,T)
        IF (MOD(2*K+H,4).NE.0) CALL TEST (26,T)
      END IF
C
C     AXIAL h00, 0k0, AND 00l
C
      IF (K.EQ.0.AND.L.EQ.0) THEN
        IF (MOD(H,2).NE.0) CALL TEST (27,T)
        IF (MOD(H,4).NE.0) CALL TEST (28,T)
      ELSE IF (H.EQ.0.AND.L.EQ.0) THEN
        IF (MOD(K,2).NE.0) CALL TEST (29,T)
        IF (MOD(K,4).NE.0) CALL TEST (30,T)
      ELSE IF (H.EQ.0.AND.K.EQ.0) THEN
        IF (MOD(L,2).NE.0) CALL TEST (31,T)
        IF (MOD(L,4).NE.0) CALL TEST (32,T)
        IF (MOD(L,3).NE.0) CALL TEST (33,T)
        IF (MOD(L,6).NE.0) CALL TEST (34,T)
      END IF
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE METRIC (A,GA,GB)
C
C     CALCULATE RECIPROCAL LATTICE PARAMETERS AND DIRECT AND RECIPROCAL
C     LATTICE METRIC TENSORS FROM DIRECT LATTICE PARAMETERS.
C
      DIMENSION A(6),B(6),GA(3,3),GB(3,3)
      PI=ACOS(-1.0)
      DEG=180/PI
      RAD=PI/180
      CA=COS(A(4)*RAD)
      CB=COS(A(5)*RAD)
      CG=COS(A(6)*RAD)
      SA=SIN(A(4)*RAD)
      SB=SIN(A(5)*RAD)
      SG=SIN(A(6)*RAD)
      V=A(1)*A(2)*A(3)*SQRT(1-CA**2-CB**2-CG**2+2*CA*CB*CG)
      B(1)=A(2)*A(3)*SA/V
      B(2)=A(1)*A(3)*SB/V
      B(3)=A(1)*A(2)*SG/V
      B(4)=(CB*CG-CA)/(SB*SG)
      B(5)=(CA*CG-CB)/(SA*SG)
      B(6)=(CA*CB-CG)/(SA*SB)
      GA(1,1)=A(1)*A(1)
      GA(1,2)=A(1)*A(2)*CG
      GA(1,3)=A(1)*A(3)*CB
      GA(2,1)=GA(1,2)
      GA(2,2)=A(2)*A(2)
      GA(2,3)=A(2)*A(3)*CA
      GA(3,1)=GA(1,3)
      GA(3,2)=GA(2,3)
      GA(3,3)=A(3)*A(3)
      GB(1,1)=B(1)*B(1)
      GB(1,2)=B(1)*B(2)*B(6)
      GB(1,3)=B(1)*B(3)*B(5)
      GB(2,1)=GB(1,2)
      GB(2,2)=B(2)*B(2)
      GB(2,3)=B(2)*B(3)*B(4)
      GB(3,1)=GB(1,3)
      GB(3,2)=GB(2,3)
      GB(3,3)=B(3)*B(3)
      B(4)=ACOS(B(4))*DEG
      B(5)=ACOS(B(5))*DEG
      B(6)=ACOS(B(6))*DEG
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE NEWABC (A,T)
      DIMENSION A(6),T(3,3),X(3),Y(3),Z(3),G(3,3),GINV(3,3)
C
C     VECTORS X, Y, AND Z EQUAL NEW A, B, AND C REFERRED TO OLD.
C
      X(1)=T(1,1)
      X(2)=T(1,2)
      X(3)=T(1,3)
      Y(1)=T(2,1)
      Y(2)=T(2,2)
      Y(3)=T(2,3)
      Z(1)=T(3,1)
      Z(2)=T(3,2)
      Z(3)=T(3,3)
C
C     NEW A, B, C FROM X, Y, Z, AND OLD METRIC MATRIX G.
C
      CALL METRIC (A,G,GINV)
      A(1)=SQRT(UTGV(X,G,X))
      A(2)=SQRT(UTGV(Y,G,Y))
      A(3)=SQRT(UTGV(Z,G,Z))
C
C     NEW ALPHA, BETA, GAMMA FROM X, Y, Z, OLD G, AND NEW A, B, C.
C
      DEG=57.2957795
      A(4)=ACOS(UTGV(Y,G,Z)/(A(2)*A(3)))*DEG
      A(5)=ACOS(UTGV(X,G,Z)/(A(1)*A(3)))*DEG
      A(6)=ACOS(UTGV(X,G,Y)/(A(1)*A(2)))*DEG
      RETURN
      END
      SUBROUTINE NEWHKL (IH,IK,IL,T)
      DIMENSION T(3,3)
      JH=T(1,1)*IH+T(1,2)*IK+T(1,3)*IL
      JK=T(2,1)*IH+T(2,2)*IK+T(2,3)*IL
      JL=T(3,1)*IH+T(3,2)*IK+T(3,3)*IL
      IH=JH
      IK=JK
      IL=JL
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE READ1 (IFILE,IEND,II,IH,IK,IL,A1,A2,A3,A4,Y,SIGY,XTIME)
      IEND=0
c      READ (IFILE,END=9,ERR=9) II,IH,IK,IL,A1,A2,A3,A4,Y,SIGY,XTIME
      read (ifile,end=9,err=9) ih,ik,il,y,sigy
      RETURN
    9 IEND=1
      RETURN
      END
      FUNCTION SINTHL(IH,IK,IL,GINV)
      DIMENSION H(3),GINV(3,3)
      H(1)=IH
      H(2)=IK
      H(3)=IL
      Q=0
      DO 1 I=1,3
      DO 1 J=1,3
      Q=Q+H(J)*GINV(I,J)*H(I)
    1 CONTINUE
      SINTHL=0.5*SQRT(Q)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE TEST (I,T)
      COMMON /STATS/ CUTOFF,NA(34),NP(34),XA(34),XP(34)
      IF (T.LT.CUTOFF) THEN
        NA(I)=NA(I)+1
        XA(I)=XA(I)+T
      ELSE
        NP(I)=NP(I)+1
        XP(I)=XP(I)+T
      END IF
      RETURN
      END
C-----------------------------------------------------------------------
      FUNCTION UTGV(U,G,V)
C
C     BILINEAR FORM = ROW VECTOR U * SQUARE MATRIX G * COLUMN VECTOR V
C
      DIMENSION U(3),G(3,3),V(3)
      UTGV=0
      DO 1 I=1,3
      DO 1 J=1,3
      UTGV=UTGV+U(J)*G(I,J)*V(I)
    1 CONTINUE
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE WRITE1 (IFILE,II,IH,IK,IL,A1,A2,A3,A4,Y,SIGY,XTIME)
c      WRITE (IFILE) II,IH,IK,IL,A1,A2,A3,A4,Y,SIGY,XTIME
      write (ifile,'(3i5,2e15.7)') ih,ik,il,y,sigy
      RETURN
      END
C-----------------------------------------------------------------------
