      PROGRAM ROGERS
C
C ROBERT H. BLESSING
C HAUPTMAN-WOODWARD INSTITUTE
C 73 HIGH STREET
C BUFFALO, NEW YORK 14203, USA
C TELEPHONE:  (716) 856-9600, EXTENSION 335
C ELECTRONIC MAIL:  Blessing@HWI.Buffalo.edu
C
C JULY 1986,..., OCTOBER 1995
C
C ESITMATES THE ABSOLUTE SCALE FACTOR AND OVERALL ANISOTROPIC MEAN-
C SQUARE ATOMIC DISPLACEMENT PARAMETERS BY MEANS OF AN ANALYSIS OF THE
C PATTERSON ORIGIN PEAK.  THE METHOD IS BASED ON AN IDEA DUE TO DONALD
C ROGERS (1965) AS DEVELOPED AND IMPLEMENTED BY R. H. BLESSING AND D. A.
C LANGS, ACTA CRYST. A44, 729-735 (1988).
C
      PARAMETER NMAX=926
      DIMENSION P(NMAX),U(NMAX),V(NMAX),W(NMAX),TEMP(NMAX)
      DIMENSION CELL(6),G(3,3),GINV(3,3),PP(3,3),B(3,3),UU(3,3)
      CHARACTER ATIME*8,ADATE*9,TITLE*80,FILE*40,FORM*11
      CHARACTER SYST*12,LATT*1,PTGP*5,LAUE*5,XRAY*2,EL*2,ELANOM*2
      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)
      CALL TIME (ATIME)
      CALL DATE (ADATE)
c      ATIME='hh:mm:ss'
c      ADATE='dd-mmm-yy'
      IO1=10
      IO2=20
      ILP=60
C
C READ CONTROL DATA.
C
      OPEN (UNIT=IO1,FILE='rogers.dat',STATUS='OLD')
      READ (IO1,'(A)') TITLE
      READ (IO1,*)     ITYPE
      READ (IO1,'(A)') FILE
      READ (IO1,'(A)') LATT
      READ (IO1,'(A)') PTGP
      READ (IO1,*)     CELL
      READ (IO1,*)     ZCELL
      READ (IO1,'(A)') XRAY
      READ (IO1,*)     NEL
      IF (ZCELL.EQ.0) ZCELL=1
      F000=0
      RMCELL=0
      SUMZSQ=0
      DO I=1,NEL
        READ (IO1,'(A2,F10.0)') EL(I),XI(I)
        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)
        SUMZSQ=SUMZSQ+(XI(I)*ZI(I))**2
      END DO
C
C READ ANOMALOUS DISPERSION CORRECTIONS FOR WAVELENGTHS OTHER THAN CU OR
C MO K-ALPHA.
C
      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,'(A2,2F10.0)') 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
      DATA CUTOFF,SMIN,SMAX,FRACTION /4*0/
      READ (IO1,*,END=6) CUTOFF,SMIN,SMAX,FRACTION
 6    CLOSE (UNIT=IO1)
      IF (CUTOFF.EQ.0) CUTOFF=-3
      IF (FRACTION.EQ.0) FRACTION=0.5
      CALL DECODE (PTGP,LAUE,SYST,LATT,MLATT,IPTGP,ILAUE)
      CALL METRIC (CELL,VCELL,G,GINV)
      DX=(RMCELL/VCELL)/0.60221
