      PROGRAM BAYES
C
C ROBERT H. BLESSING
C HAUPTMAN-WOODWARD INSTITUTE
C 73 HIGH STREET
C BUFFALO, NEW YORK 14203, USA
C TELEPHONE:  (716) 856-9600, EXTENSION 335
C E-MAIL:  blessing@hwi.buffalo.edu
C
C APRIL 1986,..., JANUARY 2001
C
C GIVEN A SET OF MEASURED INTENSITY DATA, H, K, L, FSQ, AND SIGMA(FSQ),
C THE PROGRAM USES A BAYESIAN PROBABILISTIC TREATMENT TO ESTIMATE
C STATISTICAL EXPECTATION VALUES FOR FOBS**2, SIGMA(FOBS**2), FOBS, AND
C SIGMA(FOBS).
C
C THE SUBROUTINES NEGIA AND NEGIC, WHICH ESTIMATE THE BAYESIAN POSTERIOR
C MEANS AND VARIANCES, ARE FROM SUPPLEMENTARY PUBLICATION NO. SUP 33352
C FROM:
C
C SIMON FRENCH AND KEITH WILSON (1979).  ON THE TREATMENT OF NEGATIVE
C INTENSITY OBSERVATIONS.  ACTA CRYST. A34, 517-525.
C
      CHARACTER ATIME*8,ADATE*9,TITLE*80,FILEI*80,FILEO*80,FORM*11
      CHARACTER SYST*12,LATT*1,PTGP*5,LAUE*5
      DIMENSION CELL(6),GINV(3,3)
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
      IO3=30
      ILP=60
C
C READ CONTROL DATA.                                                   
C
      OPEN (UNIT=IO1,FILE='bayes.dat',STATUS='OLD')
      READ (IO1,'(A)') TITLE
      READ (IO1,*)     ITYPE
      READ (IO1,'(A)') FILEI
      READ (IO1,'(A)') FILEO
      READ (IO1,'(A)') LATT
      READ (IO1,'(A)') PTGP
      READ (IO1,*)     CELL
      NLOCAL=0
      NOBAYES=0
      READ (IO1,*,END=5) NLOCAL
      READ (IO1,*,END=5) NOBAYES
 5    CLOSE (UNIT=IO1,STATUS='KEEP')
      CALL METRIC (CELL,GINV)
      CALL DECODE (PTGP,LAUE,SYST,LATT,MLATT,IPTGP,ILAUE)
C
C ECHO CONTROL DATA.
C
      OPEN (UNIT=ILP,FILE='bayes.lp',STATUS='NEW')
      WRITE (ILP,6000) ATIME,ADATE,TITLE
 6000 FORMAT (/1H1,130('-')/1X,'PROGRAM BAYES.  ',A,', ',A,'.  ',A)
      WRITE (ILP,'(/1X,
     &''INPUT REFLECTION FILE TYPE:   ITYPE = '',I2)') ITYPE
      WRITE (ILP,'(/1X,
     &''INPUT REFLECTION FILE NAME:   '',A//1X,
     &''OUTPUT REFLECTION FILE NAME:  '',A)')
     & FILEI,FILEO
      WRITE (ILP,'(/1X,A//1X,A,''-LATTICE''//1X,''LAUE GROUP   '',A//1X,
     &''POINT GROUP  '',A)') SYST,LATT,LAUE,PTGP
      IF (PTGP.EQ.LAUE) THEN
        WRITE (ILP,'(/1X,''CENTROSYMMETRIC STRUCTURE'')')
      ELSE
        WRITE (ILP,'(/1X,''NONCENTROSYMMETRIC STRUCTURE'')')
      END IF
      WRITE (ILP,'(/1X,
     &''LATTICE PARAMETERS''/1X,
     &''         A         B         C     ALPHA      BETA     GAMMA''/
     &1X,3F10.4,3F10.3)') CELL
      IF (NOBAYES.NE.0) WRITE (ILP,'(//1X,
     &''BAYESIAN REPLACEMENT WILL NOT BE DONE.  INSTEAD:''/1X,
     &''  FSQ = MAX(FSQ, 0.0)''/1X,
     &''  SIGMA_FSQ = MAX(SIGMA_FSQ, 0.0)''/1X,
     &''  F = SQRT(FSQ)''/1X,
     &''  IF (FSQ .GE. SIGMA_FSQ .AND. SIGMA_FSQ .GT. 0) THEN''/1X,
     &''    SIGMA_F = SIGMA_FSQ/(2*F)''/1X,
     &''  ELSE''/1X,
     &''    SIGMA_F = 0.5*SQRT(SIGMA_FSQ)''/1X,
     &''  END IF''/1X,
     &''  IF (SIGMA_FSQ .GT. 0 .AND. SIGMA_F .GT. 0) THEN''/1X,
     &''    IREJECT = 0''/1X,
     &''  ELSE''/1X,
     &''    IREJECT = 1''/1X,
     &''  END IF''/)')
C
C COPY DATA FROM THE INPUT REFLECTION DATA FILE TO A DIRECT ACCESS
C WORKING DATA FILE THAT WILL ACCOMMODATE 15-WORD, 60-BYTE RECORDS:
C
C h, k, l, Fsq, sigma(Fsq), F, sigma(F), E, sigma(E),
C icent, mult, epsilon, sin(theta)/lambda, <Fsq/epsilon>, sigma(<Fsq/e>)
C
      IF (ITYPE.LT.2) THEN
        FORM='FORMATTED'
      ELSE
        FORM='UNFORMATTED'
      END IF
      OPEN (UNIT=IO1,FILE=FILEI,STATUS='OLD',FORM=FORM)
      OPEN (UNIT=IO2,STATUS='SCRATCH',FORM='UNFORMATTED',
     & ACCESS='DIRECT',RECL=60)
      IF (IPTGP.LT.ILAUE) THEN
        JMIN=+999
        JMAX=-999
        KMIN=+999
        KMAX=-999
        LMIN=+999
        LMAX=-999
      END IF
      NDATA=0
      IEND=0
 1    CALL READ1 (ITYPE,IO1,IEND,J,K,L,FSQ,SIGFSQ)
      IF (IEND.NE.0) GO TO 9
      NDATA=NDATA+1
      WRITE (IO2,REC=NDATA) J,K,L,FSQ,SIGFSQ,
     & 0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0.0, 0.0, 0.0
      IF (IPTGP.LT.ILAUE) THEN
C
C STORE THE LAUE-UNIQUE INDEX LIMITS.
C
        CALL EQUIV (ILAUE,J,K,L)
        JMIN=MIN(JMIN,J)
        JMAX=MAX(JMAX,J)
        KMIN=MIN(KMIN,K)
        KMAX=MAX(KMAX,K)
        LMIN=MIN(LMIN,L)
        LMAX=MAX(LMAX,L)
      END IF
      GO TO 1
 9    CLOSE (UNIT=IO1,STATUS='KEEP')
      IF (IPTGP.LT.ILAUE) THEN
C
C STORE A TEMPORARY ARRAY OF THE BIJVOET- OR FRIEDEL-PAIR POPULATIONS
C (0, 1, OR 2) INDEXED BY THE PACKED LAUE-UNIQUE INDICES.
C
        NJ=JMAX-JMIN+1
        NK=KMAX-KMIN+1
        NL=LMAX-LMIN+1
        NMAX=NJ*NK*NL+NK*NL+NL
C
C ALLOCATE STORAGE FOR SORTING ARRAYS.
C
        ALLOCATE (DATA(NMAX),INDEX(NMAX))
        DO I=1,NMAX
          DATA(I)=0
        END DO
        DO I=1,NDATA
          READ (IO2,REC=I) J,K,L
          CALL EQUIV (ILAUE,J,K,L)
          JKL=(J-JMIN)*NK*NL+(K-KMIN)*NL+(L-LMIN)
          DATA(JKL)=DATA(JKL)+1
        END DO
      END IF
C
C EVALUATE ADDITIONAL REFLECTION INFORMATION AND STORE IT IN THE WORKING
C DATA FILE.
C
      SMIN=9
      SMAX=0
      DO I=1,NDATA
        READ (IO2,REC=I) J,K,L
C
C GET ACENTRIC OR CENTRIC REFLECTION FLAG.
C
        CALL CSYM (IPTGP,J,K,L,IC)
C
C GET POINT GROUP DEGENERACY AND MULTIPLICITY.
C
        CALL ESYM (IPTGP,J,K,L,IE,M)
C
C ADJUST DEGENERACY FOR LATTICE MULTIPLICITY.
C
        IE=MLATT*IE
C
C DOUBLE THE POINT GROUP MULTIPLICITY FOR ACENTRIC REFLECTIONS THAT ARE
C REPRESENTED BY ONLY ONE MEMBER OF THE BIJVOET OR FRIEDEL PAIR.
C
        IF (IC.EQ.0) THEN
          CALL EQUIV (ILAUE,J,K,L)
          JKL=(J-JMIN)*NK*NL+(K-KMIN)*NL+(L-LMIN)
          IF (DATA(JKL).EQ.1) M=2*M
        END IF
C
C GET SIN(THETA)/LAMBDA.
C
        S=SINTHL(J,K,L,GINV)
        SMIN=MIN(SMIN,S)
        SMAX=MAX(SMAX,S)
C
C RECORD THE ADDITIONAL REFLECTION INFORMATION IN THE WORKING DATA FILE.
C
        READ  (IO2,REC=I) J,K,L,FSQ,SIGFSQ
        WRITE (IO2,REC=I) J,K,L,FSQ,SIGFSQ,
     &  0.0, 0.0, 0.0, 0.0, IC,M,IE,S, 0.0, 0.0
      END DO
      WRITE (ILP,'(/1X,
     &''N    = '',I6,'' DATA FOR EVALUATING INTENSITY AVERAGES IN SHE'',
     &''LLS OF S = SIN(THETA)/LAMBDA''/1X,
     &''SMIN = '',F6.3,'' RECIPROCAL ANGSTROMS    DMAX = '',F7.2,
     &'' ANGSTROMS''/1X,
     &''SMAX = '',F6.3,''                         DMIN = '',F7.2)')
     & NDATA,SMIN,1/(2*SMIN),SMAX,1/(2*SMAX)
C
C SORT ON INCREASING SIN(THETA)/LAMBDA.
C
      DO I=1,NDATA
        READ (IO2,REC=I) (IDUMMY,J=1,3),(DUMMY,J=1,6),(IDUMMY,J=1,3),S
        DATA(I)=S
      END DO
      CALL SORT (NDATA,DATA,INDEX)
