      PROGRAM LEVY
C
C ROBERT H. BLESSING
C MEDICAL FOUNDATION OF BUFFALO
C 73 HIGH STREET
C BUFFALO, NEW YORK 14203, USA
C TELEPHONE:  (716) 856-9600, EXTENSION 335
C ELECTRONIC MAIL:  Blessing@HWI.Buffalo.Edu
C
C FEBRUARY 1995,..., JANUARY 2001
C
C A PRIORI SCALE FACTOR AND OVERALL ANISOTROPIC MEAN-SQUARE ATOMIC
C DISPLACEMENT PARAMETERS BY AN EXTENSION OF THE METHOD OF
C H.A. LEVY, W.E. THIESSEN, AND (IN PART) G.M. BROWN (1970).  A LEAST-
C SQUARES METHOD FOR CONVERTING OBSERVED INTENSITIES TO NORMALIZED
C STRUCTURE FACTORS.  AM. CRYST. ASSOC. MEETING, TULANE UNIV., NEW
C ORLEANS, MARCH 1970.  ABSRACT NO. B3.
C
C LEAST-SQUARES FIT OF WILSON EXPECTATION VALUES TO MEASURED INTENSITIES
C
C CHISQ = SUM W*DELTA**2 = SUM (DELTA/SIGMA)**2
C DELTA = YOBS - YCALC
C YOBS = FSQ(MEAS)/[EPSILON*SUM F(J)**2]
C YCAL = (SCALEK**-2)*EXP[-2*{H(J)*B(I,J)*H(I) - [H(J)*C(I,J)*H(I)]**2}]
C W = M/[SIGMA(YOBS)**2 + YCALC**2*VAR]
C M = MULTIPLICITY OF REFLECTION
C VAR = 1 OR 2 FOR, RESPECTIVELY, ACENTRIC OR CENTRIC WILSON
C DISTRIBUTIONS OF NORMALIZED INTENSITIES.
C
C YOBS IS A MEASURED SQUARED STRUCTURE FACTOR MAGNITUDE NORMALIZED TO
C THE ABSOLUTE SCALE OF SCATTERING BY ATOMS AT REST.
C
C YCALC IS A WILSON EXPECTATION VALUE ON THE RELATIVE EXPERIMENTAL
C SCALE FOR VIBRATING OR DISORDERED ATOMS.
C
C SCALEK AND B(I,J) ARE DEFINED ACCORDING TO:
C
C F(ABSOLUTE) = SCALEK*F(RELATIVE),
C
C f(T) = f(0)*EXP[-SUM(I) SUM(J) H(I)*H(J)*B(I,J)].
C
C THE C(I,J) ARE DISPERSION PARAMETERS OF THE DISTRIBUTION OVER THE
C UNIT CELL OF THE ANISOTROPIC ATOMIC DISPLACEMENT PARAMETERS B(I,J).
C
C IN THE ISOTROPIC APPROXIMATION,
C
C YCALC = (SCALEK**-2)*EXP{-2*[BISO - (CISO*S)**2]*S**2},
C
C WHERE S = SIN(THETA)/LAMBDA, AND BISO AND CISO ARE, RESPECTIVELY, THE
C ESTIMATED MEAN AND STANDARD DEVIATION OF THE B-DISTRIBUTION.
C
      CHARACTER ATIME*8,ADATE*9,TITLE*80,FILE*80,XRAY*2,SYST*12,LATT*1,
     & PTGP*5,LAUE*5,EL*2,ELANOM*2
      INTEGER ZCELL
      DIMENSION EL(20),XI(20),ZI(20),RM(20),A1(20),B1(20),A2(20),B2(20),
     & A3(20),B3(20),A4(20),B4(20),C0(20),FP(20),FPP(20)
      DIMENSION CELL(6),GA(3,3),GB(3,3),B(3,3),C(3,3),U(3,3),V(3,3),
     & COV(13,13),BIBJ(6),VAR(13),TEMP(109)
      DATA B,C,COV,TEMP /9*0,9*0,169*0,109*0/
C
C AT MOST NMAX/2 REFLECTIONS CAN BE SORTED.
C
      PARAMETER (NMAX=1000000)
      DIMENSION DATA(NMAX),INDEX(NMAX/2)
      EQUIVALENCE (INDEX(1),DATA(1+NMAX/2))
      CALL TIME (ATIME)
      CALL DATE (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='levy.dat',STATUS='OLD')
      READ (IO1,*) TITLE
      READ (IO1,*) ITYPE
      READ (IO1,*) FILE
      READ (IO1,*) LATT
      READ (IO1,*) PTGP
      READ (IO1,*) CELL
      CALL DECODE (PTGP,LAUE,SYST,LATT,MLATT,IPTGP,ILAUE)
      CALL METRIC (CELL,VCELL,GA,GB)
      ZCELL=0
      XRAY='  '
      NEL=0
      READ (IO1,*,END=4) ZCELL
      READ (IO1,*,END=4) XRAY
      READ (IO1,*,END=4) NEL
 4    CONTINUE
      IF (ZCELL.EQ.0)   ZCELL=1
      IF (XRAY.EQ.'  ') XRAY='X '
      IF (NEL.GT.0) THEN
        DO I=1,NEL
          READ (IO1,*) EL(I),XI(I)
        END DO
        ICELL=1
      ELSE
        ICELL=0
C
C DEFAULT EMPIRICAL FORMULA FOR AN AVERAGE PROTEIN, C4 H6 N O.
C
        EL(1)='C '
        EL(2)='H '
        EL(3)='N '
        EL(4)='O '
        XI(1)=4
        XI(2)=6
        XI(3)=1
        XI(4)=1
        NEL=4
      END IF
      F000=0
      RMCELL=0
      DO I=1,NEL
        CALL FTABLE (EL(I),ZI(I),RM(I),A1(I),B1(I),A2(I),B2(I),A3(I),
     &  B3(I),A4(I),B4(I),C0(I),XRAY,FP(I),FPP(I))
        XI(I)=ZCELL*XI(I)
        F000=F000+XI(I)*ZI(I)
        RMCELL=RMCELL+XI(I)*RM(I)
      END DO
      DX=(RMCELL/VCELL)/0.60221
C
C READ ANOMALOUS DISPERSION CORRECTIONS FOR WAVELENGTHS OTHER THAN CU OR
C MO K-ALPHA.
C
      NANOM=0
      READ (IO1,*,END=5) NANOM
      IF (NANOM.GT.0) THEN
        DO I=1,NEL
          FP(I)=0
          FPP(I)=0
        END DO
        DO I=1,NANOM
          READ (IO1,*) ELANOM,FPI,FPPI
          DO J=1,NEL
            IF (EL(J).EQ.ELANOM) THEN
              FP(J)=FPI
              FPP(J)=FPPI
            END IF
          END DO
        END DO
      END IF
 5    CONTINUE
      SMIN=0
      SMAX=0
      ISOLV=0
      ASOLV=0
      BSOLV=0
      READ (IO1,*,END=6) SMIN,SMAX
      READ (IO1,*,END=6) ISOLV,ASOLV,BSOLV
 6    CONTINUE
      CLOSE (UNIT=IO1,STATUS='KEEP')
      IF (ICELL.EQ.0) THEN
C
C ALONG WITH THE DEFAULT EMPIRICAL FORMULA, C4 H6 N O, FOR AN AVERAGE
C PROTEIN, ASSUME PROTEIN DENSITY 1.35 MG MM**-3 AND 50% SOLVENT VOLUME
C FRACTION.
C
        RMTEMP=1.35*0.5*VCELL*0.60221
        T=RMTEMP/RMCELL
        RMCELL=RMTEMP
        F000=F000*NINT(T)
        DO I=1,4
          XI(I)=XI(I)*NINT(T)
        END DO
      END IF
C
C CALCULATE MATTHEWS COEFFICIENT, SOLVENT VOLUME FRACTION, AND NUMBER OF
C WATER MOLECULES AT 29.9 A**3 PER H2O MOLECULE TO FILL THE SOLVENT
c VOLUME IN PROTEIN CRYSTALS.
C
      VM=VCELL/RMCELL
      VS=1-1.23*RMCELL/VCELL
      NWATER=NINT(VS*VCELL/29.9)
      IF (VS.GT.0.1.AND.NWATER.GT.0) THEN
C
C CALCULATE BULK-SOLVENT-CORRECTED UNIT-CELL SCATTERING POWER AND MASS
C DENSITY.
C
        F000=F000+NWATER*10
        DX=((RMCELL+NWATER*18.0153)/VCELL)/0.60221
      END IF