C
C PRINT CONTROL DATA.
C
      OPEN (UNIT=ILP,FILE='rogers.lp',STATUS='NEW',
     & CARRIAGECONTROL='FORTRAN')
      WRITE (ILP,600) ATIME,ADATE,TITLE
      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
      WRITE (ILP,607) INT(ZCELL),(EL(I),XI(I)/ZCELL,I=1,NEL)
      WRITE (ILP,608) NINT(F000),RMCELL/ZCELL,DX
      WRITE (ILP,609) XRAY,(EL(I),FP(I),FPP(I),I=1,NEL)
      WRITE (ILP,620) CUTOFF,FRACTION
      IF (SMAX.NE.0) WRITE (ILP,621) SMAX
      IF (SMAX.EQ.0) SMAX=9
      IF (SMIN.NE.0) WRITE (ILP,622) SMIN
 600  FORMAT ('1'/'0PROGRAM ROGERS.  ',A,1X,A,'.  ',A)
 601  FORMAT ('0INPUT REFLECTION FILE:   ',A)
 602  FORMAT ('0',A/'0',A,'-LATTICE'/'0POINT GROUP ',A/'0LAUE GROUP ',A)
 603  FORMAT ('0NONCENTROSYMMETRIC STRUCTURE')
 604  FORMAT ('0CENTROSYMMETRIC STRUCTURE')
 606  FORMAT ('0',9X,'A',9X,'B',9X,'C',5X,'ALPHA',6X,'BETA',5X,'GAMMA',
     &14X,'V'/' ',3F10.4,3F10.3,F15.2)
 607  FORMAT ('0UNIT CELL CONTENTS:'/'0  Z = ',I2,' CRYSTAL CHEMICAL UNI
     &TS PER UNIT CELL'/'0  CRYSTAL CHEMICAL UNIT:'//('   ',A2,F10.2))
 608  FORMAT ('0F000     = ',I9,' ELECTRONS'/'0MOL. WT. = ',F9.3,' DALTO
     &NS'/'0D(X-RAY) = ',F9.3,' MG MM**-3 (MILLIGRAMS PER CUBIC MILLIMET
     &ER)')
 609  FORMAT (/1X,'RADIATION:  ',A2//1X,'ANOMALOUS DISPERSION CORRECTION
     & TERMS:'/1X,'--------------------------------------'/1X,
     &'EL        FP       FPP'/1X,'--        --       ---'/
     &(1X,A2,2F10.3))
 620  FORMAT ('0THRESHOLD FOR INCLUDING REFLECTION IN PATTERSON SUMMATIO
     &N'/' CUTOFF = ',F5.2/' I.E., FSQ .LT. CUTOFF*SIGMA(FSQ) EXCLUDED'/
     &'0FRACTIONAL SHIFT MULTIPLIER FOR LEAST-SQUARES REFINEMENT'/' FRAC
     &TION = ',F5.3)
 621  FORMAT ('0MAXIMUM SIN(THETA)/LAMBDA FOR INCLUDING REFLECTION IN PA
     &TTERSON SUMMATION'/' SMAX = ',F7.4,' RECIPROCAL ANGSTROM')
 622  FORMAT ('0MINIMUM SIN(THETA)/LAMBDA FOR INCLUDING REFLECTION IN PA
     &TTERSON SUMMATION'/' SMIN = ',F7.4,' RECIPROCAL ANGSTROM')
C
C PERFORM ANALYSIS.
C
      IF (ITYPE.LT.3) THEN
        FORM='UNFORMATTED'
      ELSE
        FORM='FORMATTED'
      END IF
      OPEN (UNIT=IO1,FILE=FILE,STATUS='OLD',FORM=FORM)
      OPEN (UNIT=IO2,STATUS='SCRATCH',FORM='UNFORMATTED')
      CALL DATA (ITYPE,IO1,IO2,CUTOFF,GINV,SMIN,SMAX,MDATA,NDATA,SLIMIT,
     & NEL,XI,A1,B1,A2,B2,A3,B3,A4,B4,C0,FP,FPP,ILAUE)
      CLOSE (UNIT=IO1,STATUS='KEEP')
      WRITE (ILP,691) MDATA,NDATA,SLIMIT,1/(2*SLIMIT)
      CALL PRELIM (IO2,CELL,VCELL,GINV,N,P,PMIN,P0,PP,RADIUS)
C
C PLOT THE AXIAL PATTERSON DENSITIES.
C
      WRITE (ILP,600) ATIME,ADATE,TITLE
      CALL PLOT1 (N,P,ILP)
      WRITE (ILP,'(/1X,''P(U,V,W) ALONG [100], [010], AND [001] FOR 0 .L
     &E. U, V, W .LE. 0.5.''//1X,''P(U,V,W) = (2/VCELL)*SUM(H,K,L) [FSQ(
     &MEAS)/SUM(J) F(J)**2] COS [2*PI*(H*U + K*V + L*W)]'')')
      CALL GRID (RADIUS,GINV,U,V,W,N)
      CALL PSUM (IO2,VCELL,N,U,V,W,P,PMIN,P0)
      DO I=1,N
        TEMP(I)=P(I)
      END DO
      CALL PADJ (N,P,PMIN,P0)
      CALL PFIT (N,U,V,W,P,PMIN,P0,PP,VCELL,F000,SUMZSQ,FRACTION,NCYCLE,
     & R)
      DO I=1,N
        P(I)=TEMP(I)
      END DO
C
C PLOT THE RADIAL PROJECTION OF THE ORIGIN PEAK AND FITTED GAUSSIAN.
C
      WRITE (ILP,600) ATIME,ADATE,TITLE
      CALL PLOT2 (N,U,V,W,P,PMIN,P0,PP,RADIUS,G,ILP)
      WRITE (ILP,'(/1X,''RADIAL PROJECTION OF THE PATTERSON ORIGIN PE'',
     &''AK "O" SUPERIMPOSED ON THE FITTED GAUSSIAN "*"'')')
      CALL SOLVE (P0,PP,VCELL,G,GINV,SCALEK,B,UU,BISO,UISO)
C
C PRINT RESULTS.
C
      WRITE (ILP,600) ATIME,ADATE,TITLE
      WRITE (ILP,692) RADIUS,N,NCYCLE,R
      WRITE (ILP,693) SCALEK,B(1,1),B(2,2),B(3,3),B(1,2),B(1,3),B(2,3),
     & UU(1,1),UU(2,2),UU(3,3),UU(1,2),UU(1,3),UU(2,3),UISO,BISO
 691  FORMAT ('0NUMBER OF REFLECTIONS IN THE INPUT DATA FILE'/' MDATA ',
     &'= ',I6/'0NUMBER OF UNIQUE REFLECTION DATA USED TO CALCULATE THE P
     &ATTERSON MAP'/' NDATA = ',I6/'0DATA RESOLUTION:'/'   MAXIMUM SIN(T
     &HETA)/LAMBDA    SMAX = ',F6.3,' RECIPROCAL ANGSTROMS'/'   MINIMU',
     &'M D-SPACING            DMIN = ',F6.3,' ANGSTROMS')
 692  FORMAT ('0RESULTS:'/' -------'/'0RADIUS OF FIT OF THE TRIVARIATE G
     &AUSSIAN TO THE PATTERSON ORIGIN PEAK'/' RADIUS = ',F5.2,' ANGSTROM
     &'/'0NUMBER OF GRID POINTS CALCULATED'/' NPOINTS = ',I6/'0STATISTIC
     &S-OF-FIT:'/'0R = SQRT(SUM(WT*(POBS - PCALC)**2)/SUM(WT*POBS**2))'/
     &'0NCYCLE = ',I10/' R      = ',E10.3)
 693  FORMAT ('0SCALEK, B(I,J), AND U(I,J) ARE DEFINED ACCORDING TO'/
     &'0  F(ABSOLUTE) = SCALEK*F(RELATIVE)'/'0AND'/'0  f(T) = f(0)*EXP(-
     &SUM(I) SUM(J) H(I)*H(J)*B(I,J))'/'0WHERE'/'0  B(I,J) = 2*PI**2*AST
     &AR(I)*ASTAR(J)*U(I,J))'/'0      SCALEK'/' ',E12.4/
     &'0      B(1,1)      B(2,2)      B(2,3)      B(1,2)      B(1,3)',
     &'      B(2,3)'/' ',6E12.4/
     &'0      U(1,1)      U(2,2)      U(3,3)      U(1,2)      U(1,3)',
     &'      U(2,3)'/' ',6E12.4/
     &'0EQUIVALENT ISOTROPIC TEMPERATURE FACTORS:'/
     &'0UISO = <U**2>         = ',E12.4,' SQUARED ANGSTROMS'/
     &' BISO = 8*PI**2*<U**2> = ',E12.4,' SQUARED ANGSTROMS')
C
C WRITE OUTPUT FILE.
C
      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,705) SCALEK
      WRITE (IO1,705) B(1,1),B(1,2),B(1,3),B(2,2),B(2,3),B(3,3)
 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)
      CLOSE (UNIT=IO1,STATUS='KEEP')
      CLOSE (UNIT=IO2,STATUS='DELETE')
      STOP 'PROGRAM ROGERS FINIS!'
      END
C-----------------------------------------------------------------------
      SUBROUTINE DATA (ITYPE,IO1,IO2,CUTOFF,GINV,SMIN,SMAX,MDATA,
     & NDATA,SLIMIT,NEL,XI,A1,B1,A2,B2,A3,B3,A4,B4,C,FP,FPP,ILAUE)
