      PROGRAM EVAL
C
C ROBERT H. BLESSING
C MEDICAL FOUNDATION OF BUFFALO
C 73 HIGH STREET
C BUFFALO, NEW YORK 14203, USA
C TELEPHONE:  (716) 856-9600, EXTENSION 335
C ELECTRONIC MAIL:  Blessing@HWI.Buffalo.Edu
C
C OCTOBER 1995,..., JANUARY 2001
C
      CHARACTER ATIME*8,ADATE*9,BTIME*8,BDATE*9,TITLE*80,FILE*80,FORM*11
      CHARACTER SYST*12,LATT*1,PTGP*5,LAUE*5,XRAY*2,EL*2
      DIMENSION CELL(6),ASTAR(6),GINV(3,3),B(3,3),C(3,3),COV(13,13),
     & U(3,3),V(3,3),AIAJ(6),VAR(13),H(3)
      DIMENSION EL(20),XI(20),ZI(20),A1(20),B1(20),A2(20),B2(20),A3(20),
     & B3(20),A4(20),B4(20),C0(20),FP(20),FPP(20)
      DATA B,C,COV /9*0,9*0,169*0/
      DATA NLIST,CUTOFF,EMIN,EMAX /100,3,1.5,0.5/
      DATA XMIN,XMAX,JMIN,JMAX,KMIN,KMAX,LMIN,LMAX /9,0,+99,-99,+99,-99,
     & +99,-99/
C
C FORTRAN-90 DYNAMIC STORAGE ALLOCATION
C
      REAL,    ALLOCATABLE::DATA(:)
      INTEGER, ALLOCATABLE::INDEX(:)
C
C FORTRAN-90 TIME AND DATE TRANSLATIONS
C
      CALL VTIME (ATIME)
      CALL VDATE (ADATE)
c      ATIME='hh:mm:ss'
c      ADATE='dd-mmm-yy'
      IO1=10
      IO2=20
      ILP=60
C
C READ CONTROL DATA AND PARAMETERS FROM PROGRAM ROGERS OR PROGRAM LEVY.
C
      OPEN (UNIT=IO1,FILE='eval.dat',STATUS='OLD')
      READ (IO1,999) BTIME
      READ (IO1,999) BDATE
      READ (IO1,999) TITLE
      READ (IO1,997) ITYPE
      READ (IO1,999) FILE
      READ (IO1,999) LATT
      READ (IO1,999) PTGP
      READ (IO1,998) CELL
      READ (IO1,999) XRAY
      READ (IO1,997) NEL
      DO 5 I=1,NEL
        READ (IO1,996) EL(I),XI(I),ZI(I),A1(I),B1(I),A2(I),B2(I)
        READ (IO1,995) A3(I),B3(I),A4(I),B4(I),C0(I),FP(I),FPP(I)
 5    CONTINUE
      READ (IO1,998) SMIN,SMAX
      READ (IO1,993) ASOLV,BSOLV
      READ (IO1,994) SCALEK
      READ (IO1,994)       B(1,1),B(1,2),B(1,3),B(2,2),B(2,3),B(3,3)
      READ (IO1,994,END=9) C(1,1),C(1,2),C(1,3),C(2,2),C(2,3),C(3,3)
      READ (IO1,993,END=9) ((COV(I,J),J=1,13),I=1,13)
      READ (IO1,992,END=9) NLIST,CUTOFF,EMAX,EMIN
 999  FORMAT (1X,A)
 998  FORMAT (1X,6F10.0)
 997  FORMAT (1X,I6)
 996  FORMAT (1X,A2,F10.2,F6.0,4(1X,F10.6))
 995  FORMAT (1X,5(1X,F10.6),2X,2(1X,F10.6))
 994  FORMAT (1X,6E12.5)
 993  FORMAT (1X,7E10.3)
 992  FORMAT (I10,3F10.0)
 9    CLOSE (UNIT=IO1,STATUS='KEEP')
      CALL DECODE (PTGP,LAUE,SYST,LATT,MLATT,IPTGP)
      CALL METRIC (CELL,ASTAR,GINV)
      PI=3.14159
      IF (B(2,2).EQ.0) THEN
        BISO=B(1,1)
        CISO=C(1,1)
        UISO=BISO/(8*PI**2)
        VISO=CISO/(8*PI**2)
      ELSE
        DO 6 I=1,3
        DO 6 J=I,3
          U(I,J)=B(I,J)/(2*PI**2*ASTAR(I)*ASTAR(J))
          V(I,J)=C(I,J)/(2*PI**2*ASTAR(I)*ASTAR(J))
          B(J,I)=B(I,J)
          C(J,I)=C(I,J)
          U(J,I)=U(I,J)
          V(J,I)=V(I,J)
 6      CONTINUE
      END IF