C
C ECHO CONTROL DATA.
C
      OPEN (UNIT=ILP,FILE='levy.lp',STATUS='NEW')
      WRITE (ILP,600) ATIME,ADATE,TITLE
      WRITE (ILP,610) ITYPE
      WRITE (ILP,601) FILE
      WRITE (ILP,602) SYST,LATT,PTGP,LAUE
      IF (PTGP.EQ.LAUE) THEN
        WRITE (ILP,604)
      ELSE
        WRITE (ILP,603)
      END IF
      WRITE (ILP,606) CELL,VCELL
      IF (ICELL.EQ.0) WRITE (ILP,670)
      WRITE (ILP,607) ZCELL,(EL(I),XI(I)/ZCELL,I=1,NEL)
      WRITE (ILP,608) NINT(F000),RMCELL/ZCELL,DX,VM,100*VS,NWATER/ZCELL
      WRITE (ILP,609) XRAY,(EL(I),FP(I),FPP(I),I=1,NEL)
      IF (SMIN.NE.0) WRITE (ILP,611) SMIN,1/(2*SMIN)
      IF (SMAX.NE.0) WRITE (ILP,612) SMAX,1/(2*SMAX)
      IF (ISOLV.NE.0) THEN
        IF (ASOLV.EQ.0) ASOLV=0.75
        IF (BSOLV.EQ.0) BSOLV=250
        WRITE (ILP,613) ASOLV,BSOLV
      END IF
 600  FORMAT (/1H1,130('-')/1X,'PROGRAM LEVY.  ',A,', ',A,'.  ',A)
 610  FORMAT (/1X,'INPUT REFLECTION FILE TYPE:  ITYPE = ',I2)
 601  FORMAT (/1X,'INPUT REFLECTION FILE NAME:  ',A)
 602  FORMAT (/1X,A//1X,A,'-LATTICE'//1X,
     &'CRYSTAL CLASS POINT GROUP  ',A//1X,
     &'LAUE GROUP                 ',A)
 603  FORMAT (/1X,'NONCENTROSYMMETRIC STRUCTURE')
 604  FORMAT (/1X,'CENTROSYMMETRIC STRUCTURE')
 606  FORMAT (/1X,9X,'A',9X,'B',9X,'C',5X,'ALPHA',6X,'BETA',5X,'GAMMA',
     &14X,'V'/1X,3F10.4,3F10.3,F15.2)
 670  FORMAT (/1X,
     &'DEFAULT CRYSTAL CHEMICAL COMPOISITION FOR AN AVERAGE PROTEIN:'/
     &1X,
     &'-------------------------------------------------------------'/
     &1X,
     &'  EMPIRICAL FORMULA          C4 H6 N O'/1X,
     &'  SPECIFIC PARTIAL VOLUME    0.74 MM**3 MG-1'/1X,
     &'  MOLECULAR MASS DENSITY     1.35 MG MM**-3'/1X,
     &'  SOLVENT VOLUME FRACTION    50% H2O')
 607  FORMAT (/1X,
     &'UNIT CELL CONTENTS:'/1X,
     &'-------------------'//1X,
     &'  Z = ',I2,' CRYSTAL CHEMICAL UNITS PER UNIT CELL'//1X,
     &'  CRYSTAL CHEMICAL UNIT:'//('   ',A2,F10.2))
 608  FORMAT (/1X,
     &'F000     = ',I10,' E         (ELECTRONS)'/1X,
     &'MOL. WT. = ',F10.2,' DA        (DALTONS)'/1X,
     &'D(X-RAY) = ',E10.3,' MG MM**-3 (MILLIGRAMS PER CUBIC MILLIMETER)'
     &//1X,
     &'MATTHEWS COEFFICIENT AND SOLVENT VOLUME FRACTION:'/1X,
     &'-------------------------------------------------'/1X,
     &'  VM     = ',E10.3,' A**3 DA**-1 (CUBIC ANGSTROMS PER DALTON)'/
     &1X,
     &'  VS     = ',E10.3,' %'/1X,
     &'  N(H2O) = ',I10,' WATER MOLECULES PER CRYSTAL CHEMICAL UNIT NE',
     &'EDED TO FILL THE SOLVENT VOLUME')
 609  FORMAT (/1X,
     &'RADIATION:  ',A2//1X,
     &'ANOMALOUS DISPERSION CORRECTION TERMS:'/1X,
     &'--------------------------------------'/1X,
     &'EL        FP       FPP'/1X,
     &'--        --       ---'/
     &(1X,A2,2F10.3))
 611  FORMAT (/1X,
     &'MINIMUM SIN(THETA)/LAMBDA FOR INCLUDING REFLECTION IN LEAST-SQU',
     &'ARES FITTING'/1X,
     &'  SMIN = ',F7.4,' RECIPROCAL ANGSTROMS'/1X,
     &'  DMAX = ',F5.2,'   ANGSTROMS')
 612  FORMAT (/1X,
     &'MAXIMUM SIN(THETA)/LAMBDA FOR INCLUDING REFLECTION IN LEAST-SQU',
     &'ARES FITTING'/1X,
     &'  SMAX = ',F7.4,' RECIPROCAL ANGSTROMS'/1X,
     &'  DMIN = ',F5.2,'   ANGSTROMS')
 613  FORMAT (/1X,
     &'BULK SOLVENT COMPENSATION PARAMETERS:'/1X,
     &'-------------------------------------'/1X,
     &'  F(PROT) = F(MEAS)/(1 - ASOLV*EXP{-BSOLV*[SIN(THETA)/LAMBDA]**',
     &'2})'/1X,
     &'  ASOLV = ',E10.3/1X,
     &'  BSOLV = ',E10.3' A**2 (SQUARED ANGSTROMS)')
C
C BUILD WORKING REFLECTION DATA FILES.
C
      CALL DATAIN (FILE,ITYPE,IO1,IO2,IO3,GB,SMIN,SMAX,SLOW,SHIGH,
     & ASOLV,BSOLV,MLATT,IPTGP,ILAUE,NEL,XI,A1,B1,A2,B2,A3,B3,A4,B4,C0,
     & FP,FPP,MDATA,NDATA,DATA)
      
C
C PRINT SOME DATA SET STATISTICS.
C
      WRITE (ILP,600) ATIME,ADATE,TITLE
      WRITE (ILP,614) MDATA,NDATA,SMIN,SMAX,1/(2*SMIN),1/(2*SMAX),
     & SLOW,SHIGH,1/(2*SLOW),1/(2*SHIGH)
 614  FORMAT (/1X,
     &'NUMBER OF REFLECTIONS IN THE INPUT DATA FILE:'/1X,
     &'  MDATA = ',I6/1X,
     &'NUMBER OF REFLECTIONS USED IN THE LEAST-SQUARES FITTING:'/1X,
     &'  NDATA = ',I6//1X,
     &'RESOLUTION LIMITS OF THE INPUT DATA SET:'/1X,
     &'  S(MIN) = SIN(THETA(MIN))/LAMBDA = ',F6.3,' RECIPROCAL ANGSTRO',
     &'MS'/1X,
     &'  S(MAX)                          = ',F6.3/1X,
     &'  D(MIN) = 1/(2*S(MAX)) = ',F5.2,' ANGSTROMS'/1X,
     &'  D(MAX) = 1/(2*S(MIN)) = ',F5.2//1X,
     &'RESOLUTION LIMITS OF THE LEAST-SQUARES DATA SET:'/1X,
     &'  S(MIN) = SIN(THETA(MIN))/LAMBDA = ',F6.3,' RECIPROCAL ANGSTRO',
     &'MS'/1X,
     &'  S(MAX)                          = ',F6.3/1X,
     &'  D(MIN) = 1/(2*S(MAX)) = ',F5.2,' ANGSTROMS'/1X,
     &'  D(MAX) = 1/(2*S(MIN)) = ',F5.2)
C
C GET STARTING SCALEK, BISO, AND CISO.
C
      CALL WILSON (IO1,NDATA,SCALEK,BISO,CISO,COV)
C
C REFINE SCALEK, BISO, AND CISO.
C
      CALL KBCISO (IO1,NDATA,SCALEK,BISO,CISO,NCYCLE,R,Z,NFREE,NPAR,COV)
C
C STORE REFINED SCALEK, BISO, AND CISO.
C
      B(1,1)=BISO
      C(1,1)=CISO
      CALL STORE (NCYCLE,R,Z,NFREE,NPAR,SCALEK,B,C,COV,TEMP)
C
C IF CISO = 0, TRY REFINING AGAIN FROM PERTURBED BISO AND CISO.
C
      IF (CISO.EQ.0) THEN
        CISO=BISO/2
        BISO=2*BISO
        T=R
        NPAR=3
        CALL KBCISO (IO1,NDATA,SCALEK,BISO,CISO,NCYCLE,R,Z,NFREE,NPAR,
     &  COV)
C
C TEST FOR AN (R .LT. T) IMPROVED FIT.
C
        IF ((R-T)/T.LT.-1E-6) THEN
C
C STORE RE-REFINED SCALEK AND BISO AND CISO.
C
          B(1,1)=BISO
          C(1,1)=CISO
          CALL STORE (NCYCLE,R,Z,NFREE,NPAR,SCALEK,B,C,COV,TEMP)
        ELSE
C
C RECALL SCALEK AND BISO WITH CISO = 0.
C
          CALL RECALL (TEMP,NCYCLE,R,Z,NFREE,NPAR,SCALEK,B,C,COV)
          BISO=B(1,1)
          CISO=0
        END IF
      END IF
C
C PRINT RESULTS OF ISOTROPIC FITTING.
C
      WRITE (ILP,600) ATIME,ADATE,TITLE
      WRITE (ILP,616) NCYCLE,R,Z,NFREE,NPAR
      WRITE (ILP,617) SCALEK,BISO,CISO
      DO I=1,3
        VAR(I)=COV(I,I)
        IF (VAR(I).EQ.0) VAR(I)=1
      END DO
      WRITE (ILP,618) (SQRT(COV(I,I)),I=1,3),(I,(COV(I,J)/SQRT(VAR(I)*
     & VAR(J)),J=1,3),I=1,3)
 616  FORMAT (/1X,
     &'SCALE AND ISOTROPIC ATOMIC DISPLACEMENT AND DISTRIBUTION DISPER',
     &'SION'/1X,
     &'PARAMETERS:'//1X,
     &'  YOBS  = FSQ(MEAS)/[EPSILON*SUM F(J)**2]'/1X,
     &'  YCALC = (SCALEK**-2)*EXP{-2*[BISO - (CISO*S)**2]*S**2)}'//1X,
     &'  WHERE'/1X,
     &'    EXP{-2*[BISO - (CISO*S)**2]*S**2} = <EXP(-2*B*S**2)>'/1X,
     &'  IS THE STATISTICAL EXPECTATION VALUE OF THE DEBYE-WALLER FACT',
     &'OR FOR'/1X,
     &'  AN ASSUMED-TO-BE-NORMAL DISTRIBUTION OF ATOMIC B-VALUES WITH ',
     &'MEAN'/1X,
     &'    BISO = <B>'/1X,
     &'  AND STANDARD DEVIATION'/1X,
     &'    CISO = <(B - <B>)**2>**(1/2).'//1X,
     &'STATISTICS OF FIT:'//1X,
     &'NCYCLE = ',I11/1X,
     &'R      = ',E11.4/1X,
     &'Z      = ',E11.4/1X,
     &'NFREE  = ',I11/1X,
     &'NPAR   = ',I11)
 617  FORMAT (/1X,
     &'RESULTS:'//1X,'      SCALEK        BISO        CISO'/1X,3E12.4)
 618  FORMAT (/1X,
     &'PARAMETER E.S.D.''S'//1X,
     &'    SIGMA(K)    SIGMA(B)    SIGMA(C)'/1X,3E12.2//1X,
     &'CORRELATION MATRIX'//1X,
     &'       K  BISO  CISO'//1X,
     &'       1     2     3'//(1X,I2,3F6.2))
C
C PLOT RESULTS OF ISOTROPIC FITTING.
C
      CALL WPLOT (IO3,ILP,ATIME,ADATE,TITLE,SCALEK,BISO,CISO,MDATA,DATA,
     & INDEX)
C
C STORE SCALEK AND ISOTROPIC B AND C.
C
      B(1,1)=BISO
      C(1,1)=CISO
      CALL STORE (NCYCLE,R,Z,NFREE,NPAR,SCALEK,B,C,COV,TEMP)
C
C EVALUATE EQUIVALENT ANISOTROPIC PARAMETERS.
C
      CALL EXPAND (GB,BISO,B)
      CALL EXPAND (GB,CISO,C)
      CALL RESET (SYST,B)
      CALL RESET (SYST,C)
C
C CALCULATE U'S FROM B'S.
C
      PI=3.14159
      UISO=BISO/(8*PI**2)
      VISO=CISO/(8*PI**2)
      DO I=1,3
      DO J=I,3
        U(I,J)=B(I,J)/(2*PI**2*SQRT(GB(I,I)*GB(J,J)))
        U(J,I)=U(I,J)
        V(I,J)=C(I,J)/(2*PI**2*SQRT(GB(I,I)*GB(J,J)))
        V(J,I)=V(I,J)
      END DO
      END DO
      WRITE (ILP,600) ATIME,ADATE,TITLE
      WRITE (ILP,619) BISO,CISO,UISO,VISO,
     & B(1,1),B(2,2),B(3,3),B(1,2),B(1,3),B(2,3),
     & C(1,1),C(2,2),C(3,3),C(1,2),C(1,3),C(2,3),
     & U(1,1),U(2,2),U(3,3),U(1,2),U(1,3),U(2,3),
     & V(1,1),V(2,2),V(3,3),V(1,2),V(1,3),V(2,3)
 619  FORMAT (/1X,
     &'BISO = ',E12.4,'     CISO = ',E12.4/1X,
     &'UISO = ',E12.4,'     VISO = ',E12.4//1X,
     &'EQUIVALENT ANISOTROPIC ATOMIC DISPLACEMENT PARAMETERS'/1X,
     &'AND DISTRIBUTION DISPERSION PARAMETERS:'//1X,
     &'      B(1,1)      B(2,2)      B(3,3)      B(1,2)      B(1,3)   ',
     &'   B(2,3)'/1X,6E12.4/1X,
     &'      C(1,1)      C(2,2)      C(3,3)      C(1,2)      C(1,3)   ',
     &'   C(2,3)'/1X,6E12.4//1X,
     &'      U(1,1)      U(2,2)      U(3,3)      U(1,2)      U(1,3)   ',
     &'   U(2,3)'/1X,6E12.4/1X,
     &'      V(1,1)      V(2,2)      V(3,3)      V(1,2)      V(1,3)   ',
     &'   V(2,3)'/1X,6E12.4)
C
C REFINE SCALEK, B(I,J), AND C(I,J).
C
      T=R
      IF (CISO.EQ.0) THEN
        NPAR=7
      ELSE
        NPAR=13
      END IF
      CALL KBCIJ (IO1,NDATA,GA,GB,SCALEK,B,C,SYST,NCYCLE,R,Z,NFREE,NPAR,
     & COV)
C
C CALCULATE EQUIVALENT ISOTROPIC B AND C.
C
      CALL CONTRACT (GA,B,BISO)
      CALL CONTRACT (GA,C,CISO)
C
C CALCULATE U'S FROM B'S.
C
      UISO=BISO/(8*PI**2)
      VISO=CISO/(8*PI**2)
      DO I=1,3
      DO J=I,3
        U(I,J)=B(I,J)/(2*PI**2*SQRT(GB(I,I)*GB(J,J)))
        U(J,I)=U(I,J)
        V(I,J)=C(I,J)/(2*PI**2*SQRT(GB(I,I)*GB(J,J)))
        V(J,I)=V(I,J)
      END DO
      END DO
      WRITE (ILP,600) ATIME,ADATE,TITLE
      WRITE (ILP,624) NCYCLE,R,Z,NFREE,NPAR
      WRITE (ILP,625) SCALEK,
     & B(1,1),B(2,2),B(3,3),B(1,2),B(1,3),B(2,3),
     & C(1,1),C(2,2),C(3,3),C(1,2),C(1,3),C(2,3),
     & U(1,1),U(2,2),U(3,3),U(1,2),U(1,3),U(2,3),
     & V(1,1),V(2,2),V(3,3),V(1,2),V(1,3),V(2,3)
      WRITE (ILP,626) BISO,CISO,UISO,VISO
      WRITE (ILP,600) ATIME,ADATE,TITLE
      BIBJ(1)=GB(1,1)
      BIBJ(2)=SQRT(GB(1,1)*GB(2,2))
      BIBJ(3)=SQRT(GB(1,1)*GB(3,3))
      BIBJ(4)=GB(2,2)
      BIBJ(5)=SQRT(GB(2,2)*GB(3,3))
      BIBJ(6)=GB(3,3)
      DO I=1,13
        VAR(I)=COV(I,I)
        IF (VAR(I).EQ.0) VAR(I)=1
      END DO
      WRITE (ILP,627) (SQRT(COV(1+I,1+I))/(2*PI**2*BIBJ(I)),I=1,6),
     & (SQRT(COV(7+I,7+I))/(2*PI**2*BIBJ(I)),I=1,6),
     & (SQRT(COV(I,I)),I=1,13),
     & (I,(NINT(100*COV(I,J)/SQRT(VAR(I)*VAR(J))),J=1,13),I=1,13)
 624  FORMAT (/1X,
     &'SCALE AND ANISOTROPIC ATOMIC DISPLACEMENT AND DISTRIBUTION DISP',
     &'ERSION'/1X,
     &'PARAMETERS:'//1X,
     &'  YOBS  = FSQ(MEAS)/(EPSILON*SUM F(J)**2)'/1X,
     &'  YCALC = (SCALEK**-2)*EXP(-2*H(J)*B(I,J)*H(I) + 2*(H(J)*C(I,J)',
     &'*H(I))**2)'//1X,
     &'STATISTICS OF FIT:'//1X,
     &'NCYCLE = ',I11/1X,
     &'R      = ',E11.4/1X,
     &'Z      = ',E11.4/1X,
     &'NFREE  = ',I11/1X,
     &'NPAR   = ',I11)
 625  FORMAT (/1X,
     &'RESULTS:'//1X,
     &'SCALEK, B(I,J), C(I,J), U(I,J), AND B(I,J) ARE DEFINED ACCORDIN',
     &'G TO:'//1X,
     &'  F(ABSOLUTE) = SCALEK*F(RELATIVE)'/1X,
     &'  f(T) = f(0)*EXP(-SUM(I)SUM(J) H(I)*H(J)*B(I,J))'/1X,
     &'  f(T) = f(0)*EXP(-2*PI**2*SUM(I)SUM(J) H(I)*H(J)*ASTAR(I)*ASTA',
     &'R(J)*U(I,J))'/1X,
     &'  <EXP(-2*H(J)*B(I,J)*H(I))> = EXP(-2*H(J)*<B(I,J)>*H(I) + 2*(H',
     &'(J)*<C(I,J)>*H(I))**2)'//1X,
     &'SCALEK'/1X,E12.4//1X,
     &'OVERALL MEAN-SQUARE ATOMIC DISPLACEMENT PARAMETERS:'/1X,
     &'      B(1,1)      B(2,2)      B(3,3)      B(1,2)      B(1,3)   ',
     &'   B(2,3)'/1X,6E12.4,' (DIMENSIONLESS)'//1X,
     &'B-DISTRIBUTION DISPERSION PARAMETERS:'/1X,
     &'      C(1,1)      C(2,2)      C(3,3)      C(1,2)      C(1,3)   ',
     &'   C(2,3)'/1X,6E12.4,' (DIMENSIONLESS)'//1X,
     &'OVERALL MEAN-SQUARE ATOMIC DISPLACEMENTS:'/1X,
     &'      U(1,1)      U(2,2)      U(3,3)      U(1,2)      U(1,3)   ',
     &'   U(2,3)'/1X,6E12.4,' SQUARED ANGSTROMS'//1X,
     &'U-DISTRIBUTION DISPERSION PARAMETERS:'/1X,
     &'      V(1,1)      V(2,2)      V(3,3)      V(1,2)      V(1,3)   ',
     &'   V(2,3)'/1X,6E12.4,' SQUARED ANGSTROMS')
 626  FORMAT (/1X,
     &'EQUIVALENT ISOTROPIC VALUES (ALL IN SQUARED ANGSTROM):'//1X,
     &'BISO = ',E12.4,' = <B>'/1X,
     &'CISO = ',E12.4,' = <(B - <B>)**2>**(1/2)'//1X,
     &'UISO = ',E12.4,' = <u**2>'/1X,
     &'VISO = ',E12.4,' = <(u**2 - <u**2>)**2>**(1/2)'//1X,
     &'WHERE BISO = 8*PI**2*UISO')
 627  FORMAT (/1X,
     &'PARAMETER E.S.D.''S'//1X,
     &'NOTE THAT THE ORDER OF THE MEAN-SQUARE DISPLACEMENT PARAMETER E'',
     &RRORS IS'/1X,
     &'11, 12, 13, 22, 23, 33  NOT  11, 22, 33, 12, 13, 23.'/1X,
     &'                        ---'//!X,
     &'  SIGMA(U11)  SIGMA(U12)  SIGMA(U13)  SIGMA(U22)  SIGMA(U23)  S',
     &'IGMA(U33)'/1X,6E12.2,' SQUARED ANGSTROMS'//1X,
     &'  SIGMA(V11)  SIGMA(V12)  SIGMA(V13)  SIGMA(V22)  SIGMA(V23)  S',
     &'IGMA(V33)'/1X,6E12.2,' SQUARED ANGSTROMS'//1X,
     &'    SIGMA(K)'/1X,E12.2//1X,
     &'  SIGMA(B11)  SIGMA(B12)  SIGMA(B13)  SIGMA(B22)  SIGMA(B13)  S',
     &'IGMA(B33)'/1X,6E12.2,' (DIMENSIONLESS)'//1X,
     &'  SIGMA(C11)  SIGMA(C12)  SIGMA(C13)  SIGMA(C12)  SIGMA(C23)  S',
     &'IGMA(C33)'/1X,6E12.2,' (DIMENSIONLESS)'//1X,
     &'CORRELATION MATRIX 100*CORR(I,J):'//1X,
     &'       K B11 B12 B13 B22 B23 B33 C11 C12 C13 C22 C23 C33'//1X,
     &'       1   2   3   4   5   6   7   8   9  10  11  12  13'/
     &/(1X,I2,2X,13I4))
C
C IF THE ANISOTROPIC MODEL DID NOT IMPROVE THE FIT, RECALL THE ISOTROPIC
C MODEL 
C
      IF ((R-T)/T.LT.-1E-6) THEN
        CONTINUE
      ELSE
        CALL RECALL (TEMP,NCYCLE,R,Z,NFREE,NPAR,SCALEK,B,C,COV)
      END IF
C
C WRITE OUTPUT FILE.
C
      CLOSE (UNIT=IO1,STATUS='DELETE')
      OPEN (UNIT=IO1,FILE='eval.dat',STATUS='NEW')
      WRITE (IO1,700) ATIME
      WRITE (IO1,700) ADATE
      WRITE (IO1,700) TITLE
      WRITE (IO1,702) ITYPE
      WRITE (IO1,700) FILE
      WRITE (IO1,700) LATT
      WRITE (IO1,700) PTGP
      WRITE (IO1,701) CELL
      WRITE (IO1,700) XRAY
      WRITE (IO1,702) NEL
      DO I=1,NEL
        WRITE (IO1,703) EL(I),XI(I),ZI(I),A1(I),B1(I),A2(I),B2(I)
        WRITE (IO1,704) A3(I),B3(I),A4(I),B4(I),C0(I),FP(I),FPP(I)
      END DO
      WRITE (IO1,701) SMIN,SMAX
      WRITE (IO1,706) ASOLV,BSOLV
      WRITE (IO1,705) SCALEK
      WRITE (IO1,705) B(1,1),B(1,2),B(1,3),B(2,2),B(2,3),B(3,3)
      WRITE (IO1,705) C(1,1),C(1,2),C(1,3),C(2,2),C(2,3),C(3,3)
      WRITE (IO1,706) ((COV(I,J),J=1,13),I=1,13)
 700  FORMAT (1X,A)
 701  FORMAT (1X,3F10.4,3F10.3)
 702  FORMAT (1X,I6)
 703  FORMAT (1X,A2,F10.2,F6.0,4(1X,F10.6))
 704  FORMAT (1X,5(1X,F10.6),2X,2(1X,F10.6))
 705  FORMAT (1X,6E12.5)
 706  FORMAT (1X,7E10.3)
      CLOSE (UNIT=IO1,STATUS='KEEP')
      STOP 'PROGRAM LEVY FINIS!'
      END
C-----------------------------------------------------------------------
      SUBROUTINE AGREE (IFILE,NDATA,SCALEK,B,C,NPAR,NFREE,Z,R)
C
C AGREEMENT STATISTICS
C
      DIMENSION H(3),B(3,3),C(3,3)
      IF (NPAR.LT.7) THEN
        BISO=B(1,1)
        CISO=C(1,1)
      END IF
      NOBS=0
      CHISQ=0
      SUMWY2=0
      REWIND IFILE
      DO I=1,NDATA
        READ (IFILE) J,K,L,SSQ,Y,SIGY,MULT,VAR
        IF (NPAR.LT.7) THEN
          HTBH=BISO*SSQ
          HTCH=CISO*SSQ
        ELSE
          H(1)=J
          H(2)=K
          H(3)=L
          HTBH=VTMV(H,B)
          HTCH=VTMV(H,C)
        END IF
        YCALC=(1/SCALEK**2)*EXP(-2*(HTBH-HTCH**2))
        W=MULT/(SIGY**2+YCALC**2*VAR)
        CHISQ=CHISQ+W*(Y-YCALC)**2
        SUMWY2=SUMWY2+W*Y**2
        NOBS=NOBS+MULT
      END DO
      NFREE=NOBS-NPAR
      Z=SQRT(CHISQ/NFREE)
      R=SQRT(CHISQ/SUMWY2)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE CONTRACT (G,B,BISO)
C
C CONTRACT TENSOR BIJ'S TO EQUIVALENT SCALAR BISO.
C
      DIMENSION G(3,3),B(3,3)
      BISO=0
      DO I=1,3
      DO J=1,3
        BISO=BISO+G(I,J)*B(I,J)
      END DO
      END DO
      BISO=BISO*4/3
      RETURN
      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 DATAIN (FILE,ITYPE,IO1,IO2,IO3,GINV,S1,S2,SMIN,SMAX,
     & ASOLV,BSOLV,MLATT,IPTGP,ILAUE,NEL,X,A1,B1,A2,B2,A3,B3,A4,B4,C0,
     & FP,FPP,MDATA,NDATA,DATA)
C
C CREATE WORKING DATA FILES.
C
      CHARACTER FILE*80,FORM*11
      DIMENSION GINV(3,3),X(1),A1(1),B1(1),A2(1),B2(1),A3(1),B3(1),
     & A4(1),B4(1),C0(1),FP(1),FPP(1),DATA(1)
C
C COPY DATA FROM THE INPUT REFLECTION DATA FILE TO A TEMPORARY FILE:
C
      IF (ITYPE.LT.3) THEN
        FORM='FORMATTED'
      ELSE
        FORM='UNFORMATTED'
      END IF
      OPEN (UNIT=IO1,FILE=FILE,STATUS='OLD',FORM=FORM)
      OPEN (UNIT=IO2,STATUS='SCRATCH',FORM='UNFORMATTED')
      IF (S1.EQ.0) S1=-9
      IF (S2.EQ.0) S2=+9
      SMIN=+9
      SMAX=-9
      IF (IPTGP.LT.ILAUE) THEN
        JMIN=+999
        JMAX=-999
        KMIN=+999
        KMAX=-999
        LMIN=+999
        LMAX=-999
      END IF
      MDATA=0
      NDATA=0
      IEND=0
 1    CALL READ1 (ITYPE,IO1,IEND,IH,IK,IL,FSQ,SIGFSQ)
      IF (IEND.NE.0) GO TO 9
      S=SINTHL(IH,IK,IL,GINV)
      SMIN=MIN(SMIN,S)
      SMAX=MAX(SMAX,S)
      IF (ASOLV.NE.0) THEN
C
C COMPUTE BULK SOLVENT COMPENSATION FACTOR.
C
        Q=1/(1-ASOLV*EXP(-BSOLV*S**2))
      ELSE
        Q=1
      END IF
      WRITE (IO2) IH,IK,IL,Q**2*FSQ,Q**2*SIGFSQ
      MDATA=MDATA+1
      IF (IPTGP.LT.ILAUE) THEN
C
C STORE THE LAUE-UNIQUE INDEX LIMITS.
C
        J=IH
        K=IK
        L=IL
        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')
      ENDFILE IO2
      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
        DO I=1,NJ*NK*NL+NK*NL+NL
          DATA(I)=0
        END DO
        REWIND IO2
        DO I=1,MDATA
          READ (IO2) 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 FILES.
C
      REWIND IO2
      OPEN (UNIT=IO1,STATUS='SCRATCH',FORM='UNFORMATTED')
      OPEN (UNIT=IO3,STATUS='SCRATCH',FORM='UNFORMATTED',
     & ACCESS='DIRECT',RECL=3)
      DO I=1,MDATA
        READ (IO2) IH,IK,IL,FSQ,SIGFSQ
C
C GET ACENTRIC OR CENTRIC REFLECTION FLAG.
C
        CALL CSYM (IPTGP,IH,IK,IL,IC)
C
C GET POINT GROUP DEGENERACY AND MULTIPLICITY.
C
        CALL ESYM (IPTGP,IH,IK,IL,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
          J=IH
          K=IK
          L=IL
          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(IH,IK,IL,GINV)
C
C CALCULATE  Y = FSQ(MEAS)/(EPSILON*SUM F(J)**2)  AND  SIGMA(Y).
C
        CALL FCALC (NEL,X,A1,B1,A2,B2,A3,B3,A4,B4,C0,FP,FPP,S,SUMFSQ)
        XI=S**2
        YI=FSQ/(IE*SUMFSQ)
        SIGYI=SIGFSQ/(IE*SUMFSQ)
        IF (YI.LT.-3*SIGYI) YI=-3*SIGYI
C
C WILSON DISTRIBUTION VARIANCE = 1 OR 2 FOR ACENTRIC OR CENTRIC
C DISTRIBUTIONS, RESPECTIVELY.
C
        IF (IC.EQ.0) V=1
        IF (IC.EQ.1) V=2
C
C WRITE A DATA FILE FOR SORTING AND PLOTTING.
C
        WRITE (IO3,REC=I) XI,YI,M
C
C WRITE THE LEAST-SQUARES DATA FILE.
C
        IF (S1.LE.S.AND.S.LE.S2) THEN
          WRITE (IO1) IH,IK,IL,XI,YI,SIGYI,M,V
          NDATA=NDATA+1
        END IF
      END DO
      S1=MAX(S1,SMIN)
      S2=MIN(S2,SMAX)
      CLOSE (UNIT=IO2,STATUS='DELETE')
      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 'LEVY/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-----------------------------------------------------------------------
      FUNCTION DELTA(I,J)
C
C KRONECKER DELTA
C
      IF (I.EQ.J) THEN
        DELTA=1
      ELSE
        DELTA=0
      END IF
      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 'LEVY/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------------------------------------------------------------------------
      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-----------------------------------------------------------------------
      SUBROUTINE EXPAND (GINV,BISO,B)
C
C EXPAND SCALAR BISO TO EQUIVALENT TENSOR BIJ'S.
C
      DIMENSION GINV(3,3),B(3,3)
      DO I=1,3
      DO J=I,3
        B(I,J)=0.25*BISO*GINV(I,J)
        B(J,I)=B(I,J)
      END DO
      END DO
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE FCALC (N,X,A1,B1,A2,B2,A3,B3,A4,B4,C0,FP,FPP,S,
     & SUMFSQ)
      DIMENSION X(1),A1(1),B1(1),A2(1),B2(1),A3(1),B3(1),A4(1),B4(1),
     & C0(1),FP(1),FPP(1)
      SUMFSQ=0
      DO I=1,N
        SUMFSQ=SUMFSQ+X(I)*((A1(I)*EXP(-B1(I)*S**2)
     &                       +A2(I)*EXP(-B2(I)*S**2)
     &                        +A3(I)*EXP(-B3(I)*S**2)
     &                         +A4(I)*EXP(-B4(I)*S**2)+C0(I)+FP(I))**2
     &                          +FPP(I)**2)
      END DO
      RETURN   
      END
C-----------------------------------------------------------------------
      SUBROUTINE FTABLE (EL,Z,RM,A1,B1,A2,B2,A3,B3,A4,B4,C,XRAY,FP,FPP)
C
C NEUTRAL ATOM SCATTERING FACTOR COEFFICIENTS FROM CROMER & WABER
C (INTERNATIONAL TABLES, VOL. IV, 1974)
C
C F(S) = A1*EXP(-B1*S**2) + A2*EXP(-B2*S**2) + A3*EXP(-B3*S**2) 
C          + A4*EXP(-B4*S**2) + C
C
      CHARACTER EL*2,SYMBOL*2,XRAY*2
      DIMENSION SYMBOL(2*98),RMASS(98),AA1(98),BB1(98),AA2(98),BB2(98),
     & AA3(98),BB3(98),AA4(98),BB4(98),CC(98),FPMO(98),FPPMO(98),
     & FPCU(98),FPPCU(98),FPKISS(98)
C
C ELEMENT SYMBOLS
C
      DATA SYMBOL /
     & 'H ',                              'He',
     & 'Li','Be','B ','C ','N ','O ','F ','Ne',
     & 'Na','Mg','Al','Si','P ','S ','Cl','Ar',
     & 'K ','Ca',
     &           'Sc','Ti','V ','Cr','Mn','Fe','Co','Ni','Cu','Zn',
     &           'Ga','Ge','As','Se','Br','Kr',
     & 'Rb','Sr',
     &           'Y ','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd',
     &           'In','Sn','Sb','Te','I ','Xe',
     & 'Cs','Ba',
     &           'La',
     &                'Ce','Pr','Nd','Pm','Sm','Eu','Gd',
     &                  'Tb','Dy','Ho','Er','Tm','Yb','Lu',
     &                'Hf','Ta','W ','Re','Os','Ir','Pt','Au','Hg',
     &           'Tl','Pb','Bi','Po','At','Rn',
     & 'Fr','Ra',
     &           'Ac',
     &                'Th','Pa','U ','Np','Pu','Am','Cm','Bk','Cf',
C
     & 'H ',                              'HE',
     & 'LI','BE','B ','C ','N ','O ','F ','NE',
     & 'NA','MG','AL','SI','P ','S ','CL','AR',
     & 'K ','CA',
     &           'SC','TI','V ','CR','MN','FE','CO','NI','CU','ZN',
     &           'GA','GE','AS','SE','BR','KR',
     & 'RB','SR',
     &           'Y ','ZR','NB','MO','TC','RU','RH','PD','AG','CD',
     &           'IN','SN','SB','TE','I ','XE',
     & 'CS','BA',
     &           'LA',
     &                'CE','PR','ND','PM','SM','EU','GD',
     &                  'TB','DY','HO','ER','TM','YB','LU',
     &                'HF','TA','W ','RE','OS','IR','PT','AU','HG',
     &           'TL','PB','BI','PO','AT','RN',
     & 'FR','RA',
     &           'AC',
     &                'TH','PA','U ','NP','PU','AM','CM','BK','CF'/
C
C RELATIVE ATOMIC MASSES (NATURAL ISOTOPIC ABUNDANCES)
C
      DATA RMASS / 1.00797, 4.0026, 6.939, 9.0122, 10.811, 12.01115,
     & 14.0067, 15.9994, 18.9984, 20.183, 22.9898, 24.312, 26.9815,
     & 28.086, 30.9738, 32.064, 35.453, 39.948, 39.102, 40.08, 44.956,
     & 47.90, 50.942, 51.996, 54.938, 55.847, 58.933, 58.71, 63.54,
     & 65.37, 69.72, 72.59, 74.922, 78.96, 79.909, 83.80, 85.47, 87.62,
     & 88.905, 91.22, 92.906, 95.94, 98, 101.07, 102.905, 106.4,
     & 107.870, 112.40, 114.82, 118.69, 121.75, 127.60, 126.904, 131.30,
     & 132.905, 137.34, 138.91, 140.12, 140.907, 144.24, 147, 150.35,
     & 151.96, 157.25, 158.924, 162.50, 164.930, 167.26, 168.934,
     & 173.04, 174.97, 178.49, 180.948, 183.85, 186.2, 190.2, 192.2,
     & 195.09, 196.967, 200.59, 204.37, 207.19, 208.980, 210, 210, 222,
     & 223, 226, 227, 232.038, 231, 238.03, 237, 242, 243, 247, 247,
     & 249/
C
C SCATTERING FACTOR COEFFICIENTS
C
      DATA AA1 /   0.493002,  0.873400,  1.128200,  1.591900,  2.054500,
     &  2.310000, 12.212590,  3.048500,  3.539200,  3.955300,  4.762600,
     &  5.420400,  6.420200,  6.291500,  6.434500,  6.905300, 11.460390,
     &  7.484500,  8.218600,  8.626600,  9.189000,  9.759500, 10.297100,
     & 10.640600, 11.281900, 11.769490, 12.284090, 12.837590, 13.338000,
     & 14.074290, 15.235400, 16.081600, 16.672300, 17.000590, 17.178890,
     & 17.355490, 17.178400, 17.566290, 17.776000, 17.876490, 17.614190,
     &  3.702500, 19.130090, 19.267390, 19.295700, 19.331890, 19.280800,
     & 19.221400, 19.162390, 19.188900, 19.641790, 19.964400, 20.147200,
     & 20.293300, 20.389200, 20.336100, 20.578000, 21.167090, 22.044000,
     & 22.684490, 23.340490, 24.004190, 24.627390, 25.070900, 25.897590,
     & 26.507000, 26.904900, 27.656290, 28.181900, 28.664090, 28.947600,
     & 29.143990, 29.202390, 29.081800, 28.762100, 28.189400, 27.304900,
     & 27.005900, 16.881890, 20.680890, 27.544600, 31.061700, 33.368890,
     & 34.672600, 35.316290, 35.563090, 35.929900, 35.763000, 35.659690,
     & 35.564490, 35.884700, 36.022790, 36.187390, 36.525400, 36.670590,
     & 36.648800, 36.788100, 36.918500/
      DATA BB1 /  10.510890,  9.103700,  3.954600, 43.642700, 23.218500,
     & 20.843900,  0.005700, 13.277090, 10.282500,  8.404200,  3.285000,
     &  2.827500,  3.038700,  2.438600,  1.906700,  1.467900,  0.010400,
     &  0.907200, 12.794890, 10.442090,  9.021300,  7.850800,  6.865700,
     &  6.103800,  5.340900,  4.761100,  4.279100,  3.878500,  3.582800,
     &  3.265500,  3.066900,  2.850900,  2.634500,  2.409800,  2.172300,
     &  1.938400,  1.788800,  1.556400,  1.402900,  1.276180,  1.188650,
     &  0.277200,  0.864132,  0.808520,  0.751536,  0.698655,  0.644600,
     &  0.594600,  0.547600,  5.830300,  5.303400,  4.817420,  4.347000,
     &  3.928200,  3.569000,  3.216000,  2.948170,  2.812190,  2.773930,
     &  2.662480,  2.562700,  2.472740,  2.387900,  2.253410,  2.242560,
     &  2.180200,  2.070510,  2.073560,  2.028590,  1.988900,  1.901820,
     &  1.832620,  1.773330,  1.720290,  1.671910,  1.629030,  1.592790,
     &  1.512930,  0.461100,  0.545000,  0.655150,  0.690200,  0.704000,
     &  0.700999,  0.685870,  0.663100,  0.646453,  0.616341,  0.589092,
     &  0.563359,  0.547751,  0.529300,  0.511929,  0.499384,  0.483629,
     &  0.465154,  0.451018,  0.437533/
      DATA AA2 /   0.322912,  0.630900,  0.750800,  1.127800,  1.332600,
     &  1.020000,  3.132200,  2.286800,  2.641200,  3.112500,  3.173600,
     &  2.173500,  1.900200,  3.035300,  4.179100,  5.203400,  7.196400,
     &  6.772300,  7.439800,  7.387300,  7.367900,  7.355800,  7.351100,
     &  7.353700,  7.357300,  7.357300,  7.340900,  7.292000,  7.167600,
     &  7.031800,  6.700600,  6.374700,  6.070100,  5.819600,  5.235800,
     &  6.728600,  9.643500,  9.818400, 10.294590, 10.947990, 12.014390,
     & 17.235590, 11.094790, 12.918190, 14.350090, 15.501700, 16.688500,
     & 17.644390, 18.559600, 19.100490, 19.045500, 19.013790, 18.994900,
     & 19.029800, 19.106200, 19.296990, 19.598990, 19.769500, 19.669690,
     & 19.684690, 19.609490, 19.425790, 19.088590, 19.079800, 18.218500,
     & 17.638300, 17.294000, 16.428490, 15.885100, 15.434490, 15.220800,
     & 15.172590, 15.229290, 15.430000, 15.718890, 16.154990, 16.729590,
     & 17.763900, 18.591290, 19.041700, 19.158400, 13.063690, 12.951000,
     & 15.473290, 19.021100, 21.281600, 23.054700, 22.906400, 23.103190,
     & 23.421900, 23.294790, 23.412790, 23.596400, 23.808300, 24.099190,
     & 24.409600, 24.773600, 25.199490/
      DATA BB2 /  26.125700,  3.356800,  1.052400,  1.862300,  1.021000,
     & 10.207500,  9.893300,  5.701100,  4.294400,  3.426200,  8.842200,
     & 79.261090,  0.742600, 32.333690, 27.156990, 22.215100,  1.166200,
     & 14.840700,  0.774800,  0.659900,  0.572900,  0.500000,  0.438500,
     &  0.392000,  0.343200,  0.307200,  0.278400,  0.256500,  0.247000,
     &  0.233300,  0.241200,  0.251600,  0.264700,  0.272600, 16.579600,
     & 16.562300, 17.315090, 14.098790, 12.800600, 11.916000, 11.765990,
     &  1.095800,  8.144870,  8.434670,  8.217580,  7.989290,  7.472600,
     &  6.908900,  6.377600,  0.503100,  0.460700,  0.420885,  0.381400,
     &  0.344000,  0.310700,  0.275600,  0.244475,  0.226836,  0.222087,
     &  0.210628,  0.202088,  0.196451,  0.194200,  0.181951,  0.196143,
     &  0.202172,  0.197940,  0.223545,  0.238849,  0.257119,  9.985190,
     &  9.599900,  9.370460,  9.225900,  9.092270,  8.979480,  8.865530,
     &  8.811740,  8.621600,  8.448400,  8.707510,  2.357600,  2.923800,
     &  3.550780,  3.974580,  4.069100,  4.176190,  3.871350,  3.651550,
     &  3.462040,  3.415190,  3.325300,  3.253960,  3.263710,  3.206470,
     &  3.089970,  3.046190,  3.007750/
      DATA AA3 /   0.140191,  0.311200,  0.617500,  0.539100,  1.097900,
     &  1.588600,  2.012500,  1.546300,  1.517000,  1.454600,  1.267400,
     &  1.226900,  1.593600,  1.989100,  1.780000,  1.437900,  6.255600,
     &  0.653900,  1.051900,  1.589900,  1.640900,  1.699100,  2.070300,
     &  3.324000,  3.019300,  3.522200,  4.003400,  4.443800,  5.615800,
     &  5.165200,  4.359100,  3.706800,  3.431300,  3.973100,  5.637700,
     &  5.549300,  5.139900,  5.422000,  5.726290,  5.417320,  4.041830,
     & 12.887590,  4.649010,  4.863370,  4.734250,  5.295370,  4.804500,
     &  4.461000,  4.294800,  4.458500,  5.037100,  6.144870,  7.513800,
     &  8.976700, 10.661990, 10.887990, 11.372690, 11.851300, 12.385600,
     & 12.774000, 13.123490, 13.439590, 13.760290, 13.851790, 14.316690,
     & 14.559590, 14.558300, 14.977890, 15.154190, 15.308690, 15.100000,
     & 14.758600, 14.513500, 14.432700, 14.556400, 14.930500, 15.611490,
     & 15.713100, 25.558190, 21.657500, 15.538000, 18.442000, 16.587700,
     & 13.113800,  9.498870,  8.003700, 12.143890, 12.473890, 12.597700,
     & 12.747300, 14.189100, 14.949090, 15.640190, 16.770700, 17.341500,
     & 17.399000, 17.891900, 18.331690/
      DATA BB3 /   3.142360, 22.927590, 85.390500,103.483000, 60.349790,
     &  0.568700, 28.997490,  0.323900,  0.261500,  0.230600,  0.313600,
     &  0.380800, 31.547190,  0.678500,  0.526000,  0.253600, 18.519390,
     & 43.898300,213.186900, 85.748390,136.108000, 35.633800, 26.893790,
     & 20.262600, 17.867400, 15.353500, 13.535900, 12.176300, 11.396590,
     & 10.316300, 10.780500, 11.446800, 12.947890, 15.237190,  0.260900,
     &  0.226100,  0.274800,  0.166400,  0.125599,  0.117622,  0.204785,
     & 11.003990, 21.570690, 24.799690, 25.874890, 25.205200, 24.660500,
     & 24.700800, 25.849890, 26.890890, 27.907390, 28.528390, 27.766000,
     & 26.465890, 24.387890, 20.207300, 18.772590, 17.608300, 16.766900,
     & 15.885000, 15.100890, 14.399600, 13.754590, 12.933090, 12.664790,
     & 12.189900, 11.440690, 11.360400, 10.997500, 10.664690,  0.261033,
     &  0.275116,  0.295977,  0.321703,  0.350500,  0.382661,  0.417916,
     &  0.424593,  1.482600,  1.572900,  1.963470,  8.618000,  8.793700,
     &  9.556420, 11.382390, 14.042200, 23.105190, 19.988690, 18.598990,
     & 17.830900, 16.923490, 16.092690, 15.362190, 14.945500, 14.313590,
     & 13.434590, 12.894590, 12.404390/
      DATA AA4 /   0.040810,  0.178000,  0.465300,  0.702900,  0.706800,
     &  0.865000,  1.166300,  0.867000,  1.024300,  1.125100,  1.112800,
     &  2.307300,  1.964600,  1.541000,  1.490800,  1.586300,  1.645500,
     &  1.644200,  0.865900,  1.021100,  1.468000,  1.902100,  2.057100,
     &  1.492200,  2.244100,  2.304500,  2.348800,  2.380000,  1.673500,
     &  2.410000,  2.962300,  3.683000,  4.277900,  4.354300,  3.985100,
     &  3.537500,  1.529200,  2.669400,  3.265880,  3.657210,  3.533460,
     &  3.742900,  2.712630,  1.567560,  1.289180,  0.605844,  1.046300,
     &  1.602900,  2.039600,  2.466300,  2.682700,  2.523900,  2.273500,
     &  1.990000,  1.495300,  2.695900,  3.287190,  3.330490,  2.824280,
     &  2.851370,  2.875160,  2.896040,  2.922700,  3.545450,  2.953540,
     &  2.965770,  3.638370,  2.982330,  2.987060,  2.989630,  3.716010,
     &  4.300130,  4.764920,  5.119820,  5.441740,  5.675890,  5.833770,
     &  5.783700,  5.860000,  5.967600,  5.525930,  5.969600,  6.469200,
     &  7.025880,  7.425180,  7.443300,  2.112530,  3.210970,  4.086550,
     &  4.807030,  4.172870,  4.188000,  4.185500,  3.479470,  3.493310,
     &  4.216650,  4.232840,  4.243910/
      DATA BB4 /  57.799690,  0.982100,168.261000,  0.542000,  0.140300,
     & 51.651190,  0.582600, 32.908900, 26.147590, 21.718390,129.423900,
     &  7.193700, 85.088590, 81.693690, 68.164500, 56.171990, 47.778390,
     & 33.392890, 41.684090,178.436900, 51.353100,116.104900,102.477900,
     & 98.739890, 83.754300, 76.880490, 71.169200, 66.342100, 64.812600,
     & 58.709700, 61.413490, 54.762490, 47.797190, 43.816290, 41.432800,
     & 39.397200,164.934000,132.376000,104.354000, 87.662700, 69.795700,
     & 61.658400, 86.847190, 94.292800, 98.606200, 76.898600, 99.815590,
     & 87.482490, 92.802900, 83.957100, 75.282500, 70.840300, 66.877590,
     & 64.265790,213.904000,167.201900,133.123900,127.113000,143.643900,
     &137.902900,132.720900,128.007000,123.173900,101.397900,115.361900,
     &111.873900, 92.656600,105.703000,102.960900,100.417000, 84.329800,
     & 72.029000, 63.364390, 57.055990, 52.086100, 48.164700, 45.001090,
     & 38.610300, 36.395590, 38.324600, 45.814890, 47.257900, 48.009290,
     & 47.004500, 45.471490, 44.247290,150.645000,142.324900,117.020000,
     & 99.172190,105.251000,100.613000, 97.490790,105.979900,102.272900,
     & 88.483390, 86.003000, 83.788100/
      DATA CC /    0.003038,  0.006400,  0.037700,  0.038500, -0.193200,
     &  0.215600,-11.529000,  0.250800,  0.277600,  0.351500,  0.676000,
     &  0.858400,  1.115100,  1.140700,  1.114900,  0.866900, -9.557400,
     &  1.444500,  1.422800,  1.375100,  1.332900,  1.280700,  1.219900,
     &  1.183200,  1.089600,  1.036900,  1.011800,  1.034100,  1.191000,
     &  1.304100,  1.718900,  2.131300,  2.531000,  2.840900,  2.955700,
     &  2.825000,  3.487300,  2.506400,  1.912130,  2.069290,  3.755910,
     &  4.387500,  5.404280,  5.378740,  5.328000,  5.265930,  5.179000,
     &  5.069400,  4.939100,  4.782100,  4.590900,  4.352000,  4.071200,
     &  3.711800,  3.335200,  2.773100,  2.146780,  1.862640,  2.058300,
     &  1.984860,  2.028760,  2.209630,  2.574500,  2.419600,  3.583240,
     &  4.297280,  4.567960,  5.920460,  6.756210,  7.566720,  7.976280,
     &  8.581540,  9.243540,  9.887500, 10.472000, 11.000490, 11.472200,
     & 11.688300, 12.065790, 12.608900, 13.174590, 13.411800, 13.578200,
     & 13.677000, 13.710800, 13.690500, 13.724690, 13.621100, 13.526590,
     & 13.431400, 13.428700, 13.396590, 13.357290, 13.381190, 13.359190,
     & 13.288700, 13.275400, 13.267390/
C
C ANOMALOUS DISPERSION CORRECTIONS FOR MO AND CU X-RAYS FROM CROMER &
C LIBERMAN (INTERNATIONAL TABLES, VOL. IV, 1974)
C
      DATA FPMO /
     &  0.000,  0.000,  0.000,  0.000,  0.000,  0.002,  0.004,  0.008,
     &  0.014,  0.021,  0.030,  0.042,  0.056,  0.072,  0.090,  0.110,
     &  0.132,  0.155,  0.179,  0.203,  0.226,  0.248,  0.267,  0.284,
     &  0.295,  0.301,  0.299,  0.285,  0.263,  0.222,  0.163,  0.081,
     & -0.030, -0.178, -0.374, -0.652, -1.044, -1.657, -2.951, -2.965,
     & -2.197, -1.825, -1.590, -1.420, -1.287, -1.177, -1.085, -1.005,
     & -0.936, -0.873, -0.816, -0.772, -0.726, -0.684, -0.644, -0.613,
     & -0.588, -0.564, -0.530, -0.535, -0.530, -0.533, -0.542, -0.564,
     & -0.591, -0.619, -0.666, -0.723, -0.795, -0.884, -0.988, -1.118,
     & -1.258, -1.421, -1.598, -1.816, -2.066, -2.352, -2.688, -3.084,
     & -3.556, -4.133, -4.861, -5.924, -7.444, -8.862, -7.912, -7.620,
     & -7.725, -8.127, -8.960,-10.673,-11.158, -9.725, -8.926, -8.416,
     & -7.990, -7.683/
      DATA FPPMO /
     &  0.000,  0.000,  0.000,  0.000,  0.001,  0.002,  0.003,  0.006,
     &  0.010,  0.016,  0.025,  0.036,  0.052,  0.071,  0.095,  0.124,
     &  0.159,  0.201,  0.250,  0.306,  0.372,  0.446,  0.530,  0.624,
     &  0.729,  0.845,  0.973,  1.113,  1.266,  1.431,  1.609,  1.801,
     &  2.007,  2.223,  2.456,  2.713,  2.973,  3.264,  3.542,  0.560,
     &  0.621,  0.688,  0.759,  0.836,  0.919,  1.007,  1.101,  1.202,
     &  1.310,  1.424,  1.546,  1.675,  1.812,  1.958,  2.119,  2.282,
     &  2.452,  2.632,  2.845,  3.018,  3.225,  3.442,  3.669,  3.904,
     &  4.151,  4.410,  4.678,  4.958,  5.248,  5.548,  5.858,  6.185,
     &  6.523,  6.872,  7.232,  7.605,  7.990,  8.388,  8.798,  9.223,
     &  9.659, 10.102, 10.559, 11.042,  9.961, 10.403,  7.754,  8.105,
     &  8.472,  8.870,  9.284,  9.654,  4.148,  4.330,  4.511,  4.697,
     &  4.908,  5.107/
      DATA FPCU /
     &  0.000,  0.000,  0.001,  0.003,  0.008,  0.017,  0.029,  0.047,
     &  0.069,  0.097,  0.129,  0.165,  0.204,  0.244,  0.283,  0.319,
     &  0.348,  0.366,  0.365,  0.341,  0.285,  0.189,  0.035, -0.198,
     & -0.568, -1.179, -2.464, -2.956, -2.019, -1.612, -1.354, -1.163,
     & -1.011, -0.879, -0.767, -0.665, -0.574, -0.465, -0.386, -0.314,
     & -0.248, -0.191, -0.145, -0.105, -0.077, -0.059, -0.060, -0.079,
     & -0.126, -0.194, -0.287, -0.418, -0.579, -0.783, -1.022, -1.334,
     & -1.716, -2.170, -2.939, -3.431, -4.357, -5.696, -7.718, -9.242,
     & -9.498,-10.423,-12.255, -9.733, -8.488, -7.701, -7.713, -6.715,
     & -6.351, -6.048, -5.790, -5.581, -5.391, -5.233, -5.096, -4.990,
     & -4.883, -4.818, -4.776, -4.756, -4.772, -4.787, -4.833, -4.898,
     & -4.994, -5.091, -5.216, -5.359, -5.529, -5.712, -5.930, -6.176,
     & -6.498, -6.798/
      DATA FPPCU /
     &  0.000,  0.000,  0.000,  0.001,  0.004,  0.009,  0.018,  0.032,
     &  0.053,  0.083,  0.124,  0.177,  0.246,  0.330,  0.434,  0.557,
     &  0.702,  0.872,  1.066,  1.286,  1.533,  1.807,  2.110,  2.443,
     &  2.808,  3.204,  3.608,  0.509,  0.589,  0.678,  0.777,  0.886,
     &  1.006,  1.139,  1.283,  1.439,  1.608,  1.820,  2.025,  2.245,
     &  2.482,  2.735,  3.005,  3.296,  3.605,  3.934,  4.282,  4.653,
     &  5.045,  5.459,  5.894,  6.352,  6.835,  7.348,  7.904,  8.460,
     &  9.036,  9.648, 10.535, 10.933, 11.614, 12.320, 11.276, 11.946,
     &  9.242,  9.748,  3.704,  3.937,  4.181,  4.432,  4.693,  4.977,
     &  5.271,  5.577,  5.891,  6.221,  6.566,  6.925,  7.297,  7.686,
     &  8.089,  8.505,  8.930,  9.383,  9.843, 10.317, 10.803, 11.296,
     & 11.799, 12.330, 12.868, 13.409, 13.967, 14.536, 15.087, 15.634,
     & 16.317, 16.930/
C
C WAVELENGTH-INDEPENDENT KISSEL AND PRATT CORRECTIONS TO THE CROMER AND
C LIBERMAN FP VALUES [L. KISSEL AND R.H. PRATT (1990).  ACTA CRYST. A46,
C 170-175].
C
      DATA FPKISS /
     &  0.000,  0.000,  0.001,  0.000,  0.001,  0.001,  0.002,  0.003,
     &  0.004,  0.004,  0.006,  0.008,  0.008,  0.011,  0.012,  0.014,
     &  0.017,  0.020,  0.022,  0.025,  0.028,  0.031,  0.035,  0.039,
     &  0.042,  0.048,  0.052,  0.057,  0.061,  0.067,  0.073,  0.079,
     &  0.085,  0.092,  0.099,  0.106,  0.114,  0.122,  0.130,  0.138,
     &  0.147,  0.156,  0.166,  0.175,  0.186,  0.196,  0.207,  0.219,
     &  0.230,  0.242,  0.255,  0.267,  0.281,  0.294,  0.308,  0.323,
     &  0.338,  0.354,  0.369,  0.386,  0.402,  0.419,  0.437,  0.455,
     &  0.474,  0.493,  0.512,  0.532,  0.553,  0.574,  0.596,  0.617,
     &  0.640,  0.663,  0.687,  0.711,  0.736,  0.762,  0.788,  0.814,
     &  0.842,  0.870,  0.899,  0.928,  0.957,  0.988,  1.018,  1.050,
     &  1.083,  1.115,  1.149,  1.184,  1.219,  1.255,  1.292,  1.330,
     &  1.368,  1.407/
      IF (EL.EQ.'J ') EL='I '
      DO I=1,2*98
        IF (EL.EQ.SYMBOL(I)) GO TO 1
      END DO
      STOP 'LEVY/FTABLE:  ELEMENT SYMBOL NOT IN THE TABLE (Z .LE. 98)'
 1    CONTINUE
      IF (I.GT.98) I=I-98
      Z=I
      RM=RMASS(I)
      A1=AA1(I)
      B1=BB1(I)
      A2=AA2(I)
      B2=BB2(I)
      A3=AA3(I)
      B3=BB3(I)
      A4=AA4(I)
      B4=BB4(I)
      C=CC(I)
      IF (XRAY.EQ.'MO') THEN
        FP=FPMO(I)+FPKISS(I)
        FPP=FPPMO(I)
      ELSE IF (XRAY.EQ.'CU') THEN
        FP=FPCU(I)+FPKISS(I)
        FPP=FPPCU(I)
      ELSE
        FP=0
        FPP=0
      END IF
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE KBCISO (IFILE,NDATA,SCALEK,BISO,CISO,NCYCLE,R,Z,NFREE,
     & NPAR,COV)
C
C NON-LINEAR LEAST-SQUARES REFINEMENT OF SCALE AND ISOTROPIC
C DISPLACEMENT AND DISTRIBUTION DISPERSION PARAMETERS
C
C YOBS  = FSQ(MEAS)/[EPSILON*SUM F(J)**2]
C YCALC = SCALEK**(-2)*EXP{-2*[BISO - (CISO*S)**2]*S**2}
C
      DIMENSION AA(3,3),XX(3),BB(3),COV(13,13)
      DOUBLE PRECISION AA
      MPAR=3
      NPAR=3
C
C CALCULATE INITIAL AGREEMENT STATISTICS.
C
      CALL AGREE (IFILE,NDATA,SCALEK,BISO,CISO,NPAR,NFREE,Z,R)
C
C ITERATE THROUGH LEAST-SQUARES LOOP.
C
      NCYCLE=0
 1    NCYCLE=NCYCLE+1
C
C STORE THE AGREEMENT INDICES FROM THE PRECEDING CYCLE.
C
      ROLD=R
      ZOLD=Z
C
C REFINE DISPERSION PARAMETER ONLY IF IT IS NONZERO.
C
      IF (CISO.EQ.0) NPAR=2
C
C EVALUATE THE DERIVATIVES, AND BUILD THE NORMAL MATRIX AND VECTOR.
C
      DO J=1,NPAR
        BB(J)=0
      DO K=1,NPAR
        AA(J,K)=0
      END DO
      END DO
      REWIND IFILE
      DO I=1,NDATA
        READ (IFILE) IH,IK,IL,SSQ,YI,SIGYI,MULT,VAR
        Y=(1/SCALEK**2)*EXP(-2*(BISO*SSQ-(CISO*SSQ)**2))
        DYDK=-2*Y/SCALEK
        DYDB=-2*Y*SSQ
        DYDS=+4*Y*SSQ**2*CISO
        XX(1)=DYDK
        XX(2)=DYDB
        XX(3)=DYDS
        WI=MULT/(SIGYI**2+Y**2*VAR)
        DO J=1,NPAR
        DO K=J,NPAR
          AA(J,K)=AA(J,K)+WI*XX(J)*XX(K)
          AA(K,J)=AA(J,K)
        END DO
          BB(J)=BB(J)+WI*XX(J)*(YI-Y)
        END DO
      END DO
C
C INVERT THE NORMAL MATRIX.
C
      CALL MATINV (NPAR,MPAR,AA,DET)
      IF (DET.EQ.0) STOP 'LEVY/KBCISO:  SINGULAR NORMAL MATRIX'
C
C CALCULATE THE FULL PARAMETER SHIFTS.
C
      DO I=1,NPAR
        XX(I)=0
      DO J=1,NPAR
        XX(I)=XX(I)+AA(I,J)*BB(J)
      END DO
      END DO
      IF (NPAR.LT.MPAR) THEN
        DO I=NPAR+1,MPAR
          XX(I)=0
        END DO
      END IF
C
C INITIALIZE THE SHIFT-DAMPING FACTOR TO UNITY.
C
      F=2
 2    F=0.5*F
      IF (F.LT.1E-6) THEN
        NCYCLE=NCYCLE-1
        GO TO 9
      END IF
C
C APPLY THE DAMPED PARAMETER SHIFTS.
C
      TK=SCALEK+F*XX(1)
      TB=BISO  +F*XX(2)
      IF (NPAR.LE.2) THEN
        TC=0
      ELSE
        TC=CISO+F*XX(3)
      END IF
C
C REQUIRE POSITIVE SCALEK AND NON-NEGATIVE BISO AND CISO.
C
      IF (TK.LE.0.OR.TB.LT.0.OR.TC.LT.0) GO TO 2
C
C FIT A PARABOLA TO THE AGREEMENT INDICES FOR ZERO, HALF, AND FULL
C SHIFTS TO ESTIMATE AN OPTIMUM SHIFT-DAMPING FACTOR.
C
C SHIFT-DAMPING FACTORS, 0 .LT. F .LT. 1, MIGHT BE NECESSARY BECAUSE THE
C WEIGHTS, AS WELL AS THE RESIDUALS, DEPEND ON THE FITTED PARAMETERS.
C THIS COMPOUNDS THE NON-LINEARITY OF THE LEAST-SQUARES PROBLEM, AND
C OVER-SHIFTS CAN LEAD TO A DIVERGENT REFINEMENT.
C
      CALL AGREE (IFILE,NDATA,TK,TB,TC,NPAR,NFREE,Z,R)
      F0=0
      R0=ROLD
      F1=F
      R1=R
      X=0.5*F
      N=0
 31   N=N+1
      IF (N.LT.10) THEN
        F=0.5*(F0+F1)
      ELSE
        IF (F.LT.X) F=F0
        IF (F.GT.X) F=F1
        GO TO 39
      END IF
      TK=SCALEK+F*XX(1)
      TB=BISO  +F*XX(2)
      IF (NPAR.GT.2) TC=CISO+F*XX(3)
      CALL AGREE (IFILE,NDATA,TK,TB,TC,NPAR,NFREE,Z,R)
      RMIN=MIN(R0,R,R1)
      IF (RMIN.EQ.R0) THEN
        R1=R
        F1=F
        GO TO 31
      ELSE IF (RMIN.EQ.R1) THEN
        R0=R
        F0=F
        GO TO 31
      ELSE
        CALL PMIN (F0,R0,F,R,F1,R1,FMIN,RMIN)
        F=FMIN
      END IF
 39   TK=SCALEK+F*XX(1)
      TB=BISO  +F*XX(2)
      IF (NPAR.GT.2) TC=CISO+F*XX(3)
      CALL AGREE (IFILE,NDATA,TK,TB,TC,NPAR,NFREE,Z,R)
C
C TEST FOR CONVERGENCE.
C
      IF (R.LT.ROLD) THEN
C
C CONVERGING OR CONVERGED.  STORE SHIFTED PARAMETERS.
C
        SCALEK=TK
        BISO=TB
        IF (NPAR.GT.2) CISO=TC
C
C EVALUATE THE VARIANCE-COVARIANCE MATRIX.
C
        ZSQ=Z**2
        DO I=1,NPAR
        DO J=1,I
          COV(I,J)=ZSQ*AA(I,J)
          COV(J,I)=COV(I,J)
        END DO
        END DO
C
C TEST MAXIMUM SHIFT-TO-ERROR RATIO.
C
        X=0
        DO I=1,NPAR
          X=MAX(X,F*ABS(XX(I))/SQRT(COV(I,I)))
        END DO
        IF (X.LT.1E-6) THEN
          GO TO 9
        ELSE IF (NCYCLE.LE.20) THEN
          GO TO 1
        END IF
      ELSE
C
C CONVERGED OR DIVERGING.
C
        NCYCLE=NCYCLE-1
        R=ROLD
        Z=ZOLD
      END IF
 9    CONTINUE
      DO I=NPAR+1,13
      DO J=1,I
        COV(I,J)=0
        COV(J,I)=0
      END DO
      END DO
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE KBCIJ (IFILE,NDATA,G,GINV,SCALEK,B,C,SYSTEM,NCYCLE,R,Z,
     & NFREE,NPAR,COV)
C
C REFINE THE FIT OF THE SCALE FACTOR, OVERALL ANISOTROPIC ATOMIC
C DISPLACEMENT PARAMETERS, AND DISTRIBUTION DISPERSION PARAMETERS.
C
C YOBS = FSQ(MEAS)/[EPSILON*SUM F(J)**2]
C YCAL = SCALEK**(-2)*EXP(-2*(H(J)*B(I,J)*H(I) - (H(J)*C(I,J)*H(I))**2))
C
      CHARACTER SYSTEM*12
      DIMENSION G(3,3),GINV(3,3),H(3),B(3,3),C(3,3),AA(13,13),
     & XX(13),BB(13),COV(13,13),TB(3,3),TC(3,3)
      DOUBLE PRECISION AA
      MPAR=13
      NPAR=13
C
C CALCULATE INITIAL AGREEMENT STATISTICS.
C
      CALL AGREE (IFILE,NDATA,SCALEK,B,C,NPAR,NFREE,Z,R)
C
C ITERATE THROUGH LEAST-SQUARES LOOP.
C
      NCYCLE=0
 1    NCYCLE=NCYCLE+1
C
C STORE THE AGREEMENT INDICES FROM THE PRECEDING CYCLE.
C
      ROLD=R
      ZOLD=Z
C
C REFINE DISPERSION PARAMETERS ONLY IF THE C-TENSOR IS NONZERO.
C
      IF (C(1,1).EQ.0) NPAR=7
C
C EVALUATE THE DERIVATIVES, AND BUILD THE NORMAL MATRIX AND VECTOR.
C
      DO J=1,NPAR
        BB(J)=0
      DO K=1,NPAR
        AA(J,K)=0
      END DO
      END DO
      REWIND IFILE
      DO I=1,NDATA
        READ (IFILE) IH,IK,IL,SSQ,YI,SIGYI,MULT,VAR
        H(1)=IH
        H(2)=IK
        H(3)=IL
        HTBH=VTMV(H,B)
        HTCH=VTMV(H,C)
        Y=(1/SCALEK**2)*EXP(-2*(HTBH-HTCH**2))
        DYDK=-2*Y/SCALEK
        XX(1)=DYDK
        L=1
        DO J=1,3
        DO K=J,3
          DJK=(2-DELTA(J,K))*H(J)*H(K)
          DYDB=-2*Y*DJK
          DYDC=+4*Y*DJK*HTCH
          L=L+1
          XX(L)=DYDB
          XX(L+6)=DYDC
        END DO
        END DO
        WI=MULT/(SIGYI**2+Y**2*VAR)
        DO J=1,NPAR
        DO K=J,NPAR
          AA(J,K)=AA(J,K)+WI*XX(J)*XX(K)
          AA(K,J)=AA(J,K)
        END DO
          BB(J)=BB(J)+WI*XX(J)*(YI-Y)
        END DO
      END DO
C
C INVERT THE NORMAL MATRIX.
C
      CALL MATINV (NPAR,MPAR,AA,DET)
      IF (DET.EQ.0) STOP 'LEVY/KBCIJ:  SINGULAR NORMAL MATRIX'
C
C CALCULATE THE FULL PARAMETER SHIFTS.
C
      DO I=1,NPAR
        XX(I)=0
      DO J=1,NPAR
        XX(I)=XX(I)+AA(I,J)*BB(J)
      END DO
      END DO
      IF (NPAR.LT.MPAR) THEN
        DO I=NPAR+1,MPAR
          XX(I)=0
        END DO
      END IF
C
C INITIALIZE THE SHIFT-DAMPING FACTOR TO UNITY.
C
      F=2
 2    F=0.5*F
      IF (F.LT.1E-6) THEN
        NCYCLE=NCYCLE-1
        GO TO 9
      END IF
C
C APPLY THE DAMPED PARAMETER SHIFTS.
C
      TK=SCALEK+F*XX(1)
C
C REQUIRE POSITIVE SCALEK.
C
      IF (TK.LE.0) GO TO 2
      L=1
      DO J=1,3
      DO K=J,3
        L=L+1
        TB(J,K)=B(J,K)+F*XX(L)
        TC(J,K)=0
        IF (NPAR.GT.7) TC(J,K)=C(J,K)+F*XX(L+6)
      END DO
      END DO
      CALL RESET (SYSTEM,TB)
      CALL RESET (SYSTEM,TC)
C
C REQUIRE POSITIVE DEFINITE B- AND C-TENSORS.
C
      CALL POSDEF (TB,G,GINV,SYSTEM)
      IF (TB(1,1).LE.0) GO TO 2
      IF (NPAR.GT.7) THEN
        CALL POSDEF (TC,G,GINV,SYSTEM)
        IF (TC(1,1).LE.0) GO TO 2
      END IF
C
C FIT A PARABOLA TO THE AGREEMENT INDICES FOR ZERO, HALF, AND FULL
C SHIFTS TO ESTIMATE AN OPTIMUM SHIFT-DAMPING FACTOR.
C
C SHIFT-DAMPING FACTORS, 0 .LT. F .LT. 1, MIGHT BE NECESSARY BECAUSE THE
C WEIGHTS, AS WELL AS THE RESIDUALS, DEPEND ON THE FITTED PARAMETERS.
C THIS COMPOUNDS THE NON-LINEARITY OF THE LEAST-SQUARES PROBLEM, AND
C OVER-SHIFTS CAN LEAD TO A DIVERGENT REFINEMENT.
C
      CALL AGREE (IFILE,NDATA,TK,TB,TC,NPAR,NFREE,Z,R)
      F0=0
      R0=ROLD
      F1=F
      R1=R
      X=0.5*F
      N=0
 31   N=N+1
      IF (N.LT.10) THEN
        F=0.5*(F0+F1)
      ELSE
        IF (F.LT.X) F=F0
        IF (F.GT.X) F=F1
        GO TO 39
      END IF
      TK=SCALEK+F*XX(1)
      L=1
      DO J=1,3
      DO K=J,3
        L=L+1
        TB(J,K)=B(J,K)+F*XX(L)
        IF (NPAR.GT.7) TC(J,K)=C(J,K)+F*XX(L+6)
      END DO
      END DO
      CALL RESET (SYSTEM,TB)
      CALL RESET (SYSTEM,TC)
      CALL AGREE (IFILE,NDATA,TK,TB,TC,NPAR,NFREE,Z,R)
      RMIN=MIN(R0,R,R1)
      IF (RMIN.EQ.R0) THEN
        R1=R
        F1=F
        GO TO 31
      ELSE IF (RMIN.EQ.R1) THEN
        R0=R
        F0=F
        GO TO 31
      ELSE
        CALL PMIN (F0,R0,F,R,F1,R1,FMIN,RMIN)
        F=FMIN
      END IF
 39   TK=SCALEK+F*XX(1)
      L=1
      DO J=1,3
      DO K=J,3
        L=L+1
        TB(J,K)=B(J,K)+F*XX(L)
        IF (NPAR.GT.7) TC(J,K)=C(J,K)+F*XX(L+6)
      END DO
      END DO
      CALL RESET (SYSTEM,TB)
      CALL RESET (SYSTEM,TC)
      CALL AGREE (IFILE,NDATA,TK,TB,TC,NPAR,NFREE,Z,R)
C
C TEST FOR CONVERGENCE.
C
      IF (R.LT.ROLD) THEN
C
C CONVERGING OR CONVERGED.  STORE SHIFTED PARAMETERS.
C
        SCALEK=TK
        DO J=1,3
        DO K=1,3
          B(J,K)=TB(J,K)
          IF (NPAR.GT.7) C(J,K)=TC(J,K)
        END DO
        END DO
C
C EVALUATE THE VARIANCE-COVARIANCE MATRIX.
C
        DO I=1,NPAR
        DO J=1,I
          COV(I,J)=Z**2*AA(I,J)
          COV(J,I)=COV(I,J)
        END DO
        END DO
C
C TEST MAXIMUM SHIFT-TO-ERROR RATIO.
C
        X=0
        DO I=1,NPAR
          X=MAX(X,F*ABS(XX(I))/SQRT(COV(I,I)))
        END DO
        IF (X.LT.1E-6) THEN
          GO TO 9
        ELSE IF (NCYCLE.LE.20) THEN
          GO TO 1
        END IF
      ELSE
C
C CONVERGED OR DIVERGING.
C
        NCYCLE=NCYCLE-1
        R=ROLD
        Z=ZOLD
      END IF
 9    CONTINUE
      IF (NPAR.LT.13) THEN
        DO I=NPAR+1,13
        DO J=1,I
          COV(I,J)=0
          COV(J,I)=0
        END DO
        END DO
      END IF
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE MATINV (N,M,A,D)
C
C INVERT A SQUARE MATRIX BY GAUSS-JORDAN ELIMINATION, AND CALCULATE ITS
C DETERMINANT.
C
C A - INPUT MATRIX, WHICH IS REPLACED BY ITS INVERSE
C N - ORDER OF MATRIX
C   - LOGICAL SIZE OF ARRAY A
C M - PHYSICAL SIZE OF ARRAY A, M .GE. N
C D - DETERMINANT OF THE INPUT MATRIX
C
C RETURNS D = 0 TO FLAG A SINGULAR MATRIX.
C
C MATRIX IS SCALED TO AVOID ARITHMETIC OVERFLOW OR UNDERFLOW DURING
C INVERSION.
C
C MATRIX SCALING USES A DIAGONAL MATRIX Q WITH
C
C Q(I,I) = 1/SQRT(ABS(A(I,I))).
C
C FORM B = Q A Q,
C
C B(I,J) = Q(I,I)*A(I,J)*Q(J,J),
C
C AND INVERT B TO OBTAIN B**-1.  THEN,
C
C   B**-1   = (Q A Q)**-1
C           = (A Q)**-1 Q**-1
C   B**-1 Q = Q**-1 A**-1
C Q B**-1 Q = A**-1
C
C IN ADDITION,
C
C DET(A) = DET(B)/PROD(I=1,N) Q(I,I)**2
C
C IF N EXCEEDS MAXN, THE WORK VECTORS Q, IK, AND JK MUST BE RE-
C DIMENSIONED TO A LENGTH OF N.
C
      PARAMETER (MAXN=25)
      DIMENSION A(M,M),Q(MAXN),IK(MAXN),JK(MAXN)
      DOUBLE PRECISION A,Q,AMAX,T
C
C STORE RECIPROCAL SQUARE-ROOT DIAGONAL MAGNITUDES FOR SCALING.
C
      DO 10 I=1,N
        T=SQRT(ABS(A(I,I)))
        IF (T.EQ.0) THEN
          D=0
          RETURN
        ELSE
          Q(I)=1/T
        END IF
 10   CONTINUE
C
C SCALE MATRIX.
C
      DO 11 I=1,N
      DO 11 J=1,N
        A(I,J)=A(I,J)*Q(I)*Q(J)
 11   CONTINUE
C
C INVERT SCALED MATRIX AND CALCULATE ITS DETERMINANT.
C
      D=1
      DO 100 K=1,N
C
C FIND LARGEST ELEMENT A(I,J) IN REST OF MATRIX.
C
        AMAX=0
        DO 30 I=K,N
        DO 30 J=K,N
          IF (ABS(A(I,J)).GT.ABS(AMAX)) THEN
            AMAX=A(I,J)
            IK(K)=I
            JK(K)=J
          END IF
 30     CONTINUE
        D=D*AMAX
        IF (D.EQ.0) RETURN
C
C INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN A(K,K).
C
        I=IK(K)
        IF (I.GT.K) THEN
          DO 50 J=1,N
            T=A(K,J)
            A(K,J)=A(I,J)
            A(I,J)=-T
 50       CONTINUE
        END IF
        J=JK(K)
        IF (J.GT.K) THEN
          DO 60 I=1,N
            T=A(I,K)
            A(I,K)=A(I,J)
            A(I,J)=-T
 60       CONTINUE
        END IF
C
C ACCUMULATE ELEMENTS OF INVERSE MATRIX.
C
        DO 70 I=1,N
          IF (I.NE.K) A(I,K)=-A(I,K)/AMAX
 70     CONTINUE
        DO 80 I=1,N
        DO 80 J=1,N
          IF (I.NE.K.AND.J.NE.K) A(I,J)=A(I,J)+A(I,K)*A(K,J)
 80     CONTINUE
        DO 90 J=1,N
          IF (J.NE.K) A(K,J)=+A(K,J)/AMAX
 90     CONTINUE
        A(K,K)=1/AMAX
 100  CONTINUE
C
C RESTORE ORDERING OF MATRIX.
C
      DO 130 K=N,1,-1
        J=IK(K)
        IF (J.GT.K) THEN
          DO 110 I=1,N
            T=A(I,K)
            A(I,K)=-A(I,J)
            A(I,J)=T
 110      CONTINUE
        END IF
        I=JK(K)
        IF (I.GT.K) THEN
          DO 120 J=1,N
            T=A(K,J)
            A(K,J)=-A(I,J)
            A(I,J)=T
 120      CONTINUE
        END IF
 130  CONTINUE
C
C SCALE INVERSE MATRIX.
C
      DO 140 I=1,N
      DO 140 J=1,N
        A(I,J)=A(I,J)*Q(I)*Q(J)
 140  CONTINUE
C
C SCALE DETERMINANT, IF POSSIBLE.
C
      T=0
      DO 150 I=1,N
        T=T+LOG(Q(I))
 150  CONTINUE
      T=LOG(ABS(D))-2*T
C
C TEST AGAINST RANGE OF MACHINE ALLOWED MAGNITUDES.
C
C EIGHT-BIT EXPONENT HAS MAXIMUM MAGNITUDE 2**7 = 128.
C XMIN = 2**(-128)     = 0.29E-38      LOG(XMIN) = -88.7
C XMAX = (2**(+128))/2 = 1.70E+38      LOG(XMAX) = +88.0
C
      IF (-87.LT.T.AND.T.LT.+87) D=(D/ABS(D))*EXP(T)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE MATINV3 (A,B,D)
C
C INVERT SYMMETRICAL (3 X 3) MATRIX.
C
C A IS THE MATRIX; B IS ITS INVERSE; D IS ITS DETERMINANT.
C
      DIMENSION A(3,3),B(3,3)
      DOUBLE PRECISION A,B,D
      DO I=1,3
        K=1+MOD(I+1,3)
        M=1+MOD(I+3,3)
      DO J=1,3
        L=1+MOD(J+1,3)
        N=1+MOD(J+3,3)
        B(I,J)=A(L,K)*A(N,M)-A(L,M)*A(N,K)
      END DO
      END DO
      D=0
      DO I=1,3
        D=D+A(I,1)*B(1,I)
      END DO
      DO I=1,3
      DO J=1,3
        B(I,J)=B(I,J)/D
      END DO
      END DO
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE METRIC (A,V,GA,GB)
C
C CALCULATE RECIPROCAL LATTICE PARAMETERS AND DIRECT AND RECIPROCAL
C LATTICE METRIC TENSORS FROM DIRECT LATTICE PARAMETERS.
C
      DIMENSION A(6),B(6),GA(3,3),GB(3,3)
      PI=ACOS(-1.0)
      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 PMIN (X1,Y1,X2,Y2,X3,Y3,XMIN,YMIN)
C
C CALCULATE THE COEFFICIENTS OF THE PARABOLA,
C   Y = C0 + C1*X + C2*X**2,
C THROUGH THREE POINTS,
C   (X1,Y1), (X2,Y2), AND (X3,Y3) WITH Y2 .LT. Y1, Y3,
C AND FIND THE MINIMUM OF THE CURVE.
C
      DIMENSION A(3,3),AINV(3,3),B(3),C(3)
      DOUBLE PRECISION A,AINV,DETA
      A(1,1)=1
      A(1,2)=X1
      A(1,3)=X1*X1
      A(2,1)=1
      A(2,2)=X2
      A(2,3)=X2*X2
      A(3,1)=1
      A(3,2)=X3
      A(3,3)=X3*X3
      B(1)=Y1
      B(2)=Y2
      B(3)=Y3
      CALL MATINV3 (A,AINV,DETA)
      DO I=1,3
        C(I)=0
      DO J=1,3
        C(I)=C(I)+AINV(I,J)*B(J)
      END DO
      END DO
      XMIN=X2
      IF (C(3).GT.0) THEN
        XMIN=-0.5*C(2)/C(3)
        IF (XMIN.LT.X1) XMIN=X1
        IF (XMIN.GT.X3) XMIN=X3
        YMIN=C(1)+C(2)*XMIN+C(3)*XMIN*XMIN
      ELSE
        XMIN=X2
        YMIN=Y2
      END IF
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE POSDEF (B,G,GINV,SYSTEM)
C
C TEST NECESSARY AND SUFFICIENT CONDITION FOR POSITIVE-DEFINITE, 
C SYMMETRIC B-TENSOR.
C
      CHARACTER SYSTEM*12
      DIMENSION B(3,3),G(3,3),GINV(3,3)
      T1=B(1,1)
      T2=B(1,1)*B(2,2)-B(2,1)*B(1,2)
      T3=B(1,1)*B(2,2)*B(3,3)+B(1,2)*B(2,3)*B(3,1)+B(1,3)*B(2,1)*B(3,2)
     & -B(3,1)*B(2,2)*B(1,3)-B(3,2)*B(2,3)*B(1,1)-B(3,3)*B(2,1)*B(1,2)
      IF (MIN(T1,T2,T3).LE.0) THEN
        CALL CONTRACT (G,B,BISO)
        IF (BISO.LT.0) BISO=0
        CALL EXPAND (GINV,BISO,B)
        CALL RESET (SYSTEM,B)
      END IF
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE READ1 (ITYPE,IUNIT,IEND,IH,IK,IL,FSQ,SIGFSQ)
 1    IF (ITYPE.EQ.0) READ (IUNIT,*,END=9) IH,IK,IL,FSQ,SIGFSQ,F,SIGF
      IF (ITYPE.EQ.1) READ (IUNIT,*,END=9) IH,IK,IL,FSQ,SIGFSQ
      IF (ITYPE.EQ.2) READ (IUNIT,*,END=9) IH,IK,IL,F  ,SIGF
      IF (ITYPE.EQ.3) READ (IUNIT,  END=9) IH,IK,IL,FSQ,SIGFSQ,F,SIGF
      IF (ITYPE.EQ.4) READ (IUNIT,  END=9) IH,IK,IL,FSQ,SIGFSQ
      IF (ITYPE.EQ.5) READ (IUNIT,  END=9) IH,IK,IL,F  ,SIGF
      IF (ITYPE.EQ.2.OR.ITYPE.EQ.5) 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-----------------------------------------------------------------------
      SUBROUTINE RECALL (T,NCYCLE,R,Z,NFREE,NPAR,SCALEK,B,C,COV)
      DIMENSION T(109),B(3,3),C(3,3),COV(13,13)
      NCYCLE=T(1)
      R     =T(2)
      Z     =T(3)
      NFREE =T(4)
      NPAR  =T(5)
      SCALEK=T(6)
      I=6
      DO J=1,3
      DO K=J,3
        I=I+1
        B(J,K)=T(I)
      END DO
      END DO
      DO J=1,3
      DO K=J,3
        I=I+1
        C(J,K)=T(I)
      END DO
      END DO
      DO J=1,13
      DO K=J,13
        I=I+1
        COV(J,K)=T(I)
      END DO
      END DO
      DO I=1,3
      DO J=I,3
        B(J,I)=B(I,J)
        C(J,I)=C(I,J)
      END DO
      END DO
      DO I=1,13
      DO J=1,13
        COV(J,I)=COV(I,J)
      END DO
      END DO
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE RESET (SYSTEM,B)
      CHARACTER SYSTEM*12
      DIMENSION B(3,3)
      IF (SYSTEM.EQ.'MONOCLINIC') THEN
        B(1,2)=0
        B(2,3)=0
      ELSE IF (SYSTEM.EQ.'ORTHORHOMBIC') THEN
        B(1,2)=0
        B(1,3)=0
        B(2,3)=0
      ELSE IF (SYSTEM.EQ.'TETRAGONAL') THEN
        Q=(B(1,1)+B(2,2))/2
        B(1,1)=Q
        B(2,2)=Q
        B(1,2)=0
        B(1,3)=0
        B(2,3)=0
      ELSE IF (SYSTEM.EQ.'TRIGONAL'.OR.SYSTEM.EQ.'HEXAGONAL') THEN
        Q=(B(1,1)+B(2,2))/2
        B(1,1)=Q
        B(2,2)=Q
        B(1,2)=Q/2
        B(1,3)=0
        B(2,3)=0
      ELSE IF (SYSTEM.EQ.'RHOMBOHEDRAL') THEN
        Q=(B(1,1)+B(2,2)+B(3,3))/3
        B(1,1)=Q
        B(2,2)=Q
        B(3,3)=Q
        Q=(B(1,2)+B(1,3)+B(2,3))/3
        B(1,2)=Q
        B(1,3)=Q
        B(2,3)=Q
      ELSE IF (SYSTEM.EQ.'CUBIC') THEN
        Q=(B(1,1)+B(2,2)+B(3,3))/3
        B(1,1)=Q
        B(2,2)=Q
        B(3,3)=Q
        B(1,2)=0
        B(1,3)=0
        B(2,3)=0
      END IF
      B(2,1)=B(1,2)
      B(3,1)=B(1,3)
      B(3,2)=B(2,3)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE SORT (N,DATA,INDEX)
C
C INDEXES THE ARRAY DATA(N) AND RETURNS THE ARRAY INDEX(N) SORTED SUCH
C THAT THE VALUES DATA(INDEX(I)) ARE IN ASCENDING ORDER FOR I = 1, 2,
C ..., N.  THE INPUT VARIABLES N AND DATA(N) ARE NOT CHANGED.
C
C EMPLOYS THE "HEAPSORT" ALGORITHM.
C
C FORTRAN CODE QUOTED FROM WILLIAM H. PRESS, BRIAN P. FLANNERY, SAUL A.
C TEUKOLSKY, AND WILLIAM T. VETTERLING (1986).  NUMERICAL RECIPIES:  THE
C ART OF SCIENTIFIC COMPUTING, PP. 229-233.  CAMBRIDGE, ENGLAND:
C CAMBRIDGE UNIVERSITY PRESS.
C
      DIMENSION DATA(N),INDEX(N)
      DO 10 I=1,N
        INDEX(I)=I
 10   CONTINUE
      I=N/2+1
      J=N
 1    CONTINUE
      IF (I.GT.1) THEN
        I=I-1
        INDEXT=INDEX(I)
        T=DATA(INDEXT)
      ELSE
        INDEXT=INDEX(J)
        T=DATA(INDEXT)
        INDEX(J)=INDEX(1)
        J=J-1
        IF (J.EQ.1) THEN
          INDEX(1)=INDEXT
          RETURN
        END IF
      END IF
      K=I
      L=I+I
 2    IF (L.LE.J) THEN
        IF (L.LT.J) THEN
          IF (DATA(INDEX(L)).LT.DATA(INDEX(L+1))) L=L+1
        END IF
        IF (T.LT.DATA(INDEX(L))) THEN
          INDEX(K)=INDEX(L)
          K=L
          L=L+L
        ELSE
          L=J+1
        END IF
        GO TO 2
      END IF
      INDEX(K)=INDEXT
      GO TO 1
      END
C-----------------------------------------------------------------------
      SUBROUTINE SPLINE (N,X,Y,YPP)
C
C NATURAL CUBIC SPLINE FIT TO A TABULATED FUNCTION Y(I) = Y(X(I)), I =
C 1, 2,..., N, WITH X(1) .LT. X(2) .LT. ... .LT. X(N).  RETURNS THE
C TABULATED SECOND DERIVATIVES y'' = d2y/dx2 = YPP(I), I = 1, 2,..., N,
C WITH YPP(1) = YPP(N) = 0.
C
C FORTRAN CODE QUOTED FROM WILLIAM H. PRESS, BRIAN P. FLANNERY, SAUL A.
C TEUKOLSKY, AND WILLIAM T. VETTERLING (1986).  NUMERICAL RECIPIES:  THE
C ART OF SCIENTIFIC COMPUTING, PP. 86-89.  CAMBRIDGE, ENGLAND:
C CAMBRIDGE UNIVERSITY PRESS.
C
C MUST HAVE NMAX .GE. N.
C
      PARAMETER (NMAX=1000)
      DIMENSION X(N),Y(N),YPP(N),U(NMAX)
      YPP(1)=0
      YPP(N)=0
      U(1)=0
      U(N)=0
      DO I=2,N-1
        P=(X(I)-X(I-1))/(X(I+1)-X(I-1))
        Q=P*YPP(I-1)+2
        YPP(I)=(P-1)/Q
        U(I)=(6*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1))/
     &   (X(I)-X(I-1)))/(X(I+1)-X(I-1))-P*U(I-1))/Q
      END DO
      DO I=N-1,1,-1
        YPP(I)=YPP(I)*YPP(I+1)+U(I)
      END DO
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE SPLINT (N,X,Y,YPP,XI,YI)
C
C CUBIC SPLINE INTERPOLATION OF A TABULATED FUNCTION Y(I) = Y(X(I)), I =
C 1, 2,..., N, USING THE TABULATED SECOND DERIVATIVES YPP(I) FROM
C SUBROUTINE SPLINE.
C
C FORTRAN CODE QUOTED FROM WILLIAM H. PRESS, BRIAN P. FLANNERY, SAUL A.
C TEUKOLSKY, AND WILLIAM T. VETTERLING (1986).  NUMERICAL RECIPIES:  THE
C ART OF SCIENTIFIC COMPUTING, PP. 86-89.  CAMBRIDGE, ENGLAND:
C CAMBRIDGE UNIVERSITY PRESS.
C
      DIMENSION X(N),Y(N),YPP(N)
      IF (XI.LE.X(1)) THEN
        YI=Y(1)
        RETURN
      END IF
      IF (XI.GE.X(N)) THEN
        YI=Y(N)
        RETURN
      END IF
      I1=1
      I2=N
 1    IF (I2-I1.GT.1) THEN
        I=(I2+I1)/2
        IF (X(I).GT.XI) THEN
          I2=I
        ELSE
          I1=I
        END IF
        GO TO 1
      END IF
      H=X(I2)-X(I1)
      A=(X(I2)-XI)/H
      B=(XI-X(I1))/H
      YI=A*Y(I1)+B*Y(I2)+((A**3-A)*YPP(I1)+(B**3-B)*YPP(I2))*H**2/6
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE STORE (NCYCLE,R,Z,NFREE,NPAR,SCALEK,B,C,COV,T)
      DIMENSION B(3,3),C(3,3),COV(13,13),T(109)
      T(1)=NCYCLE
      T(2)=R
      T(3)=Z
      T(4)=NFREE
      T(5)=NPAR
      T(6)=SCALEK
      I=6
      DO J=1,3
      DO K=J,3
        I=I+1
        T(I)=B(J,K)
      END DO
      END DO
      DO J=1,3
      DO K=J,3
        I=I+1
        T(I)=C(J,K)
      END DO
      END DO
      DO J=1,13
      DO K=J,13
        I=I+1
        T(I)=COV(J,K)
      END DO
      END DO
      RETURN
      END
C-----------------------------------------------------------------------
      FUNCTION SINTHL(IH,IK,IL,GINV)
      DIMENSION H(3),GINV(3,3)
      H(1)=IH
      H(2)=IK
      H(3)=IL
      SINTHL=0.5*SQRT(VTMV(H,GINV))
      RETURN
      END
C-----------------------------------------------------------------------
      FUNCTION VTMV(A,B)
      DIMENSION A(3),B(3,3)
      VTMV=0
      DO I=1,3
      DO J=1,3
        VTMV=VTMV+A(J)*B(I,J)*A(I)
      END DO
      END DO
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE WILSON (IFILE,NDATA,SCALEK,BISO,CISO,COV)
C
C LOGARITHMICALLY LINEARIZED, WILSON-PLOT FIT TO UNAVERAGED INTENSITIES
C TO MINIMIZE
C
C                      FSQ(MEAS)
C CHISQ = SUM W*[------------------- - Q**2]**2 ,
C                EPSILON*SUM F(J)**2
C
C WHERE  Q  = Q0*EXP(Q1*S**2 + Q2*S**4) ,
C        Q0 = SCALEK**2 ,
C        Q1 = -2*<BISO> ,
C        Q2 = 2*<CISO**2> = 2*<(BISO - <BISO>)**2> ,
C AND    S  = SIN(THETA)/LAMBDA .
C
C THE LOGARITHMIC FIT MINIMIZES
C
C                           FSQ(MEAS)
C CHISQ' = SUM W'*{LN ------------------- - [2*LN(SCALEK)
C                     EPSILON*SUM F(J)**2
C                                              - 2*BISO*S**2
C                                                  + 2*CISO**2*S**4]}.
C
C INTENSITIES ARE NOT AVERAGED IN S-SHELLS BEFORE TAKING LOGARITHMS.
C LARGE INTENSITIES ARE EFFECTIVELY DOWN-WEIGHTED ON THE LOGARITHMIC
C SCALE, HENCE SCALEK IS BIASED TO A VALUE THAT IS TOO LARGE.  AND,
C SINCE THE LARGE INTENSITIES ARE AT LOW S-VALUES, BISO IS BIASED TO
C A VALUE THAT IS TOO SMALL AND CISO IS BIASED TO A VALUE THAT IS TOO
C LARGE.
C
      DIMENSION COV(13,13),A(3,3),AINV(3,3),B(3),C(3)
      REAL*8 SUMW,SUMX,SUMY,SUMX2,SUMX3,SUMX4,SUMY2,SUMXY,SUMX2Y
      N=0
      SUMW=0
      SUMX=0
      SUMY=0
      SUMX2=0
      SUMX3=0
      SUMX4=0
      SUMY2=0
      SUMXY=0
      SUMX2Y=0
      REWIND IFILE
      DO I=1,NDATA
        READ (IFILE) J,K,L,SSQ,YOBS,SIGOBS,MULT,VAR
        IF (YOBS.LT.SIGOBS) THEN
C
C THRESHOLD INTENSITY FSQMIN/2 FOR (VAR = 1) ACENTRIC REFLECTIONS, AND
C FSQMIN/3 FOR (VAR = 2) CENTRIC REFLECTIONS.
C
          IF (VAR.EQ.1) YOBS=SIGOBS/2
          IF (VAR.EQ.2) YOBS=SIGOBS/3
        END IF
C
C MINIMIZE CHISQ = SUM W*DELTA**2 = SUM (DELTA/SIGMA(DELTA))**2.
C
C ON THE DATA SCALE:
C
C   DELTA = YOBS - YCALC 
C         = YOBS - (SCALEK**-2)*EXP{-2*[BISO - (CISO*S)**2]*S**2}
C
C   SIGMA = SQRT[SIGMA(YOBS)**2 + SIGMA(YCALC)**2]
C         = SQRT[SIGMA(YOBS)**2 + YCALC**2*VAR]
C
C ON THE LOGARITHMIC FITTING SCALE:
C
C   DELTA = LN(YOBS) - LN(YCALC)
C
C   SIGMA = SQRT{[SIGMA(YOBS)/YOBS]**2 + [SIGMA(YCALC)/YCALC]**2}
C         = SQRT{[SIGMA(YOBS)/YOBS]**2 + VAR}
C
        X=SSQ
        Y=LOG(YOBS)
        W=MULT/((SIGOBS/YOBS)**2+VAR)
C
C ACCUMULATE LEAST-SQUARES SUMS.
C
        N=N+MULT
        SUMW=SUMW+W
        SUMX=SUMX+W*X
        SUMY=SUMY+W*Y
        SUMX2=SUMX2+W*X**2
        SUMX3=SUMX3+W*X**3
        SUMX4=SUMX4+W*X**4
        SUMY2=SUMY2+W*Y**2
        SUMXY=SUMXY+W*X*Y
        SUMX2Y=SUMX2Y+W*X**2*Y
      END DO
C
C STORE AND THEN INVERT THE NORMAL MATRIX
C
      A(1,1)=SUMW
      A(1,2)=SUMX
      A(1,3)=SUMX2
      A(2,1)=A(1,2)
      A(2,2)=SUMX2
      A(2,3)=SUMX3
      A(3,1)=A(1,3)
      A(3,2)=A(2,3)
      A(3,3)=SUMX4
      CALL MINV3 (A,AINV,DET)
C
C STORE THE NORMAL VECTOR AND THEN EVALUATE THE PARAMETERS AND THE
C VAR-COV MATRIX.
C
      B(1)=SUMY
      B(2)=SUMXY
      B(3)=SUMX2Y
      CALL MV (AINV,B,C)
      Q0=C(1)
      Q1=C(2)
      Q2=C(3)
C
C CHISQ = SUM W*(Y - Q)**2
C       = SUM W*Y**2 + SUM W*Q**2 - 2*SUM W*Q*Y
C
C Q     = Q0 + Q1*X + Q2*X**2
C
C Q**2  = Q0**2 + (Q1*X + Q2*X**2)**2 + 2*Q0*(Q1*X + Q2*X**2)
C       = Q0**2 + Q1**2*X**2 + Q2**2*X**4 + 2*Q1*Q2*X**3 
C                                     + 2*Q0*Q1*X + 2*Q0*Q2*X**2
C
      CHISQ=SUMY2+Q0**2*SUMW+Q1**2*SUMX2+Q2**2*SUMX4
     &              +2*(Q0*Q1*SUMX+Q0*Q2*SUMX2+Q1*Q2*SUMX3)
     &                -2*(Q0*SUMY+Q1*SUMXY+Q2*SUMX2Y)
      NPAR=3
      NFREE=N-3
      ZSQ=CHISQ/NFREE
      V0=ZSQ*AINV(1,1)
      V1=ZSQ*AINV(2,2)
      V2=ZSQ*AINV(3,3)
      COV01=ZSQ*AINV(1,2)
      COV02=ZSQ*AINV(1,3)
      COV12=ZSQ*AINV(2,3)
C
C IF THE Q2 = 2*CISO**2 COEFFICIENT OF THE TERM QUADRATIC IN S**2 IS NOT
C STATISTICALLY SIGNIFICANT, RESORT TO A SCALING FUNCTION LINEAR IN
C S**2.
C
      IF (V2.LE.0.OR.Q2.LT.SQRT(ABS(V2))) THEN
        DET=SUMW*SUMX2-SUMX*SUMX
        Q0=(SUMX2*SUMY-SUMX*SUMXY)/DET
        Q1=(-SUMX*SUMY+SUMW*SUMXY)/DET
        CHISQ=SUMY2+Q0**2*SUMW+Q1**2*SUMX2
     &                +2*Q0*Q1*SUMX
     &                  -2*(Q0*SUMY+Q1*SUMXY)
        NPAR=NPAR-1
        NFREE=NFREE+1
        ZSQ=CHISQ/NFREE
        V0=ZSQ*SUMX2/DET
        V1=ZSQ*SUMW/DET
        COV01=ZSQ*(-SUMX/DET)
        Q2=0
        V2=0
        COV02=0
        COV12=0
C
C IF THE Q1 = -2*BISO COEFFICIENT OF THE TERM LINEAR IN S**2 IS NOT
C STATISTICALLY SIGNIFICANT, RESORT TO A CONSTANT SCALING FUNCTION.
C
        IF (V1.LE.0.OR.-Q1.LT.SQRT(ABS(V1))) THEN
          Q0=SUMY/SUMW
          CHISQ=SUMY2+Q0**2*SUMW-2*Q0*SUMY
          NPAR=NPAR-1
          NFREE=NFREE+1
          ZSQ=CHISQ/NFREE
          V0=ZSQ/SUMW
          Q1=0
          V1=0
          COV01=0
        END IF
      END IF
C
C TRANSFORM PARAMETERS AND VAR-COV MATRIX FROM THE LOGARITHMIC FITTING
C SCALE TO THE DATA SCALE, AND CALCULATE AGREEMENT STATISTICS.
C
      SCALEK=EXP(-0.5*Q0)
      BISO=-0.5*Q1
      CISO=SQRT(0.5*Q2)
      DKD0=-0.5*SCALEK
      DBD1=-0.5
      IF (CISO.GT.0) THEN
        DCD2=1/CISO
      ELSE
        DCD2=0
      END IF
      COV(1,1)=DKD0**2*V0
      COV(2,2)=DBD1**2*V1
      COV(3,3)=DCD2**2*V2
      COV(1,2)=DKD0*DBD1*COV01
      COV(1,3)=DKD0*DCD2*COV02
      COV(2,3)=DBD1*DCD2*COV12
      COV(2,1)=COV(1,2)
      COV(3,1)=COV(1,3)
      COV(3,2)=COV(2,3)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE MINV3 (A,B,D)
C
C INVERSE OF A 3 X 3 MATRIX
C
      DIMENSION A(3,3),B(3,3)
C
C DETERMINANT
C
      D=A(1,1)*A(2,2)*A(3,3)+A(1,2)*A(2,3)*A(3,1)+A(1,3)*A(2,1)*A(3,2)
     &  -A(3,1)*A(2,2)*A(1,3)-A(3,2)*A(2,3)*A(1,1)-A(3,3)*A(2,1)*A(1,2)
      IF (D.EQ.0) RETURN
C
C MATRIX OF MINORS
C
      B(1,1)=A(2,2)*A(3,3)-A(3,2)*A(2,3)
      B(1,2)=A(2,1)*A(3,3)-A(3,1)*A(2,3)
      B(1,3)=A(2,1)*A(3,2)-A(3,1)*A(2,2)
      B(2,1)=A(1,2)*A(3,3)-A(3,2)*A(1,3)
      B(2,2)=A(1,1)*A(3,3)-A(3,1)*A(1,3)
      B(2,3)=A(1,1)*A(3,2)-A(3,1)*A(1,2)
      B(3,1)=A(1,2)*A(2,3)-A(2,2)*A(1,3)
      B(3,2)=A(1,1)*A(2,3)-A(2,1)*A(1,3)
      B(3,3)=A(1,1)*A(2,2)-A(2,1)*A(1,2)
C
C MATRIX OF CO-FACTORS (SIGNED MINORS)
C
      DO I=1,3
      DO J=1,3
        B(I,J)=(-1)**(I+J)*B(I,J)
      END DO
      END DO
C
C ADJOINT MATRIX (TRANSPOSED MATRIX OF CO-FACTORS)
C
      DO I=1,3-1
      DO J=I+1,3
        T=B(I,J)
        B(I,J)=B(J,I)
        B(J,I)=T
      END DO
      END DO
C
C INVERSE MATRIX
C
      DO I=1,3
      DO J=1,3
        B(I,J)=B(I,J)/D
      END DO
      END DO
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE MV (A,B,C)
C
C SQUARE MATRIX-COLUMN VECTOR MULTIPLICATION
C
C A(I,J)*B(J) = C(I)
C
      DIMENSION A(3,3),B(3),C(3)
      DO I=1,3
        C(I)=0
      DO J=1,3
        C(I)=C(I)+A(I,J)*B(J)
      END DO
      END DO
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE WPLOT (IFILE,ILP,ATIME,ADATE,TITLE,SCALEK,BISO,CISO,
     & NDATA,DATA,INDEX)
C
C WILSON-PLOT DISPLAYS OF DATA AND FITTED CURVES.
C
C X = [SIN(THETA)/LAMBDA]**2
C Y = FSQ(MEAS)/[EPSILON*SUM(J) F(J)**2]
C M = RECIPROCAL LATTICE POINT GROUP MULTIPLICITY
C
      CHARACTER ATIME*8,ADATE*9,TITLE*80,FILEO*80
      DIMENSION X(100),Y(100),YPP(100),DATA(1),INDEX(1)
C
C SORT ON [SIN(THETA)/LAMBDA]**2.
C
      DO I=1,NDATA
        READ (IFILE,REC=I) XI,YI,M
        DATA(I)=XI
      END DO
      CALL SORT (NDATA,DATA,INDEX)
C
C EVALUATE LOCAL AVERAGE INTENSITIES IN SHELLS OF SIN(THETA)/LAMBDA.
C
C OVERLAPPING MOVING WINDOW AVERAGES
C
      NSHELL=100
      NMOVE=NDATA/NSHELL
      NMOVE=MAX(NMOVE,1)
      NLOCAL=100
      NLOCAL=MIN(NLOCAL,NDATA)
      NLOCAL=MAX(NLOCAL,NMOVE)
      DO I=1,NSHELL
        X(I)=0
        Y(I)=0
        SUMM=0
        SUMX=0
        SUMY=0
        IREC=(I-1)*NMOVE
        DO N=1,NLOCAL
          IREC=IREC+1
          IF (IREC.GT.NDATA) GO TO 5
          READ (IFILE,REC=INDEX(IREC)) XI,YI,M
          SUMM=SUMM+M
          SUMX=SUMX+M*XI
          SUMY=SUMY+M*YI
        END DO
 5      CONTINUE
        IF (SUMM.GT.0) THEN
          X(I)=SUMX/SUMM
          Y(I)=SUMY/SUMM
        END IF
      END DO
      CLOSE (UNIT=IFILE,STATUS='DELETE')
C
C ELIMINATE ANY EMPTY SHELLS.
C
      N=0
      DO I=1,NSHELL
        IF (X(I).NE.0.OR.Y(I).NE.0) THEN
          N=N+1
          X(N)=X(I)
          Y(N)=Y(I)
        END IF
      END DO
      NSHELL=N
C
C ENSURE INCREASING X VALUES FOR SPLINE INTERPOLATION.
C
      N=1
      DO I=2,NSHELL
        IF (X(I).GT.X(I-1)) THEN
          N=N+1
          X(N)=X(I)
          Y(N)=Y(I)
        END IF
      END DO
      NSHELL=N
      CALL SPLINE (NSHELL,X,Y,YPP)
      WRITE (ILP,6000) ATIME,ADATE,TITLE
      CALL XYPLOT (NSHELL,X,Y,YPP,SCALEK,BISO,CISO,ILP,0)
      WRITE (ILP,6001) SCALEK,BISO,CISO
      WRITE (ILP,6000) ATIME,ADATE,TITLE
      CALL XYPLOT (NSHELL,X,Y,YPP,SCALEK,BISO,CISO,ILP,1)
      WRITE (ILP,6002) SCALEK,BISO,CISO
 6000 FORMAT (/1H1,130('-')/1X,'PROGRAM LEVY.  ',A,', ',A,'.  ',A/)
 6001 FORMAT (/1X,
     &'Y = <FMEAS**2/(EPSILON*SUM F(J)**2)>  VERSUS  X = <S**2> = <(SI',
     &'N(THETA)/LAMBDA)**2>'//1X,
     &'"*" FITTED FUNCTION:  Y = K**(-2)*EXP[-2*(BISO - CISO**2*S**2)*',
     &'S**2],  K = ',E10.3,',  BISO = ',E10.3,',  CISO = ',E10.3/1X,
     &'"O" SPLINE-INTERPOLATED DATA CURVE')
 6002 FORMAT (/1X,
     &'LOG(Y) = LOG(<FMEAS**2/(EPSILON*SUM F(J)**2)>)  VERSUS  X = <S*',
     &'*2> = <(SIN(THETA)/LAMBDA)**2>'//1X,
     &'"*" FITTED FUNCTION  Y = K**(-2)*EXP[-2*(BISO - CISO**2S**2)**2',
     &']  K = ',E10.3,'  BISO = ',E10.3,'  CISO = ',E10.3/1X,
     &'"O" SPLINE-INTERPOLATED DATA CURVE')
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE XYPLOT (N,X,Y,YPP,A0,A1,A2,ILP,ILOG)
      PARAMETER (NX=100,NY=50)
      CHARACTER AA*1
      DIMENSION AA(0:NX,0:NY),X(N),Y(N),YPP(N)
C
C BLANK OUT THE PLOT ARRAY.
C
      DO IX=0,NX
      DO IY=0,NY
        AA(IX,IY)=' '
      END DO
      END DO
C
C FILL IN THE GRID MARKS.
C
      DO IX=0,NX
        AA(IX, 0)='.'
        AA(IX,NY)='.'
      END DO
      DO IX=0,NX,10
        AA(IX, 0)='+'
        AA(IX,NY)='+'
      END DO
      DO IY=0,NY
        AA( 0,IY)='.'
        AA(NX,IY)='.'
      END DO
      DO IY=0,NY,10
        AA( 0,IY)='+'
        AA(NX,IY)='+'
      END DO
C
C SCALE PLOT AXES.
C
      XMIN=0
      XMAX=X(N)
      YMIN=+1E9
      YMAX=-1E9
      DO I=1,N
        YMIN=MIN(YMIN,Y(I))
        YMAX=MAX(YMAX,Y(I))
      END DO
      IF (A0.GT.0) YMAX=MAX(YMAX,A0**(-2))
      IF (ILOG.EQ.1) THEN
        YMIN=LOG(YMIN)
        YMAX=LOG(YMAX)
      ELSE
        YMIN=0
      END IF
C
C FILL IN FITTED FUNCTION.
C
      NI=4*NX
      DX=(XMAX-XMIN)/NI
      DO I=0,NI
        XI=XMIN+I*DX   
        IX=NINT(((XI-XMIN)/(XMAX-XMIN))*NX)
        IF (A0.GT.0.AND.A1.GE.0.AND.A2.GE.0) THEN
          YI=A0**(-2)*EXP(-2*A1*XI+2*(A2*XI)**2)
          IF (ILOG.EQ.1) YI=LOG(YI)
          IY=NINT(((YMAX-YI)/(YMAX-YMIN))*NY)
          IF (IX.GE.0.AND.IX.LE.NX.AND.IY.GE.0.AND.IY.LE.NY)
     &     AA(IX,IY)='*'
        END IF
      END DO
C
C FILL IN SPLINE-INTERPOLATED DATA CURVE.
C
      DX=(X(N)-X(1))/NI
      DO I=0,NI
        XI=X(1)+I*DX   
        IX=NINT(((XI-XMIN)/(XMAX-XMIN))*NX)
        CALL SPLINT (N,X,Y,YPP,XI,YI)
        IF (ILOG.EQ.1) THEN
          IF (YI.LE.0) THEN
            YI=YMIN
          ELSE
            YI=LOG(YI)
          END IF
        END IF
        IY=NINT(((YMAX-YI)/(YMAX-YMIN))*NY)
        IF (IX.GE.0.AND.IX.LE.NX.AND.IY.GE.0.AND.IY.LE.NY)
     &   AA(IX,IY)='O'
      END DO
C
C PRINT THE PLOT ARRAY.
C
      DY=(YMAX-YMIN)/NY
      DO IY=0,NY
        IF (MOD(IY,10).EQ.0) THEN
          WRITE (ILP,100) YMAX-IY*DY,(AA(IX,IY),IX=0,NX)
        ELSE
          WRITE (ILP,101)            (AA(IX,IY),IX=0,NX)
        END IF
      END DO
      DX=(XMAX-XMIN)/10
      WRITE (ILP,102) (XMIN+IX*DX,IX=0,10)
 100  FORMAT (1X,E10.3,101A1)
 101  FORMAT (1X,10X,  101A1)
 102  FORMAT (1X,1X,11F10.4)
      RETURN
      END
C-----------------------------------------------------------------------