C
C     EXPAND DATA TO HEMISPHERE CORRESPONDING TO -1 SYMMETRY.
C
      DIMENSION GINV(3,3),XI(1),A1(1),B1(1),A2(1),B2(1),A3(1),B3(1),
     & A4(1),B4(1),C(1),FP(1),FPP(1)
      INTEGER H
      N=IO2
      MDATA=0
      NDATA=0
      SLIMIT=0
      IEND=0
  100 CALL READ1 (ITYPE,IO1,IEND,H,K,L,FSQ,SIGFSQ)
      IF (IEND.NE.0) GO TO 99
      MDATA=MDATA+1
      IF (FSQ.LT.CUTOFF*SIGFSQ) GO TO 100
      S=SINTHL(H,K,L,GINV)
      SLIMIT=MAX(SLIMIT,S)
      IF (S.LT.SMIN.OR.S.GT.SMAX) GO TO 100
      NDATA=NDATA+1
      CALL FCALC (NEL,XI,A1,B1,A2,B2,A3,B3,A4,B4,C,FP,FPP,S,SUMFSQ)
      FSQ=FSQ/SUMFSQ
      GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13,14) ILAUE
C
C     -1
C
    1 WRITE (N) H,K,L,FSQ
      GO TO 100
C
C     2/M
C
    2 IF (K.EQ.0) FSQ=0.5*FSQ
      WRITE (N) H,K,L,FSQ
      WRITE (N) H,-K,L,FSQ
      GO TO 100
C
C     MMM
C
    3 IF (H.EQ.0) FSQ=0.5*FSQ
      IF (K.EQ.0) FSQ=0.5*FSQ
      IF (L.EQ.0) FSQ=0.5*FSQ
      WRITE (N) H,K,L,FSQ
      WRITE (N) -H,K,L,FSQ
      WRITE (N) H,-K,L,FSQ
      WRITE (N) -H,-K,L,FSQ
      GO TO 100
C
C     4/M
C
    4 IF (H.EQ.0) FSQ=0.5*FSQ
      IF (K.EQ.0) FSQ=0.5*FSQ
      IF (L.EQ.0) FSQ=0.5*FSQ
      WRITE (N) H,K,L,FSQ
      WRITE (N) K,-H,L,FSQ
      WRITE (N) -H,-K,L,FSQ
      WRITE (N) -K,H,L,FSQ
      GO TO 100
C
C     4/MMM
C
    5 IF (H.EQ.0) FSQ=0.5*FSQ
      IF (K.EQ.0) FSQ=0.5*FSQ
      IF (L.EQ.0) FSQ=0.5*FSQ
      IF (K.EQ.H) FSQ=0.5*FSQ
      WRITE (N) H,K,L,FSQ
      WRITE (N) -H,K,L,FSQ
      WRITE (N) H,-K,L,FSQ
      WRITE (N) -H,-K,L,FSQ
      WRITE (N) K,H,L,FSQ
      WRITE (N) -K,H,L,FSQ
      WRITE (N) K,-H,L,FSQ
      WRITE (N) -K,-H,L,FSQ
      GO TO 100
C
C     -3
C
    6 IF (H.EQ.0) FSQ=0.5*FSQ
      IF (K.EQ.0) FSQ=0.5*FSQ
      IF (L.EQ.0) FSQ=0.5*FSQ
      I=-H-K
      WRITE (N) H,K,L,FSQ
      WRITE (N) I,H,L,FSQ
      WRITE (N) K,I,L,FSQ
      GO TO 100
C
C     -31M
C
    7 IF (H.EQ.0) FSQ=0.5*FSQ
      IF (K.EQ.0) FSQ=0.5*FSQ
      IF (L.EQ.0) FSQ=0.5*FSQ
      IF (K.EQ.H) FSQ=0.5*FSQ
      I=-H-K
      WRITE (N) H,K,L,FSQ
      WRITE (N) I,H,L,FSQ
      WRITE (N) K,I,L,FSQ
      WRITE (N) K,H,L,FSQ
      WRITE (N) H,I,L,FSQ
      WRITE (N) I,K,L,FSQ
      GO TO 100
C
C     -3M1
C
    8 IF (H.EQ.0) FSQ=0.5*FSQ
      IF (K.EQ.0) FSQ=0.5*FSQ
      IF (L.EQ.0) FSQ=0.5*FSQ
      IF (K.EQ.H) FSQ=0.5*FSQ
      I=-H-K
      WRITE (N) H,K,L,FSQ
      WRITE (N) I,H,L,FSQ
      WRITE (N) K,I,L,FSQ
      WRITE (N) -K,-H,L,FSQ
      WRITE (N) -H,-I,L,FSQ
      WRITE (N) -I,-K,L,FSQ
      GO TO 100
C
C     6/M
C
    9 IF (H.EQ.0) FSQ=0.5*FSQ
      IF (K.EQ.0) FSQ=0.5*FSQ
      IF (L.EQ.0) FSQ=0.5*FSQ
      I=-H-K
      WRITE (N) H,K,L,FSQ
      WRITE (N) I,H,L,FSQ
      WRITE (N) K,I,L,FSQ
      WRITE (N) -H,-K,L,FSQ
      WRITE (N) -I,-H,L,FSQ
      WRITE (N) -K,-I,L,FSQ
      GO TO 100
C
C     6/MMM
C
   10 IF (H.EQ.0) FSQ=0.5*FSQ
      IF (K.EQ.0) FSQ=0.5*FSQ
      IF (L.EQ.0) FSQ=0.5*FSQ
      IF (K.EQ.H) FSQ=0.5*FSQ
      I=-H-K
      WRITE (N) H,K,L,FSQ
      WRITE (N) I,H,L,FSQ
      WRITE (N) K,I,L,FSQ
      WRITE (N) -H,-K,L,FSQ
      WRITE (N) -I,-H,L,FSQ
      WRITE (N) -K,-I,L,FSQ
      WRITE (N) -K,-H,L,FSQ
      WRITE (N) -H,-I,L,FSQ
      WRITE (N) -I,-K,L,FSQ
      WRITE (N) K,H,L,FSQ
      WRITE (N) H,I,L,FSQ
      WRITE (N) I,K,L,FSQ
      GO TO 100
C
C     R -3
C
   11 STOP 'MUST TRANSFORM FROM RHOMB. TO HEX. AXES.'
C
C     R -3M
C
   12 STOP 'MUST TRANSFORM FROM RHOMB. TO HEX. AXES.'
C
C     M3
C
   13 IF (H.EQ.0) FSQ=0.5*FSQ
      IF (K.EQ.0) FSQ=0.5*FSQ
      IF (L.EQ.0) FSQ=0.5*FSQ
      IF (H.EQ.K) FSQ=0.5*FSQ
      IF (K.EQ.L) FSQ=0.5*FSQ
      IF (L.EQ.H) FSQ=0.5*FSQ
      WRITE (N) H,K,L,FSQ
      WRITE (N) -H,K,L,FSQ
      WRITE (N) H,-K,L,FSQ
      WRITE (N) -H,-K,L,FSQ
      WRITE (N) L,H,K,FSQ
      WRITE (N) -L,H,K,FSQ
      WRITE (N) L,-H,K,FSQ
      WRITE (N) -L,-H,K,FSQ
      WRITE (N) K,L,H,FSQ
      WRITE (N) -K,L,H,FSQ
      WRITE (N) K,-L,H,FSQ
      WRITE (N) -K,-L,H,FSQ
      GO TO 100