C
C PRINT CONTROL DATA.
C
      OPEN (UNIT=ILP,FILE='eval.lp',STATUS='NEW')
      WRITE (ILP,899) ATIME,ADATE,TITLE
      WRITE (ILP,898) BTIME,BDATE
      WRITE (ILP,897) FILE
      WRITE (ILP,896) SYST,LATT,LAUE,PTGP
      IF (PTGP.EQ.LAUE) THEN
        WRITE (ILP,8961)
      ELSE
        WRITE (ILP,8960)
      END IF
      WRITE (ILP,893) CELL
      WRITE (ILP,892) (EL(I),XI(I),I=1,NEL)
      WRITE (ILP,891) XRAY,(EL(I),FP(I),FPP(I),I=1,NEL)
      IF (SMIN.GT.0.OR.SMAX.LT.9) WRITE (ILP,895) SMIN,SMAX
      IF (ASOLV.NE.0) WRITE (ILP,613) ASOLV,BSOLV
      WRITE (ILP,899) ATIME,ADATE,TITLE
      IF (B(2,2).EQ.0) THEN
        WRITE (ILP,860) SCALEK,BISO,CISO
        DO 3 I=1,3
          VAR(I)=COV(I,I)
          IF (VAR(I).EQ.0) VAR(I)=1
 3      CONTINUE
        WRITE (ILP,861) (SQRT(COV(I,I)),I=1,3),(I,(COV(I,J)/SQRT(VAR(I)*
     &   VAR(J)),J=1,3),I=1,3)
      ELSE
        WRITE (ILP,889) SCALEK,
     &   B(1,1),B(1,2),B(1,3),B(2,2),B(2,3),B(3,3),
     &   C(1,1),C(1,2),C(1,3),C(2,2),C(2,3),C(3,3),
     &   U(1,1),U(1,2),U(1,3),U(2,2),U(2,3),U(3,3),
     &   V(1,1),V(1,2),V(1,3),V(2,2),V(2,3),V(3,3)
        AIAJ(1)=ASTAR(1)**2
        AIAJ(2)=ASTAR(1)*ASTAR(2)
        AIAJ(3)=ASTAR(1)*ASTAR(3)
        AIAJ(4)=ASTAR(2)**2
        AIAJ(5)=ASTAR(2)*ASTAR(3)
        AIAJ(6)=ASTAR(3)**2
        DO 13 I=1,13
          VAR(I)=COV(I,I)
          IF (VAR(I).EQ.0) VAR(I)=1
 13     CONTINUE
        WRITE (ILP,899) ATIME,ADATE,TITLE
        WRITE (ILP,888) (SQRT(COV(1+I,1+I))/(2*PI**2*AIAJ(I)),I=1,6),
     &   (SQRT(COV(7+I,7+I))/(2*PI**2*AIAJ(I)),I=1,6),
     &   (SQRT(COV(I,I)),I=1,13),
     &   (I,(NINT(100*COV(I,J)/SQRT(VAR(I)*VAR(J))),J=1,13),I=1,13)
      END IF
 899  FORMAT (/1H1,130('-')/1X,'PROGRAM EVAL.  ',A,', ',A,'.  ',A)
 898  FORMAT (/1X,
     &'SCALE AND THERMAL PARAMETERS FROM PROGRAM ROGERS OR PROGRAM LEV',
     &'Y JOB:'/1X,A,', ',A)
 897  FORMAT (/1X,'INPUT REFLECTION FILE:   ',A)
 896  FORMAT (/1X,A//1X,A,'-LATTICE'//1X,
     &'LAUE GROUP  ',A//1X,
     &'POINT GROUP ',A)
 8960 FORMAT (/1X,'NONCENTROSYMMETRIC STRUCTURE')
 8961 FORMAT (/1X,'CENTROSYMMETRIC STRUCTURE')
 893  FORMAT (/1X,
     &'LATTICE PARAMETERS'/1X,
     &'         A         B         C     ALPHA      BETA     GAMMA'/1X,
     & 3F10.4,3F10.3)
 892  FORMAT (/1X,'UNIT CELL CONTENTS:'//(1X,A2,F10.2))
 891  FORMAT (/1X,'RADIATION:  ',A2//1X,'ANOMALOUS DISPERSION CORRECTION
     & TERMS:'/1X,'--------------------------------------'/1X,
     &'EL        FP       FPP'/1X,'--        --       ---'/
     &(1X,A2,2F10.3))
 895  FORMAT (/1X,
     &'E-VALUES ARE CALCULATED FOR ALL DATA, BUT'/1X,
     &'E-STATISTICS ARE COMPILED ONLY FOR THE DATA WITH SMIN .LE. SIN(',
     &'THETA)/LAMBDA .LE. SMAX,'/1X,
     &'SMIN = ',F6.3/1X,'SMAX = ',F6.3,' RECIPROCAL ANGSTROMS,'/1X,
     &'WHICH WERE USED TO ESTIMATE SCALEK AND B(I,J).')
 613  FORMAT (/1X,
     &'BULK SOLVENT COMPENSATION PARAMETERS:'/1X,
     &'-------------------------------------'/1X,
     &'  E(PROT) = E(MEAS)/[1 - ASOLV*EXP{-BSOLV*[SIN(THETA)/LAMBDA]**',
     &'2}'/1X,
     &'  ASOLV = ',E10.3/1X,
     &'  BSOLV = ',E10.3' A**2 (SQUARED ANGSTROMS)')
 860  FORMAT (/1X,
     &'SCALEK, BISO, AND CISO ARE DEFINED ACCORDING TO:'//1X,
     &'  F(ABSOLUTE) = SCALEK*F(RELATIVE)'/1X,
     &'  f(T) = f(0)*<EXP(-B*S**2)>, WHERE S = SIN(THETA)/LAMBDA'/1X,
     &'  <EXP(-B*S**2)> = EXP(-BISO*S**2 + CISO**2*S**4)'//1X,
     &'PARAMETERS:'//1X,
     &'      SCALEK        BISO        CISO'/1X,3E12.4)
 861  FORMAT (/1X,
     &'PARAMETER E.S.D.''S:'//1X,
     &'    SIGMA(K)    SIGMA(B)    SIGMA(C)'/1X,3E12.2//1X,
     &'CORRELATION MATRIX:'//1X,
     &'       K  BISO  CISO'//1X,
     &'       1     2     3'//(1X,I2,3F6.2))
 889  FORMAT (/1X,
     &'SCALEK, B(I,J), C(I,J), U(I,J), AND B(I,J) ARE DEFINED ACCORDIN',
     &'G TO:'//1X,
     &'  F(ABSOLUTE) = SCALEK*F(RELATIVE)'/1X,
     &'  f(T) = f(0)*EXP(-SUM(I)SUM(J) H(I)*H(J)*B(I,J))'/1X,
     &'  f(T) = f(0)*EXP(-2*PI**2*SUM(I)SUM(J) H(I)*H(J)*ASTAR(I)*ASTA',
     &'R(J)*U(I,J))'/1X,
     &'  <EXP(-2*H(J)*B(I,J)*H(I))> = EXP(-2*(H(J)*<B(I,J)>*H(I) - (H(',
     &'J)*<C(I,J)>*H(I))**2))'//1X,
     &'SCALEK'/1X,E12.4//1X,
     &'OVERALL MEAN-SQUARE ATOMIC DISPLACEMENT PARAMETERS:'/1X,
     &'      B(1,1)      B(1,2)      B(1,3)      B(2,2)      B(2,3)   ',
     &'   B(3,3)'/1X,6E12.4,' (DIMENSIONLESS)'//1X,
     &'B-DISTRIBUTION DISPERSION PARAMETERS:'/1X,
     &'      C(1,1)      C(1,2)      C(1,3)      C(2,2)      C(2,3)   ',
     &'   C(3,3)'/' ',6E12.4,' (DIMENSIONLESS)'//1X,
     &'OVERALL MEAN-SQUARE ATOMIC DISPLACEMENTS:'/1X,
     &'      U(1,1)      U(1,2)      U(1,3)      U(2,2)      U(2,3)   ',
     &'   U(3,3)'/1X,6E12.4,' SQUARED ANGSTROMS'//1X,
     &'U-DISTRIBUTION DISPERSION PARAMETERS:'/1X,
     &'      V(1,1)      V(1,2)      V(1,3)      V(2,2)      V(2,3)   ',
     &'   V(3,3)'/1X,6E12.4,' SQUARED ANGSTROMS')
 888  FORMAT (/1X,
     &'PARAMETER E.S.D.''S'//1X,
     &'NOTE THAT THE ORDER OF THE MEAN-SQUARE DISPLACEMENT PARAMETER E',
     &'RRORS IS'/1X,
     &'11, 12, 13, 22, 23, 33  NOT  11, 22, 33, 12, 13, 23.'/1X,
     &'                        ---'//1X,
     &'  SIGMA(U11)  SIGMA(U12)  SIGMA(U13)  SIGMA(U22)  SIGMA(U23)  S',
     &'IGMA(U33)'/1X,6E12.2,' SQUARED ANGSTROMS'//1X,
     &'  SIGMA(V11)  SIGMA(V12)  SIGMA(V13)  SIGMA(V22)  SIGMA(V23)  S',
     &'IGMA(V33)'/1X,6E12.2,' SQUARED ANGSTROMS'//1X,
     &'    SIGMA(K)'/1X,E12.2//1X
     &'  SIGMA(B11)  SIGMA(B12)  SIGMA(B13)  SIGMA(B22)  SIGMA(B13)  S',
     &'IGMA(B33)'/1X,6E12.2,' (DIMENSIONLESS)'//1X
     &'  SIGMA(C11)  SIGMA(C12)  SIGMA(C13)  SIGMA(C12)  SIGMA(C23)  S',
     &'IGMA(C33)'/1X,6E12.2,' (DIMENSIONLESS)'//1X,
     &'CORRELATION MATRIX 100*CORR(I,J):'//1X,
     &'       K B11 B12 B13 B22 B23 B33 C11 C12 C13 C22 C23 C33'//1X,
     &'       1   2   3   4   5   6   7   8   9  10  11  12  13'/
     &/(1X,I2,2X,13I4))
C
C CALCULATE INITIAL E-VALUES.
C
      IF (ITYPE.LT.3) THEN
        FORM='FORMATTED'
      ELSE
        FORM='UNFORMATTED'
      END IF
      OPEN (UNIT=IO1,FILE=FILE,STATUS='OLD',FORM=FORM)
      OPEN (UNIT=IO2,STATUS='SCRATCH',FORM='FORMATTED')
      SUMESQ=0
      SUMM=0
      IEND=0
 10   CALL READ1 (ITYPE,IO1,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF)
      IF (IEND.NE.0) GO TO 19
      S=SINTHL(J,K,L,GINV)
      XMIN=MIN(XMIN,S)
      XMAX=MAX(XMAX,S)
      JMIN=MIN(JMIN,J)
      JMAX=MAX(JMAX,J)
      KMIN=MIN(KMIN,K)
      KMAX=MAX(KMAX,K)
      LMIN=MIN(LMIN,L)
      LMAX=MAX(LMAX,L)
      CALL FCALC (NEL,XI,A1,B1,A2,B2,A3,B3,A4,B4,C0,FP,FPP,S,SUMFSQ)
      IF (B(2,2).EQ.0) THEN
        HTBH=BISO*S**2
        HTCH=CISO*S**2
      ELSE
        H(1)=J
        H(2)=K
        H(3)=L
        HTBH=VTMV(H,B)
        HTCH=VTMV(H,C)
      END IF
      Q=(1/SCALEK)*EXP(-HTBH+HTCH**2)
      VARQ=VARF(SCALEK,B,C,COV,S,H,HTBH,HTCH,Q)
      CALL ESYM (IPTGP,J,K,L,EPSILON,M)
      EPSILON=MLATT*EPSILON
      E=F/(Q*SQRT(EPSILON*SUMFSQ))
      SIGE=E*SQRT((SIGF/F)**2+VARQ/Q**2)
      F=F*SCALEK
      SIGF=SIGF*SCALEK
      FSQ=FSQ*SCALEK**2
      SIGFSQ=SIGFSQ*SCALEK**2
      IF (ASOLV.NE.0) THEN
        Q=1/(1-ASOLV*EXP(-BSOLV*S**2))
        E=E*Q
        SIGE=SIGE*Q
      END IF
      CALL WRITE1 (IO2,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE)
      IF (SMIN.LE.S.AND.S.LE.SMAX) THEN
        CALL CSYM (IPTGP,J,K,L,IC)
        IF (IC.EQ.0) M=2*M
        SUMESQ=SUMESQ+M*E**2
        SUMM=SUMM+M
      END IF
      GO TO 10
 19   CONTINUE
      AVESQ=SUMESQ/SUMM
      RESCL=1/SQRT(AVESQ)
      ENDFILE IO2
      CLOSE (UNIT=IO1,STATUS='KEEP')
C
C RE-SCALE E VALUES SO THAT <E**2> = 1.
C
      WRITE (ILP,899) ATIME,ADATE,TITLE
      WRITE (ILP,878) XMIN,XMAX,JMIN,JMAX,KMIN,KMAX,LMIN,LMAX,XMAX,
     &1/(2*XMAX)
 878  FORMAT (/1X,
     &'SIN(THETA)/LAMBDA AND INDEX LIMITS:'//1X,
     &F7.4,' .LE. S .LE. ',F7.4,' ANGSTROM**(-1)'//1X,
     &I4,' .LE. H .LE. ',I4/1X,
     &I4,' .LE. K .LE. ',I4/1X,
     &I4,' .LE. L .LE. ',I4//1X,
     &'DATA RESOLUTION:'//1X,
     &'SMAX = ',F6.3,' ANGSTROM**(-1)'/1X,'DMIN = ',F5.2,'  ANGSTROM')
      WRITE (ILP,880) AVESQ,RESCL
      WRITE (ILP,887)
 880  FORMAT (/1X,
     &'MULTIPLICITY-WEIGHTED AVERAGE E**2   ',F12.4//1X,
     &'OVERALL RESCALING FACTOR FOR E-VALUES',F12.4)
 887  FORMAT (/1X,
     &'   H   K   L SIN(TH)/L    D(HKL)         FSQ      SIGFSQ       ',
     &'  F      SIGF         E      SIGE'/1X,
     &'   -   -   - ---------    ------         ---      ------       ',
     &'  -      ----         -      ----'/)
 886  FORMAT (1X,3I4,F10.4,F10.2,2F12.2,2F10.2,2F10.3)
      OPEN (UNIT=IO1,FILE='edata.hkl',STATUS='NEW',FORM='FORMATTED')
      REWIND IO2
      IEND=0
      NDATA=0
 20   CALL READ2 (IO2,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE)
      IF (IEND.NE.0) GO TO 29
      NDATA=NDATA+1
      E=E*RESCL
      SIGE=SIGE*RESCL
      CALL WRITE1 (IO1,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE)
      CALL DPRINT (J,K,L,IPRINT)
      IF (IPRINT.NE.0) THEN
        S=SINTHL(J,K,L,GINV)
        D=1/(2*S)
        WRITE (ILP,886) J,K,L,S,D,FSQ,SIGFSQ,F,SIGF,E,SIGE
      END IF
      GO TO 20
 29   CLOSE (UNIT=IO2,STATUS='DELETE')
      ENDFILE IO1
      WRITE (ILP,881) NDATA
 881  FORMAT (/1X,'NDATA = ',I6,' UNIQUE REFLECTION DATA')
C
C COMPILE E-STATISTICS.
C
      WRITE (ILP,899) ATIME,ADATE,TITLE
      WRITE (ILP,800) SMIN,MIN(SMAX,XMAX)
 800  FORMAT (/1X,
     &'E-STATISTICS FOR DATA WITH SMIN .LE. SIN(THETA)/LAMBDA .LE. SMA',
     &'X, WHERE'/1X,
     &'SMIN = ',F6.3,' RECIPROCAL ANGSTROMS'/1X,
     &'SMAX = ',F6.3)
      CALL ESTATS (IO1,ILP,SMIN,MIN(SMAX,XMAX),IPTGP,MLATT,GINV)
C
C ALLOCATE STORAGE FOR SORTING ARRAYS.
C
      ALLOCATE (DATA(NDATA),INDEX(NDATA))
C
C WRITE d(hkl)-SORTED OUTPUT FILE.
C
      REWIND IO1
      OPEN (UNIT=IO2,STATUS='SCRATCH',FORM='UNFORMATTED',
     & ACCESS='DIRECT',RECL=36)
      DO I=1,NDATA
        CALL READ2 (IO1,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE)
        WRITE (IO2,REC=I)    J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE
        DATA(I)=SINTHL(J,K,L,GINV)
      END DO
      CALL SORT (NDATA,DATA,INDEX)
      CLOSE (UNIT=IO1,STATUS='KEEP')
      OPEN (UNIT=IO1,FILE='edata.ddd',STATUS='NEW',FORM='FORMATTED')
      DO I=1,NDATA
        READ (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE
        CALL WRITE1 (IO1,       J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE)
      END DO
C
C LIST LOWEST RESOLUTION E-MAGNITUDES.
C
      WRITE (ILP,899) ATIME,ADATE,TITLE
      WRITE (ILP,'(/1X,''LOWEST RESOLUTION E-MAGNITUDES:'')')
      WRITE (ILP,877)
 877  FORMAT (/1X,
     &'   H   K   L SIN(TH)/L    D(HKL)         FSQ      SIGFSQ       ',
     &'  F      SIGF         E      SIGE      RANK'/1X,
     &'   -   -   - ---------    ------         ---      ------       ',
     &'  -      ----         -      ----      ----'/)
 876  FORMAT (1X,3I4,F10.4,F10.2,2F12.2,2F10.2,2F10.3,I10)
      NLIST=MIN(NLIST,NDATA)
      REWIND IO1
      DO I=1,NLIST
        CALL READ2 (IO1,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE)
        S=SINTHL(J,K,L,GINV)
        D=1/(2*S)
        WRITE (ILP,876) J,K,L,S,D,FSQ,SIGFSQ,F,SIGF,E,SIGE,I
      END DO
C
C WRITE E-SORTED OUTPUT FILE.
C
      REWIND IO1
      DO I=1,NDATA
        CALL READ2 (IO1,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE)
        WRITE (IO2,REC=I)    J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE
        DATA(I)=E
      END DO
      CALL SORT (NDATA,DATA,INDEX)
      CLOSE (UNIT=IO1,STATUS='KEEP')
      OPEN (UNIT=IO1,FILE='edata.eee',STATUS='NEW',FORM='FORMATTED')
      DO I=NDATA,1,-1
        READ (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE
        CALL WRITE1 (IO1,       J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE)
      END DO
      CLOSE (UNIT=IO1,STATUS='KEEP')
C
C LIST LARGEST SIGNIFICANT E-MAGNITUDES.
C
      WRITE (ILP,899) ATIME,ADATE,TITLE
      WRITE (ILP,'(/1X,
     &''LARGEST SIGNIFICANT E-MAGNITUDES:''//1X,
     &''E .GE. CUTOFF*SIGMA(E),''//1X,
     &''WHERE CUTOFF = '',F5.2)') CUTOFF
      WRITE (ILP,877)
      NLIST=MIN(NLIST,NDATA)
      N=0
      DO 31 I=NDATA,1,-1
        READ (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE
        IF (E.LT.CUTOFF*SIGE) GO TO 31      
        IF (E.LT.EMIN) GO TO 41
        S=SINTHL(J,K,L,GINV)
        D=1/(2*S)
        WRITE (ILP,876) J,K,L,S,D,FSQ,SIGFSQ,F,SIGF,E,SIGE,NDATA-I+1
        N=N+1
        IF (N.EQ.NLIST) GO TO 41
 31   CONTINUE
 41   CONTINUE
C
C LIST SMALLEST SIGNIFICANT E-MAGNITUDES.
C
      WRITE (ILP,899) ATIME,ADATE,TITLE
      WRITE (ILP,'(/1X,
     &''SMALLEST SIGNIFICANT E-MAGNITUDES:''//1X,
     &''E .GE. CUTOFF*SIGMA(E),''//1X,
     &''WHERE CUTOFF = '',F5.2)') CUTOFF
      WRITE (ILP,877)
      N=0
      DO 32 I=1,NDATA,+1
        READ (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE
        IF (E.LT.CUTOFF*SIGE) GO TO 32
        IF (E.GT.EMAX) GO TO 42
        S=SINTHL(J,K,L,GINV)
        D=1/(2*S)
        WRITE (ILP,876) J,K,L,S,D,FSQ,SIGFSQ,F,SIGF,E,SIGE,-I
        N=N+1
        IF (N.EQ.NLIST) GO TO 42
 32   CONTINUE
 42   CONTINUE
      CLOSE (UNIT=IO2,STATUS='DELETE')
      WRITE (ILP,879)
 879  FORMAT (/1X,
     &'OUTPUT E-FILES ARE FORMATTED, ASCII FILES (IH,IK,IL,FSQ,SIGFSQ,F,
     &SIGF,E,SIGE)'/1X,
     &'NAMED:'/1X,
     &'  "edata.hkl"  hkl-SORTED'/1X,
     &'  "edata.ddd"  d(hkl)-SORTED'/1X,
     &'  "edata.eee"  E(hkl)-SORTED')
      WRITE (6,879)
      STOP 'PROGRAM EVAL FINIS'
      END
C-----------------------------------------------------------------------
      SUBROUTINE CSYM (PTGP,H,K,L,C)
C
C TABLE OF CENTROSYMMETRIC ROWS AND ZONES IN NON-CENTROSYMMETRIC
C RECIPROCAL LATTICE POINT GROUPS
C
C TRICL. MONOCL. ORTHO.  TETRAG.    TRIG.     RHOMB.    CUB.
C
C (1) 1  (3) 2   (6) 222  (9) 4     (17)  3   (33) R3   (38) 23
C (2) -1 (4) M   (7) MM2 (10) -4    (18)  -3  (34) R-3  (39) M3
C        (5) 2/M (8) MMM (11) 4/M
C                                   (19) 312  (35) R32  (40) 432
C                        (12) 422   (20) 321  (36) R3M  (41) -43M
C                        (13) 4MM   (21) 31M  (37) R-3M (42) M3M
C                        (14) -42M  (22) 3M1
C                        (15) -4M2  (23) -31M
C                        (16) 4/MMM (24) -3M1
C
C                                   HEXAG.
C
C                                   (25) 6
C                                   (26) -6
C                                   (27) 6/M
C
C                                   (28) 622
C                                   (29) 6MM
C                                   (30) -6M2
C                                   (31) -62M
C                                   (32) 6/MMM
C
C A TOTAL OF 42 CASES IS TABULATED BECAUSE THERE ARE TEN CASES OF
C ALTERNATIVE AXES FOR SEVEN OF THE 32 CRYSTALLOGRAPHIC POINT GROUPS:
C
C -42M  -4M2
C 3     R3
C -3    R-3
C 312   321   R32
C 31M   3M1   R3M
C -31M  -3M1  R-3M
C -62M  -6M2
C
C C = 0 FOR ACENTRIC DISTRIBUTIONS
C C = 1 FOR CENTRIC DISTRIBUTIONS
C
      INTEGER PTGP,H,C
      C=0
      GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,
     & 23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42) PTGP
 1    RETURN
 2    C=1
      RETURN
 3    IF (K.EQ.0) C=1
      RETURN
 4    IF (H.EQ.0.AND.L.EQ.0) C=1
      RETURN
 5    C=1
      RETURN
 6    IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) C=1
      RETURN
 7    IF (L.EQ.0) C=1
      RETURN
 8    C=1
      RETURN
 9    IF (L.EQ.0) C=1
      RETURN
 10   IF (L.EQ.0) C=1
      IF (H.EQ.0.AND.K.EQ.0) C=1
      RETURN
 11   C=1
      RETURN
 12   IF (L.EQ.0.OR.K.EQ.0.OR.H.EQ.0) C=1
      IF (ABS(K).EQ.ABS(H)) C=1
      RETURN
 13   IF (L.EQ.0) C=1
      RETURN
 14   IF (L.EQ.0.OR.K.EQ.0.OR.H.EQ.0) C=1
      RETURN
 15   IF (L.EQ.0) C=1
      IF (ABS(K).EQ.ABS(H)) C=1
      RETURN
 16   C=1
      RETURN
 17   RETURN
 18   C=1
      RETURN
 19   IF (K.EQ.H.OR.K.EQ.-2*H) C=1
      RETURN
 20   IF (K.EQ.0.OR.H.EQ.0) C=1
      RETURN
 21   IF ((K.EQ.0.AND.L.EQ.0).OR.(H.EQ.0.AND.L.EQ.0)) C=1
      RETURN
 22   RETURN
 23   C=1
      RETURN
 24   C=1
      RETURN
 25   IF (L.EQ.0) C=1
      RETURN
 26   IF (H.EQ.0.AND.K.EQ.0) C=1
      RETURN
 27   C=1
      RETURN
 28   IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) C=1
      IF (K.EQ.H.OR.K.EQ.-2*H) C=1
      RETURN
 29   IF (L.EQ.0) C=1
      RETURN
 30   IF (K.EQ.H.OR.K.EQ.-2*H) C=1
      RETURN
 31   IF (H.EQ.0.OR.K.EQ.0) C=1
      RETURN
 32   C=1
      RETURN
 33   RETURN
 34   C=1
      RETURN
 35   IF (H.EQ.K.OR.K.EQ.L.OR.L.EQ.H) C=1
      RETURN
 36   IF (L.EQ.0.AND.K.EQ.-H) C=1
      IF (K.EQ.0.AND.L.EQ.-H) C=1
      IF (H.EQ.0.AND.L.EQ.-K) C=1
      RETURN
 37   C=1
      RETURN
 38   IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) C=1
      RETURN
 39   C=1
      RETURN
 40   IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) C=1
      IF (ABS(H).EQ.ABS(K).OR.ABS(K).EQ.ABS(L).OR.ABS(L).EQ.ABS(H)) C=1
      RETURN
 41   IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) C=1
      RETURN
 42   C=1
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE DECODE (PTGP,LAUE,SYST,LATT,MLATT,IPTGP)
C
C CRYSTAL SYSTEM, LATTICE TYPE AND MULTIPLICITY, AND LAUE GROUP AND
C POINT GROUP SYMMETRIES
C
C THERE ARE TEN CASES OF ALTERNATIVE AXES FOR SEVEN OF THE 32
C CRYSTALLOGRAPHIC POINT GROUPS (SEE SUBROUTINE ESYM), THUS A TOTAL OF
C 42 CASES IS TABULATED.
C
C SIMILARLY, THERE ARE THREE CASES OF ALTERNATIVE AXES FOR TWO OF THE 11
C LAUE GROUPS, SO 14 CASES ARE TABULATED.
C
      CHARACTER SYST*12,LATT*1,PTGP*5,LAUE*5,SYMBOL(42)*5
      DATA SYMBOL
     &/'1    ','-1   ',
     & '2    ','M    ','2/M  ',
     & '222  ','MM2  ','MMM  ',
     & '4    ','-4   ','4/M  ',
     & '422  ','4MM  ','-42M ','-4M2 ','4/MMM',
     & '3    ','-3   ',
     & '312  ','321  ','31M  ','3M1  ','-31M ','-3M1 ',
     & '6    ','-6   ','6/M  ',
     & '622  ','6MM  ','-6M2 ','-62M ','6/MMM',
     & 'R3   ','R-3  ',
     & 'R32  ','R3M  ','R-3M ',
     & '23   ','M3   ',
     & '432  ','-43M ','M3M  '/