C
C EVALUATE LOCAL AVERAGE INTENSITIES IN SHELLS OF SIN(THETA)/LAMBDA
C CONTAINING EQUAL NUMBERS OF REFLECTIONS.
C
C AT LEAST 10, AT MOST 500, REFLECTIONS PER SHELL.
C
      IF (NLOCAL.EQ.0) THEN
        NLOCAL=MAX(NDATA/100,10)
        NLOCAL=MIN(NLOCAL,500,NDATA)
      END IF
      NHALF=NLOCAL/2
      SUMM=0
      SUM1=0
      SUM2=0
      YMAX=0
      ZMAX=0
      DO I=1,NHALF
        READ (IO2,REC=INDEX(I)) (IDUMMY,J=1,3),FSQ,(DUMMY,J=1,5),
     &  IDUMMY,M,IE
        IF (FSQ.LT.0) FSQ=0
        FSQ=FSQ/IE
        SUMM=SUMM+M
        SUM1=SUM1+M*FSQ
        SUM2=SUM2+M*FSQ**2
      END DO
      DO I=1,NDATA
        I1=I-NHALF
        I2=I+NHALF
        IF (I1.LT.1) THEN
          I1=1
          M1=0
          Y1=0
        ELSE
          READ (IO2,REC=INDEX(I1)) (IDUMMY,J=1,3),FSQ,(DUMMY,J=1,5),
     &    IDUMMY,M,IE
          IF (FSQ.LT.0) FSQ=0
          M1=M
          Y1=FSQ/IE
        END IF
        IF (I2.GT.NDATA) THEN
          I2=NDATA
          M2=0
          Y2=0
        ELSE
          READ (IO2,REC=INDEX(I2)) (IDUMMY,J=1,3),FSQ,(DUMMY,J=1,5),
     &    IDUMMY,M,IE
          IF (FSQ.LT.0) FSQ=0
          M2=M
          Y2=FSQ/IE
        END IF
        SUMM=SUMM-M1+M2
        SUM1=SUM1-M1*Y1+M2*Y2
        SUM2=SUM2-M1*Y1**2+M2*Y2**2
        N=I2-I1+1
        AVG=SUM1/SUMM
        AVGSQ=SUM2/SUMM
        RMSD=SQRT((SUMM/(SUMM-1))*(AVGSQ-AVG**2))
        AVGFSQ=AVG
        SIGAVG=RMSD/SQRT(SUMM)
        READ  (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,
     &  DUMMY,DUMMY,DUMMY,DUMMY,IC,M,IE,S
        WRITE (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,
     &  0.0,  0.0,  0.0,  0.0,  IC,M,IE,S,AVGFSQ,SIGAVG
        YMAX=MAX(YMAX,AVGFSQ)
        ZMAX=MAX(ZMAX,SIGAVG)
      END DO
C
C PRINT THE FIRST 100 REFLECTIONS SORTED ON SIN(THETA)/LAMBDA.
C
      WRITE (ILP,6000) ATIME,ADATE,TITLE
      WRITE (ILP,'(/1X,
     &''   H   K   L    FMEAS**2 SIGMA(FMEAS**2) SIN(TH)/L    D(HKL) '',
     &''   <F**2/e> SIGMA(<F**2/e>)''/1X,
     &''   -   -   -    -------- --------------- ---------    ------ '',
     &''   -------- ---------------'')')
      DO I=1,MIN(100,NDATA)
        READ (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,(DUMMY,N=1,4),
     &  (IDUMMY,N=1,3),S,AVGFSQ,SIGAVG
        D=1/(2*S)
        WRITE (ILP,'(1X,3I4,2F12.2,4X,F10.4,F10.2,2F12.2)')
     &  J,K,L,FSQ,SIGFSQ,S,D,AVGFSQ,SIGAVG
      END DO
C
C PLOT <FSQ/EPSILON> VERSUS DSTAR AND SIGMA(<FSQ/EPSILON>) VERSUS DSTAR.
C
      WRITE (ILP,6001) ATIME,ADATE,TITLE
 6001 FORMAT (/1H1,130('-')/1X,10X,'PROGRAM BAYES.  ',A,', ',A,'.  ',A/)
      CALL PLOT (ILP,IO2,NDATA,0.0,SMAX,0.0,1.1*YMAX,0)
      WRITE (ILP,'(/1X,
     &''INTENSITIES LOCALLY AVERAGED IN SHELLS OF S = SIN(THETA)/LAMB'',
     &''DA.  Y = <FMEAS**2/EPSILON> VERSUS X = DSTAR = 1/D = 2*S,''/1X,
     &''WHERE EPSILON = EPSILON(HKL) IS THE DEGENERACY OF THE RECIPRO'',
     &''CAL LATTICE POINT HKL.''/1X,
     &''NTOTAL = '',I6,'' TOTAL DATA                       SMAX  = '',
     &F10.3,'' RECIPROCAL ANGSTROMS''/1X,
     &''NLOCAL = '',I6,'' LOCAL DATA PER SHELL             DMIN  = '',
     &F9.2,''  ANGSTROMS'')')
     & NDATA,SMAX,NLOCAL,1/(2*SMAX)
      WRITE (ILP,6001) ATIME,ADATE,TITLE
      CALL PLOT (ILP,IO2,NDATA,0.0,SMAX,0.0,1.1*ZMAX,1)
      WRITE (ILP,'(/1X,
     &''STANDARD UNCERTAINTIES OF THE LOCAL INTENSITY AVERAGES'',/1X,
     &''Y = SIGMA(<FMEAS**2/EPSILON>) VERSUS X = DSTAR = 2*SIN(THETA)'',
     &''/LAMBDA'')')
C
C TEST VALIDITY OF THE WILSON DISTRIBUTIONS.
C
      CALL NZPLOT (NDATA,IO2,ILP,ATIME,ADATE,TITLE,PTGP,LAUE)
C
C OBTAIN BAYESIAN ESTIMATES OF POSTERIOR MEANS AND VARIANCES.
C
      WRITE (ILP,6000) ATIME,ADATE,TITLE
      WRITE (ILP,6010)
 6010 FORMAT (/1X,
     &'                                                           A PR',
     &'IORI                   A POSTERIORI'/1X,
     &'                   MEASUREMENT                           EXPECT',
     &'ATION              BAYESIAN  ESTIMATES'/1X,
     &'             -----------------------                     ------',
     &'-----   -------------------------------------------'/1X,
     &'   H   K   L    F**2     SIGMA(F**2) SIN(TH)/L    D(HKL)  e*<F*',
     &'*2/e>      F**2     SIGMA(F**2)     F      SIGMA(F)'/1X,
     &'   -   -   - ----------- ----------- ---------    ------ ------',
     &'-----   ----------- ----------- --------- ---------')
 6011 FORMAT (1X,3I4,2F12.2,F10.4,F10.2,F12.2,2X,2F12.2,2F10.2)
      NN=0
      NI=0
      NJ=0
      NF=0
      N1REJ=0
      N2REJ=0
      N3REJ=0
      OPEN (UNIT=IO3,STATUS='SCRATCH',FORM='UNFORMATTED',
     & ACCESS='DIRECT',RECL=44)
      IREC=0
      DO 10 I=1,NDATA
        READ  (IO2,REC=I) J,K,L,FSQ,SIGFSQ,DUMMY,DUMMY,DUMMY,DUMMY,
     &  IC,M,IE,S,AVGFSQ,SIGAVG
        XI=FSQ
        SDI=SIGFSQ
        SIGMA=IE*AVGFSQ
        IF (NOBAYES.NE.0) THEN
C
C OPTIONALLY, DO NOT DO THE BAYESIAN REPLACEMENT.
C
          IF (FSQ.LT.0) FSQ=0
          IF (SIGFSQ.LT.0) SIGFSQ=0
          F=SQRT(FSQ)
          IF (FSQ.GE.SIGFSQ.AND.SIGFSQ.GT.0) THEN
            SIGF=SIGFSQ/(2*F)
          ELSE
            SIGF=0.5*SQRT(SIGFSQ)
          END IF
          IF (SIGFSQ.GT.0.AND.SIGF.GT.0) THEN
            XJ=FSQ
            SDJ=SIGFSQ
            XF=F
            SDF=SIGF
            IREJ=0
          ELSE
            IREJ=1
          END IF
        ELSE
C
C CALL FRENCH AND WILSON SUBROUTINES TO GET BAYESIAN ESTIMATES.
C
          IREJ=0
          IF (IC.EQ.0) CALL NEGIA (XI,SDI,SIGMA,XJ,SDJ,XF,SDF,IREJ)
          IF (IC.EQ.1) CALL NEGIC (XI,SDI,SIGMA,XJ,SDJ,XF,SDF,IREJ)
        END IF
C
C ACCUMULATE STATISTICS FOR REJECTED REFLECTIONS.
C
        IF (IREJ.NE.0) THEN
          IF (IREJ.EQ.+1) N1REJ=N1REJ+1
          IF (IREJ.EQ.-1) N2REJ=N2REJ+1
          GO TO 10
        END IF
C
C REJECT REFLECTIONS WITH UNREASONABLY LARGE CHANGES IN LARGER-THAN-
C AVERAGE FSQ.
C
        IF (XJ.GT.SIGMA.AND.ABS(XJ-XI).GT.2*SDI) THEN
          N3REJ=N3REJ+1
          WRITE (IO3,REC=N3REJ) J,K,L,XI,SDI,S,SIGMA,XJ,SDJ,XF,SDF
          GO TO 10
        END IF
C
C RECORD THE ACCEPTED REFLECTIONS IN THE WORKING DATA FILE.
C
        FSQ=XJ
        SIGFSQ=SDJ
        F=XF
        SIGF=SDF
        IREC=IREC+1
        WRITE (IO2,REC=IREC) J,K,L,FSQ,SIGFSQ,F,SIGF, 0.0, 0.0,
     &  IC,M,IE,S,AVGFSQ,SIGAVG
C
C PRINT SELECTED REFELCTIONS.
C
        IF (IPRINT(J,K,L).EQ.1) WRITE (ILP,6011) J,K,L,XI,SDI,S,1/(2*S),
     &  SIGMA,XJ,SDJ,XF,SDF
C
C COMPILE SOME DATA SET STATISTICS.
C
        IF (XI.LT.0)     NN=NN+1
        IF (XI.LT.3*SDI) NI=NI+1
        IF (XJ.LT.3*SDJ) NJ=NJ+1
        IF (XF.LT.3*SDF) NF=NF+1
 10   CONTINUE
      NDATA=NDATA-N1REJ-N2REJ-N3REJ
C
C PRINT DATA SET STATISTICS.
C
      WRITE (ILP,6000) ATIME,ADATE,TITLE
      WRITE (ILP,'(/1X,
     &''N1 = '',I6 ,''  REFLECTIONS REJECTED BECAUSE MEASURED FSQ/SIG'',
     &''MA(FSQ) .LT. -4''/1X,
     &''             (UNREASONABLY NEGATIVE INTENSITY MEASUREMENTS)'')')
     & N1REJ
      WRITE (ILP,'(/1X,
     &''N2 = '',I6,''  REFLECTIONS REJECTED BECAUSE MEASURED FSQ/SIGM'',
     &''A(FSQ) - Q*SIGMA(FSQ)/<FSQ> .LT. -4''/1X,
     &''             WHERE Q = 1 FOR ACENTRIC REFLECTIONS AND Q = 0.5'',
     &'' FOR CENTRIC REFLECTIONS''/1X,
     &''             (UNREASONABLY LARGE EXPERIMENTAL ERROR ESITMATES'',
     &'')'')') N2REJ
      WRITE (ILP,'(/1X,
     &''N3 = '',I6,''  REFLECTIONS REJECTED BECAUSE ABS(DELTA(FSQ)) .'',
     &''GT. 2*SIGMA(FSQ) .AND. FSQ .GT. <FSQ>''/1X,
     &''             (UNREASONABLY LARGE CHANGES IN LARGER-THAN-LOCAL'',
     &''-AVERAGE INTENSITIES)'')') N3REJ
      WRITE (ILP,'(/1X,
     &''NP = '',I6,''  DATA PROCESSED''/1X,
     &''NN = '',I6,''  INPUT DATA WITH  FSQ .LT. 0''/1X,
     &''NI = '',I6,''  INPUT DATA WITH  FSQ .LT. 3*SIGMA(FSQ)''/1X,
     &''NJ = '',I6,''  OUTPUT DATA WITH FSQ .LT. 3*SIGMA(FSQ)''/1X,
     &''NF = '',I6,''  OUTPUT DATA WITH F   .LT. 3*SIGMA(F)'')')
     & NDATA,NN,NI,NJ,NF
      IF (N3REJ.GT.0) THEN
        DO I=1,N3REJ
          READ (IO3,REC=I) J,K,L,XI,SDI,S,SIGMA,XJ,SDJ,XF,SDF
          DATA(I)=XJ/SIGMA
        END DO
        IF (N3REJ.GT.1) THEN
          CALL SORT (N3REJ,DATA,INDEX)
        ELSE
          INDEX(1)=1
        END IF
        WRITE (ILP,'(//1X,
     &''REJECTED LARGER-THAN-LOCAL-AVERAGE INTENSITIES, LISTED IN ORD'',
     &''ER OF DECREASING E = SQRT[FSQ/(EPSILON*<FSQ/EPSILON>)]:'')')
        WRITE (ILP,6010)
        DO I=1,MIN(N3REJ,50)
          J=N3REJ-I+1
          READ (IO3,REC=INDEX(J)) J,K,L,XI,SDI,S,        SIGMA,
     &                                  XJ,SDJ,XF,SDF
          WRITE (ILP,6011)        J,K,L,XI,SDI,S,1/(2*S),SIGMA,
     &                                  XJ,SDJ,XF,SDF
        END DO
        CLOSE (UNIT=IO3,STATUS='DELETE')
      END IF
C
C ESTIMATE LOCALLY NORMALIZED STRUCTURE FACTORS,
C E(HKL) = F(HKL)/SQRT[EPSILON(HKL)*<FSQ/EPSILON>].
C
      CALL NORMAL (NDATA,IO2,IO1,ILP,ATIME,ADATE,TITLE,SMIN,SMAX,GINV,
     & DATA,INDEX,NLOCAL)
C
C WRITE OUTPUT REFLECTION FILE.
C
      OPEN (UNIT=IO1,FILE=FILEO,STATUS='NEW',FORM='FORMATTED')
      DO I=1,NDATA
        READ (IO2,REC=I) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE
        IF (E.GT.0) CALL WRITE1 (IO1,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE)
      END DO
      STOP 'PROGRAM BAYES 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,ILAUE)
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 I=1,42
        IF (PTGP.EQ.SYMBOL(I)) GO TO 1
      END DO
      STOP 'BAYES/DECODE:  ERROR IN POINT GROUP SYMBOL'
 1    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  '
      DO I=1,42
        IF (LAUE.EQ.SYMBOL(I)) GO TO 2
      END DO
 2    ILAUE=I
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-----------------------------------------------------------------------
      SUBROUTINE EQUIV (PTGP,H,K,L)
C
C EQUIVALENT, UNIQUE REFLECTION INDICES FOR THE CRYSTALLOGRAPHIC POINT
C 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 SEE INTERNATIONAL TABLES FOR CRYSTALLOGRAPHY (1983).  VOLUME A, PP.
C 750-770.
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,X,Y,Z
      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
C
C     1
C
    1 RETURN
C
C     -1
C
    2 IF (L.GT.0) RETURN
      H=-H
      K=-K
      L=-L
      IF (L.GT.0) RETURN
      IF (K.GT.0) RETURN
      H=-H
      K=-K
      IF (K.GT.0) RETURN
      IF (H.GE.0) RETURN
      H=-H
      RETURN
C
C     2
C
    3 IF (L.GT.0) RETURN
      H=-H
      L=-L
      IF (L.GT.0) RETURN
      IF (H.GE.0) RETURN
      H=-H
      RETURN
C
C     M
C
    4 K=ABS(K)
      RETURN
C
C     2/M
C
    5 K=ABS(K)
      GO TO 3
C
C     222
C
    6 IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) THEN
        H=ABS(H)
        K=ABS(K)
        L=ABS(L)
        RETURN
      END IF
      IF (H.GE.0.AND.K.GE.0) RETURN
      X=H
      Y=K
      Z=L
      H=-X
      K=-Y
      L=Z
      IF (H.GE.0.AND.K.GE.0) RETURN
      H=-X
      K=Y
      L=-Z
      IF (H.GE.0.AND.K.GE.0) RETURN
      H=X
      K=-Y
      L=-Z
      RETURN
C
C     MM2
C
    7 H=ABS(H)
      K=ABS(K)
      RETURN
C
C     MMM
C
    8 L=ABS(L)
      GO TO 7
C
C     4
C
    9 IF (H.EQ.0.AND.K.EQ.0) RETURN
      IF (H.GT.0.AND.K.GE.0) RETURN
      X=H
      Y=K
      H=-Y
      K=X
      IF (H.GT.0.AND.K.GE.0) RETURN
      H=-X
      K=-Y
      IF (H.GT.0.AND.K.GE.0) RETURN
      H=Y
      K=-X
      RETURN
C
C     -4
C
   10 IF (H.EQ.0.AND.K.EQ.0) THEN
        L=ABS(L)
        RETURN
      END IF
      IF (H.GT.0.AND.K.GE.0) RETURN
      X=H
      Y=K
      H=-X
      K=-Y
      IF (H.GT.0.AND.K.GE.0) RETURN
      H=Y
      K=-X
      L=-L
      IF (H.GT.0.AND.K.GE.0) RETURN
      H=-Y
      K=X
      RETURN
C
C     4/M
C
   11 L=ABS(L)
      GO TO 9
C
C     422
C
   12 IF (H.EQ.0.OR.K.EQ.0.OR.ABS(H).EQ.ABS(K)) L=ABS(L)
  512 IF (H.GE.K.AND.K.GE.0) RETURN
      X=H
      Y=K
      H=-X
      K=-Y
      IF (H.GE.K.AND.K.GE.0) RETURN
      H=-Y
      K=X
      IF (H.GE.K.AND.K.GE.0) RETURN
      H=Y
      K=-X
      IF (H.GE.K.AND.K.GE.0) RETURN
      H=-X
      K=Y
      L=-L
      GO TO 512
C
C     4MM
C
   13 H=ABS(H)
      K=ABS(K)
      IF (H.GE.K) RETURN
      X=H
      H=K
      K=X
      RETURN
C
C     -42M
C
   14 IF (H.EQ.0.OR.K.EQ.0) L=ABS(L)
  514 IF (H.GE.K.AND.K.GE.0) RETURN
      X=H
      Y=K
      H=-X
      K=-Y
      IF (H.GE.K.AND.K.GE.0) RETURN
      H=Y
      K=-X
      L=-L
      IF (H.GE.K.AND.K.GE.0) RETURN
      H=-Y
      K=X
      IF (H.GE.K.AND.K.GE.0) RETURN
      H=-X
      K=Y
      GO TO 514
C
C     -4M2
C
   15 H=ABS(H)
      K=ABS(K)
      IF (H.EQ.K) L=ABS(L)
      IF (H.GE.K) RETURN
      X=H
      H=K
      K=X
      L=-L
      RETURN
C
C     4/MMM
C
   16 L=ABS(L)
      GO TO 13
C
C     3
C
   17 IF (H.EQ.0.AND.K.EQ.0) RETURN
      IF (H.LT.0.AND.K.GE.0) RETURN
      X=H
      Y=K
      I=-H-K
      H=I
      K=X
      IF (H.LT.0.AND.K.GE.0) RETURN
      H=Y
      K=I
      RETURN
C
C     -3
C
   18 IF (H.EQ.0.AND.K.EQ.0) THEN
        L=ABS(L)
        RETURN
      END IF
  518 IF (H.LT.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.LT.-K))) RETURN
      X=H
      Y=K
      I=-H-K
      H=I
      K=X
      IF (H.LT.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.LT.-K))) RETURN
      H=Y
      K=I
      IF (H.LT.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.LT.-K))) RETURN
      H=-X
      K=-Y
      L=-L
      GO TO 518