C
C     M3M
C
   14 IF (H.EQ.0) FSQ=0.5*FSQ
      IF (K.EQ.0) FSQ=0.5*FSQ
      IF (L.EQ.0) FSQ=0.5*FSQ
      IF (H.EQ.K) FSQ=0.5*FSQ
      IF (K.EQ.L) FSQ=0.5*FSQ
      IF (L.EQ.H) FSQ=0.5*FSQ
      WRITE (N) H,K,L,FSQ
      WRITE (N) -H,K,L,FSQ
      WRITE (N) H,-K,L,FSQ
      WRITE (N) -H,-K,L,FSQ
      WRITE (N) L,H,K,FSQ
      WRITE (N) -L,H,K,FSQ
      WRITE (N) L,-H,K,FSQ
      WRITE (N) -L,-H,K,FSQ
      WRITE (N) K,L,H,FSQ
      WRITE (N) -K,L,H,FSQ
      WRITE (N) K,-L,H,FSQ
      WRITE (N) -K,-L,H,FSQ
      WRITE (N) L,K,H,FSQ
      WRITE (N) -L,K,H,FSQ
      WRITE (N) L,-K,H,FSQ
      WRITE (N) -L,-K,H,FSQ
      WRITE (N) H,L,K,FSQ
      WRITE (N) -H,L,K,FSQ
      WRITE (N) H,-L,K,FSQ
      WRITE (N) -H,-L,K,FSQ
      WRITE (N) K,H,L,FSQ
      WRITE (N) -K,H,L,FSQ
      WRITE (N) K,-H,L,FSQ
      WRITE (N) -K,-H,L,FSQ
      GO TO 100
   99 ENDFILE N
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE DECODE (PTGP,LAUE,SYST,LATT,MLATT,IPTGP,ILAUE)
C
C POINT GROUP, LAUE GROUP, CRYSTAL SYSTEM, AND LATTICE TYPE AND
C MULTIPLICITY
C
C THERE ARE TEN CASES OF ALTERNATIVE AXES FOR SEVEN OF THE 32
C CRYSTALLOGRAPHIC POINT GROUPS, SO A TOTAL OF 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 PTGP*5,LAUE*5,SYMBOL(42)*5,SYST*12,LATT*1
      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 '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  '
      IF (I.GE. 1) ILAUE=1
      IF (I.GE. 3) ILAUE=2
      IF (I.GE. 6) ILAUE=3
      IF (I.GE. 9) ILAUE=4
      IF (I.GE.12) ILAUE=5
      IF (I.GE.17) ILAUE=6
      IF (I.GE.19) ILAUE=7
      IF (I.EQ.20.OR.I.EQ.22.OR.I.EQ.24) ILAUE=8
      IF (I.GE.25) ILAUE=9
      IF (I.GE.28) ILAUE=10
      IF (I.GE.33) ILAUE=11
      IF (I.GE.35) ILAUE=12
      IF (I.GE.38) ILAUE=13
      IF (I.GE.40) ILAUE=14
C
C CRYSTAL SYSTEM
C
      IF (I.GE. 1) SYST='TRICLINIC'
      IF (I.GE. 3) SYST='MONOCLINIC'
      IF (I.GE. 6) SYST='ORTHORHOMBIC'
      IF (I.GE. 9) SYST='TETRAGONAL'
      IF (I.GE.17) SYST='TRIGONAL'
      IF (I.GE.25) SYST='HEXAGONAL'
      IF (I.GE.33) SYST='RHOMBOHEDRAL'
      IF (I.GE.38) SYST='CUBIC'