C
C POINT GROUP INDEX NUMBERS:
C
C TRICLINIC     1,  2;
C MONOCLINIC    3,  4,  5;
C ORTHORHOMBIC  6,  7,  8;
C TETRAGONAL    9, 10, 11,
C              12, 13, 14, 15, 16;
C TRIGONAL     17, 18,
C              19, 20, 21, 22, 23, 24;
C HEXAGONAL    25, 26, 27,
C              28, 29, 30, 31, 32;
C RHOMBOHEDRAL 33, 34,
C              35, 36, 37;
C CUBIC        38, 39,
C              40, 41, 42.
C 
      DO 1 I=1,42
        IF (PTGP.EQ.SYMBOL(I)) GO TO 2
 1    CONTINUE
      STOP 'EVAL/DECODE:  ERROR IN POINT GROUP SYMBOL'
 2    IPTGP=I
C
C LAUE GROUP
C
      IF (I.GE. 1) LAUE='-1   '
      IF (I.GE. 3) LAUE='2/M  '
      IF (I.GE. 6) LAUE='MMM  '
      IF (I.GE. 9) LAUE='4/M  '
      IF (I.GE.12) LAUE='4/MMM'
      IF (I.GE.17) LAUE='-3   '
      IF (I.GE.19) LAUE='-31M '
      IF (I.EQ.20.OR.I.EQ.22.OR.I.EQ.24) LAUE='-3M1 '
      IF (I.GE.25) LAUE='6/M  '
      IF (I.GE.28) LAUE='6/MMM'
      IF (I.GE.33) LAUE='R-3  '
      IF (I.GE.35) LAUE='R-3M '
      IF (I.GE.38) LAUE='M3   '
      IF (I.GE.40) LAUE='M3M  '
