      PROGRAM REFPK
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
C     ELECTRONIC MAIL:  Blessing@MFB.Buffalo.Edu
C
C     JULY 1993
C
      PARAMETER (NMAX=100)
      CHARACTER ATIME*8,ADATE*9,TITLE*80
      COMMON /BLOCK0/ IO1,IO2,IO3,ILP
      COMMON /BLOCK1/ ATIME,ADATE,TITLE
      COMMON /BLOCK2/ GINV(3,3),WLE,WL1,WL0,WL2,T1,T2,SIGT1,SIGT2
      COMMON /BLOCK3/ WINDOW,CUTOFF,NSMOOTH,II1,II2,XT1,XT2,SC1,SC2
      COMMON /BLOCK4/ IDIFF,JMODEL,TANTHM,RHOM,F,U(3,3),MODEL1,MODEL2,
     & Q1(3,3),Q2(3,3)
      COMMON /BLOCK5/ NSTEP,XX(NMAX),YY(NMAX),VX(NMAX),VY(NMAX)
      COMMON /BLOCK6/ I1,I2,L1,L2,X1,X2,FLL
      COMMON /BLOCK7/ VARTHE,VARTIM,TAU,VARTAU,ATT,VARATT
      COMMON /BLOCKC/ NCOND,ICOND(38)
      CALL INPUT
      CALL DATA
      CALL TFIT
      CALL UQFIT
      CALL OUTPUT
      STOP 'PROGRAM REFPK FINIS!'
      END
C-----------------------------------------------------------------------
      SUBROUTINE ACCUM1 (N,X,Y,A,B)
C
C     ACCUMULATE LEAST-SQUARES NORMAL MATRIX AND VECTOR ELEMENTS.
C
      DIMENSION X(N),A(N,N),B(N)
      DOUBLE PRECISION A
      DO 1 I=1,N
      DO 2 J=I,N
      A(I,J)=A(I,J)+X(I)*X(J)
      A(J,I)=A(I,J)
    2 CONTINUE
      B(I)=B(I)+X(I)*Y
    1 CONTINUE
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE ACCUM2 (N,X,Y1,Y2,A,B1,B2)
C
C     ACCUMULATE LEAST-SQUARES NORMAL MATRIX AND VECTOR ELEMENTS.
C
      DIMENSION X(N),A(N,N),B1(N),B2(N)
      DOUBLE PRECISION A
      DO 1 I=1,N
      DO 2 J=I,N
      A(I,J)=A(I,J)+X(I)*X(J)
      A(J,I)=A(I,J)
    2 CONTINUE
      B1(I)=B1(I)+X(I)*Y1
      B2(I)=B2(I)+X(I)*Y2
    1 CONTINUE
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE AGREE0 (IFILE,U,NOBS,NPAR,CHISQ1,CHISQ2,CHISQ3,RMSD0,
     & RMSD1,RMSDX)
C
C     STATISTICS OF FIT FOR REFLECTION PEAK POSITIONS
C
      DIMENSION U(3,3),H(3),Y(3),Z(3)
      PI=ACOS(-1.0)
      RAD=PI/180
      DEG=180/PI
      SUMDX=0
      SUMD0=0
      SUMD1=0
      CHISQ1=0
      CHISQ2=0
      CHISQ3=0
      REWIND IFILE
      DO 5 N=1,NOBS
      READ (IFILE) II,H,TWOTH,OMEGA,CHI,PHI,DSTAR,Y,Z,
     & I1,I2,XM,X0,X1,X2,L1,L2,W1,W2,TANTH
      Y1=Y(1)
      Y2=Y(2)
      Y3=Y(3)
      DO 1 J=1,3
      Y(J)=0
      DO 1 K=1,3
      Y(J)=Y(J)+H(K)*U(K,J)
    1 CONTINUE
      DO 2 J=1,3
      Y(J)=Y(J)/DSTAR
    2 CONTINUE
      IF (TWOTH.LT.0) THEN
        DO 3 J=1,3
        Y(J)=-Y(J)
    3   CONTINUE
      END IF
      CHISQ1=CHISQ1+(Y1-Y(1))**2
      CHISQ2=CHISQ2+(Y2-Y(2))**2
      CHISQ3=CHISQ3+(Y3-Y(3))**2
      COSPH=COS(PHI*RAD)
      SINPH=SIN(PHI*RAD)
      COSCH=COS(CHI*RAD)
      SINCH=SIN(CHI*RAD)
      SINOM=Y(1)*COSPH-Y(2)*SINPH
      COSOM=(Y(1)*SINPH+Y(2)*COSPH)*COSCH-Y(3)*SINCH
      DELTA=ARCTAN(SINOM,COSOM)*DEG-OMEGA
      DELTA=A180(DELTA)
      IF (TWOTH.LT.0) DELTA=-DELTA
      SUMD0=SUMD0+DELTA**2
      X0CALC=XM-DELTA
      SUMD1=SUMD1+(X0-X0CALC)**2
      DX=(W1+X2-X1+W2)/(L2-L1+1)
      SUMDX=SUMDX+DX**2
    5 CONTINUE
      NFREE=NOBS-NPAR
      CHISQ1=CHISQ1/NFREE
      CHISQ2=CHISQ2/NFREE
      CHISQ3=CHISQ3/NFREE
      RMSDX=SQRT(SUMDX/NOBS)
      RMSD0=SQRT(SUMD0/NOBS)
      RMSD1=SQRT(SUMD1/(NOBS-3*NPAR))
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE AGREE1 (MODEL,IFILE,Q1,Q2,T1,T2,TANTHM,F,WL0,
     & NOBS,NPAR,CHISQ1,CHISQ2,RMSD1,RMSD2)
      CHARACTER MODEL*4
C
C     STATISTICS OF FIT FOR REFLECTION WIDTHS
C
      DIMENSION Q1(3,3),Q2(3,3),H(3),Y(3),Z(3)
      IF (NPAR.LT.6) THEN
        U1=Q1(1,1)
        U2=Q2(1,1)
      END IF
      CHISQ1=0
      CHISQ2=0
      SUMD1=0
      SUMD2=0
      REWIND IFILE
      DO 5 N=1,NOBS
      READ (IFILE) II,H,TWOTH,OMEGA,CHI,PHI,DSTAR,Y,Z,
     & I1,I2,XM,X0,X1,X2,L1,L2,W1,W2,TANTH
      IF (NPAR.GE.6) THEN
        U1=0
        U2=0
        DO 1 I=1,3
        DO 1 J=1,3
        U1=U1+Z(J)*Q1(I,J)*Z(I)
        U2=U2+Z(J)*Q2(I,J)*Z(I)
    1   CONTINUE
      END IF
      IF (MODEL.EQ.'LRNZ') THEN
        W1CALC=SQRT(U1)+T1*TANTH
        W2CALC=SQRT(U2)+T2*TANTH
      ELSE IF (MODEL.EQ.'GAUS') THEN
        W1CALC=SQRT(U1+(T1*TANTH)**2)
        W2CALC=SQRT(U2+(T2*TANTH)**2)
      ELSE IF (MODEL.EQ.'PLMC') THEN
        W1CALC=SQRT(U1-2*F*T1*TANTH/TANTHM+T1*(TANTH/TANTHM)**2)
        W2CALC=SQRT(U2-2*F*T2*TANTH/TANTHM+T2*(TANTH/TANTHM)**2)
      END IF
      IF (MODEL.EQ.'LRNZ') THEN
        CHISQ1=CHISQ1+(W1-W1CALC)**2
        CHISQ2=CHISQ2+(W2-W2CALC)**2
      ELSE
        CHISQ1=CHISQ1+(W1**2-W1CALC**2)**2
        CHISQ2=CHISQ2+(W2**2-W2CALC**2)**2
      END IF
      SUMD1=SUMD1+(W1-W1CALC)**2
      SUMD2=SUMD2+(W2-W2CALC)**2
    5 CONTINUE
      NFREE=NOBS-NPAR
      RMSD1=SQRT(SUMD1/NFREE)
      RMSD2=SQRT(SUMD2/NFREE)
      CHISQ1=CHISQ1/NFREE
      CHISQ2=CHISQ2/NFREE
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE CENTROID (WINDOW,CUTOFF,XX,YY,VX,VY,I1,I2,X0,SIGX0)
C
C     INTENSITY-WEIGHTED CENTROID OF PEAK ABOVE LINEAR BACKGROUND
C
      DIMENSION XX(1),YY(1),VX(1),VY(1)
      PARAMETER (NMAX=100)
      DIMENSION X(NMAX),Y(NMAX),V(NMAX)
      EQUIVALENCE (V(1),X(1))
C
C     FIND THE SMALLEST X AND LARGEST X FOR WHICH THE LOCAL AVERAGE Y(X)
C     INCREASES SIGNIFICANTLY ABOVE THE LOCAL AVERAGE BACKGROUND.
C
      N=0
      DO 11 I=I1,I2,+1
      N=N+1
      Y(N)=YY(I)
      V(N)=VY(I)
   11 CONTINUE
      CALL FINDL (N,Y,V,WINDOW,CUTOFF,L)
      IF (L.EQ.0) GO TO 9
      L1=I1+L
      N=0
      DO 12 I=I2,I1,-1
      N=N+1
      Y(N)=YY(I)
      V(N)=VY(I)
   12 CONTINUE
      CALL FINDL (N,Y,V,WINDOW,CUTOFF,L)
      IF (L.EQ.0) GO TO 9
      L2=I2-L
      IF (L1.GE.L2) GO TO 9
C
C     FIT A STRAIGHT LINE TO THE BACKGROUND NEAR THE PEAK LIMITS.
C
      L=NINT(WINDOW*N)
      N=0
      DO 1 I=L1-L,L2+L
      IF (I.LT.L1.OR.I.GT.L2) THEN
        N=N+1
        X(N)=XX(I)
        Y(N)=YY(I)
      END IF
    1 CONTINUE
      CALL FITLINE (N,X,Y,C0,C1,V0,V1,COV)
C
C     FIND THE INTENSITY-WEIGHTED CENTROID BETWEEN THE PEAK LIMITS.
C
      SUMY=0
      SUMYX=0
      SUMX2V=0
      SUMXV=0
      SUMV=0
      DO 2 I=L1,L2
      XI=XX(I)
      YI=YY(I)-C0-C1*XI
      VI=VY(I)+V0+V1*XI**2+2*COV*XI
      DX=0.50*(XX(I+1)-XX(I-1))
      VD=0.25*(VX(I+1)+VX(I-1))
      VI=YI**2*VD+DX**2*VI
      YI=YI*DX
      SUMY=SUMY+YI
      SUMYX=SUMYX+YI*XI
      SUMX2V=SUMX2V+XI**2*VI
      SUMXV=SUMXV+XI*VI
      SUMV=SUMV+VI
    2 CONTINUE
C
C     TEST FOR A SIGNIFICANTLY POSITIVE PEAK ABOVE BACKGROUND BETWEEN
C     THE TRUNCATED PEAK LIMITS.
C
      IF (SUMY.LT.CUTOFF*SQRT(SUMV)) GO TO 9
      X0=SUMYX/SUMY
      SIGX0=SQRT(ABS(SUMX2V-2*X0*SUMXV+X0**2*SUMV)/SUMY**2)
      RETURN
    9 CONTINUE
C
C     RETURN X0 = 90 DEGREES TO FLAG AN ERROR CONDITION.
C
      X0=90
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE DATA
C
C     SELECT REFLECTION PROFILES FOR THE LEAST-SQUARES CALCULATIONS OF
C     PEAK POSITION AND PEAK WIDTH PARAMETERS.
C
      PARAMETER (NMAX=100)
      CHARACTER ATIME*8,ADATE*9,TITLE*80
      COMMON /BLOCK0/ IO1,IO2,IO3,ILP
      COMMON /BLOCK1/ ATIME,ADATE,TITLE
      COMMON /BLOCK2/ GINV(3,3),WLE,WL1,WL0,WL2,T1,T2,SIGT1,SIGT2
      COMMON /BLOCK3/ WINDOW,CUTOFF,NSMOOTH,II1,II2,XT1,XT2,SC1,SC2
      COMMON /BLOCK4/ IDIFF,JMODEL,TANTHM,RHOM,F,U(3,3),MODEL1,MODEL2,
     & Q1(3,3),Q2(3,3)
      COMMON /BLOCK5/ NSTEP,XX(NMAX),YY(NMAX),VX(NMAX),VY(NMAX)
      COMMON /BLOCK6/ I1,I2,L1,L2,X1,X2,FLL
      COMMON /BLOCKC/ NCOND,ICOND(38)
      DIMENSION H(3),ANGLES(4),Y(3),Z(3)
      PI=ACOS(-1.0)
      DEG=180/PI
      RAD=PI/180
C
C     INITIALIZE COUNTERS FOR REFLECTIONS THAT FAIL THE ACCEPTANCE TESTS.
C
      N0=0
      N1=0
      N2=0
C
C     INITIALIZE MINIMUM AND MAXIMUM SERIAL NUMBERS AND EXPOSURE TIMES.
C
      IIMIN=+1E6
      IIMAX=-1E6
      XTMIN=+1E5
      XTMAX=-1E5
C
C     OPEN A SCRATCH FILE FOR LEAST-SQUARES DATA.
C
      OPEN (UNIT=IO3,FILE='SCRATCH',STATUS='SCRATCH',FORM='UNFORMATTED')
C
C     LOOP THROUGH ANALYSIS OF REFLECTION PROFILE DATA.
C
      IEND=0
    1 CALL READR (IO1,IEND,II,IH,IK,IL,ANGLES,NSTEP,XX,YY,VX,VY,XTIME)
      IF (IEND.NE.0) GO TO 9
C
C     SKIP REFLECTIONS THAT ARE NOT WITHIN THE SERIAL NUMBER OR EXPOSURE
C     TIME LIMITS.
C
      IF (ABS(II).LT.II1.OR.ABS(II).GT.II2) GO TO 1
      IF (XTIME.LT.XT1.OR.XTIME.GT.XT2) GO TO 1
C
C     TEST AGAINST CONDITIONS LIMITING POSSIBLE REFLECTIONS.
C
      CALL HKLTEST (IH,IK,IL,NCOND,ICOND,IABSENT)
      IF (IABSENT.NE.0) GO TO 1
      N0=N0+1
C
C     STORE MINIMUM AND MAXIMUM SERIAL NUMBERS AND EXPOSURE TIMES.
C
      IIMIN=MIN(IIMIN,ABS(II))
      IIMAX=MAX(IIMAX,ABS(II))
      XTMIN=MIN(XTMIN,XTIME)
      XTMAX=MAX(XTMAX,XTIME)
C
C     SMOOTH THE STATISTICAL FLUCTUATIONS IN THE SCAN PROFILE.
C
      CALL SMOOTH (NSMOOTH,NSTEP,YY)
      CALL SMOOTH (NSMOOTH,NSTEP,VY)
C
C     SET THE SCAN LIMITS (DEFAULT VALUES ARE SC1 = 0, SC2 = 1).
C
      DX=XX(NSTEP)-XX(1)
      I1=NINDEX(XX(1)+SC1*DX)
      I2=NINDEX(XX(1)+SC2*DX)
C
C     FIND THE INTENSITY-WEIGHTED CENTROID OF THE PEAK.
C
      CALL CENTROID (WINDOW,CUTOFF,XX,YY,VX,VY,I1,I2,X0,SIGX0)
      IF (X0.EQ.90) THEN
        N1=N1+1
        GO TO 1
      END IF
C
C     STORE THE INDICES IN A REAL ARRAY.
C
      H(1)=IH
      H(2)=IK
      H(3)=IL
C
C     CALCULATE SIN(THETA)/LAMBDA.
C
      S=SINTHL(H,GINV)
C
C     CALCULATE SPECTRAL THETA VALUES.
C
      THE=ASIN(S*WLE)*DEG
      TH1=ASIN(S*WL1)*DEG
      TH0=ASIN(S*WL0)*DEG
      TH2=ASIN(S*WL2)*DEG
C
C     CALCULATE THE POSITIONS OF THE THE ALPHA-ONE AND ALPHA-TWO PEAKS,
C     ASSUMING THE CENTROID TO BE THE ALPHA-BAR POSITION.
C
      X1=X0+(TH1-TH0)
      X2=X0+(TH2-TH0)
      CALL XCHECK (X1,X0,X2,XX,YY,I1,I2)
C
C     ALLOW FOR THE BETA FILTER ABSORPTION EDGE.
C
      J1=I1
      IF (WLE.NE.0) THEN
        XE=X0+(THE-TH0)
        I1=MAX(I1,NINDEX(XE))
      END IF
C
C     FIND THE PEAK LIMITS.
C
      CALL LEHLAR (L)
      IF (L.EQ.0) THEN
        N2=N2+1
        GO TO 1
      END IF
      I1=J1
C
C     CALCULATE ANGULAR BASE-WIDTHS OF HALF-PEAKS BELOW THETA(ALPHA-ONE)
C     AND ABOVE THETA(ALPHA-TWO).  (L1 AND L2 ARE THE PEAK LIMITS, I.E.,
C     STEPS L1 AND L2 ARE IN THE PEAK.)
C
      W1=X1-XX(L1)
      W2=XX(L2)-X2
C
C     TRANSFORM TO SETTING ANGLES AS DEFINED BY WALTER C. HAMILTON
C     (1974).  INTERNATIONAL TABLES FOR X-RAY CRYSTALLOGRAPHY, VOL. IV,
C     PP. 273-284.
C
      CALL GEOM (IDIFF,ANGLES,TWOTH,OMEGA,CHI,PHI)