C
C LATTICE CENTERING MULTIPLICITY
C
      IF (LATT.EQ.'P') MLATT=1
      IF (LATT.EQ.'A') MLATT=2
      IF (LATT.EQ.'B') MLATT=2
      IF (LATT.EQ.'C') MLATT=2
      IF (LATT.EQ.'F') MLATT=4
      IF (LATT.EQ.'I') MLATT=2
      IF (LATT.EQ.'R') MLATT=1
      IF (LATT.EQ.'R'.AND.SYST.NE.'RHOMBOHEDRAL') MLATT=3
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE 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(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
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,98
        IF (EL.EQ.SYMBOL(I)) GO TO 1
      END DO
      STOP 'ELEMENT SYMBOL IS NOT IN THE TABLE (Z .LE. 98)'
 1    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 GRID (R,GINV,U,V,W,N)
C
C CALCULATE FRACTIONAL CRYSTAL COORDINATES OF A HEMISPHERICAL POLAR
C COORDINATE GRID.
C
C     W|   R        U = R SIN(TH) COS(PH)
C      |   /|       V = R SIN(TH) SIN(PH)
C      |TH/ |       W = R COS(TH)
C      | /  |
C      |/   |
C      _____|__V    
C     /\    |  
C    /  \   |
C   / PH \  |
C U/      \ |
C          \|
C
      DIMENSION G(3,3),GINV(3,3),T(3,3),X(3),Y(3),U(1),V(1),W(1)
      PI=ACOS(-1.0)
C
C A GRID WITH 25 POINTS ALONG RADII AT 30 DEGREE INTERVALS OF THETA AND
c PHI GIVES 1 + 25 + 25*12*3 = 926 GRID POINTS.
C
      NR=25
      NA=12
      DR=R/NR
      DA=2*PI/NA
      N=0
      DO I=0,NR
        N=N+1
        U(N)=0
        V(N)=0
        W(N)=I*DR
      END DO
      DO I=0,NA-1
        PH=I*DA
        SINPH=SIN(PH)
        COSPH=COS(PH)
      DO J=1,NA/4
        TH=J*DA
        SINTH=SIN(TH)
        COSTH=COS(TH)
      DO K=1,NR
        R=K*DR
        N=N+1
        U(N)=R*SINTH*COSPH
        V(N)=R*SINTH*SINPH
        W(N)=R*COSTH
      END DO
      END DO
      END DO
C
C TRANSFORM TO FRACTIONAL CRYSTAL COORDINATES FROM ORTHONORMAL
C COORDINATES WITH X ALONG A, Y PERPENDICULAR TO A IN THE AB PLANE, AND
c Z PERPENDICULAR TO THE AB PLANE.
C
      CALL MATINV3 (GINV,G,DETGI)
      DETG=1/DETGI
      T(1,1)=1/SQRT(G(1,1))
      T(1,2)=-G(1,2)/SQRT(DETG*G(1,1)*GINV(3,3))
      T(1,3)=GINV(3,1)/SQRT(GINV(3,3))
      T(2,1)=0
      T(2,2)=SQRT(G(1,1)/(DETG*GINV(3,3)))
      T(2,3)=GINV(3,2)/SQRT(GINV(3,3))
      T(3,1)=0
      T(3,2)=0
      T(3,3)=SQRT(GINV(3,3))
      DO I=1,N
        X(1)=U(I)
        X(2)=V(I)
        X(3)=W(I)
        DO J=1,3
          Y(J)=0
        DO K=1,3
          Y(J)=Y(J)+T(J,K)*X(K)
        END DO
        END DO
        U(I)=Y(1)
        V(I)=Y(2)
        W(I)=Y(3)
      END DO
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE MATINV (N,A,D)
C
C         PURPOSE
C           INVERT A SYMMETRIC MATRIX AND CALCULATE ITS DETERMINANT.
C
C         USAGE
C           CALL MATINV (N,A,D,DLOG)
C
C         DESCRIPTION OF PARAMETERS
C           A    - INPUT MATRIX, WHICH IS REPLACED BY ITS INVERSE
C           N    - DEGREE OF MATRIX (ORDER OF DETERMINANT)
C           D    - DETERMINANT OF EITHER THE INPUT MATRIX OR THE SCALED
C                  INPUT MATRIX
C           DLOG - THE NATURAL LOGARITHM OF THE DETERMINANT OF THE INPUT
C                  MATRIX
C
C           DLOG IS RETURNED IN CASE THE MAGNITUDE OF D = EXP(DLOG) IS
C           TOO LARGE.  THE CALLING PROGRAM SHOULD COMPARE LOG(D) AND
C           DLOG TO DETERMINE WHICH VALUE OF D WAS RETURNED.
C
C         SUBROUTINES AND FUNCTION SUBPROGRAM REQUIRED
C           NONE
C
C         COMMENT
C           IF N EXCEEDS 25, THE WORK VECTORS DIAG, IK, AND JK MUST BE
C           RE-DIMENSIONED TO A LENGTH OF N.
C
      PARAMETER (NMAX=25)
      DIMENSION A(N,N),DIAG(NMAX),IK(NMAX),JK(NMAX)
      DOUBLE PRECISION A,DIAG
C
C         SAVE SQUARE-ROOT DIAGONAL ELEMENTS OF MATRIX.
C
      DO 10 I=1,N
 10   DIAG(I)=SQRT(ABS(A(I,I)))
C
C         SCALE MATRIX TO AVOID ARITHMETIC OVERFLOW OR UNDERFLOW DURING
C         INVERSION.
C
      DO 11 I=1,N
      DO 11 J=I,N
      A(I,J)=A(I,J)/(DIAG(I)*DIAG(J))
 11   A(J,I)=A(I,J)
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
 21   DO 30 I=K,N
      DO 30 J=K,N
      IF (ABS(AMAX).GT.ABS(A(I,J))) GO TO 30
      AMAX=A(I,J)
      IK(K)=I
      JK(K)=J
 30   CONTINUE
C
C         INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN A(K,K).
C
      IF (AMAX.NE.0) GO TO 41
      D=0
      RETURN
 41   I=IK(K)
      IF (I-K) 21,51,43
 43   DO 50 J=1,N
      T=A(K,J)
      A(K,J)=A(I,J)
 50   A(I,J)=-T
 51   J=JK(K)
      IF (J-K) 21,61,53
 53   DO 60 I=1,N
      T=A(I,K)
      A(I,K)=A(I,J)
 60   A(I,J)=-T
C
C         ACCUMULATE ELEMENTS OF INVERSE MATRIX.
C
 61   DO 70 I=1,N
      IF (I.EQ.K) GO TO 70
      A(I,K)=-A(I,K)/AMAX
 70   CONTINUE
      DO 80 I=1,N
      DO 80 J=1,N
      IF (I.EQ.K) GO TO 80
      IF (J.EQ.K) GO TO 80
      A(I,J)=A(I,J)+A(I,K)*A(K,J)
 80   CONTINUE
      DO 90 J=1,N
      IF (J.EQ.K) GO TO 90
      A(K,J)=A(K,J)/AMAX
 90   CONTINUE
      A(K,K)=1/AMAX
 100  D=D*AMAX
C
C         RESTORE ORDERING OF MATRIX.
C
      DO 130 L=1,N
      K=N-L+1
      J=IK(K)
      IF (J.LE.K) GO TO 111
      DO 110 I=1,N
      T=A(I,K)
      A(I,K)=-A(I,J)
 110  A(I,J)=T
 111  I=JK(K)
      IF (I.LE.K) GO TO 130
      DO 120 J=1,N
      T=A(K,J)
      A(K,J)=-A(I,J)
 120  A(I,J)=T
 130  CONTINUE
C
C         SCALE INVERSE MATRIX.
C
      DO 140 I=1,N
      DO 140 J=I,N
      A(I,J)=A(I,J)/(DIAG(I)*DIAG(J))
 140  A(J,I)=A(I,J)
C
C         SCALE DETERMINANT (IF POSSIBLE).
C
      DLOG=0
      DO 150 I=1,N
 150  DLOG=DLOG+LOG(DIAG(I))
      DLOG=LOG(ABS(D))+2*DLOG
C
C         TEST AGAINST RANGE OF MACHINE-ALLOWED MAGNITUDES.
C
C         (QMIN = 2**(-128) = 0.29E-38 AND QMAX = 2**(+128)/2 = 1.7E+38
C         FOR VAX COMPUTERS.)
C
      QMIN=1E-37
      QMAX=1E+37
      IF (LOG(QMIN).LT.ABS(DLOG).AND.ABS(DLOG).LT.LOG(QMAX))
     1 D=(D/ABS(D))*EXP(DLOG)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE MATINV3 (A,B,D)
C
C FOR A SYMMETRIC 3 X 3 MATRIX A, B = A**(-1) AND D = DET(A).
C
      DIMENSION A(3,3),B(3,3)
      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 A  = CRYSTAL LATTICE PARAMETERS
C B  = RECIPROCAL LATTICE PARAMETERS
C GA = CRYSTAL SPACE METRIC MATRIX
C GB = RECIPROCAL SCACE METRIC MATRIX
C
      DIMENSION A(6),B(6),GA(3,3),GB(3,3)
      PI=3.141593
      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 PADJ (N,P,PMIN,P0)
C
C SCULPT A MINIMUM PROFILE ORIGIN PEAK.
C
      DIMENSION P(1)
      NR=25
      NA=12
      N=1
      Q=P(1)
      DO I=1,NR
        N=N+1
        P(N)=MIN(P(N),Q)
        Q=P(N)
      END DO
      DO I=0,NA-1
      DO J=1,NA/4
        Q=P(1)
      DO K=1,NR
        N=N+1
        P(N)=MIN(P(N),Q)
        Q=P(N)
      END DO
      END DO
      END DO
      DO I=1,N
        PMIN=MIN(PMIN,P(I))
      END DO
      P0=P(1)-PMIN
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE PFIT (N,U,V,W,P,PMIN,P0,PP,VCELL,F000,SUMZSQ,F,NCYCLE,
     & R)
C
C REFINE THE FIT OF THE TRIVARIATE GAUSSIAN DENSITY FUNCTION TO THE
C PATTERSON ORIGIN PEAK.
C
C   P(U,V,W) - PMIN = P0*EXP(-(P11*U**2 + P22*V**2 + P33*W**2
C                              + 2*P12*U*V + 2*P13*U*W + 2*P23*V*W))
C
C WHERE
C
C   PMIN = -F000**2/(SCALEK**2*VCELL*SUM(J) Z(J)**2)
C
C AND
C
C   SCALEK**2 = SQRT(DET(PIJ))/(SQRT(PI**3)*VCELL*P0)
C
      PARAMETER NMAX=926
      DIMENSION P(1),U(1),V(1),W(1),WT(NMAX),PP(3,3)
      DIMENSION AA(8,8),BB(8),XX(8),QQ(8,7),QQT(7,8),CC(7,7),DD(7),UU(7)
      DOUBLE PRECISION AA,BB,CC,DD
      NOBS=N
      NPAR=8
      NCON=1
      PI=3.141593
      CONST1=VCELL*SQRT(PI**3)
      CONST2=VCELL*SUMZSQ
C
C CALCULATE INITIAL GOODNESS-OF-FIT STATISTICS.
C
      CHISQ=0
      SUMSQ=0
      DO I=1,N
        CHISQ=CHISQ+(P(I)-PMIN-PUVW(U(I),V(I),W(I),P0,PP))**2
        SUMSQ=SUMSQ+P(I)**2
      END DO
      Z=SQRT(CHISQ/(NOBS-NPAR))
      R=SQRT(CHISQ/SUMSQ)
C
C ITERATE THROUGH LEAST-SQUARES LOOP.
C
      NCYCLE=0
 1    NCYCLE=NCYCLE+1
C
C EVALUATE THE DERIVATIVES, AND BUILD THE NORMAL MATRIX AND VECTOR.
C
      NOBS=0
      DO J=1,NPAR
      DO K=1,NPAR
        AA(J,K)=0
      END DO
        BB(J)=0
      END DO
      DO I=1,N
        CALL WEIGHT (P(I)-PMIN,P0,PP,U(I),V(I),W(I),Z,WT(I))
        IF (WT(I).NE.0) THEN
          UI=U(I)
          VI=V(I)
          WI=W(I)
          YI=PUVW(UI,VI,WI,P0,PP)
          XX(1)=-UI**2*YI
          XX(2)=-VI**2*YI
          XX(3)=-WI**2*YI
          XX(4)=-2*UI*VI*YI
          XX(5)=-2*UI*WI*YI
          XX(6)=-2*VI*WI*YI
          XX(7)=YI/P0
          XX(8)=1
          YI=P(I)-PMIN-YI
          NOBS=NOBS+1
          DO J=1,NPAR
          DO K=J,NPAR
            AA(J,K)=AA(J,K)+WT(I)*XX(J)*XX(K)
            AA(K,J)=AA(J,K)
          END DO
            BB(J)=BB(J)+WT(I)*XX(J)*YI
          END DO
        END IF
      END DO
C
C BUILD PMIN CONSTRAINT MATRIX QQ(8,7).
C
C CONSTRAINT MATRIX ELEMENT (I,J) IS THE DERIVATIVE OF UNCONSTRAINED
C PARAMETER (I) WITH RESPECT TO CONSTRAINED PARAMETER (J).
C
      DO I=1,7
      DO J=1,7
        QQ(I,J)=0
      END DO
        QQ(I,I)=1
      END DO
      DETP=PP(1,1)*PP(2,2)*PP(3,3)+PP(1,2)*PP(2,3)*PP(3,1)
     &      +PP(1,3)*PP(2,1)*PP(3,2)-PP(3,1)*PP(2,2)*PP(1,3)
     &        -PP(3,2)*PP(2,3)*PP(1,1)-PP(3,3)*PP(2,1)*PP(1,2)
      D1=(PP(2,2)*PP(3,3)-PP(2,3)**2)/(2*DETP)
      D2=(PP(1,1)*PP(3,3)-PP(1,3)**2)/(2*DETP)
      D3=(PP(1,1)*PP(2,2)-PP(1,2)**2)/(2*DETP)
      D4=(2*PP(1,3)*PP(2,3)-PP(1,2)*PP(3,3))/(2*DETP)
      D5=(2*PP(1,2)*PP(3,1)-PP(1,3)*PP(2,2))/(2*DETP)
      D6=(2*PP(2,1)*PP(3,1)-PP(2,3)*PP(1,1))/(2*DETP)
      T=PMIN
      SCLSQ=SQRT(DETP)/(CONST1*P0)
      PMIN=-F000**2/(SCLSQ*CONST2)
      QQ(8,1)=-0.5*PMIN*D1/DETP
      QQ(8,2)=-0.5*PMIN*D2/DETP
      QQ(8,3)=-0.5*PMIN*D3/DETP
      QQ(8,4)=-0.5*PMIN*D4/DETP
      QQ(8,5)=-0.5*PMIN*D5/DETP
      QQ(8,6)=-0.5*PMIN*D6/DETP
      QQ(8,7)=PMIN/P0
      PMIN=T
C
C APPLY CONSTRAINTS TO THE NORMAL MATRIX AND VECTOR.
C
      DO I=1,NPAR-NCON
      DO J=1,NPAR
        QQT(I,J)=QQ(J,I)
      END DO
      END DO
      DO I=1,NPAR-NCON
      DO J=1,NPAR-NCON
        CC(I,J)=0
        DO K=1,NPAR
        DO L=1,NPAR
          CC(I,J)=CC(I,J)+QQT(I,K)*AA(K,L)*QQ(L,J)
        END DO
        END DO
      END DO
      END DO
      DO I=1,NPAR-NCON
        DD(I)=0
        DO J=1,NPAR
          DD(I)=DD(I)+QQT(I,J)*BB(J)
        END DO
      END DO
C
C INVERT THE CONSTRAINED NORMAL MATRIX.
C
      CALL MATINV (NPAR-NCON,CC,DET)
      IF (DET.EQ.0) STOP 'SINGULAR NORMAL MATRIX'
C
C CALCULATE THE CONSTRAINED PARAMETER SHIFTS.
C
      DO I=1,NPAR-NCON
        UU(I)=0
      DO J=1,NPAR-NCON
        UU(I)=UU(I)+CC(I,J)*DD(J)
      END DO
      END DO
C
C CALCULATE AND APPLY THE UNCONSTRAINED PARAMETER SHIFTS.
C
      DO I=1,NPAR
        XX(I)=0
      DO J=1,NPAR-NCON
        XX(I)=XX(I)+QQ(I,J)*UU(J)
      END DO
      END DO
      PP(1,1)=PP(1,1)+F*XX(1)
      PP(2,2)=PP(2,2)+F*XX(2)
      PP(3,3)=PP(3,3)+F*XX(3)
      PP(1,2)=PP(1,2)+F*XX(4)
      PP(1,3)=PP(1,3)+F*XX(5)
      PP(2,3)=PP(2,3)+F*XX(6)
      PP(2,1)=PP(1,2)
      PP(3,1)=PP(1,3)
      PP(3,2)=PP(2,3)
      P0=P0+F*XX(7)
      PMIN=PMIN+F*XX(8)
C
C CALCULATE GOODNESS-OF-FIT STATISTICS.
C
      CHISQ=0
      SUMSQ=0
      SUMWT=0
      DO I=1,N
        CHISQ=CHISQ+WT(I)*(P(I)-PMIN-PUVW(U(I),V(I),W(I),P0,PP))**2
        SUMSQ=SUMSQ+WT(I)*P(I)**2
        SUMWT=SUMWT+WT(I)
      END DO
      Z1=Z
      R1=R
      Z=SQRT((CHISQ/SUMWT)*NOBS/(NOBS-(NPAR-NCON)))
      R=SQRT(CHISQ/SUMSQ)
C
C TEST FOR CONVERGENCE.
C
      T=(R-R1)/R1
      IF (T.GT.0) GO TO 9
      IF (T.GT.-1E-6) RETURN
      IF (NCYCLE.LE.50) GO TO 1
 9    CONTINUE
C
C DIVERGENT REFINEMENT.  BACK-SHIFT THE PARAMETERS.
C
      NCYCLE=NCYCLE-1
      Z=Z1
      R=R1
      PP(1,1)=PP(1,1)-F*XX(1)
      PP(2,2)=PP(2,2)-F*XX(2)
      PP(3,3)=PP(3,3)-F*XX(3)
      PP(1,2)=PP(1,2)-F*XX(4)
      PP(1,3)=PP(1,3)-F*XX(5)
      PP(2,3)=PP(2,3)-F*XX(6)
      PP(2,1)=PP(1,2)
      PP(3,1)=PP(1,3)
      PP(3,2)=PP(2,3)
      P0=P0-F*XX(7)
      PMIN=PMIN-F*XX(8)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE PLOT1 (N,P,ILP)
      PARAMETER (NX=100,NY=50)
      CHARACTER AA*1
      DIMENSION AA(0:NX,0:NY),P(1),Y(3)
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=0.5
      YMIN=+1E9
      YMAX=-1E9
      DO I=1,3*N
        YMIN=MIN(YMIN,P(I))
        YMAX=MAX(YMAX,P(I))
      END DO
C
C FILL IN DATA CURVE.
C
      DX=(XMAX-XMIN)/N
      DO I=1,N
        XI=XMIN+(I-1)*DX   
        IX=NINT(((XI-XMIN)/(XMAX-XMIN))*NX)
        Y(3)=P(I)
        Y(2)=P(I+N)
        Y(1)=P(I+2*N)
        DO J=1,3
          YI=Y(J)
          IY=NINT(((YMAX-YI)/(YMAX-YMIN))*NY)
          IF (IX.GE.0.AND.IX.LE.NX.AND.IY.GE.0.AND.IY.LE.NY) THEN
            IF (J.EQ.1) AA(IX,IY)='W'
            IF (J.EQ.2) AA(IX,IY)='V'
            IF (J.EQ.3) AA(IX,IY)='U'
          END IF
        END DO
      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.2)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE PLOT2 (N,U,V,W,P,PMIN,P0,PP,RADIUS,G,ILP)
      PARAMETER (NX=100,NY=50)
      CHARACTER AA*1
      DIMENSION AA(0:NX,0:NY)
      DIMENSION U(1),V(1),W(1),P(1),PP(3,3),G(3,3),R(3)
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=RADIUS
      YMIN=0
      YMAX=P(1)-PMIN