C
C CRYSTAL SYSTEM
C
      IF (I.GE. 1) SYST='TRICLINIC'
      IF (I.GE. 3) SYST='MONOCLINIC'
      IF (I.GE. 6) SYST='ORTHORHOMBIC'
      IF (I.GE. 9) SYST='TETRAGONAL'
      IF (I.GE.17) SYST='TRIGONAL'
      IF (I.GE.25) SYST='HEXAGONAL'
      IF (I.GE.33) SYST='RHOMBOHEDRAL'
      IF (I.GE.38) SYST='CUBIC'
C
C LATTICE CENTERING MULTIPLICITY
C
      IF (LATT.EQ.'P') MLATT=1
      IF (LATT.EQ.'A') MLATT=2
      IF (LATT.EQ.'B') MLATT=2
      IF (LATT.EQ.'C') MLATT=2
      IF (LATT.EQ.'F') MLATT=4
      IF (LATT.EQ.'I') MLATT=2
      IF (LATT.EQ.'R') MLATT=1
      IF (LATT.EQ.'R'.AND.SYST.NE.'RHOMBOHEDRAL') MLATT=3
      RETURN
      END
C-----------------------------------------------------------------------
      FUNCTION DELTA(I,J)
C
C KRONECKER DELTA
C
      DELTA=0
      IF (I.EQ.J) DELTA=1
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE DPRINT (IH,IK,IL,IPRINT)
C
C DECIDE TO PRINT OR NOT.
C
C IPRINT = 0  DO NOT
C        = 1  DO PRINT REFLECIOTN
C
C PRINT SPECIAL REFLECTIONS OF THE TYPE H00, 0K0, 00L, HH0, H0H,0KK, AND
C HHH.
C
      IPRINT=0
      J=ABS(IH)
      K=ABS(IK)
      L=ABS(IL)
      IF (J.EQ.0.AND.K.EQ.0) GO TO 1
      IF (J.EQ.0.AND.L.EQ.0) GO TO 1
      IF (K.EQ.0.AND.L.EQ.0) GO TO 1
      IF (J.EQ.K.AND.L.EQ.0) GO TO 1
      IF (J.EQ.L.AND.K.EQ.0) GO TO 1
      IF (K.EQ.L.AND.J.EQ.0) GO TO 1
      IF (J.EQ.K.AND.K.EQ.L) GO TO 1
      RETURN
 1    IPRINT=1
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE ESTATS (IFILE,ILP,SMIN,SMAX,IPTGP,MLATT,GINV)
C
C E-STATISTICS                  THEORETICAL MOMENTS
C                           ACENTRIC CENTRIC HYPERCENTRIC
C
C A1 = <E>                     0.866   0.798   0.718
C A2 = <E**2>                  1.000   1.000   1.000 
C A3 = <E**3>                  1.329   1.596   1.916 
C A4 = <E**4>                  2.000   3.000   4.500 
C A5 = <E**5>                  3.323   6.383  12.260 
C A6 = <E**6>                  6.000  15.000  37.500 
C B1 = <ABS(E**2 - 1)>         0.736   0.968   1.145 
C B2 = <(E**2 - 1)**2>         1.000   2.000   3.500 
C B3 = <(E**2 - 1)**3>         2.000   8.000  26.000 
C C3 = <ABS(E**2 - 1)**3>      2.415   8.691  26.903
C
C F1  = FRACTION E .GE. 1      0.368   0.320
C F2                    2      0.018   0.050
C F3                    3      0.0001  0.003
C
C F05 = FRACTION E .LT. 0.5    0.22    0.38
C F15 =            .GT. 1.5    0.105   0.13
C
      DIMENSION F1(3),F2(3),F3(3),NI(7),A1(7),A2(7),A3(7),A4(7),A5(7),
     & A6(7),B1(7),B2(7),B3(7),C3(7),GINV(3,3),H(3),SI(10),NJ(10),
     & SJ(10),E1(10),E2(10),F05(10),F15(10)
      DATA F1,F2,F3,NI,A1,A2,A3,A4,A5,A6,B1,B2,B3,C3,NJ,SJ,E1,E2,F05,F15
     & /146*0/
      DIMENSION W(50),X(50),Y(50),Z(50)
C
C MOMENTS OF THE E-DISTRIBUTION
C
      REWIND IFILE
      IEND=0
      NREC=0
      NA=0
      NC=0
 1    CALL READ2 (IFILE,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE)
      IF (IEND.NE.0) GO TO 9
      NREC=NREC+1
      S=SINTHL(J,K,L,GINV)
      IF (SMIN.LE.S.AND.S.LE.SMAX) THEN
        CALL ESYM (IPTGP,J,K,L,EPSILON,M)
        CALL CSYM (IPTGP,J,K,L,IC)
        IF (IC.EQ.0) M=2*M
        DO 10 I=1,7
          IF (I.EQ.2.AND.IC.NE.0) GO TO 10
          IF (I.EQ.3.AND.IC.NE.1) GO TO 10
          IF (I.EQ.4.AND.(J.EQ.0.OR.K.EQ.0.OR.L.EQ.0)) GO TO 10
          IF (I.EQ.5.AND.J.NE.0) GO TO 10
          IF (I.EQ.6.AND.K.NE.0) GO TO 10
          IF (I.EQ.7.AND.L.NE.0) GO TO 10
          NI(I)=NI(I)+M
          A1(I)=A1(I)+M*E
          A2(I)=A2(I)+M*E**2
          A3(I)=A3(I)+M*E**3
          A4(I)=A4(I)+M*E**4
          A5(I)=A5(I)+M*E**5
          A6(I)=A6(I)+M*E**6
          B1(I)=B1(I)+M*ABS(E**2-1)
          B2(I)=B2(I)+M*(E**2-1)**2
          B3(I)=B3(I)+M*(E**2-1)**3
          C3(I)=C3(I)+M*ABS(E**2-1)**3
          IF (I.LE.3) THEN
            IF (E.GE.1) F1(I)=F1(I)+M
            IF (E.GE.2) F2(I)=F2(I)+M
            IF (E.GE.3) F3(I)=F3(I)+M
            IF (E.LT.0.5) F05(I)=F05(I)+M
            IF (E.GT.1.5) F15(I)=F15(I)+M
          END IF
 10     CONTINUE
        IF (IC.EQ.0) NA=NA+1
        IF (IC.EQ.1) NC=NC+1
      END IF
      GO TO 1
 9    CONTINUE
      DO I=1,7
        IF (NI(I).GT.0) THEN
          A1(I)=A1(I)/NI(I)
          A2(I)=A2(I)/NI(I)
          A3(I)=A3(I)/NI(I)
          A4(I)=A4(I)/NI(I)
          A5(I)=A5(I)/NI(I)
          A6(I)=A6(I)/NI(I)
          B1(I)=B1(I)/NI(I)
          B2(I)=B2(I)/NI(I)
          B3(I)=B3(I)/NI(I)
          C3(I)=C3(I)/NI(I)
        END IF
      END DO
      DO I=1,3
        IF (NI(I).GT.0) THEN
          F1(I)=F1(I)/NI(I)
          F2(I)=F2(I)/NI(I)
          F3(I)=F3(I)/NI(I)
          F05(I)=F05(I)/NI(I)
          F15(I)=F15(I)/NI(I)
        END IF
      END DO
      WRITE (ILP,900)
      WRITE (ILP,901) A1, 0.866,  0.798,  0.718
      WRITE (ILP,902) A2, 1.000,  1.000,  1.000
      WRITE (ILP,903) A3, 1.329,  1.596,  1.916
      WRITE (ILP,904) A4, 2.000,  3.000,  4.500
      WRITE (ILP,905) A5, 3.323,  6.383, 12.260
      WRITE (ILP,906) A6, 6.000, 15.000, 37.500
      WRITE (ILP,907) B1, 0.736,  0.968,  1.145
      WRITE (ILP,908) B2, 1.000,  2.000,  3.500
      WRITE (ILP,909) B3, 2.000,  8.000, 26.000
      WRITE (ILP,910) C3, 2.415,  8.691, 26.903
      WRITE (ILP,911) NI
      WRITE (ILP,912) F1, 0.368,  0.320,
     &                F2, 0.018,  0.050,
     &                F3, 0.0001, 0.003,
     &                (F05(I),I=1,3), 0.22,  0.38,
     &                (F15(1),I=1,3), 0.105, 0.13
 900  FORMAT (/1X,
     &'                                           EXPERIMENTAL DATA AV',
     &'ERAGES                            THEORETICAL MOMENTS'/1X,
     &'                      -----------------------------------------',
     &'---------------------------  -------------------------------'/1X,
     &'AVERAGE               ALL DATA  ACENTRIC   CENTRIC       HKL   ',
     &'    0KL       H0L       HK0  ACENTRIC   CENTRIC HYPERCENTRIC'/1X,
     &'                      --------  --------   -------       ---   ',
     &'    ---       ---       ---  --------   ------- ------------')
 901  FORMAT (' <E>                 ',10F10.3)
 902  FORMAT (' <E**2>              ',10F10.3)
 903  FORMAT (' <E**3>              ',10F10.3)
 904  FORMAT (' <E**4>              ',10F10.3)
 905  FORMAT (' <E**5>              ',10F10.3)
 906  FORMAT (' <E**6>              ',10F10.3)
 907  FORMAT (' <ABS(E**2 - 1)>     ',10F10.3)
 908  FORMAT (' <(E**2 - 1)**2>     ',10F10.3)
 909  FORMAT (' <(E**2 - 1)**3>     ',10F10.3)
 910  FORMAT (' <ABS(E**2 -1)**3>   ',10F10.3)
 911  FORMAT ('                       --------------------------------',
     &'------------------------------------  -------------------------',
     &'------'/1X,'MULT.-WTD. NO.      ',7I10)
 912  FORMAT (/1X,
     &'                                        EXPERIMENTAL           ',
     &'   THEORETICAL'/1X,
     &'                                ----------------------------  -',
     &'-----------------'/1X,
     &'                                ALL DATA  ACENTRIC   CENTRIC  A',
     &'CENTRIC   CENTRIC'/1X,
     &'                                --------  --------   -------  -',
     &'-------   -------'/1X,
     &'FRACTION OF DATA WITH E .GE. 1',5F10.3/1X,
     &'                             2',5F10.3/1X,
     &'                             3',F11.4,F10.4,F9.3,F11.4,F9.3/1X,
     &'                                ----------------------------  ',
     &'------------------'/1X,
     &'FRACTION OF DATA WITH E .LT. 0.5',F7.2,1X,4(F9.2,1X)/1X,
     &'                      E .GT. 1.5',F7.2,1X,2(F9.2,1X),F10.3,F9.2/
     &1X,
     &'                                ----------------------------  -',
     &'-----------------')