C
C     312
C
   19 IF (H.EQ.0.OR.K.EQ.0.OR.-H.EQ.K) L=ABS(L)
  519 IF (H.GE.0.AND.K.GE.0) RETURN
      X=H
      Y=K
      I=-H-K
      H=I
      K=X
      IF (H.GE.0.AND.K.GE.0) RETURN
      H=Y
      K=I
      IF (H.GE.0.AND.K.GE.0) RETURN
      H=-Y
      K=-X
      L=-L
      GO TO 519
C
C     321
C
   20 IF (H.EQ.K.OR.-2*H.EQ.K.OR.-H.EQ.2*K) L=ABS(L)
  520 IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN
      X=H
      Y=K
      I=-H-K
      H=I
      K=X
      IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN
      H=Y
      K=I
      IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN
      H=Y
      K=X
      L=-L
      GO TO 520
C
C     31M
C
   21 IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN
      X=H
      Y=K
      I=-H-K
      H=I
      K=X
      IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN
      H=Y
      K=I
      IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN
      H=Y
      K=X
      GO TO 21
C
C     3M1
C
   22 IF (H.GE.0.AND.K.GE.0) RETURN
      X=H
      Y=K
      I=-H-K
      H=I
      K=X
      IF (H.GE.0.AND.K.GE.0) RETURN
      H=Y
      K=I
      IF (H.GE.0.AND.K.GE.0) RETURN
      H=-Y
      K=-X
      GO TO 22
C
C     -31M
C
   23 IF (H.EQ.0.OR.K.EQ.0.OR.-H.EQ.K) L=ABS(L)
  523 IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H.AND.(L.GT.0.OR.(L.EQ.0.AND.
     & K.GE.0))) RETURN
      X=H
      Y=K
      I=-H-K
      H=I
      K=X
      IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H.AND.(L.GT.0.OR.(L.EQ.0.AND.
     & K.GE.0))) RETURN
      H=Y
      K=I
      IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H.AND.(L.GT.0.OR.(L.EQ.0.AND.
     & K.GE.0))) RETURN
      H=Y
      K=X
      IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H.AND.(L.GT.0.OR.(L.EQ.0.AND.
     & K.GE.0))) RETURN
      H=X
      K=I
      IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H.AND.(L.GT.0.OR.(L.EQ.0.AND.
     & K.GE.0))) RETURN
      H=I
      K=Y
      IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H.AND.(L.GT.0.OR.(L.EQ.0.AND.
     & K.GE.0))) RETURN
      H=-X
      K=-Y
      L=-L
      GO TO 523
C
C     -3M1
C
   24 IF (H.EQ.K.OR.-2*H.EQ.K.OR.-H.EQ.2*K) L=ABS(L)
  524 IF (H.GE.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.GE.K))) RETURN
      X=H
      Y=K
      I=-H-K
      H=I
      K=X
      IF (H.GE.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.GE.K))) RETURN
      H=Y
      K=I
      IF (H.GE.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.GE.K))) RETURN
      H=Y
      K=X
      L=-L
      IF (H.GE.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.GE.K))) RETURN
      H=X
      K=I
      IF (H.GE.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.GE.K))) RETURN
      H=I
      K=Y
      IF (H.GE.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.GE.K))) RETURN
      H=-X
      K=-Y
      GO TO 524
C
C     6
C
   25 IF (H.EQ.0.AND.K.EQ.0) RETURN
      IF (H.GT.0.AND.K.GE.0) RETURN
      X=H
      Y=K
      I=-H-K
      H=I
      K=X
      IF (H.GT.0.AND.K.GE.0) RETURN
      H=Y
      K=I
      IF (H.GT.0.AND.K.GE.0) RETURN
      H=-X
      K=-Y
      IF (H.GT.0.AND.K.GE.0) RETURN
      H=-I
      K=-X
      IF (H.GT.0.AND.K.GE.0) RETURN
      H=-Y
      K=-I
      RETURN