C
C     CALCULATE THE DISPLACEMENT OF THE PEAK CENTROID FROM THE MIDPOINT
C     OF THE SCAN RANGE.
C
      XM=0.5*(XX(I1)+XX(I2))
      DELTA=X0-XM
C
C     GET UNIT DIRECTION VECTORS FOR THE DIFFRACTOMETER Y- AND Z-AXES.
C
      CALL YZ (DELTA,TWOTH,OMEGA,CHI,PHI,Y,Z)
C
C     CALCULATE DSTAR AND TAN(THETA).
C
      DSTAR=2*S
      TANTH=TAN(ASIN(S*WL0))
C
C     WRITE DATA TO A FILE OF LEAST-SQUARES OBSERVATIONS.
C
      WRITE (IO3) II,H,TWOTH,OMEGA,CHI,PHI,DSTAR,Y,Z,
     & I1,I2,XM,X0,X1,X2,L1,L2,W1,W2,TANTH
C
C     LOOP BACK FOR NEXT REFLECTION.
C
      GO TO 1
C
C     END LOOP THROUGH REFLECTION DATA.
C
    9 ENDFILE IO3
      CLOSE (UNIT=IO1)
C
C     WRITE DATA SET LIMITS TO BGLP.DAT FILE.
C
      WRITE (IO2,601) IIMIN,IIMAX,XTMIN-0.005,XTMAX+0.005
  601 FORMAT (1X,2I10,2F10.3)