C
C E-STATISTICS IN EQUAL-VOLUME SHELLS OF SIN(THETA)/LAMBDA
C
      N=10
      DO I=1,N
        SI(I)=(FLOAT(I)/N)**(1/3.0)*SMAX
      END DO
      REWIND IFILE
      DO 2 IREC=1,NREC
        CALL READ2 (IFILE,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE)
        S=SINTHL(J,K,L,GINV)
        IF (SMIN.LE.S.AND.S.LE.SMAX) THEN
          CALL ESYM (IPTGP,J,K,L,EPSILON,M)
          CALL CSYM (IPTGP,J,K,L,IC)
          IF (IC.EQ.0) M=2*M
          DO I=1,N
            IF (S.LE.SI(I)) THEN
              NJ(I)=NJ(I)+M
              SJ(I)=SJ(I)+M*S
              E1(I)=E1(I)+M*E
              E2(I)=E2(I)+M*E**2
              IF (E.LT.0.5) F05(I)=F05(I)+M
              IF (E.GT.1.5) F15(I)=F15(I)+M
              GO TO 2
            END IF
          END DO
        END IF
 2    CONTINUE
      DO I=1,N
        IF (NJ(I).GT.0) THEN
          SJ(I)=SJ(I)/NJ(I)
          E1(I)=E1(I)/NJ(I)
          E2(I)=E2(I)/NJ(I)
          F05(I)=F05(I)/NJ(I)
          F15(I)=F15(I)/NJ(I)
        END IF
      END DO
      WRITE (ILP,920) (2*SJ(I),E2(I),E1(I),F05(I),F15(I),NJ(I),I=1,N)
 920  FORMAT (/1X,
     &'E-STATISTICS IN EQUAL-VOLUME SHELLS OF DSTAR = 2*SIN(THETA)/LAM',
     &'BDA'//1X,
     &'   <DSTAR>    <E**2>       <E> F(E < 0.5) F(E > 1.5) MULT.-WTD.',
     &' NO.'/1X,
     &'   -------    ------       --- ---------- ---------- ----------',
     &'----'/(1X,3F10.3,2F10.2,I10))
C
C E-STATISTICS IN INDEX PARITY GROUPS
C
      DO I=1,8
        NJ(I)=0
        E1(I)=0
        E2(I)=0
        F05(I)=0
        F15(I)=0
      END DO
      REWIND IFILE
      DO IREC=1,NREC
        CALL READ2 (IFILE,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE)
        S=SINTHL(J,K,L,GINV)
        IF (SMIN.LE.S.AND.S.LE.SMAX) THEN
          CALL ESYM (IPTGP,J,K,L,EPSILON,M)
          CALL CSYM (IPTGP,J,K,L,IC)
          IF (IC.EQ.0) M=2*M
          CALL PARITY (I,J,K,L)
          NJ(I)=NJ(I)+M
          E1(I)=E1(I)+M*E
          E2(I)=E2(I)+M*E**2
          IF (E.LT.0.5) F05(I)=F05(I)+M
          IF (E.GT.1.5) F15(I)=F15(I)+M
        END IF
      END DO
      DO I=1,8
        IF (NJ(I).GT.0) THEN
          E1(I)=E1(I)/NJ(I)
          E2(I)=E2(I)/NJ(I)
          F05(I)=F05(I)/NJ(I)
          F15(I)=F15(I)/NJ(I)
        END IF
      END DO
      WRITE (ILP,931) (E2(I),E1(I),F05(I),F15(I),NJ(I),I=1,8)
 931  FORMAT (/1X,
     &'E-STATISTICS FOR INDEX PARITY GROUPS'//1X,
     &'    PARITY    <E**2>       <E> F(E < 0.5) F(E > 1.5) MULT.-WTD.',
     &' NO.'/1X,
     &'    ------    ------       --- ---------- ---------- ----------',
     &'----'/
     &'        EEE',2F10.3,2F10.2,I10/
     &'        EEO',2F10.3,2F10.2,I10/
     &'        EOE',2F10.3,2F10.2,I10/
     &'        OEE',2F10.3,2F10.2,I10/
     &'        EOO',2F10.3,2F10.2,I10/
     &'        OEO',2F10.3,2F10.2,I10/
     &'        OOE',2F10.3,2F10.2,I10/
     &'        OOO',2F10.3,2F10.2,I10)
C
C PLOT THE E-DISTRIBUTIONS P(E) VS. E, WHERE 0 < P < 1 AND 0 < E < 5.0.
C
      IF (NA.GT.0) THEN
C
C ACENTRIC DISTRIBUTION
C
        N=MIN(50,NINT(0.1*NA))
        DO I=1,N
          X(I)=0
          Y(I)=0
        END DO
        DX=5.0/N
        REWIND IFILE
        DO 20 IREC=1,NREC
          CALL READ2 (IFILE,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE)
          S=SINTHL(J,K,L,GINV)
          IF (S.LT.SMIN.OR.S.GT.SMAX) GO TO 20
          CALL CSYM (IPTGP,J,K,L,IC)
          IF (IC.EQ.1) GO TO 20
          CALL ESYM (IPTGP,J,K,L,EPSILON,M)
          M=2*M
          DO I=1,N
            IF (E.LE.I*DX) THEN
              X(I)=X(I)+M*E
              Y(I)=Y(I)+M
              GO TO 20
            END IF
          END DO
 20     CONTINUE
C
C EVALUATE BIN AVERAGES, AND ELIMINATE ANY EMPTY BINS.
C
        M=N
        N=0
        DO I=1,M
          IF (Y(I).GT.0) THEN
            N=N+1
            X(N)=X(I)/Y(I)
            Y(N)=Y(I)
          END IF
        END DO
C
C ENSURE INCREASING X VALUES.
C
        M=N
        N=1
        DO I=2,M
          IF (X(I).GT.X(I-1)) THEN
            N=N+1
            X(N)=X(I)
            Y(N)=Y(I)
          END IF
        END DO
C
C NORMALIZE BIN AREAS.
C
        A=Y(1)*(X(2)-X(1))+Y(N)*(X(N)-X(N-1))
        DO I=2,N-1
          A=A+0.5*Y(I)*(X(I+1)-X(I-1))
        END DO
        DO I=1,N
          Y(I)=Y(I)/A
        END DO
        WRITE (ILP,600)
        CALL PLOTA (N,X,Y,ILP)
      END IF
 600  FORMAT (/1H1,130('-')/1X,10X,
     &'PROGRAM EVAL.  ACENTRIC E-DISTRIBUTION'/)
      IF (NC.GT.0) THEN
C
C CENTRIC DISTRIBUTION
C
        N=MIN(50,NINT(0.1*NC))
        DO I=1,N
          X(I)=0
          Y(I)=0
        END DO
        DX=5.0/N
        REWIND IFILE
        DO 21 IREC=1,NREC
          CALL READ2 (IFILE,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE)
          S=SINTHL(J,K,L,GINV)
          IF (S.LT.SMIN.OR.S.GT.SMAX) GO TO 21
          CALL ESYM (IPTGP,J,K,L,EPSILON,M)
          CALL CSYM (IPTGP,J,K,L,IC)
          IF (IC.EQ.0) GO TO 21
          DO I=1,N
            IF (E.LE.I*DX) THEN
              X(I)=X(I)+M*E
              Y(I)=Y(I)+M
              GO TO 21
            END IF
          END DO
 21     CONTINUE
C
C EVALUATE BIN AVERAGES, AND ELIMINATE ANY EMPTY BINS.
C
        M=N
        N=0
        DO I=1,M
          IF (Y(I).GT.0) THEN
            N=N+1
            X(N)=X(I)/Y(I)
            Y(N)=Y(I)
          END IF
        END DO
C
C ENSURE INCREASING X VALUES.
C
        M=N
        N=1
        DO I=2,M
          IF (X(I).GT.X(I-1)) THEN
            N=N+1
            X(N)=X(I)
            Y(N)=Y(I)
          END IF
        END DO
C
C NORMALIZE BIN AREAS.
C
        A=Y(1)*(X(2)-X(1))+Y(N)*(X(N)-X(N-1))
        DO I=2,N-1
          A=A+0.5*Y(I)*(X(I+1)-X(I-1))
        END DO
        DO I=1,N
          Y(I)=Y(I)/A
        END DO
        WRITE (ILP,610)
        CALL PLOTC (N,X,Y,ILP)
      END IF
 610  FORMAT (/1H1,130('-')/1X,10X,
     &'PROGRAM EVAL.  CENTRIC E-DISTRIBUTION'/)
C
C PLOT <E**2> VERSUS <DSTAR> = 2*<SIN(THETA)/LAMBDA>.
C
      N=MIN(50,NINT(FLOAT(NREC)/50))
      DO I=1,N
        W(I)=(FLOAT(I)/N)**(1/3.0)*SMAX
      END DO
      DO I=1,N
        X(I)=0
        Y(I)=0
        Z(I)=0
      END DO
      REWIND IFILE
      DO 3 IREC=1,NREC
        CALL READ2 (IFILE,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE)
        S=SINTHL(J,K,L,GINV)
        IF (S.LT.SMIN.OR.S.GT.SMAX) GO TO 3
        CALL ESYM (IPTGP,J,K,L,EPSILON,M)
        CALL CSYM (IPTGP,J,K,L,IC)
        IF (IC.EQ.0) M=2*M
        DO I=1,N
          IF (S.LE.W(I)) THEN
            X(I)=X(I)+M*S
            Y(I)=Y(I)+M*E**2
            Z(I)=Z(I)+M
            GO TO 3
          END IF
        END DO
 3    CONTINUE
C
C EVALUATE BIN AVERAGES, AND ELIMINATE ANY EMPTY BINS.
C
      M=N
      N=0
      DO I=1,M
        IF (Z(I).GT.0) THEN
          N=N+1
          X(N)=X(I)/Z(I)
          Y(N)=Y(I)/Z(I)
        END IF
      END DO
C
C ENSURE INCREASING X VALUES.
C
      M=N
      N=1
      DO I=2,M
        IF (X(I).GT.X(I-1)) THEN
          N=N+1
          X(N)=X(I)
          Y(N)=Y(I)
        END IF
      END DO
      WRITE (ILP,630)
      CALL PLOTS (N,X,Y,SMAX,ILP)
 630  FORMAT (/1H1,130('-')/1X,10X,
     &'PROGRAM EVAL.  <E**2> VERSUS <DSTAR> = 2*<SIN(THETA)/LAMBDA> DI',
     &'STRIBUTION'/)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE ESYM (PTGP,H,K,L,E,M)