C
C     -6
C
   26 L=ABS(L)
      GO TO 17
C
C     6/M
C
   27 L=ABS(L)
      GO TO 25
C
C     622
C
   28 IF (H.EQ.0.OR.K.EQ.0.OR.ABS(H).EQ.ABS(K).OR.ABS(2*H).EQ.ABS(K).OR.
     & ABS(H).EQ.ABS(2*K)) L=ABS(L)
  528 IF (H.GE.K.AND.K.GE.0) RETURN
      X=H
      Y=K
      I=-H-K
      H=I
      K=X
      IF (H.GE.K.AND.K.GE.0) RETURN
      H=Y
      K=I
      IF (H.GE.K.AND.K.GE.0) RETURN
      H=-X
      K=-Y
      IF (H.GE.K.AND.K.GE.0) RETURN
      H=-I
      K=-X
      IF (H.GE.K.AND.K.GE.0) RETURN
      H=-Y
      K=-I
      IF (H.GE.K.AND.K.GE.0) RETURN
      H=Y
      K=X
      L=-L
      GO TO 528
C
C     6MM
C
   29 IF (H.GE.K.AND.K.GE.0) RETURN
      X=H
      Y=K
      I=-H-K
      H=I
      K=X
      IF (H.GE.K.AND.K.GE.0) RETURN
      H=Y
      K=I
      IF (H.GE.K.AND.K.GE.0) RETURN
      H=-X
      K=-Y
      IF (H.GE.K.AND.K.GE.0) RETURN
      H=-I
      K=-X
      IF (H.GE.K.AND.K.GE.0) RETURN
      H=-Y
      K=-I
      IF (H.GE.K.AND.K.GE.0) RETURN
      H=Y
      K=X
      GO TO 29
C
C     -6M2
C
   30 L=ABS(L)
      IF (H.EQ.0.AND.K.EQ.0) RETURN
  530 IF (H.GT.0.AND.K.GE.0) RETURN
      X=H
      Y=K
      I=-H-K
      H=I
      K=X
      IF (H.GT.0.AND.K.GE.0) RETURN
      H=Y
      K=I
      IF (H.GT.0.AND.K.GE.0) RETURN
      H=-Y
      K=-X
      GO TO 530