C
C FILL IN THE ORIGIN PEAK PROFILE AND FITTED GAUSSIAN.
C
      DO I=1,N
        R(1)=U(I)
        R(2)=V(I)
        R(3)=W(I)
        XI=0
        DO J=1,3
        DO K=1,3
          XI=XI+R(J)*R(K)*G(J,K)
        END DO
        END DO
        XI=SQRT(XI)
        IX=NINT(((XI-XMIN)/(XMAX-XMIN))*NX)
C
C FITTED GAUSSIAN
C
        YI=PUVW(U(I),V(I),W(I),P0,PP)
        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)='*'
C
C PEAK PROFILE
C
        YI=P(I)-PMIN
        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.3)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE PRELIM (IFILE,CELL,VCELL,GINV,N,P,PMIN,P0,PP,R)
C
C CALCULATE THE PATTERSON DENSITY ALONG THE A, B, AND C AXES OUT TO
C U, V, AND W = 0.5, AND FORM A FIRST APPROXIMATION TO A TRIVARIATE
C GAUSSIAN ORIGIN PEAK FUNCTION.
C
      DIMENSION CELL(6),GINV(3,3),P(1),SIGMA(3),QQ(3,3),PP(3,3)
C
C SUM THE AXIAL PATTERSON DENSITIES FROM 0 TO 0.5.
C
      N=200
      DO I=1,3*N
        P(I)=0
      END DO
      DELTA=0.5/N
      PI=3.141593
      CONST=2*PI*DELTA
      REWIND IFILE
 1    READ (IFILE,END=9) IH,IK,IL,FSQ
      FH=CONST*IH
      FK=CONST*IK
      FL=CONST*IL
      DO I=0,N-1
        U=FH*I
        V=FK*I
        W=FL*I
        IU=I+1
        IV=I+N+1
        IW=I+2*N+1
        P(IU)=P(IU)+FSQ*COS(U)
        P(IV)=P(IV)+FSQ*COS(V)
        P(IW)=P(IW)+FSQ*COS(W)
      END DO
      GO TO 1
 9    DO I=1,3*N
        P(I)=(2/VCELL)*P(I)
      END DO