C
C RECIPROCAL LATTICE POINT DEGENERACIES, E = EPSILON(HKL), AND
C MULTIPLICITIES, M = M(HKL), FOR THE CRYSTALLOGRAPHIC POINT GROUPS
C
C AS TABULATED BY:
C
C HITOSHI IWASAKI AND TETSUZO ITO (1977).  VALUES OF EPSILON FOR
C OBTAINING NORMALIZED STRUCTURE FACTORS.  ACTA CRYST. A33, 227-229.
C
C TRICL. MONOCL. ORTHO.  TETRAG.    TRIG.     RHOMB.    CUB.
C
C (1) 1  (3) 2   (6) 222 (9)  4     (17) 3    (33) R3   (38) 23
C (2) -1 (4) M   (7) MM2 (10) -4    (18) -3   (34) R-3  (39) M3
C        (5) 2/M (8) MMM (11) 4/M   
C                                   (19) 312  (35) R32  (40) 432
C                        (12) 422   (20) 321  (36) R3M  (41) -43M
C                        (13) 4MM   (21) 31M  (37) R-3M (42) M3M
C                        (14) -42M  (22) 3M1
C                        (15) -4M2  (23) -31M
C                        (16) 4/MMM (24) -3M1
C 
C                                    HEXAG.
C
C                                   (25) 6
C                                   (26) -6
C                                   (27) 6/M
C
C                                   (28) 622
C                                   (29) 6MM
C                                   (30) -6M2
C                                   (31) -62M
C                                   (32) 6/MMM
C
C A TOTAL OF 42 CASES IS TABULATED BECAUSE THERE ARE TEN CASES OF
C ALTERNATIVE AXES FOR SEVEN OF THE 32 CRYSTALLOGRAPHIC POINT GROUPS:
C
C -42M  -4M2
C 3     R3
C -3    R-3
C 312   321   R32
C 31M   3M1   R3M
C -31M  -3M1  R-3M
C -62M  -6M2
C
      INTEGER PTGP,H
      E=1
      GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,
     & 23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42) PTGP
 1    M=1
      GO TO 99
 2    M=2
      GO TO 99
 3    M=2
      IF (H.EQ.0.AND.L.EQ.0) E=2
      GO TO 99
 4    M=2
      IF (K.EQ.0) E=2
      GO TO 99
 5    M=4
      IF (K.EQ.0) E=2
      IF (H.EQ.0.AND.L.EQ.0) E=2
      GO TO 99
 6    M=4
      IF (K.EQ.0.AND.L.EQ.0) E=2
      IF (H.EQ.0.AND.L.EQ.0) E=2
      IF (H.EQ.0.AND.K.EQ.0) E=2
      GO TO 99
 7    M=4
      IF (H.EQ.0) E=2
      IF (K.EQ.0) E=2
      IF (H.EQ.0.AND.K.EQ.0) E=4
      GO TO 99
 8    M=8
      IF (H.EQ.0) E=2
      IF (K.EQ.0) E=2
      IF (L.EQ.0) E=2
      IF (K.EQ.0.AND.L.EQ.0) E=4
      IF (H.EQ.0.AND.L.EQ.0) E=4
      IF (H.EQ.0.AND.K.EQ.0) E=4
      GO TO 99
 9    M=4
      IF (H.EQ.0.AND.K.EQ.0) E=4
      GO TO 99
 10   M=4
      IF (H.EQ.0.AND.K.EQ.0) E=2
      GO TO 99
 11   M=8
      IF (L.EQ.0) E=2
      IF (H.EQ.0.AND.K.EQ.0) E=4
      GO TO 99
 12   M=8
      IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-H)) E=2
      IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0)) E=2
      IF (H.EQ.0.AND.K.EQ.0) E=4
      GO TO 99
 13   M=8
      IF (K.EQ.0.OR.H.EQ.0) E=2
      IF (K.EQ.H.OR.K.EQ.-H) E=2
      IF (H.EQ.0.AND.K.EQ.0) E=8
      GO TO 99
 14   M=8
      IF (K.EQ.H.OR.K.EQ.-H) E=2
      IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0)) E=2
      IF (H.EQ.0.AND.K.EQ.0) E=4
      GO TO 99
 15   M=8
      IF (K.EQ.0.OR.H.EQ.0) E=2
      IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-H)) E=2
      IF (H.EQ.0.AND.K.EQ.0) E=4
      GO TO 99
 16   M=16
      IF (K.EQ.0.OR.H.EQ.0.OR.L.EQ.0) E=2
      IF (K.EQ.H.OR.K.EQ.-H) E=2
      IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-H)) E=4
      IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0)) E=4
      IF (H.EQ.0.AND.K.EQ.0) E=8
      GO TO 99
 17   M=3
      IF (H.EQ.0.AND.K.EQ.0) E=3
      GO TO 99
 18   M=6
      IF (H.EQ.0.AND.K.EQ.0) E=3
      GO TO 99
 19   M=6
      IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H)) E=2
      IF (H.EQ.0.AND.K.EQ.0) E=3
      GO TO 99
 20   M=6
      IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K)) E=2
      IF (H.EQ.0.AND.K.EQ.0) E=3
      GO TO 99
 21   M=6
      IF (K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K) E=2
      IF (H.EQ.0.AND.K.EQ.0) E=6
      GO TO 99
 22   M=6
      IF (K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H) E=2
      IF (H.EQ.0.AND.K.EQ.0) E=6
      GO TO 99
 23   M=12
      IF (K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K) E=2
      IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H)) E=2
      IF (H.EQ.0.AND.K.EQ.0) E=6
      GO TO 99
 24   M=12
      IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K)) E=2
      IF (K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H) E=2
      IF (H.EQ.0.AND.K.EQ.0) E=6
      GO TO 99
 25   M=6
      IF (H.EQ.0.AND.K.EQ.0) E=6
      GO TO 99
 26   M=6
      IF (L.EQ.0) E=2
      IF (H.EQ.0.AND.K.EQ.0) E=3
      GO TO 99
 27   M=12
      IF (L.EQ.0) E=2
      IF (H.EQ.0.AND.K.EQ.0) E=6
      GO TO 99
 28   M=12
      IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K)) E=2
      IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H)) E=2
      IF (H.EQ.0.AND.K.EQ.0) E=6
      GO TO 99
 29   M=12
      IF (K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K) E=2
      IF (K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H) E=2
      IF (H.EQ.0.AND.K.EQ.0) E=12
      GO TO 99
 30   M=12
      IF (L.EQ.0) E=2
      IF (K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H) E=2
      IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H)) E=4
      IF (H.EQ.0.AND.K.EQ.0) E=6
      GO TO 99
 31   M=12
      IF (L.EQ.0) E=2
      IF (K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K) E=2
      IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K)) E=4
      IF (H.EQ.0.AND.K.EQ.0) E=6
      GO TO 99
 32   M=24
      IF (L.EQ.0) E=2
      IF (K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K) E=2
      IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K)) E=4
      IF (K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H) E=2
      IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H)) E=4
      IF (H.EQ.0.AND.K.EQ.0) E=12
      GO TO 99
 33   M=3
      IF (K.EQ.H.AND.L.EQ.H) E=3
      GO TO 99
 34   M=6
      IF (K.EQ.H.AND.L.EQ.H) E=3
      GO TO 99
 35   M=6
      IF (L.EQ.0.AND.K.EQ.-H) E=2
      IF (K.EQ.0.AND.L.EQ.-H) E=2
      IF (H.EQ.0.AND.K.EQ.-L) E=2
      IF (K.EQ.H.AND.L.EQ.H) E=3
      GO TO 99
 36   M=6
      IF (K.EQ.H.OR.L.EQ.H.OR.L.EQ.K) E=2
      IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=2
      IF (K.EQ.H.AND.L.EQ.H) E=6
      GO TO 99
 37   M=12
      IF (K.EQ.H.OR.L.EQ.H.OR.K.EQ.L) E=2
      IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=2
      IF (K.EQ.-H.AND.L.EQ.0) E=2
      IF (L.EQ.-H.AND.K.EQ.0) E=2
      IF (L.EQ.-K.AND.H.EQ.0) E=2
      IF (K.EQ.H.AND.L.EQ.H) E=6
      GO TO 99
 38   M=12
      IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=3
      IF (K.EQ.0.AND.L.EQ.0) E=2
      IF (H.EQ.0.AND.L.EQ.0) E=2
      IF (H.EQ.0.AND.K.EQ.0) E=2
      GO TO 99
 39   M=24
      IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=3
      IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) E=2
      IF (K.EQ.0.AND.L.EQ.0) E=4
      IF (H.EQ.0.AND.L.EQ.0) E=4
      IF (H.EQ.0.AND.K.EQ.0) E=4
      GO TO 99
 40   M=24
      IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=3
      IF (L.EQ.0.AND.ABS(K).EQ.ABS(H)) E=2
      IF (K.EQ.0.AND.ABS(L).EQ.ABS(H)) E=2
      IF (H.EQ.0.AND.ABS(K).EQ.ABS(L)) E=2
      IF (K.EQ.0.AND.L.EQ.0) E=4
      IF (H.EQ.0.AND.L.EQ.0) E=4
      IF (H.EQ.0.AND.K.EQ.0) E=4
      GO TO 99
 41   M=24
      IF (ABS(K).EQ.ABS(H).OR.ABS(L).EQ.ABS(H).OR.ABS(K).EQ.ABS(L)) E=2
      IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=6
      IF (K.EQ.0.AND.L.EQ.0) E=4
      IF (H.EQ.0.AND.L.EQ.0) E=4
      IF (H.EQ.0.AND.K.EQ.0) E=4
      GO TO 99
 42   M=48
      IF (ABS(K).EQ.ABS(H).OR.ABS(L).EQ.ABS(H).OR.ABS(K).EQ.ABS(L)) E=2
      IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=6
      IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) E=2
      IF (L.EQ.0.AND.ABS(K).EQ.ABS(H)) E=4
      IF (K.EQ.0.AND.ABS(L).EQ.ABS(H)) E=4
      IF (H.EQ.0.AND.ABS(K).EQ.ABS(L)) E=4
      IF (K.EQ.0.AND.L.EQ.0) E=8
      IF (H.EQ.0.AND.L.EQ.0) E=8
      IF (H.EQ.0.AND.K.EQ.0) E=8
 99   M=M/E
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE FCALC (N,X,A1,B1,A2,B2,A3,B3,A4,B4,C0,FP,FPP,S,
     & SUMFSQ)
      DIMENSION X(1),A1(1),B1(1),A2(1),B2(1),A3(1),B3(1),A4(1),B4(1),
     & C0(1),FP(1),FPP(1)
      SUMFSQ=0
      DO 1 I=1,N
        SUMFSQ=SUMFSQ+X(I)*((A1(I)*EXP(-B1(I)*S**2)
     &                       +A2(I)*EXP(-B2(I)*S**2)
     &                        +A3(I)*EXP(-B3(I)*S**2)
     &                         +A4(I)*EXP(-B4(I)*S**2)+C0(I)+FP(I))**2
     &                          +FPP(I)**2)
 1    CONTINUE
      RETURN   
      END
C-----------------------------------------------------------------------
      SUBROUTINE FSQTOF (FSQ,SIGFSQ,F,SIGF)
      IF (FSQ.LE.0) THEN
        F=0
      ELSE
        F=SQRT(FSQ)
      END IF
      IF (FSQ.LE.SIGFSQ) THEN
        SIGF=0.5*SQRT(SIGFSQ)
      ELSE
        SIGF=0.5*SIGFSQ/F
      END IF
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE FTOFSQ (F,SIGF,FSQ,SIGFSQ)
      IF (F.LT.0) F=0
      IF (SIGF.LT.0) SIGF=0
      FSQ=F**2
      SIGFSQ=MAX(2*F*SIGF,4*SIGF**2)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE METRIC (A,B,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)
      RAD=PI/180
      DEG=180/PI
      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 PARITY (I,J,K,L)