C
C     -62M
C
   31 L=ABS(L)
      IF (H.EQ.0.AND.K.EQ.0) RETURN
  531 IF (H.GT.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN
      X=H
      Y=K
      I=-H-K
      H=I
      K=X
      IF (H.GT.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN
      H=Y
      K=I
      IF (H.GT.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN
      H=Y
      K=X
      GO TO 531
C
C     6/MMM
C
   32 L=ABS(L)
      GO TO 29
C
C     R3
C
   33 STOP 'BAYES/EQUIV:  RE-INDEX ON HEXAGONAL AXES'
C
C     R-3
C
   34 GO TO 33
C
C     R32
C
   35 GO TO 33
C
C     R3M
C
   36 GO TO 33
C
C     R-3M
C
   37 GO TO 33
C
C     23
C
   38 IF (H.EQ.K.AND.K.EQ.ABS(L)) RETURN
  538 IF (H.GT.ABS(L).AND.K.GE.ABS(L)) RETURN
      X=H
      Y=K
      Z=L
      H=-X
      K=-Y
      L=Z
      IF (H.GT.ABS(L).AND.K.GE.ABS(L)) RETURN
      H=-X
      K=Y
      L=-Z
      IF (H.GT.ABS(L).AND.K.GE.ABS(L)) RETURN
      H=X
      K=-Y
      L=-Z
      IF (H.GT.ABS(L).AND.K.GE.ABS(L)) RETURN
      H=Z
      K=X
      L=Y
      GO TO 538
C
C     M3
C
   39 H=ABS(H)
      K=ABS(K)
      L=ABS(L)
      IF (H.EQ.K.AND.K.EQ.L) RETURN
  539 IF (H.GT.L.AND.K.GE.L) RETURN
      X=H
      Y=K
      Z=L
      H=Z
      K=X
      L=Y
      IF (H.GT.L.AND.K.GE.L) RETURN
      H=Y
      K=Z
      L=X
      GO TO 539
C
C     432
C
   40 IF (ABS(H).EQ.ABS(K).AND.ABS(K).EQ.ABS(L)) THEN
        H=ABS(H)
        K=ABS(K)
        L=ABS(L)
        RETURN
      END IF
  540 IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN
      X=H
      Y=K
      Z=L
      H=Z
      K=X
      L=Y
      IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN
      H=Y
      K=Z
      L=X
      IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN
      H=-X
      K=-Y
      L=Z
      IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN
      H=Z
      K=-X
      L=-Y
      IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN
      H=-Y
      K=Z
      L=-X
      IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN
      H=-X
      K=Y
      L=-Z
      IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN
      H=-Z
      K=-X
      L=Y
      IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN
      H=Y
      K=-Z
      L=-X
      IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN
      H=X
      K=-Y
      L=-Z
      IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN
      H=-Z
      K=X
      L=-Y
      IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN
      H=-Y
      K=-Z
      L=X
      IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN
      H=Y
      K=X
      L=-Z
      GO TO 540
C
C     -43M
C
   41 IF (H.GE.K.AND.K.GE.ABS(L)) RETURN
      X=H
      Y=K
      Z=L
      H=Z
      K=X
      L=Y
      IF (H.GE.K.AND.K.GE.ABS(L)) RETURN
      H=Y
      K=Z
      L=X
      IF (H.GE.K.AND.K.GE.ABS(L)) RETURN
      H=-X
      K=-Y
      L=Z
      IF (H.GE.K.AND.K.GE.ABS(L)) RETURN
      H=Z
      K=-X
      L=-Y
      IF (H.GE.K.AND.K.GE.ABS(L)) RETURN
      H=-Y
      K=Z
      L=-X
      IF (H.GE.K.AND.K.GE.ABS(L)) RETURN
      H=-X
      K=Y
      L=-Z
      IF (H.GE.K.AND.K.GE.ABS(L)) RETURN
      H=-Z
      K=-X
      L=Y
      IF (H.GE.K.AND.K.GE.ABS(L)) RETURN
      H=Y
      K=-Z
      L=-X
      IF (H.GE.K.AND.K.GE.ABS(L)) RETURN
      H=X
      K=-Y
      L=-Z
      IF (H.GE.K.AND.K.GE.ABS(L)) RETURN
      H=-Z
      K=X
      L=-Y
      IF (H.GE.K.AND.K.GE.ABS(L)) RETURN
      H=-Y
      K=-Z
      L=X
      IF (H.GE.K.AND.K.GE.ABS(L)) RETURN
      H=Y
      K=X
      L=Z
      GO TO 41
C
C     M3M
C
   42 H=ABS(H)
      K=ABS(K)
      L=ABS(L)
      IF (H.GE.K.AND.K.GE.L) RETURN
      X=H
      Y=K
      Z=L
      H=Z
      K=X
      L=Y
      IF (H.GE.K.AND.K.GE.L) RETURN
      H=Y
      K=Z
      L=X
      IF (H.GE.K.AND.K.GE.L) RETURN
      H=Z
      K=Y
      L=X
      IF (H.GE.K.AND.K.GE.L) RETURN
      H=X
      K=Z
      L=Y
      IF (H.GE.K.AND.K.GE.L) RETURN
      H=Y
      K=X
      L=Z
      RETURN
      END
C------------------------------------------------------------------------
      FUNCTION ERRF(X)
C
C RETURNS THE ERROR FUNCTION
C
C ERF(X) = [2/SQRT(PI)]*INTEGRAL[0,ABS(X)] EXP(-T**2) DT
C
C RATIONAL APPROXIMATION FORMULA NO. 7.1.26 FROM:
C
C MILTON ABRAMOWITZ AND IRENE A. STEGUN (1965).  HANDBOOK OF
C MATHEMATICAL FUNCTIONS, P. 299.  NEW YORK:  DOVER PUBLICATIONS, INC.
C
      REAL*8 Z,ONE,T,P,A1,A2,A3,A4,A5
      DATA ONE,P,A1,A2,A3,A4,A5 /1.D0, 0.3275911D0, 0.254829592D0,
     & -0.284496736D0, 1.421413741D0, -1.43152027D0, 1.061405429D0/
      Z=ABS(X)
      T=ONE/(ONE+P*Z)
      T=T*(A1+T*(A2+T*(A3+T*(A4+T*A5))))
      ERRF=ONE-T*EXP(-Z**2)
      IF (X.LT.0) ERRF=-ERRF
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE ESTATS (NREC,IFILE,ILP,ATIME,ADATE,TITLE,SMIN,SMAX)
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)
      CHARACTER ATIME*8,ADATE*9,TITLE*80
      WRITE (ILP,'(/1H1,130(''-'')/1X,
     &''SUBPROGRAM NORMAL/BAYES.  '',A,'', '',A,''.  '',A)')
     & ATIME,ADATE,TITLE
      WRITE (ILP,'(/1X,
     &''E-STATISTICS FOR NORMALIZED STRUCTURE FACTORS ESTIMATED AS E('',
     &''HKL) = F(HKL)/SQRT[EPSILON(HKL)*<FSQ/EPSILON>]'')')
C
C MOMENTS OF THE E-DISTRIBUTION
C
      NA=0
      NC=0
      DO IREC=1,NREC
        READ (IFILE,REC=IREC) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE,
     &  IC,M,IE,S,AVGFSQ
        IF (SMIN.LE.S.AND.S.LE.SMAX.AND.E.GT.0) THEN
          DO 1 I=1,7
            IF (I.EQ.2.AND.IC.NE.0) GO TO 1
            IF (I.EQ.3.AND.IC.NE.1) GO TO 1
            IF (I.EQ.4.AND.(J.EQ.0.OR.K.EQ.0.OR.L.EQ.0)) GO TO 1
            IF (I.EQ.5.AND.J.NE.0) GO TO 1
            IF (I.EQ.6.AND.K.NE.0) GO TO 1
            IF (I.EQ.7.AND.L.NE.0) GO TO 1
            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
 1        CONTINUE
          IF (IC.EQ.0) NA=NA+1
          IF (IC.EQ.1) NC=NC+1
        END IF
      END DO
      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,
     &'                      -----------------------------------------',
     &'---------------------------  -------------------------------'/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 DSTAR = 2*SIN(THETA)/LAMBDA
C
      N=10
      DO I=1,N
        SI(I)=(FLOAT(I)/N)**(1/3.0)*SMAX
      END DO
      DO 2 IREC=1,NREC
        READ (IFILE,REC=IREC) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE,
     &  IC,M,IE,S,AVGFSQ
        IF (S.LT.SMIN.OR.S.GT.SMAX.OR.E.LE.0) GO TO 2
        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
 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
      DO IREC=1,NREC
        READ (IFILE,REC=IREC) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE,
     &  IC,M,IE,S,AVGFSQ
        IF (SMIN.LE.S.AND.S.LE.SMAX.AND.E.GT.0) THEN
          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,
     &'    ------    ------       --- ---------- ---------- ----------',
     &'----'/1X,
     &'       EEE',2F10.3,2F10.2,I10/1X,
     &'       EEO',2F10.3,2F10.2,I10/1X,
     &'       EOE',2F10.3,2F10.2,I10/1X,
     &'       OEE',2F10.3,2F10.2,I10/1X,
     &'       EOO',2F10.3,2F10.2,I10/1X,
     &'       OEO',2F10.3,2F10.2,I10/1X,
     &'       OOE',2F10.3,2F10.2,I10/1X,
     &'       OOO',2F10.3,2F10.2,I10)
C
C PLOT 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
        DO 10 IREC=1,NREC
          READ (IFILE,REC=IREC) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE,
     &    IC,M,IE,S,AVGFSQ
          IF (S.LT.SMIN.OR.S.GT.SMAX.OR.E.LE.0) GO TO 10
          IF (IC.EQ.1) GO TO 10
          DO I=1,N
            IF (E.LE.I*DX) THEN
              X(I)=X(I)+M*E
              Y(I)=Y(I)+M
              GO TO 10
            END IF
          END DO
 10     CONTINUE
        DO I=1,N
          IF (Y(I).GT.0) X(I)=X(I)/Y(I)
        END DO
        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,6000) ATIME,ADATE,TITLE
        CALL PLOTA (N,X,Y,ILP)
      END IF
 6000 FORMAT (/1H1,130('-')/1X,10X,
     &'SUBPROGRAM NORMAL/BAYES.  ',A,', ',A,'.  ',A/)
      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
        DO 11 IREC=1,NREC
          READ (IFILE,REC=IREC) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE,
     &    IC,M,IE,S,AVGFSQ
          IF (S.LT.SMIN.OR.S.GT.SMAX.OR.E.LE.0) GO TO 11
          IF (IC.EQ.0) GO TO 11
          DO I=1,N
            IF (E.LE.I*DX) THEN
              X(I)=X(I)+M*E
              Y(I)=Y(I)+M
              GO TO 11
            END IF
          END DO
 11     CONTINUE
        DO I=1,N
          IF (Y(I).GT.0) X(I)=X(I)/Y(I)
        END DO
        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,6000) ATIME,ADATE,TITLE
        CALL PLOTC (N,X,Y,ILP)
      END IF
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
      DO 3 IREC=1,NREC
        READ (IFILE,REC=IREC) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE,
     &  IC,M,IE,S,AVGFSQ
        IF (S.LT.SMIN.OR.S.GT.SMAX.OR.E.LE.0) GO TO 3
        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 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
      WRITE (ILP,6000) ATIME,ADATE,TITLE
      CALL PLOTS (N,X,Y,SMAX,ILP)
      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
      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-----------------------------------------------------------------------
      FUNCTION IPRINT(IH,IK,IL)
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 METRIC (A,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 NEGIA (XI,SDI,SIGMA,XJ,SDJ,XF,SDF,IREJ)
C
C     FROM SUPPLEMENTARY PUBLICATION NO. SUP 33352 FROM:
C
C     SIMON FRENCH AND KEITH WILSON (1979).  ON THE TREATMENT OF
C     NEGATIVE INTENSITY OBSERVATIONS.  ACTA CRYST. A34, 517-525.
C
C     GIVEN AN OBSERVED INTENSITY AND ITS ASSOCIATED STANDARD DEVIATION,
C     THE SUBROUTINE CALCULATES AND RETURNS THE POSTERIOR MEAN AND
C     STANDARD DEVIATION OF BOTH THE TRUE INTENSITY AND THE TRUE
C     STRUCTURE FACTOR MODULUS.
C
C     THE PRIOR DISTRIBUTION USED IS WILSON'S DISTRIBUTION FOR AN
C     ACENTRIC REFLECTION.  FOR EXACT FORM SEE FRENCH AND WILSON (1977),
C     EQUATION (2.2).
C
C     ACCURACY - PROBABLY BETTER THAN 5%
C
C     ARGUMENTS:
C
C     XI    - OBSERVED INTENSITY
C     SDI   - S.D. OF OBSERVATION
C     SIGMA - PARAMETER FOR WILSON'S DISTRIBUTION
C     XJ    - POSTERIOR MEAN FOR THE TRUE INTENSITY
C     SDJ   - POSTERIOR S.D. FOR THE TRUE INTENSITY
C     XF    - POSTERIOR MEAN FOR THE TRUE STRUCTURE FACTOR MODULUS
C     SDF   - POSTERIOR S.D. FOR THE TRUE STRUCTURE FACTOR MODULUS
C     IREJ  - ON ENTRY: = -1  IF ONLY XF AND SDF REQUIRED
C                       =  0  IF ALL XJ, SDJ, XF, AND SDJ REQUIRED
C                       = +1  IF ONLY XJ AND SDJ REQUIRED
C             ON EXIT:  = -1  IF XI/SDI - SDI/SIGMA < -4
C                       =  0  IF ALL O.K.
C                       = +1  IF XI/SDI < -4, I.E., OBSERVATION IS
C                             ESSENTIALLY IMPOSSIBLE
C
C     ARRAYS ZJ, ZJSD, ZF, ZFSD CONTAIN THE TABULATED VALUES OF XJ, SDJ,
C     XF, AND SDF FOR XI/SDI - SDI/SIGMA = -4.0, (0.1), 3.0.  LINEAR
C     INTERPOLATION IS USED IN THESE TABLES.
C
      DIMENSION ZJ(71),ZJSD(71),ZF(71),ZFSD(71)
      DATA ZJ   /0.226, 0.230, 0.235, 0.240, 0.246, 0.251, 0.257,
     &           0.263, 0.270, 0.276, 0.283, 0.290, 0.298, 0.306,
     &           0.314, 0.323, 0.332, 0.341, 0.351, 0.362, 0.373,
     &           0.385, 0.397, 0.410, 0.424, 0.439, 0.454, 0.470,
     &           0.487, 0.505, 0.525, 0.545, 0.567, 0.590, 0.615,
     &           0.641, 0.668, 0.698, 0.729, 0.762, 0.798, 0.835,
     &           0.875, 0.917, 0.962, 1.009, 1.059, 1.112, 1.167,
     &           1.226, 1.287, 1.352, 1.419, 1.490, 1.563, 1.639,
     &           1.717, 1.798, 1.882, 1.967, 2.055, 2.145, 2.236,
     &           2.329, 2.422, 2.518, 2.614, 2.710, 2.808, 2.906,
     &           3.004/
      DATA ZJSD /0.217, 0.221, 0.226, 0.230, 0.235, 0.240, 0.245,
     &           0.250, 0.255, 0.261, 0.267, 0.273, 0.279, 0.286,
     &           0.292, 0.299, 0.307, 0.314, 0.322, 0.330, 0.339,
     &           0.348, 0.357, 0.367, 0.377, 0.387, 0.398, 0.409,
     &           0.421, 0.433, 0.446, 0.459, 0.473, 0.488, 0.503,
     &           0.518, 0.535, 0.551, 0.568, 0.586, 0.604, 0.622,
     &           0.641, 0.660, 0.679, 0.698, 0.718, 0.737, 0.757,
     &           0.776, 0.795, 0.813, 0.831, 0.848, 0.865, 0.881,
     &           0.895, 0.909, 0.921, 0.933, 0.943, 0.953, 0.961,
     &           0.968, 0.974, 0.980, 0.984, 0.988, 0.991, 0.994,
     &           0.996/
      DATA ZF   /0.423, 0.428, 0.432, 0.437, 0.442, 0.447, 0.453,
     &           0.458, 0.464, 0.469, 0.475, 0.482, 0.488, 0.495,
     &           0.502, 0.509, 0.516, 0.524, 0.532, 0.540, 0.549,
     &           0.557, 0.567, 0.576, 0.586, 0.597, 0.608, 0.619,
     &           0.631, 0.643, 0.656, 0.670, 0.684, 0.699, 0.714,
     &           0.730, 0.747, 0.765, 0.783, 0.802, 0.822, 0.843,
     &           0.865, 0.887, 0.911, 0.935, 0.960, 0.987, 1.014,
     &           1.042, 1.070, 1.100, 1.130, 1.161, 1.192, 1.224,
     &           1.257, 1.289, 1.322, 1.355, 1.388, 1.421, 1.454,
     &           1.487, 1.519, 1.551, 1.583, 1.615, 1.646, 1.676,
     &           1.706/
      DATA ZFSD /0.216, 0.218, 0.220, 0.222, 0.224, 0.226, 0.229,
     &           0.231, 0.234, 0.236, 0.239, 0.241, 0.244, 0.247,
     &           0.250, 0.253, 0.256, 0.259, 0.262, 0.266, 0.269,
     &           0.272, 0.276, 0.279, 0.283, 0.287, 0.291, 0.295,
     &           0.298, 0.302, 0.307, 0.311, 0.315, 0.319, 0.324,
     &           0.328, 0.332, 0.337, 0.341, 0.345, 0.349, 0.353,
     &           0.357, 0.360, 0.364, 0.367, 0.369, 0.372, 0.374,
     &           0.375, 0.376, 0.377, 0.377, 0.377, 0.376, 0.374,
     &           0.372, 0.369, 0.366, 0.362, 0.358, 0.353, 0.348,
     &           0.343, 0.338, 0.332, 0.327, 0.321, 0.315, 0.310,
     &           0.304/
C
      IOP=IREJ
      IREJ=0
C
C     IS OBSERVATION POSSIBLE?
C
      H=XI/SDI
      IF (H+4.0) 5,8,8
C
    5 CONTINUE
      IREJ=1
      GO TO 15
C
C     LOCATE RANGE OF DISTRIBUTION.
C
    8 CONTINUE
      H=H-SDI/SIGMA
      IF (H+4.0) 10,20,20
C
C     DISTRIBUTION IS ESSENTIALLY NEGATIVE.
C
   10 CONTINUE
      IREJ=-1
   15 CONTINUE
      XJ=0.0
      SDJ=1.0
      XF=0.0
      SDF=1.0
      RETURN
C
   20 CONTINUE
      IF (H-3.0) 60,30,30
C
C     DISTRIBUTION IS ESSENTIALLY ALL POSTIVE.
C
   30 CONTINUE
      H=H*SDI
      IF (IOP) 40,35,35
   35 CONTINUE
      XJ=H
      SDJ=SDI
   40 CONTINUE
      IF (IOP) 45,45,50
   45 CONTINUE
      XF=SQRT(H)
      SDF=0.5*SDI/XF
   50 CONTINUE
      RETURN
C
C     LOCATE XI IN THE TABLES.
C
   60 CONTINUE
      H=H+4.0
      IH=H*10.0
      A=(H-IH*0.1)*10.0
      B=1.0-A
      IH=IH+1
      IHP1=IH+1
C
C     INTERPOLATE IN THE TABLES.
C
      IF (IOP) 70,65,65
   65 CONTINUE
      XJ=(A*ZJ(IHP1)+B*ZJ(IH))*SDI
      SDJ=(A*ZJSD(IHP1)+B*ZJSD(IH))*SDI
   70 CONTINUE
      IF (IOP) 75,75,80
   75 CONTINUE
      SDF=SQRT(SDI)
      XF=(A*ZF(IHP1)+B*ZF(IH))*SDF
      SDF=(A*ZFSD(IHP1)+B*ZFSD(IH))*SDF
   80 CONTINUE
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE NEGIC (XI,SDI,SIGMA,XJ,SDJ,XF,SDF,IREJ)
C
C     FROM SUPPLEMENTARY PUBLICATION NO. SUP 33352 FROM:
C
C     SIMON FRENCH AND KEITH WILSON (1979).  ON THE TREATMENT OF
C     NEGATIVE INTENSITY OBSERVATIONS.  ACTA CRYST. A34, 517-525.
C
C     GIVEN AN OBSERVED INTENSITY AND ITS ASSOCIATED STANDARD DEVIATION,
C     THE SUBROUTINE CALCULATES AND RETURNS THE POSTERIOR MEAN AND
C     STANDARD DEVIATION OF BOTH THE TRUE INTENSITY AND THE TRUE
C     STRUCTURE FACTOR MODULUS.
C
C     THE PRIOR DISTRIBUTION USED IS WILSON'S DISTRIBUTION FOR A
C     CENTRIC REFLECTION.  FOR EXACT FORM SEE FRENCH AND WILSON (1977),
C     EQUATION (2.3).
C
C     ACCURACY - PROBABLY BETTER THAN 5%
C
C     ARGUMENTS:
C
C     XI    - OBSERVED INTENSITY
C     SDI   - S.D. OF OBSERVATION
C     SIGMA - PARAMETER FOR WILSON'S DISTRIBUTION
C     XJ    - POSTERIOR MEAN FOR THE TRUE INTENSITY
C     SDJ   - POSTERIOR S.D. FOR THE TRUE INTENSITY
C     XF    - POSTERIOR MEAN FOR THE TRUE STRUCTURE FACTOR MODULUS
C     SDF   - POSTERIOR S.D. FOR THE TRUE STRUCTURE FACTOR MODULUS
C     IREJ  - ON ENTRY: = -1  IF ONLY XF AND SDF REQUIRED
C                       =  0  IF ALL XJ, SDJ, XF, AND SDJ REQUIRED
C                       = +1  IF ONLY XJ AND SDJ REQUIRED
C             ON EXIT:  = -1  IF XI/SDI - SDI/(2*SIGMA) < -4
C                       =  0  IF ALL O.K.
C                       = +1  IF XI/SDI < -4, I.E., OBSERVATION IS
C                             ESSENTIALLY IMPOSSIBLE
C
C     ARRAYS ZJ, ZJSD, ZF, ZFSD CONTAIN THE TABULATED VALUES OF XJ, SDJ,
C     XF, AND SDF FOR XI/SDI - SDI/(2*SIGMA) = -4.0, (0.1), 4.0.  LINEAR
C     INTERPOLATION IS USED IN THESE TABLES.
C
      DIMENSION ZJ(81),ZJSD(81),ZF(81),ZFSD(81)
      DATA ZJ   /0.114, 0.116, 0.119, 0.122, 0.124, 0.127, 0.130, 0.134,
     &           0.137, 0.141, 0.145, 0.148, 0.153, 0.157, 0.162, 0.166,
     &           0.172, 0.177, 0.183, 0.189, 0.195, 0.202, 0.209, 0.217,
     &           0.225, 0.234, 0.243, 0.253, 0.263, 0.275, 0.287, 0.300,
     &           0.314, 0.329, 0.345, 0.363, 0.382, 0.402, 0.425, 0.449,
     &           0.475, 0.503, 0.534, 0.567, 0.603, 0.642, 0.684, 0.730,
     &           0.779, 0.833, 0.890, 0.952, 1.018, 1.089, 1.164, 1.244,
     &           1.327, 1.416, 1.508, 1.603, 1.703, 1.805, 1.909, 2.015,
     &           2.123, 2.233, 2.343, 2.453, 2.564, 2.674, 2.784, 2.894,
     &           3.003, 3.112, 3.220, 3.328, 3.435, 3.541, 3.647, 3.753,
     &           3.962/
      DATA ZJSD/ 0.158, 0.161, 0.165, 0.168, 0.172, 0.176, 0.179, 0.184,
     &           0.188, 0.192, 0.197, 0.202, 0.207, 0.212, 0.218, 0.224,
     &           0.230, 0.236, 0.243, 0.250, 0.257, 0.265, 0.273, 0.282,
     &           0.291, 0.300, 0.310, 0.321, 0.332, 0.343, 0.355, 0.368,
     &           0.382, 0.397, 0.412, 0.428, 0.445, 0.463, 0.481, 0.501,
     &           0.521, 0.543, 0.565, 0.589, 0.613, 0.638, 0.664, 0.691,
     &           0.718, 0.745, 0.773, 0.801, 0.828, 0.855, 0.881, 0.906,
     &           0.929, 0.951, 0.971, 0.989, 1.004, 1.018, 1.029, 1.038,
     &           1.044, 1.049, 1.052, 1.054, 1.054, 1.053, 1.051, 1.049,
     &           1.047, 1.044, 1.041, 1.039, 1.036, 1.034, 1.031, 1.029,
     &           1.028/
      DATA ZF   /0.269, 0.272, 0.276, 0.279, 0.282, 0.286, 0.289, 0.293,
     &           0.297, 0.301, 0.305, 0.309, 0.314, 0.318, 0.323, 0.328,
     &           0.333, 0.339, 0.344, 0.350, 0.356, 0.363, 0.370, 0.377,
     &           0.384, 0.392, 0.400, 0.409, 0.418, 0.427, 0.438, 0.448,
     &           0.460, 0.471, 0.484, 0.498, 0.512, 0.527, 0.543, 0.560,
     &           0.578, 0.597, 0.618, 0.639, 0.662, 0.687, 0.713, 0.740,
     &           0.769, 0.800, 0.832, 0.866, 0.901, 0.938, 0.976, 1.016,
     &           1.057, 1.098, 1.140, 1.183, 1.227, 1.270, 1.313, 1.356,
     &           1.398, 1.439, 1.480, 1.519, 1.558, 1.595, 1.632, 1.667,
     &           1.701, 1.735, 1.767, 1.799, 1.829, 1.859, 1.889, 1.917,
     &           1.945/
      DATA ZFSD /0.203, 0.205, 0.207, 0.209, 0.211, 0.214, 0.216, 0.219,
     &           0.222, 0.224, 0.227, 0.230, 0.233, 0.236, 0.239, 0.243,
     &           0.246, 0.250, 0.253, 0.257, 0.261, 0.265, 0.269, 0.273,
     &           0.278, 0.283, 0.288, 0.293, 0.298, 0.303, 0.309, 0.314,
     &           0.320, 0.327, 0.333, 0.340, 0.346, 0.353, 0.361, 0.368,
     &           0.375, 0.383, 0.390, 0.398, 0.405, 0.413, 0.420, 0.427,
     &           0.433, 0.440, 0.445, 0.450, 0.454, 0.457, 0.459, 0.460,
     &           0.460, 0.458, 0.455, 0.451, 0.445, 0.438, 0.431, 0.422,
     &           0.412, 0.402, 0.392, 0.381, 0.370, 0.360, 0.349, 0.339,
     &           0.330, 0.321, 0.312, 0.304, 0.297, 0.290, 0.284, 0.278,
     &           0.272/
C
      IOP=IREJ
      IREJ=0
C
C     IS OBSERVATION POSSIBLE?
C
      H=XI/SDI
      IF (H+4.0) 5,8,8
C
    5 CONTINUE
      IREJ=1
      GO TO 15
C
    8 CONTINUE
      H=H-0.5*SDI/SIGMA
      IF (H+4.0) 10,20,20
C
   10 CONTINUE
      IREJ=-1
   15 CONTINUE
      XJ=0.0
      SDJ=1.0
      XF=0.0
      SDF=1.0
      RETURN
C
   20 CONTINUE
      IF (H-4.0) 60,30,30
C
C     DISTRIBUTION IS ESSENTIALLY ALL POSTIVE.
C
   30 CONTINUE
      XI2=1.0/(H*H)
      XI4=XI2*XI2
      XXF=SQRT(H)*(1.0-0.375*XI2-0.6796875*XI4)
      XSDF=SQRT(H*(0.25*XI2+0.46875*XI4))
      IF (IOP) 40,35,35
   35 CONTINUE
      XJ=H*(1.0-0.5*XI2-0.75*XI4)
      SDJ=SDI*XSDF*2.0*XXF
      XJ=XJ*SDI
   40 CONTINUE
      IF (IOP) 45,45,50
   45 CONTINUE
      SDF=SQRT(SDI)
      XF=XXF*SDF
      SDF=XSDF*SDF
   50 CONTINUE
      RETURN
C
C     LOCATE XI IN THE TABLES.
C
   60 CONTINUE
      H=H+4.0
      IH=H*10.0
      A=(H-IH*0.1)*10.0
      B=1.0-A
      IH=IH+1
      IHP1=IH+1
C
C     INTERPOLATE IN THE TABLES.
C
      IF (IOP) 70,65,65
   65 CONTINUE
      XJ=(A*ZJ(IHP1)+B*ZJ(IH))*SDI
      SDJ=(A*ZJSD(IHP1)+B*ZJSD(IH))*SDI
   70 CONTINUE
      IF (IOP) 75,75,80
   75 CONTINUE
      SDF=SQRT(SDI)
      XF=(A*ZF(IHP1)+B*ZF(IH))*SDF
      SDF=(A*ZFSD(IHP1)+B*ZFSD(IH))*SDF
   80 CONTINUE
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE NORMAL (NDATA,IO2,IO1,ILP,ATIME,ADATE,TITLE,SMIN,SMAX,
     & GINV,DATA,INDEX,NLOCAL)
C
C LOCALLY NORMALIZED STRUCTURE FACTOR MAGNITUDES
C
C E(HKL) = F(HKL)/SQRT[EPSILON(HKL)*<FSQ/EPSILON>]
C
      CHARACTER ATIME*8,ADATE*9,TITLE*80,FILEO*80
      DIMENSION GINV(3,3),DATA(1),INDEX(1)
C
C SORT REFLECTIONS ON INCREASING SIN(THETA)/LAMBDA.
C
      DO I=1,NDATA
        READ (IO2,REC=I) (IDUMMY,J=1,3),(DUMMY,J=1,6),(IDUMMY,J=1,3),S
        DATA(I)=S
      END DO
      CALL SORT (NDATA,DATA,INDEX)
C
C EVALUATE LOCAL AVERAGE INTENSITIES IN SHELLS OF SIN(THETA)/LAMBDA
C CONTAINING EQUAL NUMBERS OF REFLECTIONS.
C
      IF (NLOCAL.EQ.0) THEN
        NLOCAL=MAX(NDATA/100,10)
        NLOCAL=MIN(NLOCAL,500,NDATA)
      END IF
      NHALF=NLOCAL/2
      SUMM=0
      SUM1=0
      SUM2=0
      DO I=1,NHALF
        READ (IO2,REC=INDEX(I)) (IDUMMY,J=1,3),FSQ,(DUMMY,J=1,5),
     &  IDUMMY,M,IE
        IF (FSQ.LT.0) FSQ=0
        FSQ=FSQ/IE
        SUMM=SUMM+M
        SUM1=SUM1+M*FSQ
        SUM2=SUM2+M*FSQ**2
      END DO
      DO I=1,NDATA
        I1=I-NHALF
        I2=I+NHALF
        IF (I1.LT.1) THEN
          I1=1
          M1=0
          Y1=0
        ELSE
          READ (IO2,REC=INDEX(I1)) (IDUMMY,J=1,3),FSQ,(DUMMY,J=1,5),
     &    IDUMMY,M,IE
          IF (FSQ.LT.0) FSQ=0
          M1=M
          Y1=FSQ/IE
        END IF
        IF (I2.GT.NDATA) THEN
          I2=NDATA
          M2=0
          Y2=0
        ELSE
          READ (IO2,REC=INDEX(I2)) (IDUMMY,J=1,3),FSQ,(DUMMY,J=1,5),
     &    IDUMMY,M,IE
          IF (FSQ.LT.0) FSQ=0
          M2=M
          Y2=FSQ/IE
        END IF
        SUMM=SUMM-M1+M2
        SUM1=SUM1-M1*Y1+M2*Y2
        SUM2=SUM2-M1*Y1**2+M2*Y2**2
        N=I2-I1+1
        AVG=SUM1/SUMM
        AVGSQ=SUM2/SUMM
        RMSD=SQRT((FLOAT(N)/(N-1))*(AVGSQ-AVG**2))
        AVGFSQ=AVG
        SIGAVG=RMSD/SQRT(FLOAT(N))
        READ  (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,F,SIGF,DUMMY,DUMMY,
     &  IC,M,IE,S
        WRITE (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,F,SIGF,  0.0,  0.0,
     &  IC,M,IE,S,AVGFSQ,SIGAVG
      END DO
C
C EVALUATE LOCALLY NORMALIZED STRUCTURE FACTOR MAGNITUDES, AND FIND RE-
C SCALING FACTOR TO FORCE <E**2> = 1.
C
      SUMM=0
      SUM2=0
      DO I=1,NDATA
        READ (IO2,REC=I) J,K,L,FSQ,SIGFSQ,F,SIGF,DUMMY,DUMMY,
     &  IC,M,IE,S,AVGFSQ
        IF (F.GT.0) THEN
          E=F/SQRT(IE*AVGFSQ)
          SUMM=SUMM+M
          SUM2=SUM2+M*E**2
        END IF
      END DO
      AVESQ=SUM2/SUMM
      RESCL=1/SQRT(AVESQ)
      WRITE (ILP,6000) ATIME,ADATE,TITLE
      WRITE (ILP,6100)
      WRITE (ILP,6200) AVESQ,RESCL
      WRITE (ILP,6001)
 6000 FORMAT (/1H1,130('-')/1X,
     &'SUBPROGRAM NORMAL/BAYES.  ',A,', ',A,'.  ',A)
 6100 FORMAT (/1X,
     &'NORMALIZED AND RE-SCALED STRUCTURE FACTORS E(HKL) = F(HKL)/SQRT',
     &'[EPSILON(HKL)*<FSQ/EPSILON>)')
 6200 FORMAT (/1X,
     &'MULTIPLICITY-WEIGHTED AVERAGE ESQ BEFORE RE-SCALING:  <E**2> = ',
     &'                 ',F6.3/1X,
     &'E-MAGNITUDE RE-SCALING FACTOR TO FORCE <E**2> = 1:    RESCl  = ',
     &'1/SQRT(<E**2>) = ',F6.3)
 6001 FORMAT (/1X,
     &'   H   K   L SIN(TH)/L    D(HKL)         FSQ      SIGFSQ       ',
     &'  F      SIGF         E      SIGE'/1X,
     &'   -   -   - ---------    ------         ---      ------       ',
     &'  -      ----         -      ----')
 6002 FORMAT (1X,3I4,F10.4,F10.2,2F12.2,2F10.2,2F10.3)
C
C STORE RE-SCALED E-MAGNITUDES IN THE WORKING DATA FILE.
C
      DO I=1,NDATA
        READ (IO2,REC=I) J,K,L,FSQ,SIGFSQ,F,SIGF,DUMMY,DUMMY,
     &  IC,M,IE,S,AVGFSQ,SIGAVG
        IF (F.GT.0) THEN
          E=F/SQRT(IE*AVGFSQ)
          SIGE=E*SQRT((SIGF**2/(IE*F**2))+(0.5*SIGAVG/AVGFSQ)**2)
          E=RESCL*E
          SIGE=RESCL*SIGE
          WRITE (IO2,REC=I) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE,
     &    IC,M,IE,S,AVGFSQ,SIGAVG
          IF (IPRINT(J,K,L).EQ.1) THEN
            D=1/(2*S)
            WRITE (ILP,6002) J,K,L,S,D,FSQ,SIGFSQ,F,SIGF,E,SIGE
          END IF
        END IF
      END DO
C
C COMPILE AND PRINT E-STATISTICS.
C
      CALL ESTATS (NDATA,IO2,ILP,ATIME,ADATE,TITLE,SMIN,SMAX)
C
C SORT ON d(hkl).
C
      DO I=1,NDATA
        READ (IO2,REC=I) J,K,L
        DATA(I)=SINTHL(J,K,L,GINV)
      END DO
      CALL SORT (NDATA,DATA,INDEX)
C
C LIST LOWEST RESOLUTION E-MAGNITUDES.
C
      WRITE (ILP,6000) ATIME,ADATE,TITLE
      WRITE (ILP,'(/1X,''LOWEST RESOLUTION E-MAGNITUDES:'')')
      WRITE (ILP,6005)
      DO I=1,MIN(100,NDATA)
        READ (IO2,REC=INDEX(I)) J,K,L,    FSQ,SIGFSQ,F,SIGF,E,SIGE
        S=SINTHL(J,K,L,GINV)
        D=1/(2*S)
        WRITE (ILP,6006)        J,K,L,S,D,FSQ,SIGFSQ,F,SIGF,E,SIGE,I
      END DO
C
C WRITE E-SORTED OUTPUT FILE.
C
      DO I=1,NDATA
        READ  (IO2,REC=I) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE
        DATA(I)=E
      END DO
      CALL SORT (NDATA,DATA,INDEX)
      OPEN (UNIT=IO1,FILE='edata.bayes',STATUS='NEW',FORM='FORMATTED')
      DO I=NDATA,1,-1
        READ (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE
        IF (E.GT.0) CALL WRITE1 (IO1,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE)
      END DO
      CLOSE (UNIT=IO1,STATUS='KEEP')
C
C LIST LARGEST AND SMALLEST SIGNIFICANT E-MAGNITUDES.
C
      WRITE (ILP,6000) ATIME,ADATE,TITLE
      WRITE (ILP,6051)
      WRITE (ILP,6005)
      I=NDATA
      N=0
      DO WHILE (I.GE.1.AND.N.LT.100)
        READ  (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE,
     &  IC,M,IE,S,AVGFSQ
        IF (E.GT.3*SIGE) THEN
          D=1/(2*S)
          WRITE (ILP,6006) J,K,L,S,D,FSQ,SIGFSQ,F,SIGF,E,SIGE,NDATA-I+1
          N=N+1
        END IF
        I=I-1
      END DO
      WRITE (ILP,6000) ATIME,ADATE,TITLE
      WRITE (ILP,6052)
      WRITE (ILP,6005)
      I=1
      N=0
      DO WHILE (I.LE.NDATA.AND.N.LT.100)
        READ  (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE,
     &  IC,M,IE,S,AVGFSQ
        IF (E.GT.3*SIGE) THEN
          D=1/(2*S)
          WRITE (ILP,6006) J,K,L,S,D,FSQ,SIGFSQ,F,SIGF,E,SIGE,-I
          N=N+1
        END IF
        I=I+1
      END DO
      WRITE (ILP,6009)
 6051 FORMAT (/1X,'LARGEST E-MAGNITUDES WITH  E .GT. 3*SIGMA(E)')
 6052 FORMAT (/1X,'SMALLEST E-MAGNITUDES WITH  E .GT. 3*SIGMA(E)')
 6005 FORMAT (/1X,
     &'   H   K   L SIN(TH)/L    D(HKL)         FSQ      SIGFSQ       ',
     &'  F      SIGF         E      SIGE      RANK'/1X,
     &'   -   -   - ---------    ------         ---      ------       ',
     &'  -      ----         -      ----      ----')
 6006 FORMAT (1X,3I4,F10.4,F10.2,2F12.2,2F10.2,2F10.3,I10)
 6009 FORMAT (/1X,
     &'OUTPUT FILE OF E-VALUES ESTIMATED AS'/1X,'  E(HKL) = F(HKL)/SQR',
     &'T[EPSILON(HKL)*<FSQ(MEAS)/EPSILON>]'/1X,
     &'IS AN E-SORTED, FORMATTED, ASCII FILE (IH,IK,IL,FSQ,SIGFSQ,F,SI',
     &'GF,E,SIGE)'/1X,
     &'NAMED:  "edata.bayes"')
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE NZPLOT (NDATA,IFILE,ILP,ATIME,ADATE,TITLE,PTGP,LAUE)
C
C COMPILE STATISTICS ON AND PLOT THE CUMULATIVE WISLON DISTRIBUTIONS.
C
      PARAMETER (NX=100,NY=50)
      CHARACTER AA*1
      DIMENSION AA(0:NX,0:NY)
      CHARACTER ATIME*8,ADATE*9,TITLE*80,PTGP*5,LAUE*5
      PARAMETER (NZ=10)
      DIMENSION ZI(NZ),FA(NZ),FC(NZ)
      WRITE (ILP,'(/1H1,130(''-'')/1X,
     &''PROGRAM BAYES.  '',A,'', '',A,''.  '',A)') ATIME,ADATE,TITLE
      WRITE (ILP,'(/1X,
     &''TESTS OF THE VALIDITY OF THE WILSON INTENSITY DISTRIBUTIONS''//
     &1X,
     &''POINT GROUP  '',A/1X,''LAUE GROUP   '',A)') PTGP,LAUE
      SUMMA=0
      SUMFA=0
      SUMZA=0
      SUMASQ=0
      SUMMC=0
      SUMFC=0
      SUMZC=0
      SUMCSQ=0
      DO I=1,NZ
        ZI(I)=FLOAT(I)/NZ
        FA(I)=0
        FC(I)=0
      END DO
      DO I=1,NDATA
        READ (IFILE,REC=I) J,K,L,FSQ,SIGFSQ,DUMMY,DUMMY,DUMMY,DUMMY,
     &  IC,M,IE,S,AVGFSQ
        IF (FSQ.LT.-4*SIGFSQ) FSQ=-4*SIGFSQ
        Z=FSQ/(IE*AVGFSQ)
        F=SQRT(MAX(0.0,Z))
        IF (IC.EQ.0) THEN
C
C ACENTRIC REFLECTION
C
          SUMMA=SUMMA+M
          SUMFA=SUMFA+M*F
          SUMZA=SUMZA+M*Z
          SUMASQ=SUMASQ+M*Z**2
          DO J=1,NZ
            IF (Z.LE.ZI(J)) FA(J)=FA(J)+M 
          END DO
        ELSE
C
C CENTRIC REFLECTION
C
          SUMMC=SUMMC+M
          SUMFC=SUMFC+M*F
          SUMZC=SUMZC+M*Z
          SUMCSQ=SUMCSQ+M*Z**2
          DO J=1,NZ
            IF (Z.LE.ZI(J)) FC(J)=FC(J)+M 
          END DO
        END IF
      END DO
      IF (SUMMA.EQ.0) THEN
        WRITE (ILP,'(/1X,
     &''NO ACENTRIC REFLECTIONS IN THE GIVEN POINT GROUP'')')
      ELSE
        AVGZA=SUMZA/SUMMA
        VARZA=SUMASQ/SUMMA-(SUMZA/SUMMA)**2
        RHOA=(SUMFA/SUMMA)**2/AVGZA
        DO I=1,NZ
          FA(I)=FA(I)/SUMMA
        END DO
        WRITE (ILP,'(//1X,
     &''ACENTRIC WILSON DISTRIBUTION STATISTICS ON Z = FSQ/(EPSILON*<'',
     &''FSQ/EPSILON>''/1X,
     &''--------''/1X,
     &''                           EXPERIMENTAL  THEORETICAL''/1X,
     &''                           ------------  -----------''/1X,
     &''MULT.-WTD. NO. OF REFLNS.  '',I9/1X,
     &''<Z>                        '',F9.3,5X,'' 1.000''/1X,
     &''<(Z - <Z>)**2>             '',F9.3,5X,'' 1.000''/1X,
     &''RHO = <SQRT(Z)>**2/<Z>     '',F9.3,5X,'' 0.785 = PI/4'')')
     & NINT(SUMMA),AVGZA,VARZA,RHOA
      END IF
      IF (SUMMC.EQ.0) THEN
        WRITE (ILP,'(/1X,
     &''NO CENTRIC REFLECTIONS IN THE GIVEN POINT GROUP'')')
      ELSE
        AVGZC=SUMZC/SUMMC
        VARZC=SUMCSQ/SUMMC-(SUMZC/SUMMC)**2
        RHOC=(SUMFC/SUMMC)**2/AVGZC
        DO I=1,NZ
          FC(I)=FC(I)/SUMMC
        END DO
        WRITE (ILP,'(//1X,
     &''CENTRIC WILSON DISTRIBUTION STATISTICS ON Z = FSQ/(EPSILON*<F'',
     &''SQ/EPSILON>''/1X,
     &''-------''/1X,
     &''                           EXPERIMENTAL  THEORETICAL''/1X,
     &''                           ------------  -----------''/1X,
     &''MULT.-WTD. NO. OF REFLNS.  '',I9/1X,
     &''<Z>                        '',F9.3,5X,'' 1.000''/1X,
     &''<(Z - <Z>)**2>             '',F9.3,5X,'' 2.000''/1X,
     &''RHO = <SQRT(Z)>**2/<Z>     '',F9.3,5X,'' 0.637 = 2/PI'')')
     & NINT(SUMMC),AVGZC,VARZC,RHOC
      END IF
C
C ERASE THE PLOT ARRAY.
C
      DO  IX=0,NX
      DO  IY=0,NY
        AA(IX,IY)=' '
      END DO
      END DO
C
C PLOT 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 PLOT THE THEORETICAL ACENTRIC AND CENTRIC CUMULATIVE DISTRIBUTION
C CURVES, N(Z) VS. Z FOR 0 .LE. Z .LE. 1.
C
      XMIN=0
      XMAX=1
      YMIN=0
      YMAX=1
      XDEL=XMAX-XMIN
      YDEL=YMAX-YMIN
      DX=(XMAX-XMIN)/(2*NX)
      DO I=0,(2*NX)
        XI=-XMIN+I*DX
        IX=NINT(NX*(XI-XMIN)/XDEL)
        IX=MAX(IX,0)
        IX=MIN(IX,NX)
C
C ACENTRIC DISTRIBUTION N(Z) = 1 - EXP(-Z)
C
        YI=1-EXP(-XI)
        IY=NINT(NY*(YMAX-YI)/YDEL)
        IY=MAX(IY,0)
        IY=MIN(IY,NY)
        AA(IX,IY)='0'
C
C CENTRIC DISTRIBUTION N(Z) = ERF[SQRT(Z/2)]
C
        YI=ERRF(SQRT(0.5*XI))
        IY=NINT(NY*(YMAX-YI)/YDEL)
        IY=MAX(IY,0)
        IY=MIN(IY,NY)
        AA(IX,IY)='1'
      END DO
C
C PLOT THE EXPERIMENTAL DATA POINTS.
C
      DO I=1,NZ
        XI=ZI(I)
        IX=NINT(NX*(XI-XMIN)/XDEL)
        IX=MAX(IX,0)
        IX=MIN(IX,NX)
C
C ACENTRIC POINTS
C
        IF (SUMMA.GT.0) THEN
          YI=FA(I)
          IY=NINT(NY*(YMAX-YI)/YDEL)
          IY=MAX(IY,0)
          IY=MIN(IY,NY)
          AA(IX,IY)='A'
        END IF
C
C CENTRIC POINTS
C
        IF (SUMMC.GT.0) THEN
          YI=FC(I)
          IY=NINT(NY*(YMAX-YI)/YDEL)
          IY=MAX(IY,0)
          IY=MIN(IY,NY)
          AA(IX,IY)='C'
        END IF
      END DO
C
C PRINT THE PLOT ARRAY.
C
      WRITE (ILP,'(/1H1,130(''-'')/1X,10X,
     &''PROGRAM BAYES.  '',A,'', '',A,''.  '',A/)') ATIME,ADATE,TITLE
      DY=(YMAX-YMIN)/NY
      DO IY=0,NY
        IF (MOD(IY,10).EQ.0) THEN
          WRITE (ILP,'(1X,F10.1,101A1)') YMAX-IY*DY,(AA(IX,IY),IX=0,NX)
        ELSE
          WRITE (ILP,'(1X,10X,101A1)')              (AA(IX,IY),IX=0,NX)
        END IF
      END DO
      DX=(XMAX-XMIN)/10
      WRITE (ILP,'(2X,11F10.1)') (XMIN+IX*DX,IX=0,10)
      WRITE (ILP,'(/1X,10X,
     &''CUMULATIVE WILSON DISTRIBUTIONS N(Z) VS. Z = FSQ/(EPSILON*<FS'',
     &''Q/EPSILON>) FOR 0 .LE. Z .LE. 1''//1X,10X,
     &''ACENTRIC REFLECTIONS:  THEORETICAL    N(Z) = 1 - EXP(-Z)     '',
     &''  "0",    EXPERIMENTAL    "A"''/1X,10X,
     &''CENTRIC       "     :       "         N(Z) = ERF[SQRT(Z/2)]  '',
     &''  "1",         "          "C"'')')
      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 PLOT (ILP,IFILE,NDATA,XMIN,XMAX,YMIN,YMAX,IFLAG)
      PARAMETER (NX=100,NY=50)
      CHARACTER AA*1
      DIMENSION AA(0:NX,0:NY)
C
C ERASE THE PLOT ARRAY.
C
      DO  IX=0,NX
      DO  IY=0,NY
        AA(IX,IY)=' '
      END DO
      END DO
C
C PLOT 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 PLOT DATA POINTS.
C
      XDEL=XMAX-XMIN
      YDEL=YMAX-YMIN
      DO I=1,NDATA
        READ (IFILE,REC=I) (IDUMMY,J=1,3),(DUMMY,J=1,6),(IDUMMY,J=1,3),
     &   S,YI,ZI
        XI=S
        IF (IFLAG.NE.0) YI=ZI
        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 PRINT THE PLOT ARRAY.
C
      DY=(YMAX-YMIN)/NY
      DO IY=0,NY
        IF (MOD(IY,10).EQ.0) THEN
          WRITE (ILP,'(1X,E10.3,101A1)') YMAX-IY*DY,(AA(IX,IY),IX=0,NX)
        ELSE
          WRITE (ILP,'(1X,10X,  101A1)')            (AA(IX,IY),IX=0,NX)
        END IF
      END DO
      DX=(XMAX-XMIN)/10
      WRITE (ILP,'(2X,11F10.3)') (2*(XMIN+IX*DX),IX=0,10)
      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 ERASE THE PLOT ARRAY.
C
      DO  IX=0,NX
      DO  IY=0,NY
        AA(IX,IY)=' '
      END DO
      END DO
C
C PLOT 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 PLOT 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 PLOT 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,'(1X,F10.1,101A1)') YMAX-IY*DY,(AA(IX,IY),IX=0,NX)
        ELSE
          WRITE (ILP,'(1X,10X,  101A1)')            (AA(IX,IY),IX=0,NX)
        END IF
      END DO
      DX=(XMAX-XMIN)/10
      WRITE (ILP,'(2X,11F10.1)') (XMIN+IX*DX,IX=0,10)
      WRITE (ILP,'(/1X,10X,
     &''ACENTRIC DISTRIBUTION:  P(E) = 2*E*EXP(-E**2),  "*" THEORETIC'',
     &''AL,  "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 ERASE THE PLOT ARRAY.
C
      DO  IX=0,NX
      DO  IY=0,NY
        AA(IX,IY)=' '
      END DO
      END DO
C
C PLOT 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 PLOT 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 PLOT 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,'(1X,F10.1,101A1)') YMAX-IY*DY,(AA(IX,IY),IX=0,NX)
        ELSE
          WRITE (ILP,'(1X,10X,  101A1)')            (AA(IX,IY),IX=0,NX)
        END IF
      END DO
      DX=(XMAX-XMIN)/10
      WRITE (ILP,'(2X,11F10.1)') (XMIN+IX*DX,IX=0,10)
      WRITE (ILP,'(/1X,10X,
     &''CENTRIC DISTRIBUTION:  P(E) = SQRT(2/PI)*EXP(-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 ERASE THE PLOT ARRAY.
C
      DO  IX=0,NX
      DO  IY=0,NY
        AA(IX,IY)=' '
      END DO
      END DO
C
C PLOT 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(TH)/L < 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 PLOT <E**2> = 1.
C
      DO IX=0,NX
        AA(IX,30)='.'
      END DO
C
C PLOT 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 PLOT 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,'(1X,F10.1,101A1)') YMAX-IY*DY,(AA(IX,IY),IX=0,NX)
        ELSE
          WRITE (ILP,'(1X,10X,  101A1)')            (AA(IX,IY),IX=0,NX)
        END IF
      END DO
      DX=(XMAX-XMIN)/10
      WRITE (ILP,'(2X,11F10.3)') (2*(XMIN+IX*DX),IX=0,10)
      WRITE (ILP,'(/1X,10X,
     &''<E**2> VERSUS <DSTAR> = 2*<SIN(THETA)/LAMBDA> DISTRIBUTION'')')
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE READ1 (ITYPE,IFILE,IEND,IH,IK,IL,FSQ,SIGFSQ)
 1    IF (ITYPE.EQ.0) READ (IFILE,*,END=9) IH,IK,IL,FSQ,SIGFSQ
      IF (ITYPE.EQ.1) READ (IFILE,*,END=9) IH,IK,IL,F,  SIGF
      IF (ITYPE.EQ.2) READ (IFILE,  END=9) IH,IK,IL,FSQ,SIGFSQ
      IF (ITYPE.EQ.3) READ (IFILE,  END=9) IH,IK,IL,F,  SIGF
      IF (ITYPE.EQ.1.OR.ITYPE.EQ.3) THEN
        IF (F.LT.0) F=0
        IF (SIGF.LT.0) SIGF=0
        FSQ=F**2
        SIGFSQ=MAX(2*F*SIGF,4*SIGF**2)
      END IF
      IF (SIGFSQ.LE.0) GO TO 1
      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
      Q=0
      DO I=1,3
      DO J=1,3
        Q=Q+H(J)*GINV(I,J)*H(I)
      END DO
      END DO
      SINTHL=0.5*SQRT(Q)
      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 ADAPTED 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. 229-233.  CAMBRIDGE, ENGLAND:
C CAMBRIDGE UNIVERSITY PRESS.
C
      DIMENSION DATA(N),INDEX(N)
      DO I=1,N
        INDEX(I)=I
      END DO
      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
      DO WHILE (L.LE.J)
        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
      END DO
      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 ADAPTED 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=100)
      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 ADAPTED 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
      DO WHILE (I2-I1.GT.1)
        I=(I2+I1)/2
        IF (X(I).GT.XI) THEN
          I2=I
        ELSE
          I1=I
        END IF
      END DO
      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-----------------------------------------------------------------------
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=====================================================================
      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-----------------------------------------------------------------------