C
C FIND THE MINIMUM DENSITY.
C
      PMIN=0
      DO I=1,3*N
        PMIN=MIN(PMIN,P(I))
      END DO
C
C FIND THE HALF-HEIGHT RADII (ANGSTROMS) OF THE ORIGIN PEAK.
C
      Q=0.5*(P(1)+PMIN)
      U=0
      V=0
      W=0
      DO I=2,N
        IU=I
        IV=I+N
        IW=I+2*N
        IF (U.EQ.0.AND.P(IU).LE.Q) U=IU
        IF (V.EQ.0.AND.P(IV).LE.Q) V=IV
        IF (W.EQ.0.AND.P(IW).LE.Q) W=IW
        IF (U*V*W.NE.0) GO TO 5
      END DO
 5    IU=U
      IV=V
      IW=W
      U=IU-(P(IU)-Q)/(P(IU)-P(IU-1))
      V=IV-(P(IV)-Q)/(P(IV)-P(IV-1))-N
      W=IW-(P(IW)-Q)/(P(IW)-P(IW-1))-2*N
      U=U*DELTA*CELL(1)
      V=V*DELTA*CELL(2)
      W=W*DELTA*CELL(3)
C
C APPROXIMATE THE ORIGIN PEAK FUNCTION, P(U) - PMIN = P0*EXP(-UT*PP*U).
C
      P0=P(1)-PMIN