C
C 1 EEE
C 2 EEO
C 3 EOE
C 4 OEE
C 5 EOO
C 6 OEO
C 7 OOE
C 8 OOO
C
      P=MOD(ABS(J),2)
      Q=MOD(ABS(K),2)
      R=MOD(ABS(L),2)
      IF (P.EQ.0.AND.Q.EQ.0.AND.R.EQ.0) I=1
      IF (P.EQ.0.AND.Q.EQ.0.AND.R.EQ.1) I=2
      IF (P.EQ.0.AND.Q.EQ.1.AND.R.EQ.0) I=3
      IF (P.EQ.1.AND.Q.EQ.0.AND.R.EQ.0) I=4
      IF (P.EQ.0.AND.Q.EQ.1.AND.R.EQ.1) I=5
      IF (P.EQ.1.AND.Q.EQ.0.AND.R.EQ.1) I=6
      IF (P.EQ.1.AND.Q.EQ.1.AND.R.EQ.0) I=7
      IF (P.EQ.1.AND.Q.EQ.1.AND.R.EQ.1) I=8
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE PLOTA (N,X,Y,ILP)
      PARAMETER (NX=100,NY=50)
      CHARACTER AA*1
      DIMENSION AA(0:NX,0:NY),X(N),Y(N)
C
C BLANK OUT THE PLOT ARRAY.
C
      DO  IX=0,NX
      DO  IY=0,NY
        AA(IX,IY)=' '
      END DO
      END DO
C
C FILL IN THE GRID MARKS.
C
      DO IX=0,NX
        AA(IX, 0)='.'
        AA(IX,NY)='.'
      END DO
      DO IX=0,NX,10
        AA(IX, 0)='+'
        AA(IX,NY)='+'
      END DO
      DO IY=0,NY
        AA( 0,IY)='.'
        AA(NX,IY)='.'
      END DO
      DO IY=0,NY,10
        AA( 0,IY)='+'
        AA(NX,IY)='+'
      END DO
C
C SET THE PLOT LIMITS, 0 < E < 5.0 AND 0 < P(E) < 1.
C
      XMIN=0
      XMAX=5
      YMIN=0
      YMAX=1
      XDEL=XMAX-XMIN
      YDEL=YMAX-YMIN
C
C FILL IN THE DISTRIBUTION FUNCTION.
C
      DX=XDEL/(2*NX)
      DO I=0,(2*NX)
        XI=I*DX
        YI=2*XI*EXP(-XI**2)
        IX=NINT(NX*(XI-XMIN)/XDEL)
        IY=NINT(NY*(YMAX-YI)/YDEL)
        IX=MAX(IX,0)
        IX=MIN(IX,NX)
        IY=MAX(IY,0)
        IY=MIN(IY,NY)
        AA(IX,IY)='*'
      END DO
C
C FILL IN THE EXPERIMENTAL DATA POINTS.
C
      DO I=1,N
        IX=NINT(NX*(X(I)-XMIN)/XDEL)
        IY=NINT(NY*(YMAX-Y(I))/YDEL)
        IX=MAX(IX,0)
        IX=MIN(IX,NX)
        IY=MAX(IY,0)
        IY=MIN(IY,NY)
        AA(IX,IY)='O'
      END DO
C
C PRINT THE PLOT ARRAY.
C
      DY=(YMAX-YMIN)/NY
      DO IY=0,NY
        IF (MOD(IY,10).EQ.0) THEN
          WRITE (ILP,100) YMAX-IY*DY,(AA(IX,IY),IX=0,NX)
        ELSE
          WRITE (ILP,101)            (AA(IX,IY),IX=0,NX)
        END IF
      END DO
      DX=(XMAX-XMIN)/10
      WRITE (ILP,102) (XMIN+IX*DX,IX=0,10)
  100 FORMAT (' ',F10.1,101A1)
  101 FORMAT (' ',10X,  101A1)
  102 FORMAT (' ',1X,11F10.1)
      WRITE (ILP,'(/1X,10X,''ACENTRIC DISTRIBUTION:  P(E) = 2*E*EXP(-E**
     &2),  "*" THEORETICAL,  "O" EXPERIMENTAL'')')
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE PLOTC (N,X,Y,ILP)
      PARAMETER (NX=100,NY=50)
      CHARACTER AA*1
      DIMENSION AA(0:NX,0:NY),X(N),Y(N)
C
C BLANK OUT THE PLOT ARRAY.
C
      DO  IX=0,NX
      DO  IY=0,NY
        AA(IX,IY)=' '
      END DO
      END DO
C
C FILL IN THE GRID MARKS.
C
      DO IX=0,NX
        AA(IX, 0)='.'
        AA(IX,NY)='.'
      END DO
      DO IX=0,NX,10
        AA(IX, 0)='+'
        AA(IX,NY)='+'
      END DO
      DO IY=0,NY
        AA( 0,IY)='.'
        AA(NX,IY)='.'
      END DO
      DO IY=0,NY,10
        AA( 0,IY)='+'
        AA(NX,IY)='+'
      END DO
C
C SET THE PLOT LIMITS, 0 < E < 5.0 AND 0 < P(E) < 1.
C
      XMIN=0
      XMAX=5
      YMIN=0
      YMAX=1
      XDEL=XMAX-XMIN
      YDEL=YMAX-YMIN
C
C FILL IN THE DISTRIBUTION FUNCTION.
C
      DX=XDEL/(2*NX)
      C=SQRT(2/3.141593)
      DO I=0,(2*NX)
        XI=I*DX
        YI=C*EXP(-0.5*XI**2)
        IX=NINT(NX*(XI-XMIN)/XDEL)
        IY=NINT(NY*(YMAX-YI)/YDEL)
        IX=MAX(IX,0)
        IX=MIN(IX,NX)
        IY=MAX(IY,0)
        IY=MIN(IY,NY)
        AA(IX,IY)='*'
      END DO
C
C FILL IN THE EXPERIMENTAL DATA POINTS.
C
      DO I=1,N
        IX=NINT(NX*(X(I)-XMIN)/XDEL)
        IY=NINT(NY*(YMAX-Y(I))/YDEL)
        IX=MAX(IX,0)
        IX=MIN(IX,NX)
        IY=MAX(IY,0)
        IY=MIN(IY,NY)
        AA(IX,IY)='O'
      END DO
C
C PRINT THE PLOT ARRAY.
C
      DY=(YMAX-YMIN)/NY
      DO IY=0,NY
        IF (MOD(IY,10).EQ.0) THEN
          WRITE (ILP,100) YMAX-IY*DY,(AA(IX,IY),IX=0,NX)
        ELSE
          WRITE (ILP,101)            (AA(IX,IY),IX=0,NX)
        END IF
      END DO
      DX=(XMAX-XMIN)/10
      WRITE (ILP,102) (XMIN+IX*DX,IX=0,10)
  100 FORMAT (' ',F10.1,101A1)
  101 FORMAT (' ',10X,  101A1)
  102 FORMAT (' ',1X,11F10.1)
      WRITE (ILP,'(/1X,10X,''CENTRIC DISTRIBUTION:  P(E) = SQRT(2/PI)*EX
     &P(-0.5*E**2),  "*" THEORETICAL,  "O" EXPERIMENTAL'')')
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE PLOTS (N,X,Y,SMAX,ILP)
      PARAMETER (NX=100,NY=50)
      CHARACTER AA*1
      DIMENSION AA(0:NX,0:NY),X(N),Y(N),YPP(50)
C
C BLANK OUT THE PLOT ARRAY.
C
      DO  IX=0,NX
      DO  IY=0,NY
        AA(IX,IY)=' '
      END DO
      END DO
C
C FILL IN THE GRID MARKS.
C
      DO IX=0,NX
        AA(IX, 0)='.'
        AA(IX,NY)='.'
      END DO
      DO IX=0,NX,10
        AA(IX, 0)='+'
        AA(IX,NY)='+'
      END DO
      DO IY=0,NY
        AA( 0,IY)='.'
        AA(NX,IY)='.'
      END DO
      DO IY=0,NY,10
        AA( 0,IY)='+'
        AA(NX,IY)='+'
      END DO
C
C SET THE PLOT LIMITS, 0 < SIN(THETA)/LAMBDA < SMAX AND 0 < E(S) < 2.5.
C
      XMIN=0
      XMAX=SMAX
      YMIN=0
      YMAX=2.5
      XDEL=XMAX-XMIN
      YDEL=YMAX-YMIN
C
C MARK <E**2> = 1.
C
      DO IX=0,NX
        AA(IX,30)='.'
      END DO
C
C FILL IN THE SPLINE-INTERPOLATED DATA POINTS.
C
      CALL SPLINE (N,X,Y,YPP)
      DX=(X(N)-X(1))/(2*NX)
      DO I=0,(2*NX)
        XI=X(1)+I*DX
        CALL SPLINT (N,X,Y,YPP,XI,YI)
        IX=NINT(NX*(XI-XMIN)/XDEL)
        IY=NINT(NY*(YMAX-YI)/YDEL)
        IX=MAX(IX,0)
        IX=MIN(IX,NX)
        IY=MAX(IY,0)
        IY=MIN(IY,NY)
        AA(IX,IY)='*'
      END DO
C
C FILL IN THE EXPERIMENTAL DATA POINTS.
C
      DO I=1,N
        IX=NINT(NX*(X(I)-XMIN)/XDEL)
        IY=NINT(NY*(YMAX-Y(I))/YDEL)
        IX=MAX(IX,0)
        IX=MIN(IX,NX)
        IY=MAX(IY,0)
        IY=MIN(IY,NY)
        AA(IX,IY)='O'
      END DO
C
C PRINT THE PLOT ARRAY.
C
      DY=(YMAX-YMIN)/NY
      DO IY=0,NY
        IF (MOD(IY,10).EQ.0) THEN
          WRITE (ILP,100) YMAX-IY*DY,(AA(IX,IY),IX=0,NX)
        ELSE
          WRITE (ILP,101)            (AA(IX,IY),IX=0,NX)
        END IF
      END DO
      DX=(XMAX-XMIN)/10
      WRITE (ILP,102) (2*(XMIN+IX*DX),IX=0,10)
  100 FORMAT (' ',F10.1,101A1)
  101 FORMAT (' ',10X,  101A1)
  102 FORMAT (' ',1X,11F10.4)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE READ1 (ITYPE,IUNIT,IEND,IH,IK,IL,FSQ,SIGFSQ,F,SIGF)
      FSQ=0
      SIGFSQ=0
      F=0
      SIGF=0
 1    IF (ITYPE.EQ.0) READ (IUNIT,*,END=9) IH,IK,IL,FSQ,SIGFSQ,F,SIGF
      IF (ITYPE.EQ.1) READ (IUNIT,*,END=9) IH,IK,IL,FSQ,SIGFSQ
      IF (ITYPE.EQ.2) READ (IUNIT,*,END=9) IH,IK,IL,F,  SIGF
      IF (ITYPE.EQ.3) READ (IUNIT,  END=9) IH,IK,IL,FSQ,SIGFSQ,F,SIGF
      IF (ITYPE.EQ.4) READ (IUNIT,  END=9) IH,IK,IL,FSQ,SIGFSQ
      IF (ITYPE.EQ.5) READ (IUNIT,  END=9) IH,IK,IL,F,  SIGF
      IF (SIGFSQ.LE.0.AND.SIGF.LE.0) GO TO 1
      IF (FSQ.LE.0.AND.F.LE.0) GO TO 1
      IF ((ITYPE.EQ.1.OR.ITYPE.EQ.4).AND.SIGFSQ.GT.0)
     & CALL FSQTOF (FSQ,SIGFSQ,F,SIGF)
      IF ((ITYPE.EQ.2.OR.ITYPE.EQ.5).AND.SIGF.GT.0)
     & CALL FTOFSQ (F,SIGF,FSQ,SIGFSQ)
      IF (F*SIGF.LE.0) GO TO 1
      RETURN
 9    IEND=1
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE READ2 (IFILE,IEND,IH,IK,IL,FSQ,SIGFSQ,F,SIGF,E,SIGE)
      READ (IFILE,*,END=9) IH,IK,IL,FSQ,SIGFSQ,F,SIGF,E,SIGE
      RETURN
 9    IEND=1
      RETURN
      END