C
C     PRINT SOME DATA SET STATISTICS.
C
      WRITE (ILP,610) ATIME,ADATE,TITLE
      WRITE (ILP,611) IIMIN,IIMAX,XTMIN-0.005,XTMAX+0.005
      WRITE (ILP,612) N0,N1,N2,N0-(N1+N2)
  610 FORMAT ('1'/'0PROGRAM REFPK.  ',A,1X,A,'.  ',A)
  611 FORMAT ('0BEGINNING AND ENDING REFLECTION SERIAL NUMBERS AND X-RAY
     & EXPOSURE TIMES (HOURS)'/'0     IIMIN     IIMAX               XTMI
     &N     XTMAX'/' '2I10,10X,2F10.2)
  612 FORMAT ('0STATISTICS ON REFLECTIONS PROCESSED:'/'0N0 = ',I5,'  TOT
     &AL REFLECTIONS (EXCLUDING SYMMETRY FORBIDDEN REFLECTIONS)'/' N1 ',
     &'= ',I5,'  REFLECTIONS DID NOT HAVE A SIGNIFICANT PEAK ABOVE LINEA
     &R BACKGROUND'/' N2 = ',I5,'  ADDITIONAL REFLECTIONS DID NOT GIVE W
     &ELL-DEFINED PEAK LIMITS'/' N0 - (N1 + N2) = ',I5,'  REFLECTIONS US
     &ED FOR LEAST-SQUARES CALCULATIONS OF PEAK POSITION AND WIDTH PARAM
     &ETERS')
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE DPRINT (II,IH,IK,IL,IPRINT)
C
C     IPRINT = 1 FOR H00, 0K0, 00L, HH0, 0KK, H0H, AND HHH.
C     IPRINT = 0 OTHERWISE.
C
      IF (II.LT.0) GO TO 10
      J=ABS(IH)
      K=ABS(IK)
      L=ABS(IL)
      IF (J.EQ.0.AND.K.EQ.0) GO TO 11
      IF (J.EQ.0.AND.L.EQ.0) GO TO 11
      IF (K.EQ.0.AND.L.EQ.0) GO TO 11
      IF (J.EQ.K.AND.L.EQ.0) GO TO 11
      IF (J.EQ.L.AND.K.EQ.0) GO TO 11
      IF (K.EQ.L.AND.J.EQ.0) GO TO 11
      IF (J.EQ.K.AND.K.EQ.L) GO TO 11
   10 IPRINT=0
      RETURN
   11 IPRINT=1
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE EIGEN (N,M,A,R)
C
C     COMPUTE EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC MATRIX
C
C     A - ORIGINAL MATRIX, DESTROYED IN COMPUTATION.
C         RESULTANT EIGENVALUES ARE DEVELOPED IN DIAGONAL OF MATRIX A IN
C         DESCENDING ORDER.
C     R - RESULTANT MATRIX OF EIGENVECTORS (STORED COLUMNWISE, IN SAME
C         SEQUENCE AS EIGENVALUES)
C     N - ORDER OF MATRICES A AND R
C       - LOGICAL SIZE OF ARRAYS A AND R
C     M - PHYSICAL SIZE OF ARRAYS A AND R, M .GE. N
C
C     E.G., FOR N = 3 SUBROUTINE EIGEN RETURNS
C
C         (W1 0  0 )
C     A = (0  W2 0 ) , WITH W1 .GE. W2 .GE. W3, AND
C         (0  0  W3)
C
C         (X1 X2 X3)
C     R = (Y1 Y2 Y3) .
C         (Z1 Z2 Z3)
C
C     THE DIAGONAL ELEMENTS OF A ARE THE EIGENVALUES, AND THE COLUMNS OF
C     R ARE THE CORRESPONDING EIGENVECTORS.
C
C     DIAGONALIZATION METHOD ORIGINATED BY JACOBI AND ADAPTED BY VON
C     NEUMANN FOR LARGE COMPUTERS AS FOUND IN 'MATHEMATICAL METHODS FOR
C     DIGITAL COMPUTERS', VOL. I, EDITED BY A. RALSTON AND H. S. WILF,
C     JOHN WILEY AND SONS, NEW YORK, 1962, CHAPTER 7.
C
      DIMENSION A(M,M),R(M,M)
C
C     GENERATE IDENTITY MATRIX AS FIRST APPROXIMATION TO THE MATRIX OF
C     EIGENVECTORS.
C
      DO 10 I=1,N
      R(I,I)=1
   10 CONTINUE
      DO 15 I=1,N-1
      DO 15 J=I+1,N
      R(I,J)=0
      R(J,I)=0
   15 CONTINUE
C
C     CALCULATE INITIAL AND FINAL OFF-DIAGONAL NORMS FOR SETTING AND
C     TESTING THRESHOLD.
C
      Q=0
      DO 20 I=1,N-1
      DO 20 J=I+1,N
      Q=Q+A(I,J)**2
   20 CONTINUE
      Q=SQRT(2*Q)
      E=1E-10
      IF (Q.LT.E) GO TO 41
      QF=Q*E/N
C
C     LOWER THRESHOLD AT EACH ITERATION.
C
      DO 40 ICYCLE=1,50
      Q=Q/N
C
C     LOOP THROUGH OFF-DIAGONAL PIVOT ELEMENTS.
C
      DO 30 I=1,N-1
      DO 30 J=I+1,N
      IF (ABS(A(I,J)).GT.Q) THEN
C
C       CALCULATE ROTATION TO ANNHILATE OFF-DIAGONAL ELEMENT A(I,J).
C
        U=-A(I,J)
        V=(A(I,I)-A(J,J))/2
        W=SIGN(1.0,V)*U/SQRT(U**2+V**2)
C
C       SINE AND COSINE OF ROTATION ANGLE
C
        S=W/SQRT(2*(1+SQRT(1-W**2)))
        C=SQRT(1-S**2)
        C2=C**2
        S2=S**2
        CS=C*S
C
C       TRANSFORM AND REPLACE THE ELEMENTS IN THE ITH AND JTH COLUMNS OF
C       THE A MATRIX AND OF THE R MATRIX.
C
        DO 35 K=1,N
C
C       SKIP ELEMENTS OF THE PIVOT SET A(I,J), (A(J,I)), A(I,I), A(J,J).
C
        IF (K.NE.I.AND.K.NE.J) THEN
          BKI=A(K,I)*C-A(K,J)*S
          BKJ=A(K,I)*S+A(K,J)*C
          A(K,I)=BKI
          A(K,J)=BKJ
          A(I,K)=A(K,I)
          A(J,K)=A(K,J)
        END IF
        BKI=R(K,I)*C-R(K,J)*S
        BKJ=R(K,I)*S+R(K,J)*C
        R(K,I)=BKI
        R(K,J)=BKJ
   35   CONTINUE
C
C       COMPUTE THE NEW OFF-DIAGONAL PIVOT ELEMENT AND ITS ASSOCIATED
C       DIAGONAL ELEMENTS.
C
        BII=A(I,I)*C2+A(J,J)*S2-2*A(I,J)*CS
        BJJ=A(I,I)*S2+A(J,J)*C2+2*A(I,J)*CS
        BIJ=(A(I,I)-A(J,J))*CS+A(I,J)*(C2-S2)
        A(I,I)=BII
        A(J,J)=BJJ
        A(I,J)=BIJ
        A(J,I)=A(I,J)
      END IF
   30 CONTINUE
      IF (Q.LT.QF) GO TO 41  
   40 CONTINUE
   41 CONTINUE
C
C     SORT ON DECREASING EIGENVALUE.
C
      DO 45 I=1,N-1
      DO 45 J=I+1,N
      IF (A(I,I).LT.A(J,J)) THEN
        T=A(I,I)
        A(I,I)=A(J,J)
        A(J,J)=T
        DO 50 K=1,N
        T=R(K,I)
        R(K,I)=R(K,J)
        R(K,J)=T
   50   CONTINUE
      END IF
   45 CONTINUE
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE EXTREM (IXMIN,IXMAX,IX)
      PARAMETER (N1=4,N2=16,N3=50)
      DIMENSION IXMIN(N1,N2,N3),IXMAX(N1,N2,N3),ITYPE(N1),IX(N2)
      DATA ITYPE /7, 9, 15, 16/
C
C         TABULATE IN RANKED ORDER THE N3 SMALLEST AND N3 LARGEST
C         VALUES OF THE VARIABLE OF TYPE I1.
C
C           I1 TYPE OF VARIABLE
C           I2 VALUE OF VARIABLE
C           I3 RANK IN LIST
C
      DO 10 I1=1,N1
            I2=ITYPE(I1)
      DO 11 I3=1,N3
      IF (IX(I2).GE.IXMIN(I1,I2,I3)) GO TO 11
      DO 12 J3=N3,I3+1,-1
      DO 12 I2=1,N2
 12   IXMIN(I1,I2,J3)=IXMIN(I1,I2,J3-1)
      DO 13 I2=1,N2
 13   IXMIN(I1,I2,I3)=IX(I2)
      GO TO 20
 11   CONTINUE
 20   CONTINUE
            I2=ITYPE(I1)
      DO 21 I3=1,N3
      IF (IX(I2).LE.IXMAX(I1,I2,I3)) GO TO 21
      DO 22 J3=N3,I3+1,-1
      DO 22 I2=1,N2
 22   IXMAX(I1,I2,J3)=IXMAX(I1,I2,J3-1)
      DO 23 I2=1,N2
 23   IXMAX(I1,I2,I3)=IX(I2)
      GO TO 10
 21   CONTINUE
 10   CONTINUE
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE FINDL (N,Y,V,WINDOW,CUTOFF,L)
C
C     FIND THE X AT WHICH THE LOCAL AVERAGE Y(X) FIRST INCREASES
C     SIGNIFICANTLY ABOVE THE LOCAL AVERAGE BACKGROUND.
C
      DIMENSION Y(N),V(N)
      L=NINT(WINDOW*N)
      DO 1 I=L+1,N-L
      YB=0
      VB=0
      YI=0
      VI=0
      DO 2 J=0,L
      YB=YB+Y(I-J)
      VB=VB+V(I-J)
      YI=YI+Y(I+J)
      VI=VI+V(I+J)
    2 CONTINUE
      IF (YI-YB.GT.CUTOFF*SQRT(VI+VB)) THEN
        L=I
        RETURN
      END IF
    1 CONTINUE
      L=0
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE FITLINE (N,X,Y,C0,C1,V0,V1,COV)
C
C     LEAST-SQUARES STRAIGHT LINE BACKGROUND
C
C     GIVE EQUAL WEIGHT TO ALL BACKGROUND POINTS BECAUSE:  (1) IT IS
C     EXPECTED THAT THE BACKGROUND SHOULD BE APPROXIMATELY CONSTANT, AND
C     THUS CONSTANT WEIGHTS SHOULD BE APPROPRIATE; AND (2) WEIGHTS
C     PROPORTIONAL TO THE RECIPROCALS OF THE VARIANCES OF THE
C     EXPERIMENTAL MEASUREMENTS, ALTHOUGH APPROPRIATE FOR NORMALLY
C     DISTRIBUTED DATA, ARE NOT APPROPRIATE FOR POISSON DISTRIBUTED
C     DATA.  WEIGHT = 1/VARIANCE WOULD BIAS THE FIT TOWARD THE
C     BACKGROUND POINTS WITH THE LOWEST COUNT RATES.  (SEE P. F. PRICE
C     (1979).  ACTA CRYST. A35, 57-60.)
C
      DIMENSION X(N),Y(N)
      SUMX=0
      SUMY=0
      SUMX2=0
      SUMXY=0
      DO 1 I=1,N
      XI=X(I)
      YI=Y(I)
      SUMX=SUMX+XI
      SUMY=SUMY+YI
      SUMX2=SUMX2+XI*XI
      SUMXY=SUMXY+XI*YI
    1 CONTINUE
      DET=N*SUMX2-SUMX*SUMX
      C0=(SUMX2*SUMY-SUMX*SUMXY)/DET
      C1=(-SUMX*SUMY+N*SUMXY)/DET
      V0=SUMX2/DET
      V1=N/DET
      COV=-SUMX/DET
      CHISQ=0
      DO 2 I=1,N
      CHISQ=CHISQ+(Y(I)-C0-C1*X(I))**2
    2 CONTINUE
      CHISQ=CHISQ/(N-2)
      V0=CHISQ*V0
      V1=CHISQ*V1
      COV=CHISQ*COV
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE GEOM (IDIFF,ANGLES,TWOTH,OMEGA,CHI,PHI)
C
C DIFFRACTOMETER GEOMETRIES
C
      DIMENSION ANGLES(4)
      REAL KAPPA
C
C KAPPA AXIS INCLINATION ANGLE, ALPHA = PI/3.6 = 50 DEGREES
C
      DATA SINA, COSA /0.766044, 0.642788/
      DATA DEG, RAD /57.2957795, 0.017453293/
      IF (IDIFF.EQ.1) THEN
C
C HAMILTON'S DIFFRACTOMETER AXES
C
C WALTER C. HAMILTON (1974).  INTERNATIONAL TABLES FOR X-RAY
C CRYSTALLOGRAPHY, VOL. IV, PP. 273-284.
C
        TWOTH = ANGLES(1)
        OMEGA = ANGLES(2)
        CHI   = ANGLES(3)
        PHI   = ANGLES(4)
      ELSE IF (IDIFF.EQ.2) THEN
C
C BUSING'S AND LEVY'S DIFFRACTOMETER AXES
C
C W. R. BUSING AND H. A. LEVY (1967).  ACTA CRYST. 22, 457-464.
C
        TWOTH = ANGLES(1)
        OMEGA = ANGLES(2)
        CHI   = ANGLES(3)
        PHI   = ANGLES(4)
C
C TRANSFORM TO SETTING ANGLES AS DEFINED BY HAMILTON.
C
        OMEGA = -OMEGA
        CHI   = -CHI
        PHI   = -PHI
      ELSE IF (IDIFF.EQ.3) THEN
C
C SIEMENS (NEE NICOLET, NEE SYNTEX) P3 DIFFRACTOMETER
C
        TWOTH = ANGLES(1)
        OMEGA = ANGLES(2)
        PHI   = ANGLES(3)
        CHI   = ANGLES(4)
C
C TRANSFORM TO SETTING ANGLES AS DEFINED BY HAMILTON.
C
        THETA =  TWOTH/2
        OMEGA =  OMEGA - THETA
        OMEGA = -OMEGA
      ELSE IF (IDIFF.EQ.4) THEN
C
C ENRAF-NONIUS CAD4 DIFFRACTOMETER KAPPA ANGLES
C
        THETA = ANGLES(1)
        PHI   = ANGLES(2)
        OMEGA = ANGLES(3)
        KAPPA = ANGLES(4)
C
C TRANSFORM FROM KAPPA TO EULERIAN SETTING ANGLES.
C
        KAPPA = KAPPA*RAD
        SINX  = SINA*SIN(KAPPA/2)
        COSX  = SQRT(COSA**2 + SINA**2*COS(KAPPA/2)**2)
        CHI   = 2*ARCTAN(SINX,COSX)
C
C WITH ALPHA = 50 DEGREES, KAPPA IS MECHANICALLY LIMITED TO THE RANGE 0
C TO 180 DEGREES.  THIS LIMITS CHI TO THE RANGE -100 TO +100 DEGREES.
C
        SINX  = COSA*SIN(KAPPA/2)/COS(CHI/2)
        COSX  = COS(KAPPA/2)/COS(CHI/2)
        DELTA = ARCTAN(SINX,COSX)*DEG
        OMEGA = OMEGA + DELTA
        PHI   = PHI   + DELTA
        CHI   = CHI*DEG
C
C TRANSFORM TO SETTING ANGLES AS DEFINED BY HAMILTON.
C
        TWOTH =  2*THETA
        OMEGA =  OMEGA - THETA
        OMEGA = -OMEGA
        CHI   = -CHI
        PHI   = -PHI
      END IF
C
C RETURN ANGLES -180 .LE. A .LE. +180.
C
      TWOTH = A180(TWOTH)
      OMEGA = A180(OMEGA)
      PHI   = A180(PHI)
      CHI   = A180(CHI)
      RETURN
      END
C-----------------------------------------------------------------------
      FUNCTION ARCTAN(SINX,COSX)
      IF (SINX.LT.-1) SINX=-1
      IF (SINX.GT.+1) SINX=+1
      IF (COSX.LT.-1) COSX=-1
      IF (COSX.GT.+1) COSX=+1
      ARCTAN=ATAN2(SINX,COSX)
      RETURN
      END
C-----------------------------------------------------------------------
      FUNCTION A180(A)
      A=MOD(A,360.0)
      IF (A.LT.-180) A=A+360
      IF (A.GT.+180) A=A-360
      A180=A
      RETURN
      END
C-----------------------------------------------------------------------
      FUNCTION A360(A)
      A=MOD(A,360.0)
      IF (A.LT.0) A=A+360
      A360=A
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE HKLCOND (H,I)
C
C     CONDITIONS LIMITING POSSIBLE REFLECTIONS
C
      CHARACTER*18 H,A
      DIMENSION A(38)
      DATA A
     &/'HKL,H+K=2N        ','HKL,K+L=2N        ','HKL,L+H=2N        ',
     & 'HKL,H+K,K+L,L+H=2N','HKL,H+K+L=2N      ','HKL,-H+K+L=3N     ',
     & 'HKL,H-K+L=3N      ','HKL,H-K=3N        ','HK0,H=2N          ',
     & 'HK0,K=2N          ','HK0,H+K=2N        ','HK0,H+K=4N        ',
     & '0KL,K=2N          ','0KL,L=2N          ','0KL,K+L=2N        ',
     & '0KL,K+L=4N        ','H0L,L=2N          ','H0L,H=2N          ',
     & 'H0L,L+H=2N        ','H0L,L+H=4N        ','HH(-2H)L,L=2N     ',
     & 'H(-H)0L,L=2N      ','HHL,L=2N(R-AXES)  ','HHL,L=2N          ',
     & 'HKH,K=2N          ','HKK,H=2N          ','HHL,2H+L=4N       ',
     & 'HKH,2H+K=4N       ','HKK,2K+H=4N       ','H00,H=2N          ',
     & 'H00,H=4N          ','0K0,K=2N          ','0K0,K=4N          ',
     & '00L,L=2N          ','00L,L=4N          ','000L,L=2N         ',
     & '000L,L=3N         ','000L,L=6N         '/
      IF (I.EQ.0) THEN
        DO 1 I=1,38
        IF (H.EQ.A(I)) RETURN
    1   CONTINUE
        STOP 'INPUT ERROR.  HKL CONDITION IS NOT IN THE TABLE.'
      ELSE
        H=A(I)
        RETURN
      END IF
      END
C-----------------------------------------------------------------------
      SUBROUTINE HKLTEST (H,K,L,NCOND,ICOND,IABSENT)
C
C     CONDITIONS LIMITING POSSIBLE REFLECTIONS
C
C     ( 1) HKL,H+K=2N
C     ( 2) HKL,K+L=2N
C     ( 3) HKL,L+H=2N
C     ( 4) HKL,H+K,K+L=2N
C     ( 5) HKL,H+K+L=2N
C     ( 6) HKL,-H+K+L=3N
C     ( 7) HKL,H-K+L=3N
C     ( 8) HKL,H-K=3N
C     ( 9) HK0,H=2N
C     (10) HK0,K=2N
C     (11) HK0,H+K=2N
C     (12) HK0,H+K=4N
C     (13) 0KL,K=2N
C     (14) 0KL,L=2N
C     (15) 0KL,K+L=2N
C     (16) 0KL,K+L=4N
C     (17) H0L,L=2N
C     (18) H0L,H=2N
C     (19) H0L,L+H=2N
C     (20) H0L,L+H=4N
C     (21) HH(-2H)L,L=2N
C     (22) H(-H)0L,L=2N
C     (23) HHL,L=2N (RHOMBOHEDRAL AXES)
C     (24) HHL,L=2N
C     (25) HKH,K=2N
C     (26) HKK,H=2N
C     (27) HHL,2H+L=4N
C     (28) HKH,2H+K=4N
C     (29) HKK,2K+H=4N
C     (30) H00,H=2N
C     (31) H00,H=4N
C     (32) 0K0,K=2N
C     (33) 0K0,K=4N
C     (34) 00L,L=2N
C     (35) 00L,L=4N
C     (36) 000L,L=2N
C     (37) 000L,L=3N
C     (38) 000L,L=6N
C
C     SET IABSENT = 0 FOR AN ALLOWED REFLECTION,
C                 = 1 FOR A FORBIDDEN REFLECTION.
C
      INTEGER H
      DIMENSION ICOND(38)
      IF (NCOND.EQ.0) THEN
        IABSENT=0
        RETURN
      END IF
      DO 90 I=1,NCOND
      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) ICOND(I)
    1 IF (MOD(H+K,2).NE.0) GO TO 91
      GO TO 90
    2 IF (MOD(K+L,2).NE.0) GO TO 91
      GO TO 90
    3 IF (MOD(L+H,2).NE.0) GO TO 91
      GO TO 90
    4 IF (MOD(H+K,2).NE.0.OR.MOD(K+L,2).NE.0) GO TO 91
      GO TO 90
    5 IF (MOD(H+K+L,2).NE.0) GO TO 91
      GO TO 90
    6 IF (MOD(-H+K+L,3).NE.0) GO TO 91
      GO TO 90
    7 IF (MOD(H-K+L,3).NE.0) GO TO 91
      GO TO 90
    8 IF (MOD(H-K,3).NE.0) GO TO 91
      GO TO 90
    9 IF (L.EQ.0.AND.MOD(H,2).NE.0) GO TO 91
      GO TO 90
   10 IF (L.EQ.0.AND.MOD(K,2).NE.0) GO TO 91
      GO TO 90
   11 IF (L.EQ.0.AND.MOD(H+K,2).NE.0) GO TO 91
      GO TO 90
   12 IF (L.EQ.0.AND.MOD(H+K,4).NE.0) GO TO 91
      GO TO 90
   13 IF (H.EQ.0.AND.MOD(K,2).NE.0) GO TO 91
      GO TO 90
   14 IF (H.EQ.0.AND.MOD(L,2).NE.0) GO TO 91
      GO TO 90
   15 IF (H.EQ.0.AND.MOD(K+L,2).NE.0) GO TO 91
      GO TO 90
   16 IF (H.EQ.0.AND.MOD(K+L,4).NE.0) GO TO 91
      GO TO 90
   17 IF (K.EQ.0.AND.MOD(L,2).NE.0) GO TO 91
      GO TO 90
   18 IF (K.EQ.0.AND.MOD(H,2).NE.0) GO TO 91
      GO TO 90
   19 IF (K.EQ.0.AND.MOD(L+H,2).NE.0) GO TO 91
      GO TO 90
   20 IF (K.EQ.0.AND.MOD(L+H,4).NE.0) GO TO 91
      GO TO 90
   21 IF (K.EQ.H.AND.MOD(L,2).NE.0) GO TO 91
      GO TO 90
   22 IF (K.EQ.-H.AND.MOD(L,2).NE.0) GO TO 91
      GO TO 90
   23 IF (K.EQ.H.AND.MOD(L,2).NE.0) GO TO 91
      GO TO 90
   24 IF (ABS(K).EQ.ABS(H).AND.MOD(L,2).NE.0) GO TO 91
      GO TO 90
   25 IF (ABS(L).EQ.ABS(H).AND.MOD(K,2).NE.0) GO TO 91
      GO TO 90
   26 IF (ABS(L).EQ.ABS(K).AND.MOD(H,2).NE.0) GO TO 91
      GO TO 90
   27 IF (ABS(K).EQ.ABS(H).AND.MOD(2*H+L,4).NE.0) GO TO 91
      GO TO 90
   28 IF (ABS(L).EQ.ABS(H).AND.MOD(2*H+K,4).NE.0) GO TO 91
      GO TO 90
   29 IF (ABS(L).EQ.ABS(K).AND.MOD(2*K+H,4).NE.0) GO TO 91
      GO TO 90
   30 IF (K.EQ.0.AND.L.EQ.0.AND.MOD(H,2).NE.0) GO TO 91
      GO TO 90
   31 IF (K.EQ.0.AND.L.EQ.0.AND.MOD(H,4).NE.0) GO TO 91
      GO TO 90
   32 IF (H.EQ.0.AND.L.EQ.0.AND.MOD(K,2).NE.0) GO TO 91
      GO TO 90
   33 IF (H.EQ.0.AND.L.EQ.0.AND.MOD(K,4).NE.0) GO TO 91
      GO TO 90
   34 IF (H.EQ.0.AND.K.EQ.0.AND.MOD(L,2).NE.0) GO TO 91
      GO TO 90
   35 IF (H.EQ.0.AND.K.EQ.0.AND.MOD(L,4).NE.0) GO TO 91
      GO TO 90
   36 IF (H.EQ.0.AND.K.EQ.0.AND.MOD(L,2).NE.0) GO TO 91
      GO TO 90
   37 IF (H.EQ.0.AND.K.EQ.0.AND.MOD(L,3).NE.0) GO TO 91
      GO TO 90
   38 IF (H.EQ.0.AND.K.EQ.0.AND.MOD(L,6).NE.0) GO TO 91
   90 CONTINUE
      IABSENT=0
      RETURN
   91 IABSENT=1
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE INPUT
C
C     READ CONTROL DATA, AND SET DEFAULT VALUES.
C
      PARAMETER (NMAX=100)
      CHARACTER ATIME*8,ADATE*9,TITLE*80,FILE1*80,FILE2*80,HCOND*18
      CHARACTER DIFF*4,SOURCE*2,TARGET*2,BETFIL*2,FILTER*2,MONOCR*2
      COMMON /BLOCK0/ IO1,IO2,IO3,ILP
      COMMON /BLOCK1/ ATIME,ADATE,TITLE
      COMMON /BLOCK2/ GINV(3,3),WLE,WL1,WL0,WL2,T1,T2,SIGT1,SIGT2
      COMMON /BLOCK3/ WINDOW,CUTOFF,NSMOOTH,II1,II2,XT1,XT2,SC1,SC2
      COMMON /BLOCK4/ IDIFF,JMODEL,TANTHM,RHOM,F,U(3,3),MODEL1,MODEL2,
     & Q1(3,3),Q2(3,3)
      COMMON /BLOCK5/ NSTEP,XX(NMAX),YY(NMAX),VX(NMAX),VY(NMAX)
      COMMON /BLOCK6/ I1,I2,L1,L2,X1,X2,FLL
      COMMON /BLOCK7/ VARTHE,VARTIM,TAU,VARTAU,ATT,VARATT
      COMMON /BLOCKC/ NCOND,ICOND(38)
      DIMENSION CELL(6)
      DIMENSION TARGET(3),WLA1(3),WLA2(3),DWLA1(3),DWLA2(3)
      DIMENSION BETFIL(5),ABSEDG(5)
      DIMENSION THMC(2,3),RHOMC(4)
      DIMENSION ER(2,4),CC(4,3,4)
      PI=ACOS(-1.0)
C
C     K-ALPHA X-RAY WAVELENGTHS
C
      DATA TARGET /   'CU',    'MO',    'AG'/
      DATA WLA1   /1.54056, 0.70930, 0.55941/
      DATA WLA2   /1.54439, 0.71359, 0.56380/
C
C     K-ALPHA SPECTRAL LINE-WIDTHS (FULL-WIDTHS AT HALF-HEIGHT)
C
      DATA DWLA1   /0.00058, 0.00029, 0.00028/
      DATA DWLA2   /0.00077, 0.00032, 0.00029/
C
C     K-BETA FILTER ABSORPTION EDGE WAVELENGTHS
C
      DATA BETFIL /   'NI',    'NB',    'ZR',    'PD',    'RH'/
      DATA ABSEDG /1.48807, 0.65298, 0.68883, 0.50920, 0.53395/
C
C     MONOCHROMATOR BRAGG ANGLES
C
C     GRAPHITE (0 0 0 2) REFLECTION
C     QUARTZ   (1 0-1 1) REFLECTION
C
C     GRAPHITE, P 6(3)/M M C, Z = 4, A = 2.455(2), C = 6.69 A,
C     D(0 0 0 2) = 3.345 A
C
C     QUARTZ, P 3(1,2) 2 1, Z = 3, A = 4.910(10), C = 5.394(10) A,
C     D(1 0-1 1) = 3.339 A
C
C            GRAPHITE QUARTZ
C
C     CU         TH11   TH21
C     MO         TH12   TH22
C     AG         TH13   TH23
C
      DATA THMC /13.3,  13.3,
     &            6.1,   6.1,
     &            4.8,   4.8/
C
C     MONOCHROMATOR-DIFFRACTOMETER GEOMETRY ANGLE
C
C     PARALLEL       RHO =  0 OR
C     PERPENDICULAR  RHO = 90 DEGREES
C
C                  INT.TAB.  BUS.LEV.       P3       CAD4
C
      DATA RHOMC /       0,        0,        0,        90/
C 
C     ERROR ESTIMATES FOR DIFFRACTOMETER SCAN ANGLE (DEGREES THETA) AND
C     SCAN-ANGLE DRIVE-PULSE TIMING (MICORSECONDS)
C
C             SIGTHE SIGTIM
C
C     INT.TAB.  ER11   ER21
C     BUS.LEV.  ER12   ER22
C     P3        ER13   ER23
C     CAD4      ER14   ER24
C
      DATA ER / 0.001,  1.0,
     &          0.001,  1.0,
     &          0.005, 10.0,
     &          0.005, 10.0/
C
C     COINCIDENCE CORRECTION DATA (DEAD TIME (MICROSECONDS) AND
C     ATTENUATOR FACTOR)
C
C                    TAU SIGTAU    ATT SIGATT
C
C     INT.TAB. CU  CC111  CC211  CC311  CC411
C              MO  CC121  CC221  CC321  CC421
C              AG  CC131  CC231  CC331  CC431
C     BUS.LEV. CU  CC112  CC212  CC312  CC412
C              MO  CC122  CC222  CC322  CC422
C              AG  CC132  CC232  CC332  CC432
C     P3       CU  CC113  CC213  CC313  CC413
C              MO  CC123  CC223  CC323  CC423
C              AG  CC133  CC233  CC333  CC433
C     CAD4     CU  CC114  CC214  CC314  CC414
C              MO  CC124  CC224  CC324  CC424
C              AG  CC134  CC234  CC334  CC434
C
      DATA CC /      0.0,   0.0,   1.0,   0.0,
     &               0.0,   0.0,   1.0,   0.0,
     &               0.0,   0.0,   1.0,   0.0,
     &               0.0,   0.0,   1.0,   0.0,
     &               0.0,   0.0,   1.0,   0.0,
     &               0.0,   0.0,   1.0,   0.0,
     &               2.5,   0.2,  31.4,   0.5,
     &               2.48,  0.12, 27.01,  0.41,
     &               0.0,   0.0,   1.0,   0.0,
     &               0.0,   0.0,  25.6,   0.3,
     &               0.0,   0.0,  21.5,   0.3,
     &               0.0,   0.0,   1.0,   0.0/
C
C     DUMMY CONTROL VARIABLES
C
      FILTER='  '
      MONOCR='  '
      THM=0
      RHOM=0
      FRACTD=0
      SIGFD=0
      F=0
      WL1=0
      WL2=0
      DWL1=0
      DWL2=0
      WLE=0
      SIGTHE=0
      SIGTIM=0
      TAU=0
      SIGTAU=0
      ATT=0
      SIGATT=0
      II1=0
      II2=0
      XT1=0
      XT2=0
      SC1=0
      SC2=0
      WINDOW=0
      CUTOFF=0
      FLL=0
      WT2=0
      NSMOOTH=0
      JMODEL=0
C
C     START.
C
      CALL TIME (ATIME)
      CALL DATE (ADATE)
C
C     I/O DEVICES.  LOGICAL UNIT NUMBERS
C
      IO1=1
      IO2=2
      IO3=3
      ILP=6
C
C     READ CONTROL DATA.
C
      OPEN (UNIT=IO1,FILE='refpk.dat',STATUS='OLD')
      READ (IO1,500) TITLE
      READ (IO1,500) FILE1
      READ (IO1,500) FILE2
      READ (IO1,501) CELL
      READ (IO1,505) DIFF
      READ (IO1,502) SOURCE
      READ (IO1,502,END=9) FILTER
      READ (IO1,502,END=9) MONOCR
      READ (IO1,501,END=9) THM,RHOM,FRACTD,SIGFD
      READ (IO1,501,END=9) WL1,WL2,DWL1,DWL2,WLE
      READ (IO1,501,END=9) SIGTHE,SIGTIM
      READ (IO1,501,END=9) TAU,SIGTAU,ATT,SIGATT
      READ (IO1,503,END=9) II1,II2,XT1,XT2
      READ (IO1,501,END=9) SC1,SC2
      READ (IO1,501,END=9) WINDOW
      READ (IO1,501,END=9) CUTOFF
      READ (IO1,501,END=9) FLL
      READ (IO1,501,END=9) WT2
      READ (IO1,504,END=9) NSMOOTH
      READ (IO1,504,END=9) JMODEL
  500 FORMAT (A)
  501 FORMAT (8F10.0)
  505 FORMAT (A4)
  502 FORMAT (A2)
  503 FORMAT (2I10,2F10.0)
  504 FORMAT (8I10)
    9 CLOSE (UNIT=IO1)
C
C     RECIPROCAL LATTICE METRIC TENSOR
C
      CALL METRIC (CELL,GINV)
C
C     DIFFRACTOMETER TYPE
C
      IDIFF=0
      IF (DIFF.EQ.'H   ') IDIFF=1
      IF (DIFF.EQ.'BL  ') IDIFF=2
      IF (DIFF.EQ.'P3  ') IDIFF=3
      IF (DIFF.EQ.'CAD4') IDIFF=4
      IF (IDIFF.EQ.0) STOP 'UNKNOWN DIFFRACTOMETER TYPE'
C
C     WAVELENGTHS AND SPECTRAL WIDTHS
C
      IXRAY=0
      INEUTRON=0
      DO 10 I=1,3
      IF (SOURCE.EQ.TARGET(I)) IXRAY=I
   10 CONTINUE
      IF (SOURCE.EQ.'N '.OR.SOURCE.EQ.'NE') INEUTRON=1
      IF (IXRAY.NE.0) THEN
        I=IXRAY
        IF (WL1.EQ.0) WL1=WLA1(I)
        IF (WL2.EQ.0) WL2=WLA2(I)
        IF (DWL1.EQ.0) DWL1=DWLA1(I)
        IF (DWL2.EQ.0) DWL2=DWLA2(I)
      END IF
      IF (WT2.EQ.0) WT2=0.5
      WL0=(WL1+WL2*WT2)/(1+WT2)
C
C     BETA FILTER ABSORPTION EDGE
C
      I=0
      DO 20 J=1,5
      IF (FILTER.EQ.BETFIL(J)) I=J
   20 CONTINUE
      IF (I.NE.0.AND.WLE.EQ.0) WLE=ABSEDG(I)
C
C     MONOCHROMATOR VARIABLES
C
      IF (MONOCR.NE.'  '.OR.THM.NE.0) THEN
C
C       BRAGG ANGLE
C
        IF (MONOCR.EQ.'GR') ICRYST=1
        IF (MONOCR.EQ.'QU') ICRYST=2
        IF (THM.EQ.0) THM=THMC(ICRYST,IXRAY)
        TANTHM=TAN(THM*PI/180)
C
C       RHO ANGLE (PARALLEL OR PERPENDICULAR GEOMETRY)
C
        IF (RHOM.EQ.0) RHOM=RHOMC(IDIFF)
C
C       DYNAMIC DIFFRACTION (PERFECT CRYSTAL) FRACTION
C
        IF (FRACTD.NE.0) SIGFD=0.05*FRACTD
      END IF
C
C     ESTIMATED ERRORS FOR DIFFRACTOMETER SCAN ANGLE AND DRIVE PULSE
C     TIMING
C
      J=IDIFF
      IF (SIGTHE.EQ.0) SIGTHE=ER(1,J)
      IF (SIGTIM.EQ.0) SIGTIM=ER(2,J)
C
C     CONSTANTS FOR COINCIDENCE CORRECTIONS
C
      J=IXRAY
      K=IDIFF
      IF    (TAU.EQ.0)    TAU=CC(1,J,K)
      IF (SIGTAU.EQ.0) SIGTAU=CC(2,J,K)
      IF    (ATT.EQ.0)    ATT=CC(3,J,K)
      IF (SIGATT.EQ.0) SIGATT=CC(4,J,K)
C
C     BEGINNING AND ENDING REFLECTION SERIAL NUMBERS AND X-RAY
C     EXPOSURE TIMES
C
      IF (II1.EQ.0) II1=-1E6
      IF (II2.EQ.0) II2=+1E6
      IF (XT1.EQ.0) XT1=-1E5
      IF (XT2.EQ.0) XT2=+1E5
C
C     INITIAL SCAN LIMITS (DECIMAL FRACTIONS OF SCAN WIDTH)
C
      IF (SC1.LT.0) SC1=0
      IF (SC2.LE.0.OR.SC2.GT.1) SC2=1
C
C     WINDOW WIDTH FOR INITIAL SCAN FOR PEAK LIMITS
C
      IF (WINDOW.EQ.0) WINDOW=1/12.0
C
C     E.S.D. MULTIPLIER FOR STATISTICAL SIGNIFICANCE TESTS
C
      IF (CUTOFF.EQ.0) CUTOFF=3
C
C     NUMBER OF DATA-SMOOTHING PASSES FOR PEAK LOCATION
C
      IF (NSMOOTH.EQ.0) NSMOOTH=2
      IF (NSMOOTH.LT.0) NSMOOTH=0
C
C     MULTIPLIER FOR LEHMANN-LARSEN BIAS CORRECTION
C
      IF (FLL.EQ.0) FLL=2
C
C     CALCULATE SPECTRAL DISPERSION COEFFICIENTS FOR ALPHA-ONE AND
C     ALPHA-TWO LINES.
C
C     DELTA(THETA) = (180/PI)*(DELTA(LAMBDA)/LAMBDA)*TAN(THETA)
C
C     LORENTZIAN LINE SHAPE
C
C       Y/YMAX = 1/(1 + (X - X0)**2/(W/2)**2)
C
C       AT (X - X0) = 5*(W/2),
C       WHERE W/2 = HALF-WIDTH AT HALF-HEIGHT,
C       (X - X0) = 2.5*W
C       Y/YMAX = 1/(1 + 5**2) = 1/26 = 0.03846
C
C     GAUSSIAN LINE SHAPE
C
C       Y/YMAX = EXP(-(X - X0)**2/(2*SIGMA**2))
C
C       AT Y/YMAX = 1/2,
C       1/2 = EXP(-(W/2)**2/(2*SIGMA**2))
C       LN(1/2) = -(W/2)**2/(2*SIGMA**2)
C       (W/2)/SIGMA = SQRT(-2*LN(1/2)) = SQRT(2*LN(2)) = 1.1774
C
C       W/2   = 1.1774*SIGMA
C       SIGMA = 0.8493*(W/2)
C
C       AT Y/YMAX = 1/26,
C       (X - X0) = SIGMA*SQRT(-2*LN(1/26))      = 2.5527*SIGMA
C       (X - X0) = (W/2)*SQRT(LN(1/26)/LN(1/2)) = 2.1680*(W/2)
C       (X - X0) = 1.084*W
C
      CL=2.500*180/PI
      CG=1.084*180/PI
      T1=DWL1/WL1
      T2=DWL2/WL2
      T1L=CL*T1
      T2L=CL*T2
      T1G=CG*T1
      T2G=CG*T2
C
C     THE CHARACTERISTIC EMISSION LINE SHAPES SHOULD BE APPROXIMATELY
C     LORENTZIAN.
C
      T1=T1L
      T2=T2L
      SIGT1=0
      SIGT2=0
C
C     SET RATIO OF DISPERSION COEFFICIENTS FOR PARALLEL-GEOMETRY
C     MONOCHROMATOR.
C
C     THE WIDTH OF THE TWICE-REFLECTED BEAM IS GIVEN BY
C
C     W = SQRT(A - 2*B*TANTH/TANTHM + C*(TANTH/TANTHM)**2),
C
C     WHERE
C
C     B = ALPHA**2 + 2*ETA**2,
C     C = ALPHA**2 + ETA**2,
C
C     ALPHA IS THE DIVERGENCE OF THE BEAM TO THE MONOCHROMATOR CRYSTAL,
C     AND ETA IS THE MOSAIC SPREAD OF THE MONOCHROMATOR CRYSTAL.
C
C     TH AND THM HAVE THE SAME SIGN FOR PARALLEL GEOMETRY, OPPOSITE
C     SIGNS FOR FOR ANTI-PARALLEL GEOMETRY.
C
C     W = SQRT(A - 2*F*C*TANTH/TANTHM + C*(TANTH/TANTHM)**2),
C
C     WHERE F = B/C = (ALPHA**2 + 2*ETA**2)/(ALPHA**2 + ETA**2), AND
C     1 .LE. F .LE. 2.
C
C     FOR FIXED- OR ROTATING-ANODE DIVERGENT X-RAY SOURCES, EXPECT
C     EFFECTIVELY ALPHA = 2*ETA AND THEREFORE F = 6/5.
C
      IF (RHOM.EQ.0.AND.F.EQ.0) F=1.2
C
C     WRITE CONTROL DATA TO BGLP.DAT FILE.
C
      OPEN (UNIT=IO2,FILE='bglp.dat',STATUS='NEW')
      WRITE (IO2,600) ATIME
      WRITE (IO2,600) ADATE
      WRITE (IO2,600) TITLE
      WRITE (IO2,600) FILE1
      WRITE (IO2,600) FILE2
      WRITE (IO2,601) (CELL(I),I=1,6)
      WRITE (IO2,602) IDIFF,INEUTRON
      WRITE (IO2,601) THM,RHOM,FRACTD,SIGFD,F
      WRITE (IO2,601) WLE,WL1,WL0,WL2
      WRITE (IO2,601) SIGTHE,SIGTIM
      WRITE (IO2,601) TAU,SIGTAU,ATT,SIGATT
      WRITE (IO2,601) WINDOW
      WRITE (IO2,601) CUTOFF
      WRITE (IO2,602) NSMOOTH
      WRITE (IO2,601) SC1,SC2
  600 FORMAT (1X,A)
  601 FORMAT (1X,6F10.5)
  602 FORMAT (1X,2I10)
C
C     PRINT CONTROL DATA.
C
      OPEN (UNIT=ILP,FILE='refpk.lp',STATUS='NEW',CARRIAGECONTROL=
     & 'FORTRAN')
      WRITE (ILP,610) ATIME,ADATE,TITLE
      WRITE (ILP,611) FILE1
      WRITE (ILP,612) FILE2
      WRITE (ILP,613) CELL
      IF (IDIFF.EQ.1) WRITE (ILP,614)
      IF (IDIFF.EQ.2) WRITE (ILP,619)
      IF (IDIFF.EQ.3) WRITE (ILP,615)
      IF (IDIFF.EQ.4) WRITE (ILP,616)
      WRITE (ILP,629) SOURCE
      IF (FILTER.NE.'  ') WRITE (ILP,630) FILTER
      IF (THM.NE.0) WRITE (ILP,618) MONOCR,THM,RHOM,FRACTD,SIGFD
      WRITE (ILP,631) WLE,WL1,WL0,WL2,DWL1,DWL2
      WRITE (ILP,632) T1,T2
      IF (WT2.NE.0.5) WRITE (ILP,622) WT2
      WRITE (ILP,617) SIGTHE,SIGTIM
      WRITE (ILP,641) TAU,SIGTAU,ATT,SIGATT
      IF (II1.GT.-1E6.OR.II2.LT.+1E6.OR.XT1.GT.-1E5.OR.XT2.LT.+1E5)
     & WRITE (ILP,640) II1,II2,XT1,XT2
      IF (SC1.GT.0.OR.SC2.LT.1) WRITE (ILP,642) SC1,SC2
      WRITE (ILP,648) WINDOW
      WRITE (ILP,646) CUTOFF
      IF (FLL.NE.2) WRITE (ILP,621) FLL
      IF (NSMOOTH.NE.2) WRITE (ILP,647) NSMOOTH
C
C     READ AND PRINT CONDITIONS LIMITING POSSIBLE REFLECTIONS.
C
      OPEN (UNIT=IO1,FILE='hklcond.dat',STATUS='UNKNOWN')
      WRITE (ILP,610) ATIME,ADATE,TITLE
      WRITE (ILP,643)
      NCOND=0
    7 READ (IO1,500,END=8) HCOND
      WRITE (ILP,644) HCOND
      I=0
      CALL HKLCOND (HCOND,I)
      NCOND=NCOND+1
      ICOND(NCOND)=I
      GO TO 7
    8 IF (NCOND.EQ.0) THEN
        WRITE (ILP,645)
        CLOSE (UNIT=IO1,DISPOSE='DELETE')
      ELSE
        CLOSE (UNIT=IO1)
      END IF
C
C     CONVERT FROM MICROSECONDS TO SECONDS.
C
      TAU=TAU*1E-6
      SIGTAU=SIGTAU*1E-6
      SIGTIM=SIGTIM*1E-6
C
C     CALCULATE VARIANCE VALUES.
C
      VARTAU=SIGTAU**2
      VARATT=SIGATT**2
      VARTHE=SIGTHE**2
      VARTIM=SIGTIM**2
C
C     OPEN INPUT REFLECTION DATA FILE.
C
      OPEN (UNIT=IO1,FILE=FILE1,STATUS='OLD',FORM='UNFORMATTED',
     & READONLY)
  610 FORMAT ('1'/'0PROGRAM REFPK.  ',A,1X,A,'.  ',A)
  611 FORMAT ('0INPUT FILE:   ',A)
  612 FORMAT ('0OUTPUT FILE:  ',A)
  613 FORMAT ('0LATTICE PARAMETERS'/'0',9X,'A',9X,'B',9X,'C',5X,'ALPHA',
     &6X,'BETA',5X,'GAMMA'/' ',3F10.4,3F10.3)
  614 FORMAT ('0INTERNATIONAL TABLES VOL. IV DIFFRACTOMETER AXES')
  619 FORMAT ('0BUSING AND LEVY DIFFRACTOMETER AXES')
  615 FORMAT ('0SIEMENS (NEE NICOLET, NEE SYNTEX) P3 DIFFRACTOMETER')
  616 FORMAT ('0ENRAF-NONIUS CAD4 DIFFRACTOMETER')
  618 FORMAT ('0CRYSTAL MONOCHROMATOR, BRAGG ANGLE, GEOMETRY ANGLE, AN',
     &'D DYNAMIC DIFFRACTION (PERFECT CRYSTAL) FRACTION'/'0CRYSTAL     T
     &HETA(M)    RHO(M) FRAC(DYN) SIG(FRAC)'/' ',A2,8X,2F10.2,2E10.3)
  629 FORMAT (' '/'$RADIATION:  ',A)
  630 FORMAT ('+  FILTER:  ',A)
  631 FORMAT ('0WAVELENGTHS'/'0FILT. EDGE ALPHA-ONE ALPHA-BAR ALPHA-TWO'
     &/' ',4F10.5/'0SPECTRAL LINE WIDTHS'/'0           ALPHA-ONE      ',
     &'     ALPHA-TWO'/' ',2(10X,F10.5))
  632 FORMAT ('0EXPECTED, LOWER LIMIT VALUES FOR T*TAN(THETA) SPECTRAL L
     &INE BROADENING COEFFICIENTS'/'0                  T1             ',
     &'     T2'/' ',2(10X,F10.4))
  622 FORMAT ('0INTENSITY OF ALPHA-TWO RELATIVE TO ALPHA-ONE'/'0       W
     &T2'/' ',F10.2)
  617 FORMAT ('0ERROR ESTIMATES FOR DIFFRACTOMETER SCAN ANGLE (DEGREES T
     &HETA) AND SCAN-ANGLE DRIVE-PULSE TIMING (MICROSECONDS)'/'0    SIGT
     &HE    SIGTIM'/' ',F10.4,F10.1)
  641 FORMAT ('0COUNTER DEAD TIME (MICROSECONDS) AND BEAM ATTENUATOR FAC
     &TOR'/'0       TAU    SIGTAU       ATT    SIGATT'/' ',4F10.2)
  640 FORMAT ('0BEGINNING AND ENDING REFLECTION SERIAL NUMBERS AND X-RAY
     & EXPOSURE TIMES (HOURS)'/'0     IIMIN     IIMAX               XTMI
     &N     XTMAX'/' ',2I10,10X,2F10.2)
  642 FORMAT ('0SCAN LIMITS (DECIMAL FRACTION OF SCAN WIDTH)'/'0      ',
     &'  I1        I2'/' ',2F10.4)
  648 FORMAT ('0WINDOW WIDTH FOR INITIAL SCAN FOR PEAK LIMITS (DECIMAL F
     &RACTION OF SCAN WIDTH)'/'0    WINDOW'/' ',F10.4)
  646 FORMAT ('0E.S.D. MULTIPLIER FOR STATISTICAL SIGNIFICANCE TESTS'/'0
     &    CUTOFF'/' ',F10.2)
  621 FORMAT ('0PROPORTIONALITY FACTOR FOR THE LEHMANN-LARSEN BIAS CORRE
     &CTION'/'0       FLL'/' ',F10.2)
  647 FORMAT ('0NUMBER OF DATA-SMOOTHING PASSES FOR PEAK LOCATION'/'0 ',
     &'  NSMOOTH'/' ',I10)
  643 FORMAT ('0CONDITIONS LIMITING POSSIBLE REFLECTIONS'/)
  644 FORMAT (' ',2X,A)
  645 FORMAT (' NONE')
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE LEHLAR (L)
C
C     LOCATE THE PEAK LIMITS BY THE METHOD OF LEHMANN AND LARSEN, WHICH
C     SEEKS THE TWO POINTS IN THE PROFILE BETWEEN WHICH INTEGRATION OF
C     THE NET INTENSITY GIVES MAXIMUM I/SIGMA(I).
C
C     M. S. LEHMANN AND F. K. LARSEN (1974).  ACTA CRYST. A30, 580-584.
C
C     R. H. BLESSING, P. COPPENS, AND P. J. BECKER (1974).  J. APPL.
C     CRYST. 7, 448-492.
C
C     RETURN L = 0 TO FLAG AN ERROR CONDITION.
C
      PARAMETER (NMAX=100)
      COMMON /BLOCK5/ NSTEP,XX(NMAX),YY(NMAX),VX(NMAX),VY(NMAX)
      COMMON /BLOCK6/ I1,I2,L1,L2,X1,X2,FLL
      DIMENSION Y(NMAX),V(NMAX)
C
C     ANALYZE LOWER ANGLE ALPHA-ONE HALF-PEAK.
C
      N=0
      DO 1 I=I1,NINDEX(X1),+1
      N=N+1
      Y(N)=YY(I)
      V(N)=VY(I)
    1 CONTINUE
      CALL LLIMIT (N,Y,V,FLL,L)
      IF (L.EQ.0) RETURN
      L1=I1+L-1
C
C     ANALYZE HIGHER ANGLE ALPHA-TWO HALF-PEAK.
C
      N=0
      DO 2 I=I2,NINDEX(X2),-1
      N=N+1
      Y(N)=YY(I)
      V(N)=VY(I)
    2 CONTINUE
      CALL LLIMIT (N,Y,V,FLL,L)
      IF (L.EQ.0) RETURN
      L2=I2-L+1
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE LLCORR (N,Y,F,L)
C
C     CORRECT FOR THE MAXIMUM I/SIGMA(I) PEAK LIMIT BIAS.
C
C     THE METHOD FINDS PEAK LIMITS THAT ARE SLIGHTLY WITHIN THE PEAKS.
C     LEHMANN AND LARSEN HAVE SHOWN THAT THE BIAS OF THE PEAK LIMITS IS
C     PROPORTIONAL TO 1/(2 + R), WHERE R IS THE PEAK-TO-BACKGROUND
C     RATIO.
C
C     RETURN L = 0 TO FLAG AN ERROR CONDITION.
C
      DIMENSION Y(N)
      P=0
      B=0
      DO 1 I=1,L-1
      B=B+Y(I)
    1 CONTINUE
      DO 2 I=L,N
      P=P+Y(I)
    2 CONTINUE
      NP=N-L+1
      NB=L-1
      B=B*NP/NB
      R=(P-B)/B
      IF (R.LE.0) THEN
        L=0
        RETURN
      END IF
      Q=F/(2+R)
C
C     RMIN = 0, QMAX = 0.5*F.  DEFAULT:  F = 2.0, QMAX = 1.
C
      L=L-NINT(Q*NP)
      IF (L.LT.0) L=0
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE LLIMIT (N,Y,V,F,L)
C
C     FIND MAXIMUM I/SIGMA(I) PEAK LIMIT.
C
      DIMENSION Y(N),V(N)
      IF (N.LT.4) GO TO 9
      P=0
      B=0
      VP=0
      VB=0
      XP=N
      XB=0
      QMAX=0
      DO 1 I=1,N
      P=P+Y(I)
      VP=VP+V(I)
    1 CONTINUE
      DO 2 I=1,N-1
      P=P-Y(I)
      B=B+Y(I)
      VP=VP-V(I)
      VB=VB+V(I)
      XP=XP-1
      XB=XB+1
      Q=XP/XB
      Q=(P-Q*B)/SQRT(VP+Q**2*VB)
      IF (Q.GT.QMAX) THEN
        QMAX=Q
        L=I
      END IF
C
C     STORE THE Q = I/SIGMA(I) DATA IN THE NO-LONGER NEEDED V-ARRAY.
C
      V(I)=Q
    2 CONTINUE
C
C     FIT A THREE- TO NINE-POINT PARABOLA TO THE I/SIGMA(I) DATA, AND
C     FIND ITS MAXIMUM.
C 
      M=MIN(L-1,N-L-1,4)
      IF (M.LT.1) GO TO 9
      J=0
      DO 11 I=L-M,L+M
      J=J+1
      V(J)=V(I)
   11 CONTINUE
      K=J
      DO 12 I=L-M,L+M
      K=K+1
      V(K)=I
   12 CONTINUE
      CALL PFIT (J,V(J+1),V(1))
      C1=V(J+2)
      C2=V(J+3)
      IF (C2.GE.0) GO TO 9
      L=NINT(-0.5*C1/C2)
      IF (L.LE.1.OR.L.GE.N-1) GO TO 9
      CALL LLCORR (N,Y,F,L)
      RETURN
    9 CONTINUE
C
C     RETURN L = 0 TO FLAG AN ERROR CONDITION.
C
      L=0
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE MATINV (N,M,A,D)
C
C     INVERT A SQUARE MATRIX BY GAUSS-JORDAN ELIMINATION, AND CALCULATE
C     ITS 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 A SYMMETRIC 3 X 3 MATRIX AND CALCULATE ITS DETERMINANT.
C
C     A IS MATRIX, B IS INVERSE, D IS DETERMINANT OF A.
C
      DIMENSION A(3,3),B(3,3)
      DO 1 I=1,3
      K=1+MOD(I+1,3)
      M=1+MOD(I+3,3)
      DO 1 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)
    1 CONTINUE
      D=0
      DO 2 I=1,3
      D=D+A(I,1)*B(1,I)
    2 CONTINUE
      DO 3 I=1,3
      DO 3 J=1,3
      B(I,J)=B(I,J)/D
    3 CONTINUE
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE METRIC (A,GB)
C
C     CALCULATE RECIPROCAL LATTICE PARAMETERS AND DIRECT AND RECIPROCAL
C     LATTICE METRIC TENSORS FROM DIRECT LATTICE PARAMETERS.
C
      DIMENSION A(6),B(6),GA(3,3),GB(3,3)
      PI=ACOS(-1.0)
      DEG=180/PI
      RAD=PI/180
      CA=COS(A(4)*RAD)
      CB=COS(A(5)*RAD)
      CG=COS(A(6)*RAD)
      SA=SIN(A(4)*RAD)
      SB=SIN(A(5)*RAD)
      SG=SIN(A(6)*RAD)
      V=A(1)*A(2)*A(3)*SQRT(1-CA**2-CB**2-CG**2+2*CA*CB*CG)
      B(1)=A(2)*A(3)*SA/V
      B(2)=A(1)*A(3)*SB/V
      B(3)=A(1)*A(2)*SG/V
      B(4)=(CB*CG-CA)/(SB*SG)
      B(5)=(CA*CG-CB)/(SA*SG)
      B(6)=(CA*CB-CG)/(SA*SB)
      GA(1,1)=A(1)*A(1)
      GA(1,2)=A(1)*A(2)*CG
      GA(1,3)=A(1)*A(3)*CB
      GA(2,1)=GA(1,2)
      GA(2,2)=A(2)*A(2)
      GA(2,3)=A(2)*A(3)*CA
      GA(3,1)=GA(1,3)
      GA(3,2)=GA(2,3)
      GA(3,3)=A(3)*A(3)
      GB(1,1)=B(1)*B(1)
      GB(1,2)=B(1)*B(2)*B(6)
      GB(1,3)=B(1)*B(3)*B(5)
      GB(2,1)=GB(1,2)
      GB(2,2)=B(2)*B(2)
      GB(2,3)=B(2)*B(3)*B(4)
      GB(3,1)=GB(1,3)
      GB(3,2)=GB(2,3)
      GB(3,3)=B(3)*B(3)
      B(4)=ACOS(B(4))*DEG
      B(5)=ACOS(B(5))*DEG
      B(6)=ACOS(B(6))*DEG
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE MV (N,A,B,X)
C
C     MATRIX A TIMES COLUMN VECTOR B
C
      DIMENSION A(N,N),B(N),X(N)
      DOUBLE PRECISION A
      DO 1 I=1,N
      X(I)=0
      DO 1 J=1,N
      X(I)=X(I)+A(I,J)*B(J)
    1 CONTINUE
      RETURN
      END
C-----------------------------------------------------------------------
      FUNCTION NINDEX(X)
C
C     FIND THE INDEX I OF THE ABSCISSA XX(I) NEAREST THE VALUE X.
C
C     THE ABSCISSAE IN THE ARRAY XX(NSTEP) MUST BE ORDERED SUCH THAT
C     XX(I) INCREASES WITH I.
C
      PARAMETER (NMAX=100)
      COMMON /BLOCK5/ NSTEP,XX(NMAX),YY(NMAX),VX(NMAX),VY(NMAX)
      IF (X.LT.XX(1)) THEN
        NINDEX=0
        RETURN
      END IF
      IF (X.GT.XX(NSTEP)) THEN
        NINDEX=NSTEP+1
        RETURN
      END IF
      DO 1 I=2,NSTEP
      IF (XX(I).GE.X) THEN
        IF (X-XX(I-1).LT.XX(I)-X) THEN
          NINDEX=I-1
        ELSE
          NINDEX=I
        END IF
        RETURN
      END IF
    1 CONTINUE
      END
C-----------------------------------------------------------------------
      SUBROUTINE OUTPUT
C
C         PRINTED OUTPUT COMPARING OBSERVED AND CALCULATED PEAK
C         POSITIONS AND PEAK LIMITS
C
      PARAMETER (NMAX=100,M1=4,M2=16,M3=50)
      CHARACTER ATIME*8,ADATE*9,TITLE*80,HEADING*120
      DIMENSION HEADING(M1),JXMIN(M1,M2,M3),JXMAX(M1,M2,M3),JX(M2)
      COMMON /BLOCK0/ IO1,IO2,IO3,ILP
      COMMON /BLOCK1/ ATIME,ADATE,TITLE
      COMMON /BLOCK2/ GINV(3,3),WLE,WL1,WL0,WL2,T1,T2,SIGT1,SIGT2
      COMMON /BLOCK4/ IDIFF,JMODEL,TANTHM,RHOM,F,U(3,3),MODEL1,MODEL2,
     & Q1(3,3),Q2(3,3)
      DIMENSION H(3),Y(3),Z(3)
      DATA DEG, RAD /57.2957795, 0.017453293/
      CHARACTER TYPE1*4,TYPE2*4
      IF (MODEL1.EQ.1) TYPE1='LRNZ'
      IF (MODEL1.EQ.2) TYPE1='GAUS'
      IF (MODEL1.EQ.3) TYPE1='PLMC'
      IF (MODEL2.EQ.1) TYPE2='LRNZ'
      IF (MODEL2.EQ.2) TYPE2='GAUS'
      IF (MODEL2.EQ.3) TYPE2='PLMC'
C
C         INITIALIZE ARRAYS FOR CATALOGING PROFILES WITH EXTREME
C         FEATURES.
C
      DO 10 J1=1,M1
      DO 10 J2=1,M2
      DO 10 J3=1,M3
      JXMIN(J1,J2,J3)=+1E6
 10   JXMAX(J1,J2,J3)=-1E6
C
C         LOOP THROUGH DATA FILE OF LEAST-SQUARES OBSERVATIONS.
C
      WRITE (ILP,610) ATIME,ADATE,TITLE
      WRITE (ILP,635)
      NOBS=0
      REWIND IO3
 1    READ (IO3,END=9) II,H,TWOTH,OMEGA,CHI,PHI,DSTAR,Y,Z,
     & I1,I2,XM,X0,X1,X2,L1,L2,W1,W2,TANTH
      IF (II.LT.0) GO TO 1
      NOBS=NOBS+1
C
C         CALCULATE PEAK CENTROID.
C
      DO 11 I=1,3
      Y(I)=0
      DO 11 J=1,3
 11   Y(I)=Y(I)+H(J)*U(J,I)
      DO 12 I=1,3
 12   Y(I)=Y(I)/DSTAR
      IF (TWOTH.LT.0) THEN
        DO 13 I=1,3
 13     Y(I)=-Y(I)
      END IF
      COSPH=COS(PHI*RAD)
      SINPH=SIN(PHI*RAD)
      COSCH=COS(CHI*RAD)
      SINCH=SIN(CHI*RAD)
      SINOM=Y(1)*COSPH-Y(2)*SINPH
      COSOM=(Y(1)*SINPH+Y(2)*COSPH)*COSCH-Y(3)*SINCH
      DELTA=ARCTAN(SINOM,COSOM)*DEG-OMEGA
      DELTA=A180(DELTA)
      IF (TWOTH.LT.0) DELTA=-DELTA
      X0CALC=XM-DELTA
C
C         CALCULATE PEAK LIMITS.
C
      U1=0
      U2=0
      Q=0
      DO 15 I=1,3
      DO 15 J=1,3
      U1=U1+Z(J)*Q1(I,J)*Z(I)
 15   U2=U2+Z(J)*Q2(I,J)*Z(I)
      IF (TYPE1.EQ.'LRNZ') W1CALC=SQRT(U1)+T1*TANTH
      IF (TYPE1.EQ.'GAUS') W1CALC=SQRT(U1+(T1*TANTH)**2)
      IF (TYPE1.EQ.'PLMC') W1CALC=SQRT(U1-2*F*T1*TANTH/TANTHM
     &                                   +T1*(TANTH/TANTHM)**2)
      IF (TYPE2.EQ.'LRNZ') W2CALC=SQRT(U2)+T2*TANTH
      IF (TYPE2.EQ.'GAUS') W2CALC=SQRT(U2+(T2*TANTH)**2)
      IF (TYPE2.EQ.'PLMC') W2CALC=SQRT(U2-2*F*T2*TANTH/TANTHM
     &                                   +T2*(TANTH/TANTHM)**2)
      L1CALC=NINDEX(X1-W1CALC)
      L2CALC=NINDEX(X2+W2CALC)
C
C         PRINT SELECTED REFLECTION RECORDS.
C
      I=II
      J=H(1)
      K=H(2)
      L=H(3)
      CALL DPRINT (I,J,K,L,IPRINT)
      IF (IPRINT.EQ.1) THEN
        I0=NINDEX(X0)
        I0CALC=NINDEX(X0CALC)
        WRITE (ILP,636) I,J,K,L,I0,I0CALC,L1,L2,L1CALC,L2CALC
      END IF
C
C         CATALOG PROFILES WITH EXTREME FEATURES.
C
      JX(1)=I
      JX(2)=J
      JX(3)=K
      JX(4)=L
      JX(5)=NINDEX(X0)
      JX(6)=NINDEX(X0CALC)
      JX(7)=1000*(X0-XM)
      JX(8)=1000*(X0CALC-XM)
      JX(9)=1000*(X0-X0CALC)
      JX(10)=L1
      JX(11)=L2
      JX(12)=L1CALC
      JX(13)=L2CALC
      W=W1+X2-X1+W2
      WCALC=W1CALC+X2-X1+W2CALC
      JX(14)=1000*W
      JX(15)=1000*WCALC
      JX(16)=1000*(W-WCALC)
      CALL EXTREM (JXMIN,JXMAX,JX)
C
C         LOOP BACK FOR NEXT REFLECTION.
C
      GO TO 1
C
C         END LOOP THROUGH REFLECTION DATA.
C
 9    CLOSE (UNIT=IO3)
C
C         PRINT CATALOGS OF EXTREME CASES.
C
      WRITE (ILP,610) ATIME,ADATE,TITLE
      WRITE (ILP,649)
      DO 29 J1=1,M1
      WRITE (ILP,610) ATIME,ADATE,TITLE
      WRITE (ILP,645) HEADING(J1)
      DO 39 J3=1,MIN(M3,NOBS)
      I1=JXMIN(J1,1,J3)
      I2=JXMIN(J1,2,J3)
      I3=JXMIN(J1,3,J3)
      I4=JXMIN(J1,4,J3)
      I5=JXMIN(J1,5,J3)
      I6=JXMIN(J1,6,J3)
      F7=JXMIN(J1,7,J3)*0.001
      F8=JXMIN(J1,8,J3)*0.001
      F9=JXMIN(J1,9,J3)*0.001
      I10=JXMIN(J1,10,J3)
      I11=JXMIN(J1,11,J3)
      I12=JXMIN(J1,12,J3)
      I13=JXMIN(J1,13,J3)
      F14=JXMIN(J1,14,J3)*0.001
      F15=JXMIN(J1,15,J3)*0.001
      F16=JXMIN(J1,16,J3)*0.001
      IF (F8.LT.-999.99) F8=-999.99
      IF (F8.GT.+999.99) F8=+999.99
      IF (F9.LT.-999.99) F9=-999.99
      IF (F9.GT.+999.99) F9=+999.99
 39   WRITE (ILP,646) I1,I2,I3,I4,I5,I6,F7,F8,F9,I10,I11,I12,I13,
     & F14,F15,F16
      WRITE (ILP,610) ATIME,ADATE,TITLE
      WRITE (ILP,645) HEADING(J1)
      DO 49 J3=1,MIN(M3,NOBS)
      I1=JXMAX(J1,1,J3)
      I2=JXMAX(J1,2,J3)
      I3=JXMAX(J1,3,J3)
      I4=JXMAX(J1,4,J3)
      I5=JXMAX(J1,5,J3)
      I6=JXMAX(J1,6,J3)
      F7=JXMAX(J1,7,J3)*0.001
      F8=JXMAX(J1,8,J3)*0.001
      F9=JXMAX(J1,9,J3)*0.001
      I10=JXMAX(J1,10,J3)
      I11=JXMAX(J1,11,J3)
      I12=JXMAX(J1,12,J3)
      I13=JXMAX(J1,13,J3)
      F14=JXMAX(J1,14,J3)*0.001
      F15=JXMAX(J1,15,J3)*0.001
      F16=JXMAX(J1,16,J3)*0.001
      IF (F8.LT.-999.99) F8=-999.99
      IF (F8.GT.+999.99) F8=+999.99
      IF (F9.LT.-999.99) F9=-999.99
      IF (F9.GT.+999.99) F9=+999.99
 49   WRITE (ILP,646) I1,I2,I3,I4,I5,I6,F7,F8,F9,I10,I11,I12,I13,
     & F14,F15,F16
 29   CONTINUE
 610  FORMAT ('1'/'0PROGRAM REFPK.  ',A,1X,A,'.  ',A)
 635  FORMAT ('0     I     H   K   L     X0O X0C     L1O L2O     L1C L2C
     &'/)
 636  FORMAT (' ',I6,2X,3I4,3(4X,2I4))
 649  FORMAT (
     &'0THE FOLLOWING TABLES LIST REFLECTIONS WITH EXTREME VALUES OF'/
     &' SEVERAL CHARACTERISTIC PROPERTIES.  THE TABLES ARE PRESENTED'/
     &' IN PAIRS IN WHICH THE FIRST TABLE LISTS THE REFLECTIONS WITH'/
     &' THE SMALLEST VALUES OF A GIVEN PROPERTY AND THE SECOND THOSE'/
     &' WITH THE LARGEST VALUES.  THE PROPERTY ON WHICH THE TABLE IS'/
     &' BASED IS STATED IN THE HEADING ABOVE THE TABLE COLUMN'/
     &' HEADINGS, AND THE ENTRIES IN THE TABLE ARE SORTED IN ORDER OF'/
     &' INCREASING VALUE OF THE PROPERTY IN THE TABLES OF SMALLEST'/
     &' VALUES, AND DECREASING VALUE IN THE TABLES OF LARGEST VALUES.')
      DATA HEADING /
     &'DX0O, DISPLACEMENT OF OBSERVED PEAK CENTROID FROM MIDPOINT OF SCA
     &N (DEGREES THETA)',
     &'D0=X0O-X0C, DIFFERENCE BETWEEN OBSERVED AND CALCULATED CENTROIDS 
     &(DEGREES THETA)',
     &'WC, CALCULATED BASE WIDTH (DEGREES THETA)',
     &'WO-WC, DIFFERENCE BETWEEN OBSERVED AND CALCULATED FULL-PEAK BASE 
     &WIDTH (DEGREES THETA)'/
 645  FORMAT ('0',A/
     &'0     I     H   K   L     X0O X0C      DX0O    DX0C X0O-X0C',
     &'     L1O L2O   L1C L2C   L2O-L1O L2C-L1C   WO-WC'
     &/
     &'                                                         D0',
     &'                              WO      WC'/)
 646  FORMAT (' ',I6,2X,3I4,4X,2I4,2X,3F8.2,4X,2I4,2X,2I4,2X,3F8.2)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE PFIT (N,X,Y)
C
C     FIT A PARABOLA, Y = C0 + C1*X + C2*X**2.
C
C     RETURN C0, C1, AND C2 IN X(1), X(2), AND X(3), RESPECTIVELY.
C
      DIMENSION X(N),Y(N),A(3,3),AINV(3,3),B(3)
      IF (N.LT.3) GO TO 9
      X0=0.5*(X(1)+X(N))
      IF (N.GT.3) GO TO 5
C
C     FOR THREE POINTS, FIND THE PARABOLA THROUGH THE POINTS.
C
      X1=X(1)-X0
      X2=X(2)-X0
      X3=X(3)-X0
      Y1=Y(1)
      Y2=Y(2)
      Y3=Y(3)
      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
      GO TO 2
    5 CONTINUE
C
C     FOR MORE THAN THREE POINTS, FIND THE LEAST-SQUARES PARABOLA.
C
      SUMX=0
      SUMX2=0
      SUMX3=0
      SUMX4=0
      SUMY=0
      SUMXY=0
      SUMX2Y=0
      DO 1 I=1,N
      XI=X(I)-X0
      YI=Y(I)
      SUMX=SUMX+XI
      SUMX2=SUMX2+XI*XI
      SUMX3=SUMX3+XI*XI*XI
      SUMX4=SUMX4+XI*XI*XI*XI
      SUMY=SUMY+YI
      SUMXY=SUMXY+XI*YI
      SUMX2Y=SUMX2Y+XI*XI*YI
    1 CONTINUE
      A(1,1)=N
      A(1,2)=SUMX
      A(1,3)=SUMX2
      A(2,1)=SUMX
      A(2,2)=SUMX2
      A(2,3)=SUMX3
      A(3,1)=SUMX2
      A(3,2)=SUMX3
      A(3,3)=SUMX4
      B(1)=SUMY
      B(2)=SUMXY
      B(3)=SUMX2Y
    2 CALL MATINV3 (A,AINV,DETA)
      IF (DETA.EQ.0) GO TO 9
      DO 3 I=1,3
      X(I)=0
      DO 3 J=1,3
      X(I)=X(I)+AINV(I,J)*B(J)
    3 CONTINUE
C
C     ADJUST FOR SHIFTED ORIGIN.
C
      X(1)=X(1)-X(2)*X0+X(3)*X0**2
      X(2)=X(2)-2*X(3)*X0
      RETURN
    9 CONTINUE
C
C     RETURN C2 = 0 TO FLAG AN ERROR CONDITION.
C
      X(3)=0
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE POSDEF (Q,ILP)
      DIMENSION Q(3,3)
      R=Q(1,1)
      S=Q(1,1)*Q(2,2)-Q(2,1)*Q(1,2)
      T=Q(1,1)*Q(2,2)*Q(3,3)+Q(1,2)*Q(2,3)*Q(3,1)+Q(1,3)*Q(2,1)*Q(3,2)
     &   -Q(3,1)*Q(2,2)*Q(1,3)-Q(3,2)*Q(2,3)*Q(1,1)-Q(3,3)*Q(2,1)*Q(1,2)
      IF (MIN(R,S,T).GT.0) RETURN
      WRITE (ILP,6000) ((Q(I,J),J=1,3),I=1,3)
      T=(Q(1,1)+Q(2,2)+Q(3,3))/3
      DO 1 I=1,3
      DO 1 J=1,3
      IF (I.EQ.J) THEN
        Q(I,J)=T
      ELSE
        Q(I,J)=0
      END IF
    1 CONTINUE
      WRITE (ILP,6001) ((Q(I,J),J=1,3),I=1,3)
 6000 FORMAT (/'0ATTENTION!  FITTED PEAK-WIDTH ANISOTROPY TENSOR WAS NOT
     & POSITIVE-DEFINITE,'/' ----------'/'0Q11 Q12 Q13   ',3E10.3/' Q2',
     &'1 Q22 Q23 = ',3E10.3/' Q31 Q32 Q33   ',3E10.3/'0AND HAS BEEN REPL
     &ACED BY ITS ISOTROPIC EQUIVALENT')
 6001 FORMAT ('0Q11 Q12 Q13   ',3E10.3/' Q21 Q22 Q23 = ',3E10.3/' Q31 Q3
     &2 Q33   ',3E10.3)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE READ1 (IFILE,IEND,II,IH,IK,IL,ANGLES,WIDTH,SPEED,XTIME,
     & NSTEP,YY)
      DIMENSION ANGLES(4),YY(1)
      INTEGER*2 JH,JK,JL,JY(96)
      READ (IFILE,END=9) II,JH,JK,JL,ANGLES,WIDTH,SPEED,(X,I=1,5),XTIME,
     & JY
      IH=JH
      IK=JK
      IL=JL
      NSTEP=96
      DO 1 I=1,NSTEP
      YY(I)=JY(I)
    1 CONTINUE
      IEND=0
      RETURN
    9 IEND=1
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE READR (IFILE,IEND,II,IH,IK,IL,ANGLES,NSTEP,XX,YY,VX,VY,
     & XTIME)
C
C     SUBROUTINE READR MUST RETURN:
C
C       INTEGER*4 VARIABLES:
C
C         II        - MEASUREMENT SERIAL NUMBER
C         IH,IK,IL  - REFLECTION INDICES
C         NSTEP     - NUMBER OF SCAN STEPS
C
C       REAL*4 ARRAYS:
C
C         ANGLES(4) - DIFFRACTOMETER SETTING ANGLES
C         XX(NSTEP) - ROCKING-ANGLE STEPS
C         YY(NSTEP) - COUNT-RATE STEPS
C         VX(NSTEP) - ROCKING-ANGLE VARIANCES
C         VY(NSTEP) - COUNT-RATE VARIANCES
C
C       REAL*4 VARIABLE:
C
C         XTIME     - RADIATION EXPOSURE TIME
C
C     THE DIFFRACTOMETER ANGLES MUST BE IN DEGREES, AND THEIR ORDER IS
C     IMPORTANT:
C
C         DIFF.     ANGLES(1)  (2)    (3)    (4)
C
C         INT.TAB.  TWO-THETA  OMEGA  CHI    PHI
C         BUS.LEV.  TWO-THETA  OMEGA  CHI    PHI
C         P3        TWO-THETA  OMEGA  PHI    CHI
C         CAD4      THETA      PHI    OMEGA  KAPPA
C
C     THE ROCKING-ANGLE STEPS MUST INCREASE (FROM NEGATIVE TO POSITIVE)
C     IN ORDER OF INCREASING ABSOLUTE VALUE OF THE BRAGG ANGLE, THETA.
C
C     RECOMMENDED UNITS:
C
C         ROCKING ANGLE  - DEGREES (THETA)
C         COUNT RATES    - COUNTS PER SECOND
C         EXPOSURE TIME  - HOURS
C
      DIMENSION ANGLES(4),XX(1),YY(1),VX(1),VY(1)
      IEND=0
      CALL READ1 (IFILE,IEND,II,IH,IK,IL,ANGLES,WIDTH,SPEED,XTIME,NSTEP,
     & YY)
      IF (IEND.NE.0) RETURN
C
C     CONVERT FROM STEP-SCAN COUNT PROFILE TO COUNT-RATE VERSUS ROCKING-
C     ANGLE PROFILE.
C
      CALL XYDATA (WIDTH,SPEED,NSTEP,XX,YY,VX,VY)
      RETURN
      END
C-----------------------------------------------------------------------
      FUNCTION SINTHL(H,GINV)
      DIMENSION H(3),GINV(3,3)
      Q=0
      DO 1 I=1,3
      DO 1 J=1,3
      Q=Q+H(J)*GINV(I,J)*H(I)
    1 CONTINUE
      SINTHL=0.5*SQRT(Q)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE SMOOTH (NSMOOTH,N,Y)
C
C     SMOOTHING OF EQUALLY SPACED DATA BY FOURTH DIFFERENCES
C
C     EQUIVALENT TO A FIVE-POINT LEAST-SQUARES FIT OF A PARABOLA,
C
C     Y = A0 + A1*X + A2*X**2,
C
C     TO THE (I - 2)ND, (I - 1)ST, ITH, (I + 1)ST, AND (I + 2)ND
C     POINTS, I.E., TO X = -2, -1, 0, 1, AND 2.
C
C     C. LANCZOS (1956).  APPLIED ANALYSIS, PP. 316-320.  PRENTICE-
C     HALL, ENGLEWOOD CLIFFS, NEW JERSEY.
C
      DIMENSION Y(N),T(-2:+2)
      IF (NSMOOTH.LE.0.OR.N.LT.6) RETURN
      DO 9 I=1,NSMOOTH
      DO 1 J=-2,+2
      T(J)=Y(3+J)
    1 CONTINUE
      D3=T(1)-3*T(0)+3*T(-1)-T(-2)
      D4=T(2)-4*T(1)+6*T(0)-4*T(-1)+T(-2)
      Y(1)=Y(1)+D3/5+3*D4/35
      Y(2)=Y(2)-2*D3/5-D4/7
      Y(3)=Y(3)-3*D4/35
      DO 2 J=4,N-2
      DO 3 K=-2,+1
      T(K)=T(K+1)
    3 CONTINUE
      T(2)=Y(J+2)
      D4=T(2)-4*T(1)+6*T(0)-4*T(-1)+T(-2)
      Y(J)=Y(J)-3*D4/35
    2 CONTINUE
      D3=T(2)-3*T(1)+3*T(0)-T(-1)
      Y(N-1)=Y(N-1)+2*D3/5-D4/7
      Y(N)=Y(N)-D3/5+3*D4/35
      DO 9 J=1,N
      IF (Y(J).LT.0) Y(J)=0
    9 CONTINUE
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE TFIT
C
C     LEAST-SQUARES FIT OF ISOTROPIC, SCALAR COEFFICIENTS OF TAN(THETA)
C     SPECTRAL DISPERSION TERMS
C
      COMMON /BLOCK0/ IO1,IO2,IO3,ILP
      COMMON /BLOCK2/ GINV(3,3),WLE,WL1,WL0,WL2,T1,T2,SIGT1,SIGT2
      COMMON /BLOCK4/ IDIFF,JMODEL,TANTHM,RHOM,F,U(3,3),MODEL1,MODEL2,
     & Q1(3,3),Q2(3,3)
      CHARACTER MODEL*4,TYPE1*4,TYPE2*4
      DIMENSION H(3),Y(3),Z(3),XX(3)
      DIMENSION AL(2,2),B1L(2),B2L(2),X1L(2),X2L(2)
      DIMENSION AG(2,2),B1G(2),B2G(2),X1G(2),X2G(2)
      DIMENSION AM(3,3),B1M(3),B2M(3),X1M(3),X2M(3)
      DOUBLE PRECISION AG,AL,AM
C
C     ZERO NORMAL MATRICES AND VECTORS.
C
      NPAR=2
      DO 20 I=1,NPAR
      DO 21 J=I,NPAR
      AL(I,J)=0
      AG(I,J)=0
   21 CONTINUE
      B1L(I)=0
      B2L(I)=0
      B1G(I)=0
      B2G(I)=0
   20 CONTINUE
      NPAR=3
      DO 30 I=1,NPAR
      DO 31 J=I,NPAR
      AM(I,J)=0
   31 CONTINUE
      B1M(I)=0
      B2M(I)=0
   30 CONTINUE
C
C     LOOP THROUGH DATA FILE OF LEAST-SQUARES OBSERVATIONS TO BUILD
C     NORMAL EQUATIONS.
C
      NOBS=0
      REWIND IO3
    1 READ (IO3,END=9) II,H,TWOTH,OMEGA,CHI,PHI,DSTAR,Y,Z,
     & I1,I2,XM,X0,X1,X2,L1,L2,W1,W2,TANTH
C
C     ACCUMULATE NORMAL MATRICES AND VECTORS.
C
      NOBS=NOBS+1
      IF (TANTHM.EQ.0.OR.RHOM.EQ.90) THEN
C
C       BETA-FILTER OR PERPENDICULAR MONOCHROMATOR GEOMETRY
C
        NPAR=2
C
C       LORENTZIAN CONVOLUTION MODEL, W = SQRT(Q) + T*TANTH
C
        XX(1)=1
        XX(2)=TANTH
        Y1=W1
        Y2=W2
        CALL ACCUM2 (NPAR,XX,Y1,Y2,AL,B1L,B2L)
C
C       GAUSSIAN CONVOLUTION MODEL, W**2 = Q + T**2*(TANTH)**2
C
        XX(1)=1
        XX(2)=TANTH**2
        Y1=W1**2
        Y2=W2**2
        CALL ACCUM2 (NPAR,XX,Y1,Y2,AG,B1G,B2G)
      ELSE
C
C       PARALLEL MONOCHROMATOR GEOMETRY.  GAUSSIAN CONVOLUTION MODEL,
C       W**2 = Q - 2*F*T*TANTH/TANTHM + T*(TANTH/TANTHM)**2 ,  1 < F < 2
C
        NPAR=3
        XX(1)=1
        XX(2)=-2*TANTH/TANTHM
        XX(3)=(TANTH/TANTHM)**2
        Y1=W1**2
        Y2=W2**2
        CALL ACCUM2 (NPAR,XX,Y1,Y2,AM,B1M,B2M)
      END IF
C
C     LOOP BACK FOR NEXT REFLECTION.
C
      GO TO 1
C
C     END LOOP THROUGH REFLECTION DATA.
C
    9 CONTINUE
C
C     SOLVE LEAST-SQUARES NORMAL EQUATIONS, AND CALCULATE STATISTICS OF
C     FIT AND PARAMETER ERRORS.
C
      IF (TANTHM.EQ.0.OR.RHOM.EQ.90) THEN
C
C       BETA-FILTER OR PERPENDICULAR MONOCHROMATOR
C
        F=0
        NPAR=2
        IF (NOBS.LE.NPAR) STOP 'NOT ENOUGH DATA IN SUBROUTINE TFIT'
C
C       LORENTZIAN MODEL.  Q = X(1)**2, T = X(2)
C
        CALL MATINV (NPAR,NPAR,AL,DET)
        IF (DET.EQ.0) STOP 'ILL-CONDITIONED MATRIX IN SUBROUTINE TFIT'
        CALL MV (NPAR,AL,B1L,X1L)
        CALL MV (NPAR,AL,B2L,X2L)
        MODEL='LRNZ'
        Q1L=X1L(1)**2
        Q2L=X2L(1)**2
        T1L=X1L(2)
        T2L=X2L(2)
        CALL AGREE1 (MODEL,IO3,Q1L,Q2L,T1L,T2L,0.0,0.0,WL0,
     &   NOBS,NPAR,CHSQ1L,CHSQ2L,RMSD1L,RMSD2L)
        SIGT1L=SQRT(CHSQ1L*AL(2,2))
        SIGT2L=SQRT(CHSQ2L*AL(2,2))
C
C       GAUSSIAN MODEL.  Q = X(1), T = SQRT(X(2))
C
        CALL MATINV (NPAR,NPAR,AG,DET)
        IF (DET.EQ.0) STOP 'ILL-CONDITIONED MATRIX IN SUBROUTINE TFIT'
        CALL MV (NPAR,AG,B1G,X1G)
        CALL MV (NPAR,AG,B2G,X2G)
        MODEL='GAUS'
        Q1G=X1G(1)
        Q2G=X2G(1)
        T1G=0
        T2G=0
        IF (X1G(2).GT.0) T1G=SQRT(X1G(2))
        IF (X2G(2).GT.0) T2G=SQRT(X2G(2))
        CALL AGREE1 (MODEL,IO3,Q1G,Q2G,T1G,T2G,0.0,0.0,WL0,
     &   NOBS,NPAR,CHSQ1G,CHSQ2G,RMSD1G,RMSD2G)
        SIGT1G=0
        SIGT2G=0
        IF (T1G.NE.0) SIGT1G=SQRT(CHSQ1G*AG(2,2))/(2*T1G)
        IF (T2G.NE.0) SIGT2G=SQRT(CHSQ2G*AG(2,2))/(2*T2G)
      ELSE
C
C       PARALLEL MONOCHROMATOR GEOMETRY.  Q = X(1), F*T = X(2), T = X(3)
C
        NPAR=3
        IF (NOBS.LE.NPAR) STOP 'NOT ENOUGH DATA IN SUBROUTINE TFIT'
        CALL MATINV (NPAR,NPAR,AM,DET)
        IF (DET.EQ.0) STOP 'ILL-CONDITIONED MATRIX IN SUBROUTINE TFIT'
        CALL MV (NPAR,AM,B1M,X1M)
        CALL MV (NPAR,AM,B2M,X2M)
        MODEL='PLMC'
        Q1M=X1M(1)
        Q2M=X2M(1)
        F1=X1M(2)/X1M(3)
        F2=X2M(2)/X2M(3)
        F=(F1+F2)/2
        IF (F.LT.1) F=1
        IF (F.GT.2) F=2
        T1=X1M(3)
        T2=X2M(3)
        CALL AGREE1 (MODEL,IO3,Q1M,Q2M,T1,T2,TANTHM,F,WL0,
     &   NOBS,NPAR,CHISQ1,CHISQ2,RMSD1,RMSD2)
        SIGF1=F1*SQRT(CHISQ1*((AM(2,2)/(F1*T1))**2+(AM(3,3)/T1)**2))
        SIGF2=F2*SQRT(CHISQ2*((AM(2,2)/(F2*T2))**2+(AM(3,3)/T2)**2))
        SIGF=SQRT(SIGF1**2+SIGF2**2)/2
        SIGT1=SQRT(CHISQ1*AM(3,3))
        SIGT2=SQRT(CHISQ2*AM(3,3))
      END IF
C
C     CHOOSE MODELS.
C
      IF (TANTHM.EQ.0.OR.RHOM.EQ.90) THEN
        IF (JMODEL.EQ.0) THEN
C
C         FOR EACH HALF-WIDTH, CHOOSE THE MODEL WITH SMALLER ROOT-MEAN-
C         SQUARE ERROR-OF-FIT.
C
          RMSD1=MIN(RMSD1L,RMSD1G)
          RMSD2=MIN(RMSD2L,RMSD2G)
          IF (RMSD1.EQ.RMSD1L) TYPE1='LRNZ'
          IF (RMSD1.EQ.RMSD1G) TYPE1='GAUS'
          IF (RMSD2.EQ.RMSD2L) TYPE2='LRNZ'
          IF (RMSD2.EQ.RMSD2G) TYPE2='GAUS'
        ELSE IF (JMODEL.EQ.1.) THEN
          TYPE1='LRNZ'
          TYPE2='LRNZ'
        ELSE IF (JMODEL.EQ.2.) THEN
          TYPE1='GAUS'
          TYPE2='GAUS'
        END IF
      ELSE
        TYPE1='PLMC'
        TYPE2='PLMC'
      END IF
      IF (TYPE1.EQ.'LRNZ') MODEL1=1
      IF (TYPE1.EQ.'GAUS') MODEL1=2
      IF (TYPE1.EQ.'PLMC') MODEL1=3
      IF (TYPE2.EQ.'LRNZ') MODEL2=1
      IF (TYPE2.EQ.'GAUS') MODEL2=2
      IF (TYPE2.EQ.'PLMC') MODEL2=3
      IF (TYPE1.EQ.'PLMC') RETURN
C
C     EVALUATE AND STORE T AND SIGMA(T) VALUES.
C
C     EXPECTED, LOWER LIMIT VALUES T1 AND T2 WERE CALCULATED FOR
C     TRUNCATED LORENTZIAN SHAPED SPECTRAL LINES IN SUBROUTINE INPUT.
C
      T1E=T1
      T2E=T2
      IF (TYPE1.EQ.'LRNZ') THEN          
        T1=T1L
        SIGT1=SIGT1L
      ELSE IF (TYPE1.EQ.'GAUS') THEN
        T1=T1G
        SIGT1=SIGT1G
      END IF
      IF (TYPE2.EQ.'LRNZ') THEN          
        T2=T2L
        SIGT2=SIGT2L
      ELSE IF (TYPE2.EQ.'GAUS') THEN
        T2=T2G
        SIGT2=SIGT2G
      END IF
      IF (ABS(T1-T1E).GT.1.96*SIGT1.OR.ABS(T2-T2E).GT.1.96*SIGT2)
     & WRITE (ILP,6000) T1E,T1,SIGT1,T2E,T2,SIGT2
      IF (T1.LT.0.OR.T2.LT.0) WRITE (ILP,6001)
      IF (T1.LT.0) T1=0
      IF (T2.LT.0) T2=0
 6000 FORMAT (/'0ATTENTION!  FITTED SPECTRAL LINE-BROADENING COEFFICIENT
     &(S) OF T*TAN(THETA) DIFFER SIGNIFICANTLY FROM EXPECTED VALUE(S).'/
     &' ----------'/'0EXPECTED:  T1 = ',F8.4,'    FITTED:  T1 = ',E10.3,
     &'  SIGMA(T1) = ',E10.3/'            T2 = ',F8.4,'             T2 =
     & ',E10.3,'  SIGMA(T2) = ',E10.3)
 6001 FORMAT (/'0NEGATIVE T*TAN(THETA) COEFFIFIENT(S) HAVE BEEN RESET TO
     & T = 0.')
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE UQFIT
C
C     LEAST-SQUARES FIT OF ADJUSTED DIFFRACTOMETER ORIENTATION MATRIX,
C     U, AND ANISOTROPIC PEAK WIDTH TENSORS, Q.
C
      CHARACTER ATIME*8,ADATE*9,TITLE*80
      COMMON /BLOCK0/ IO1,IO2,IO3,ILP
      COMMON /BLOCK1/ ATIME,ADATE,TITLE
      COMMON /BLOCK2/ GINV(3,3),WLE,WL1,WL0,WL2,T1,T2,SIGT1,SIGT2
      COMMON /BLOCK4/ IDIFF,JTYPE,TANTHM,RHOM,F,U(3,3),MODEL1,MODEL2,
     & Q1(3,3),Q2(3,3)
      DIMENSION H(3),Y(3),Z(3)
      DIMENSION AU1(3,3),BU1(3),AU2(3,3),BU2(3),AU3(3,3),BU3(3)
      DIMENSION U1(3),U2(3),U3(3)
      DIMENSION XX(6),AA(6,6),B1(6),B2(6),XX1(6),XX2(6)
      DIMENSION SIGU(3,3),SIGQ1(3,3),SIGQ2(3,3)
      DIMENSION AU(3,3),QCRYST(3,3),UCRYST(3,3)
      DOUBLE PRECISION AU1,AU2,AU3,AA
      CHARACTER TYPE1*4,TYPE2*4
      IF (MODEL1.EQ.1) TYPE1='LRNZ'
      IF (MODEL1.EQ.2) TYPE1='GAUS'
      IF (MODEL1.EQ.3) TYPE1='PLMC'
      IF (MODEL2.EQ.1) TYPE2='LRNZ'
      IF (MODEL2.EQ.2) TYPE2='GAUS'
      IF (MODEL2.EQ.3) TYPE2='PLMC'
C
C     ZERO LEAST-SQUARES NORMAL MATRICES AND VECTORS.
C
      NPAR=3
      DO 10 I=1,NPAR
      DO 11 J=1,NPAR
      AU1(I,J)=0
      AU2(I,J)=0
      AU3(I,J)=0
   11 CONTINUE
      BU1(I)=0
      BU2(I)=0
      BU3(I)=0
   10 CONTINUE
      NPAR=6
      DO 20 I=1,NPAR
      DO 21 J=I,NPAR
      AA(I,J)=0
   21 CONTINUE
      B1(I)=0
      B2(I)=0
   20 CONTINUE
C
C     LOOP THROUGH DATA FILE OF LEAST-SQUARES OBSERVATIONS TO BUILD
C     NORMAL EQUATIONS.
C
      NOBS=0
      REWIND IO3
    1 READ (IO3,END=9) II,H,TWOTH,OMEGA,CHI,PHI,DSTAR,Y,Z,
     & I1,I2,XM,X0,X1,X2,L1,L2,W1,W2,TANTH
C
C     ACCUMULATE LEAST-SQUARES NORMAL MATRICES AND VECTORS.
C
      NOBS=NOBS+1
C
C     U-MATRIX
C
      NPAR=3
      CALL ACCUM1 (NPAR,H,DSTAR*Y(1),AU1,BU1)
      CALL ACCUM1 (NPAR,H,DSTAR*Y(2),AU2,BU2)
      CALL ACCUM1 (NPAR,H,DSTAR*Y(3),AU3,BU3)
C
C     Q-TENSORS
C
      NPAR=6
      K=0
      DO 5 I=1,3
      DO 5 J=I,3
      K=K+1
      Q=Z(I)*Z(J)
      IF (I.NE.J) Q=2*Q
      XX(K)=Q
    5 CONTINUE
C
C     LORENTZIAN OR GAUSSIAN CONVOLUTION MODEL FOR BETA-FILTERED
C     CHARACTERISTIC X-RAYS OR PERPENDICULAR MONOCHROMATOR GEOMETRY,
C     GAUSSIAN CONVOLUTION MODEL FOR PARALLEL MONOCHROMATOR GEOMETRY
C
      IF (TYPE1.EQ.'LRNZ') Y1=(W1-T1*TANTH)**2
      IF (TYPE1.EQ.'GAUS') Y1=W1**2-(T1*TANTH)**2
      IF (TYPE1.EQ.'PLMC') Y1=W1**2-(-2*F*T1*TANTH/TANTHM+T1*(TANTH/
     & TANTHM)**2)
      IF (TYPE2.EQ.'LRNZ') Y2=(W2-T2*TANTH)**2
      IF (TYPE2.EQ.'GAUS') Y2=W2**2-(T2*TANTH)**2
      IF (TYPE2.EQ.'PLMC') Y2=W2**2-(-2*F*T2*TANTH/TANTHM+T2*(TANTH/
     & TANTHM)**2)
      CALL ACCUM2 (NPAR,XX,Y1,Y2,AA,B1,B2)
C
C     LOOP BACK FOR NEXT REFLECTION.
C
      GO TO 1
C
C     END LOOP THROUGH REFLECTION DATA.
C
    9 CONTINUE
C
C     SOLVE LEAST-SQUARES NORMAL EQUATIONS, CALCULATE STATISTICS OF FIT,
C     AND ESTIMATE PARAMETER ERRORS, ASSUMING COVARIANCES ARE
C     NEGLIGIBLE.
C
C     DIFFRACTOMETER ORIENTATION MATRIX
C
      NPAR=3
      IF (NOBS.LE.NPAR) STOP 'NOT ENOUGH DATA IN SUBROUTINE UQFIT'
      CALL MATINV (NPAR,NPAR,AU1,DET1)
      CALL MATINV (NPAR,NPAR,AU2,DET2)
      CALL MATINV (NPAR,NPAR,AU3,DET3)
      IF (DET1.EQ.0.OR.DET2.EQ.0.OR.DET3.EQ.0)
     & STOP 'ILL-CONDITIONED MATRIX IN SUBROUTINE UQFIT'
      CALL MV (NPAR,AU1,BU1,U1)
      CALL MV (NPAR,AU2,BU2,U2)
      CALL MV (NPAR,AU3,BU3,U3)
      DO 40 I=1,3
      U(I,1)=U1(I)
      U(I,2)=U2(I)
      U(I,3)=U3(I)
   40 CONTINUE
      CALL AGREE0 (IO3,U,NOBS,NPAR,CHISQ1,CHISQ2,CHISQ3,RMSDX0,RMSDX1,
     & RMSDTH)
      DO 41 I=1,3
      SIGU(I,1)=SQRT(AU1(I,I)*CHISQ1)
      SIGU(I,2)=SQRT(AU2(I,I)*CHISQ2)
      SIGU(I,3)=SQRT(AU3(I,I)*CHISQ3)
   41 CONTINUE
C
C     ANISOTROPIC PEAK WIDTH TENSORS
C
      NPAR=6
      IF (NOBS.LE.NPAR) STOP 'NOT ENOUGH DATA IN SUBROUTINE UQFIT'
      CALL MATINV (NPAR,NPAR,AA,DET)
      IF (DET.EQ.0) STOP 'ILL-CONDITIONED MATRIX IN SUBROUTINE UQFIT'
      CALL MV (NPAR,AA,B1,XX1)
      CALL MV (NPAR,AA,B2,XX2)
      K=0
      DO 60 I=1,3
      DO 60 J=I,3
      K=K+1
      Q1(I,J)=XX1(K)
      Q2(I,J)=XX2(K)
      Q1(J,I)=Q1(I,J)
      Q2(J,I)=Q2(I,J)
   60 CONTINUE
      CALL POSDEF (Q1,ILP)
      CALL POSDEF (Q2,ILP)
      CALL AGREE1 (TYPE1,IO3,Q1,Q2,T1,T2,TANTHM,F,WL0,
     & NOBS,NPAR,CHISQ1,DUMMY1,RMSD1,DUMMY2)
      CALL AGREE1 (TYPE2,IO3,Q1,Q2,T1,T2,TANTHM,F,WL0,
     & NOBS,NPAR,DUMMY1,CHISQ2,DUMMY2,RMSD2)
      K=0
      DO 61 I=1,3
      DO 61 J=I,3
      K=K+1
      SIGQ1(I,J)=SQRT(AA(K,K)*CHISQ1)
      SIGQ2(I,J)=SQRT(AA(K,K)*CHISQ2)
      SIGQ1(J,I)=SIGQ1(I,J)
      SIGQ2(J,I)=SIGQ2(I,J)
   61 CONTINUE
C
C     WRITE RESULTS TO BGLP.DAT FILE.
C
      WRITE (IO2,601)    ((U(I,J),J=1,3),I=1,3)
      WRITE (IO2,601) ((SIGU(I,J),J=1,3),I=1,3)
      WRITE (IO2,602) TYPE1
      WRITE (IO2,601)    ((Q1(I,J),J=1,3),I=1,3)
      WRITE (IO2,601) ((SIGQ1(I,J),J=1,3),I=1,3)
      WRITE (IO2,601)    T1
      WRITE (IO2,601) SIGT1
      WRITE (IO2,602) TYPE2
      WRITE (IO2,601)    ((Q2(I,J),J=1,3),I=1,3)
      WRITE (IO2,601) ((SIGQ2(I,J),J=1,3),I=1,3)
      WRITE (IO2,601)    T2
      WRITE (IO2,601) SIGT2
      WRITE (IO2,601) F
C
C     IF R.M.S. ERROR IN CALCULATED PEAK CENTROID POSITIONS EXCEEDS
C     HALF THE R.M.S. PEAK DISPLACEMENT, SET DEFAULT PEAK POSITION TO
C     SCAN MIDPOINT.
C
      IF (RMSDX1.GT.0.5*RMSDX0) WRITE (IO2,603) 0.5,0.0,0.0
  601 FORMAT (3(1X,3E15.8/))
  602 FORMAT (1X,A)
  603 FORMAT (1X,3F10.2)
      CLOSE (UNIT=IO2,DISPOSE='SAVE')
C
C     PRINT RESULTS.
C
      WRITE (ILP,610) ATIME,ADATE,TITLE
      NPAR=9
      WRITE (ILP,620) NOBS,NPAR,RMSDX0,RMSDX1,RMSDTH
      WRITE (ILP,610) ATIME,ADATE,TITLE
      NPAR=7
      IF (TYPE1.EQ.'PLMC') NPAR=8
      WRITE (ILP,634) NOBS,NPAR,TYPE1,RMSD1,TYPE2,RMSD2,RMSDTH
      WRITE (ILP,610) ATIME,ADATE,TITLE
      WRITE (ILP,621) ((U(I,J),J=1,3),I=1,3),((SIGU(I,J),J=1,3),I=1,3)
      WRITE (ILP,610) ATIME,ADATE,TITLE
      WRITE (ILP,637)
      IF (TYPE1.EQ.'LRNZ') WRITE (ILP,647)
      IF (TYPE1.EQ.'GAUS') WRITE (ILP,648)
      IF (TYPE1.EQ.'PLMC') WRITE (ILP,658) F
      WRITE (ILP,632)    ((Q1(I,J),J=1,3),I=1,3)
      WRITE (ILP,633) ((SIGQ1(I,J),J=1,3),I=1,3)
      WRITE (ILP,643)    T1
      WRITE (ILP,644) SIGT1
      WRITE (ILP,638)
      IF (TYPE2.EQ.'LRNZ') WRITE (ILP,647)
      IF (TYPE2.EQ.'GAUS') WRITE (ILP,648)
      IF (TYPE2.EQ.'PLMC') WRITE (ILP,658) F
      WRITE (ILP,632)    ((Q2(I,J),J=1,3),I=1,3)
      WRITE (ILP,633) ((SIGQ2(I,J),J=1,3),I=1,3)
      WRITE (ILP,643)    T2
      WRITE (ILP,644) SIGT2
      WRITE (ILP,639)
C
C     TRANSFORM PEAK WIDTH TENSORS TO CRYSTALLOGRAPHIC RECIPROCAL
C     LATTICE AXES, AND FIND EIGENVALUES AND EIGENVECTORS.
C
      WRITE (ILP,610) ATIME,ADATE,TITLE
      DO 19 I=1,3
      AI=1/SQRT(GINV(I,I))
      DO 19 J=1,3
      AU(I,J)=AI*U(I,J)
   19 CONTINUE
      CALL UQUT (AU,Q1,QCRYST)
      WRITE (ILP,641) ((QCRYST(I,J),J=1,3),I=1,3)
      CALL EIGEN (3,3,QCRYST,UCRYST)
      WRITE (ILP,642) (QCRYST(I,I),(UCRYST(J,I),J=1,3),I=1,3)
      CALL UQUT (AU,Q2,QCRYST)
      WRITE (ILP,653) ((QCRYST(I,J),J=1,3),I=1,3)
      CALL EIGEN (3,3,QCRYST,UCRYST)
      WRITE (ILP,642) (QCRYST(I,I),(UCRYST(J,I),J=1,3),I=1,3)
      WRITE (ILP,654)
  610 FORMAT ('1'/'0PROGRAM REFPK.  ',A,1X,A,'.  ',A)
  620 FORMAT ('0STATISTICS OF FIT FOR PEAK POSITIONS:'/
     &'0RMSD0 = SQRT(SUM((X0 - XM)**2)/NOBS)'/
     &' RMSD1 = SQRT(SUM((X0 - X0CALC)**2)/(NOBS - NPAR))'/
     &'0XM IS THE SCAN MIDPOINT.'/
     &' X0 IS THE INTENSITY-WEIGHTED CENTROID OF THE PEAK ABOVE BACKGROU
     &ND.'/
     &' X0CALC IS THE PEAK POSITION CALCULATED FROM THE REFLECTION INDIC
     &ES AND THE RECALCULATED ORIENTATION MATRIX.'/
     &'0NOBS  = ',I6/
     &' NPAR  = ',I6/
     &'0RMSD0 = ',F6.3,' DEGREES THETA'/
     &' RMSD1 = ',F6.3,' DEGREES THETA'/
     &'0ROOT-MEAN-SQUARE STEP WIDTH:'/
     &'0RMSDX = ',F6.3,' DEGREES THETA')
  634 FORMAT ('0STATISTICS OF FIT FOR PEAK WIDTHS:'/
     &'0RMSD(I) = SQRT[SUM([WOBS(I) - WCALC(I)]**2)/(NOBS - NPAR)]'/
     &'0I = 1 OR 2 FOR THE LOWER ANGLE OR HIGHER ANGLE HALF-PEAK, RESPEC
     &TIVELY.'/
     &'0THE WOBS ARE FROM THE LEHMANN-LARSEN PEAK LIMITS.'/
     &' THE WCALC ARE FROM THE FITTED PEAK WIDTH TENSORS.'/
     &'0NOBS = ',I6/
     &' NPAR = ',I6/
     &'0LOWER ANGLE HALF-PEAK, W(1):'/
     &'0MODEL(1) = ',1X,A/
     &'  RMSD(1) = ',F6.3,' DEGREES THETA'/
     &'0HIGHER ANGLE HALF-PEAK, W(2):'/
     &'0MODEL(2) = ',1X,A/
     &'  RMSD(2) = ',F6.3,' DEGREES THETA'/
     &'0ROOT-MEAN-SQUARE STEP WIDTH:'/
     &'0RMSDX    = ',F6.3,' DEGREES THETA')
  621 FORMAT ('0PEAK POSITION MATRIX'/
     &'0DIFFRACTOMETER ORIENTATION MATRIX RECALCULATED BY LEAST-SQUARES'
     &/' FIT TO OMEGA VALUES FOR THE INTENSITY-WEIGHTED PEAK CENTROIDS'/
     &' ALONG WITH THE DIFFRACTOMETER ANGLES CHI AND PHI.'/'0DIFFRACTOME
     &TER AXES AS DEFINED BY WALTER  C. HAMILTON (1974).'/' INTERNATIONA
     &L TABLES FOR X-RAY CRYSTALLOGRAPHY, VOL. IV, PP. 273-284.'/
     &'0U(1,1)    U(1,2)    U(1,3)      ',3F12.8/
     &' U(2,1)    U(2,2)    U(2,3)    = ',3F12.8/
     &' U(3,1)    U(3,2)    U(3,3)      ',3F12.8/
     &'0SIG(1,1)  SIG(1,2)  SIG(1,3)    ',3F12.8/
     &' SIG(2,1)  SIG(2,2)  SIG(2,3)  = ',3F12.8/
     &' SIG(3,1)  SIG(3,2)  SIG(3,3)    ',3F12.8)
  637 FORMAT ('0PEAK WIDTH TENSORS'/'0LOWER ANGLE HALF-PEAK, W1')
  638 FORMAT ('0HIGHER ANGLE HALF-PEAK, W2')
  632 FORMAT (
     &'0Q(1,1)    Q(1,2)    Q(1,3)      ',3F8.4/
     &' Q(2,1)    Q(2,2)    Q(2,3)    = ',3F8.4/
     &' Q(3,1)    Q(3,2)    Q(3,3)      ',3F8.4)
  633 FORMAT (
     &'0SIG(1,1)  SIG(1,2)  SIG(1,3)    ',3F8.4/
     &' SIG(2,1)  SIG(2,2)  SIG(2,3)  = ',3F8.4/
     &' SIG(3,1)  SIG(3,2)  SIG(3,3)    ',3F8.4)
  643 FORMAT ('0                    T         = ',F8.4)
  644 FORMAT ('0                    SIG       = ',F8.4)
  647 FORMAT ('0LORENTZIAN MODEL'/
     &'0W = (Z(K)*Q(J,K)*Z(J))**(1/2) + T*TAN(THETA)')
  648 FORMAT ('0GAUSSIAN MODEL'/
     &'0W = (Z(K)*Q(J,K)*Z(J) + (T*TAN(THETA))**2)**(1/2)')
  658 FORMAT ('0GAUSSIAN MODEL FOR PARALLEL MONOCHROMATOR GEOMETRY'/
     &'0W = (Z(K)*Q(J,K)*Z(J) - 2*F*T*TAN(THETA)/TAN(THETA(M)) + T*(TAN(
     &THETA)/TAN(THETA(M)))**2)**(1/2),  F = ',F4.1)
  639 FORMAT ('0W1 AND W2 ARE THE BASE WIDTHS OF THE HALF-PEAKS BELOW TH
     &ETA(ALPHA-1)'/' AND ABOVE THETA(ALPHA-2), RESPECTIVELY, IN UNITS O
     &F DEGREES OF CRYSTAL'/' ROTATION (THETA).'/'0THE Z(3) ARE THE COMP
     &ONENTS OF A UNIT VECTOR NORMAL TO THE INCIDENT'/' AND DIFFRACTED B
     &EAM VECTORS.'/'0THE Q(3,3) ARE THE ELEMENTS OF SYMMETRIC TENSORS T
     &HAT ARE RELATED TO'/' THE ANISOTROPY OF THE SIZE AND MOSAICITY O',
     &'F THE SPECIMEN CRYSTAL.'/'0THE T ARE (SCALAR) COEFFICIENTS PROPOR
     &TIONAL TO THE SPECTRAL BAND'/' WIDTH OF THE CHARACTERISTIC ALPHA-1
     & AND ALPHA-2 LINES.'/'0THE VECTORS Z AND TENSORS Q ARE REFERRED TO
     & CRYSTAL-FIXED CARTESIAN AXES'/' THAT ARE COINCIDENT WITH THE DIFF
     &RACTOMETER AXES AT OMEGA = PHI = CHI = 0.')
  641 FORMAT ('0ANISOTROPIC PEAK WIDTH TENSORS TRANSFORMED TO NORMALIZED
     &'/' CRYSTALLOGRAPHIC RECIPROCAL LATTICE AXES'/'0  Q(CRYST) = (A*U)
     &*Q*(A*U)-TRANSPOSE'/'0  Z(CRYST) = Z*(A*U)-INVERSE,'/'0WHERE THE M
     &ATRIX'/'0      (1/ASTAR 0       0      )'/'   A = (0       1/BSTAR
     & 0      )'/'       (0       0       1/CSTAR)'/'0NORMALIZES TO BASI
     &S VECTORS OF UNIT LENGTH.'//'0LOWER ANGLE HALF-PEAK, Q1'/'0QCRYST(
     &I,J)'//3(3F8.4/))
  642 FORMAT (' EIGENVALUE(I)  EIGENVECTOR(I,J), (I DOWN, J ACROSS)'//
     &3(' 'F8.4,6X,3F8.4/))
  653 FORMAT ('0HIGHER ANGLE HALF-PEAK, Q2'/'0QCRYST(I,J)'//3(3F8.4/))
  654 FORMAT ('0NARROWEST REFLECTION PEAK OCCURS WHEN THE EIGENVECTOR OF
     & THE SMALLEST'/' EIGENVALUE IS PERPENDICULAR TO THE EQUATORIAL PLA
     &NE, AND WIDEST PEAK'/' WHEN THE EIGENVECTOR OF THE LARGEST EIGENVA
     &LUE IS PERPENDICULAR.'/'0WITH EIGENVECTOR(I,J) PERPENDICULAR TO TH
     &E EQUATORIAL PLANE, THE WIDTH'/' OF THE REFLECTION HALF-PEAK IS SQ
     &RT(EIGENVALUE(I)) (DEGREES THETA).')
  651 FORMAT ('0FAILURE.  NOBS = ',I2,' .LE. NPAR = ',I2)
  652 FORMAT ('0FAILURE.  DET  =  0.  NORMAL MATRIX SINGULAR.')
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE UQUT (U,Q,T)
C
C     MATRIX TRIPLE PRODUCT, T = U*Q*U-TRANSPOSE
C
      DIMENSION U(3,3),Q(3,3),T(3,3)
      DO 1 I=1,3
      DO 1 J=1,3
      T(I,J)=0
      DO 1 K=1,3
      DO 1 L=1,3
      T(I,J)=T(I,J)+U(I,K)*Q(K,L)*U(J,L)
    1 CONTINUE
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE XCHECK (X1,X0,X2,X,Y,I1,I2)
C
C     FOR WEAK, HIGH-ANGLE REFLECTIONS, SUBROUTINE CENTROID SOMETIMES
C     FINDS X0 TOO CLOSE TO THE ALPHA-ONE PEAK.  ENSURE THAT X1 IS AT
C     THE ALPHA-ONE PEAK.
C
      DIMENSION X(1),Y(1)
      D=(X(I2)-X(I1))/(I2-I1+1)
      IF (X0-X1.LE.D) RETURN
      J1=NINDEX(X1)
      J0=NINDEX(X0)
      J1=MAX(J1,I1+1)
      J0=MIN(J0,I2-1)
      IF (J0-J1.LE.1) RETURN
      YMAX=0
      DO 1 I=J1,J0
      YI=Y(I-1)+Y(I)+Y(I+1)
      IF (YI.GT.YMAX) THEN
        YMAX=YI
        IMAX=I
      END IF
    1 CONTINUE
      D=X(IMAX)-X1
      IF (D.GT.0) THEN
        X1=X1+D
        X0=X0+D
        X2=X2+D
      END IF
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE XYDATA (WIDTH,SPEED,NSTEP,XX,YY,VX,VY)
C
C     CONVERT FROM STEP-SCAN COUNT PROFILE TO COUNT-RATE VERSUS ROCKING-
C     ANGLE PROFILE, AND APPLY COINCIDENCE CORRECTIONS FOR COUNTING
C     LOSSES.
C
C     RETURN:
C
C         ROCKING-ANGLE STEPS XX(NSTEP) IN DEGREES (THETA) INCREASING
C         (FROM NEGATIVE TO POSITIVE) IN ORDER OF INCREASING ABSOLUTE
C         VALUE OF THE BRAGG ANGLE, THETA;
C
C         COUNT-RATE STEPS YY(NSTEP) IN COUNTS PER SECOND;
C
C         VARIANCE STEPS VX(NSTEP) AND VY(NSTEP).
C
C     THIS VERSION OF SUBROUTINE XYDATA TREATS OMEGA, OMEGA/THETA, OR
C     OMEGA/TWO-THETA STEP SCANS WITH:
C
C         NSTEP EQUALLY SPACED STEPS,
C         -----
C         TOTAL SCAN WIDTH IN DEGREES OF OMEGA ROTATION,
C                    -----
C         SCAN SPEED IN DEGREES (OMEGA) PER MINUTE.
C              -----
C
C     NEGATIVE SCAN WIDTH IS A FLAG THAT THE ATTENUATOR WAS USED.
C
      DIMENSION XX(1),YY(1),VX(1),VY(1)
      COMMON /BLOCK7/ VARTHE,VARTIM,TAU,VARTAU,ATT,VARATT
C
C     ROCKING-ANGLE ABSCISSAE
C
      X0=0.5*(1+NSTEP)
      DX=ABS(WIDTH)/NSTEP
      VARD=2*VARTHE/NSTEP**2
      DO 1 I=1,NSTEP
      XX(I)=(I-X0)*DX
      VX(I)=(I-X0)**2*VARD
    1 CONTINUE
C
C     COUNT-RATE ORDINATES
C
      SPEED=SPEED/60
      T=(ABS(WIDTH)/NSTEP)/SPEED
      VART=2*(2*VARTHE/SPEED**2+VARTIM)/NSTEP**2
      DO 2 I=1,NSTEP
      Y=YY(I)
C
C     ALLOW FOR STEPS WITH VERY FEW OR NO COUNTS.
C
      Y=Y+0.5
      VARY=Y
      X=Y/T
      VARX=(VARY+X**2*VART)/T**2
C
C     DEAD TIME CORRECTION
C
      Y=X/(1-TAU*X)
      VARY=Y**4*(VARX/X**4+VARTAU)
C
C     ATTENUATOR CORRECTION
C
      IF (WIDTH.LT.0) THEN
        X=Y
        VARX=VARY
        Y=ATT*X
        VARY=X**2*VARATT+ATT**2*VARX
      END IF
      YY(I)=Y
      VY(I)=VARY
    2 CONTINUE
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE YZ (DELTA,TWOTH,OMEGA,CHI,PHI,Y,Z)
C
C     GET THE COMPONENTS ALONG CARTESIAN AXES FIXED IN THE CRYSTAL OF A
C     UNIT VECTOR, Y, PARALLEL TO THE RECIPROCAL LATTICE DIFFRACTION
C     VECTOR, DSTAR, AND OF A UNIT VECTOR, Z, NORMAL TO THE EQUATORIAL
C     PLANE OF THE INCIDENT AND DIFFRACTED BEAM VECTORS.
C
C     AXES AS DEFINED BY WALTER C. HAMILTON (1974).  INTERNATIONAL
C     TABLES FOR X-RAY CRYSTALLOGRAPHY, VOL. IV, PP. 273-284.
C
C     THE CRYSTAL-FIXED CARTESIAN AXES ARE COINCIDENT WITH THE
C     DIFFRACTOMETER AXES AT OMEGA = PHI = CHI = 0.
C
C     DELTA IS THE DISPLACEMENT OF THE PEAK CENTROID FROM THE SCAN
C     MIDPOINT.
C
C     THE POSITIVE SENSE OF THE OMEGA ROTATION IS DEFINED SUCH THAT
C     DELTA(OMEGA) = -DELTA(X0).
C
      DIMENSION Y(3),Z(3)
      DATA RAD /0.017453293/
      IF (TWOTH.LT.0) DELTA=-DELTA
      COSOM=COS((OMEGA-DELTA)*RAD)
      SINOM=SIN((OMEGA-DELTA)*RAD)
      COSPH=COS(PHI*RAD)
      SINPH=SIN(PHI*RAD)
      COSCH=COS(CHI*RAD)
      SINCH=SIN(CHI*RAD)
      Y(1)= COSPH*SINOM+SINPH*COSCH*COSOM
      Y(2)=-SINPH*SINOM+COSPH*COSCH*COSOM
      Y(3)=-SINCH*COSOM
      IF (TWOTH.LT.0) THEN
        DO 1 I=1,3
        Y(I)=-Y(I)
    1   CONTINUE
      END IF
      Z(1)=SINPH*SINCH
      Z(2)=COSPH*SINCH
      Z(3)=COSCH
      RETURN
      END
C-----------------------------------------------------------------------