C
C THE TRIVARIATE GAUSSIAN EXPONENTIAL ARGUMENT -U(J)*PP(I,J)*U(I)
C CORRESPONDS TO THE EXPONENTIAL AGRUMENT OF A UNIVARIATE NORMAL
C DISTRIBUTION FUNCTION,
C
C P(X) = 1/(SIGMA*SQRT(2*PI))*EXP(-(X - MU)**2/(2*SIGMA**2)),
C
C AND FROM THE HALF-HEIGHT RADIUS R,
C
C SIGMA = R/SQRT(-2*LN(0.5)).
C
      SIGMA(1)=U/1.177410
      SIGMA(2)=V/1.177410
      SIGMA(3)=W/1.177410
      DO I=1,3
      DO J=I,3
        QQ(I,J)=2*SIGMA(I)*SIGMA(J)*GINV(I,J)
        QQ(J,I)=QQ(I,J)
      END DO
      END DO
      CALL MATINV3 (QQ,PP,DETQQ)
C
C FOR A NORMAL GAUSSIAN, P(X) < 0.003 IF ABS(X - MU) > 3*SIGMA.
C
      R=3*MAX(SIGMA(1),SIGMA(2),SIGMA(3))
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE PSUM (IFILE,VCELL,N,U,V,W,P,PMIN,P0)
C
C CALCULATE THE PATTERSON DENSITY ON A PRE-CALCULATED HEMISPHERICAL
C POLAR COORDINATE GRID AROUND THE ORIGIN.
C
      DIMENSION P(1),U(1),V(1),W(1)
      TWOPI=2*3.141593
      DO I=1,N
        U(I)=TWOPI*U(I)
        V(I)=TWOPI*V(I)
        W(I)=TWOPI*W(I)
        P(I)=0
      END DO
      REWIND IFILE
 1    READ (IFILE,END=9) IH,IK,IL,FSQ
      DO I=1,N
        P(I)=P(I)+FSQ*COS(IH*U(I)+IK*V(I)+IL*W(I))
      END DO
      GO TO 1
 9    DO I=1,N
        U(I)=U(I)/TWOPI
        V(I)=V(I)/TWOPI
        W(I)=W(I)/TWOPI
        P(I)=(2/VCELL)*P(I)
      END DO
      RETURN
      END
C-----------------------------------------------------------------------
      FUNCTION PUVW(U,V,W,P0,PP)
      DIMENSION PP(3,3),UU(3)
      UU(1)=U
      UU(2)=V
      UU(3)=W
      Q=0
      DO I=1,3
      DO J=1,3
        Q=Q+UU(J)*PP(I,J)*UU(I)
      END DO
      END DO
      PUVW=P0*EXP(-Q)
      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-----------------------------------------------------------------------
      FUNCTION SINTHL (IH,IK,IL,GINV)
      DIMENSION GINV(3,3),H(3)
      H(1)=IH
      H(2)=IK
      H(3)=IL
      Q=0
      DO I=1,3
      DO J=1,3
        Q=Q+H(J)*GINV(I,J)*H(I)
      END DO
      END DO
      SINTHL=0.5*SQRT(Q)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE SOLVE (P0,P,VCELL,G,GINV,SCALEK,B,U,BISO,UISO)
C
C OBTAIN THE SCALE FACTOR, SCALEK, AND OVERALL ANISOTROPIC THERMAL
C VIBRATION TENSORS, B(I,J) AND U(I,J), DEFINED BY
C
C   F(ABSOLUTE) = SCALEK*F(RELATIVE)
C
C AND
C
C   F(T) = F(0)*EXP(-SUM(I) SUM(J) H(I)*H(J)*B(I,J))
C
C WHERE
C
C   B(I,J) = -2*PI**2*ASTAR(I)*ASTAR(J)*U(I,J))
C
      DIMENSION P(3,3),PINV(3,3),G(3,3),GINV(3,3),B(3,3),U(3,3)
      PI=3.141593
      CALL MATINV3 (P,PINV,DETP)
      SCALEK=SQRT(SQRT(DETP)/(SQRT(PI**3)*VCELL*P0))
      DO I=1,3
      DO J=I,3
        B(I,J)=PI**2*PINV(I,J)/2
        B(J,I)=B(I,J)
        U(I,J)=B(I,J)/(2*PI**2*SQRT(GINV(I,I)*GINV(J,J)))
        U(J,I)=U(I,J)
      END DO
      END DO
      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
      UISO=BISO/(8*PI**2)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE WEIGHT (P,P0,PP,U,V,W,Z,WT)
C
C DOWN-WEIGHT OFF-ORIGIN PEAKS.
C
      DIMENSION PP(3,3)
      IF (P.LT.0.5*P0.AND.P-PUVW(U,V,W,P0,PP).GT.Z) THEN
        WT=0.1
      ELSE
        WT=1
      END IF
      RETURN
      END
C-----------------------------------------------------------------------