C-----------------------------------------------------------------------
      FUNCTION SINTHL(IH,IK,IL,GINV)
      DIMENSION H(3),GINV(3,3)
      H(1)=IH
      H(2)=IK
      H(3)=IL
      SINTHL=0.5*SQRT(VTMV(H,GINV))
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE SORT (N,DATA,INDEX)
C
C INDEXES THE ARRAY DATA(N) AND RETURNS THE ARRAY INDEX(N) SORTED SUCH
C THAT THE VALUES DATA(INDEX(I)) ARE IN ASCENDING ORDER FOR I = 1, 2,
C ..., N.  THE INPUT VARIABLES N AND DATA(N) ARE NOT CHANGED.
C
C EMPLOYS THE "HEAPSORT" ALGORITHM.
C
C FORTRAN CODE QUOTED FROM WILLIAM H. PRESS, BRIAN P. FLANNERY, SAUL A.
C TEUKOLSKY, AND WILLIAM T. VETTERLING (1986).  NUMERICAL RECIPIES:  THE
C                                               --------- --------   ---
C ART OF SCIENTIFIC COMPUTING, PP. 229-233.  CAMBRIDGE, ENGLAND:
C --- -- ---------- ---------
C CAMBRIDGE UNIVERSITY PRESS.
C
      DIMENSION DATA(N),INDEX(N)
      DO 10 I=1,N
        INDEX(I)=I
 10   CONTINUE
      I=N/2+1
      J=N
 1    CONTINUE
      IF (I.GT.1) THEN
        I=I-1
        INDEXT=INDEX(I)
        T=DATA(INDEXT)
      ELSE
        INDEXT=INDEX(J)
        T=DATA(INDEXT)
        INDEX(J)=INDEX(1)
        J=J-1
        IF (J.EQ.1) THEN
          INDEX(1)=INDEXT
          RETURN
        END IF
      END IF
      K=I
      L=I+I
 2    IF (L.LE.J) THEN
        IF (L.LT.J) THEN
          IF (DATA(INDEX(L)).LT.DATA(INDEX(L+1))) L=L+1
        END IF
        IF (T.LT.DATA(INDEX(L))) THEN
          INDEX(K)=INDEX(L)
          K=L
          L=L+L
        ELSE
          L=J+1
        END IF
        GO TO 2
      END IF
      INDEX(K)=INDEXT
      GO TO 1
      END
C-----------------------------------------------------------------------
      SUBROUTINE SPLINE (N,X,Y,YPP)
C
C NATURAL CUBIC SPLINE FIT TO A TABULATED FUNCTION Y(I) = Y(X(I)), I =
C 1, 2,..., N, WITH X(1) .LT. X(2) .LT. ... .LT. X(N).  RETURNS THE
C TABULATED SECOND DERIVATIVES y'' = d2y/dx2 = YPP(I), I = 1, 2,..., N,
C WITH YPP(1) = YPP(N) = 0.
C
C FORTRAN CODE QUOTED FROM WILLIAM H. PRESS, BRIAN P. FLANNERY, SAUL A.
C TEUKOLSKY, AND WILLIAM T. VETTERLING (1986).  NUMERICAL RECIPIES:  THE
C ART OF SCIENTIFIC COMPUTING, PP. 86-89.  CAMBRIDGE, ENGLAND:
C CAMBRIDGE UNIVERSITY PRESS.
C
C MUST HAVE NMAX .GE. N.
C
      PARAMETER (NMAX=1000)
      DIMENSION X(N),Y(N),YPP(N),U(NMAX)
      YPP(1)=0
      YPP(N)=0
      U(1)=0
      U(N)=0
      DO I=2,N-1
        P=(X(I)-X(I-1))/(X(I+1)-X(I-1))
        Q=P*YPP(I-1)+2
        YPP(I)=(P-1)/Q
        U(I)=(6*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1))/
     &   (X(I)-X(I-1)))/(X(I+1)-X(I-1))-P*U(I-1))/Q
      END DO
      DO I=N-1,1,-1
        YPP(I)=YPP(I)*YPP(I+1)+U(I)
      END DO
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE SPLINT (N,X,Y,YPP,XI,YI)
C
C CUBIC SPLINE INTERPOLATION OF A TABULATED FUNCTION Y(I) = Y(X(I)), I =
C 1, 2,..., N, USING THE TABULATED SECOND DERIVATIVES YPP(I) FROM
C SUBROUTINE SPLINE.
C
C FORTRAN CODE QUOTED FROM WILLIAM H. PRESS, BRIAN P. FLANNERY, SAUL A.
C TEUKOLSKY, AND WILLIAM T. VETTERLING (1986).  NUMERICAL RECIPIES:  THE
C ART OF SCIENTIFIC COMPUTING, PP. 86-89.  CAMBRIDGE, ENGLAND:
C CAMBRIDGE UNIVERSITY PRESS.
C
      DIMENSION X(N),Y(N),YPP(N)
      IF (XI.LE.X(1)) THEN
        YI=Y(1)
        RETURN
      END IF
      IF (XI.GE.X(N)) THEN
        YI=Y(N)
        RETURN
      END IF
      I1=1
      I2=N
 1    IF (I2-I1.GT.1) THEN
        I=(I2+I1)/2
        IF (X(I).GT.XI) THEN
          I2=I
        ELSE
          I1=I
        END IF
        GO TO 1
      END IF
      H=X(I2)-X(I1)
      A=(X(I2)-XI)/H
      B=(XI-X(I1))/H
      YI=A*Y(I1)+B*Y(I2)+((A**3-A)*YPP(I1)+(B**3-B)*YPP(I2))*H**2/6
      RETURN
      END
C-----------------------------------------------------------------------
      FUNCTION VARF(SCALEK,B,C,COV,S,H,HTBH,HTCH,Q)
      DIMENSION B(3,3),C(3,3),COV(13,13),H(3),DQDP(13)
      DQDP(1)=-Q/SCALEK
      IF (B(2,2).EQ.0) THEN
        DQDP(2)=-Q*S**2
        DQDP(3)=2*Q*C(1,1)*S**4
        N=3
      ELSE
        N=1
        DO 1 I=1,3
        DO 1 J=I,3
          N=N+1
          DQDP(N)=-Q*(2-DELTA(I,J))*H(I)*H(J)
 1      CONTINUE
        DO 2 I=1,3
        DO 2 J=I,3
          N=N+1
          DQDP(N)=2*Q*HTCH*(2-DELTA(I,J))*H(I)*H(J)
 2      CONTINUE
      END IF
      VARF=0
      DO 3 I=1,N
      DO 3 J=I,N
        VARF=VARF+(2-DELTA(I,J))*DQDP(I)*DQDP(J)*COV(I,J)
 3    CONTINUE
      RETURN
      END
C-----------------------------------------------------------------------
      FUNCTION VTMV(A,B)
C
C QUADRATIC FORM, VTMV, EQUALS ROW VECTOR, A-TRANSPOSE, TIMES SQUARE
C MATRIX, B, TIMES COLUMN VECTOR, A.
C
      DIMENSION A(3),B(3,3)
      VTMV=0
      DO 1 I=1,3
      DO 1 J=1,3
        VTMV=VTMV+A(J)*B(I,J)*A(I)
 1    CONTINUE
      RETURN
      END
C-----------------------------------------------------------------------
C=====================================================================
C Machine-dependent system-service routines (such as TIME and DATE)
C=====================================================================
      SUBROUTINE VFDATE(DAYTIME)
C
C This routine returns the date/time as a 20 charater string.
C Truncation occurs if the length of the input string is less than 20.
C The format of the returned string is "dd mmm yyyy hh:mm:ss".
C Note: this routine is Y2K compliant.
C
C +++++++++++++++++++++++++++++++++++++++
C Machine dependent intel/Win32 version.
C +++++++++++++++++++++++++++++++++++++++
C
      IMPLICIT NONE
C input/output
      CHARACTER*(*) DAYTIME
C local
      CHARACTER DATE*8, TIME*10, ZONE*5, MONTH(12)*3, TEMP*20
      INTEGER VALUES(8)
      
      DATA MONTH/'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun',
     &           'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'/
C
      CALL DATE_AND_TIME(DATE, TIME, ZONE, VALUES)
      TEMP = DATE(7:8) // ' ' // MONTH(VALUES(2)) // ' ' //
     &       DATE(1:4) // ' ' // TIME(1:2) // ':' // TIME(3:4) //
     &       ':' // TIME(5:6)
      DAYTIME = TEMP
      RETURN
      END
C=====================================================================
      SUBROUTINE VDATE(ADATE)
C
C This routine returns the date as a 9 charater string.  Truncation
C occurs if the length of the input string is less than 9.  The format
C of the returned string is "dd-mmm-yy".
C
C Note: this routine is not Y2K compliant.
C
C +++++++++++++++++++++++++++++++++++++++
C Machine dependent intel/Win32 version.
C +++++++++++++++++++++++++++++++++++++++
C
      IMPLICIT NONE
C input/output
      CHARACTER*(*) ADATE
C local
      CHARACTER TEMP*20
C
      CALL VFDATE(TEMP)
      ADATE = TEMP(1:2) // '-' // TEMP(4:6) // '-' // TEMP(10:11)     
      RETURN
      END
C=====================================================================
      SUBROUTINE VTIME(ATIME)
C
C This routine returns the time of day as an 8 character string.
C Truncation occurs if the length of the input string is less than 8.
C The format of the returned string is "hh:mm:ss".
C
C +++++++++++++++++++++++++++++++++++++++
C Machine dependent intel/Win32 version.
C +++++++++++++++++++++++++++++++++++++++
C
      IMPLICIT NONE
C input/output
      CHARACTER*(*) ATIME
C local
      CHARACTER TEMP*20
C
      CALL VFDATE(TEMP)
      ATIME = TEMP(13:20)
      RETURN
      END
C=====================================================================
C-----------------------------------------------------------------------
      SUBROUTINE WRITE1 (IFILE,IH,IK,IL,FSQ,SIGFSQ,F,SIGF,E,SIGE)
      WRITE (IFILE,'(1X,3I5,4(E12.4,E10.2))')
     & IH,IK,IL,FSQ,SIGFSQ,F,SIGF,E,SIGE
      RETURN
      END
C-----------------------------------------------------------------------
