      PROGRAM BGLP
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 /BLOCKC/ NCOND,ICOND(38)
      COMMON /BLOCK1/ ATIME,ADATE,TITLE
      COMMON /BLOCK2/ GINV(3,3),WLE,WL1,WL0,WL2
      COMMON /BLOCK3/ WINDOW,CUTOFF,NSMOOTH,II1,II2,XT1,XT2,SC1,SC2,
     & XX0,XW1,XW2
      COMMON /BLOCK4/ C0KIN,C1KIN,C2KIN,C0DYN,C1DYN,C2DYN,FRACTD,VARFD,
     & INEUTRON
      COMMON /BLOCK5/ IDIFF,U(3,3),MODEL1,MODEL2,Q1(3,3),Q2(3,3),T1,T2,
     & F,SIGU(3,3),SIGQ1(3,3),SIGQ2(3,3),SIGT1,SIGT2,TANTHM
      COMMON /BLOCK6/ VARTHE,VARTIM,TAU,VARTAU,ATT,VARATT
      COMMON /BLOCK7/ NSTEP,XX(NMAX),YY(NMAX),VX(NMAX),VY(NMAX)
      COMMON /BLOCK8/ TH0,X0,SIGX0,XE,X1,X2,W1,W2,I1,I2,L1,L2,R,SIGR,EPS
      CALL INPUT
      CALL PROCESS
      CALL OUTPUT
      STOP 'PROGRAM BGLP FINIS!'
      END
C-----------------------------------------------------------------------
      SUBROUTINE BGLP1
C
C     LINEAR BACKGROUND
C
      PARAMETER (NMAX=100)
      DIMENSION X(NMAX),Y(NMAX)
      COMMON /BLOCK7/ NSTEP,XX(NMAX),YY(NMAX),VX(NMAX),VY(NMAX)
      COMMON /BLOCK8/ TH0,X0,SIGX0,XE,X1,X2,W1,W2,I1,I2,L1,L2,R,SIGR,EPS
      DATA RAD /0.017453293/
C
C     FIT A STRAIGHT LINE TO THE PROFILE IN THE TWO BACKGROUND REGIONS,
C     REQUIRING THAT, FOR THE PURPOSE OF FITTING THE LINE, EACH
C     BACKGROUND SAMPLE INCLUDES AT LEAST TWO POINTS.
C
      J1=MAX(L1-1,I1+1)
      J2=MIN(L2+1,I2-1)
      N=0
      DO 1 I=I1,J1
      N=N+1
      X(N)=XX(I)
      Y(N)=YY(I)
    1 CONTINUE
      DO 2 I=J2,I2
      N=N+1
      X(N)=XX(I)
      Y(N)=YY(I)
    2 CONTINUE
      CALL FITLINE (N,X,Y,C0,C1,V0,V1,COV)
C
C     SUBTRACT LINEAR BACKGROUND POINT-BY-POINT, DIVIDE OUT LORENTZ AND
C     POLARIZATION FACTORS POINT-BY-POINT, AND INTEGRATE THE PEAK
C     PROFILE BETWEEN PEAK LIMITS.
C
      R=0
      VR=0
      DO 9 I=MAX(I1+1,L1),MIN(L2,I2-1)
C
C     BACKGROUND SUBTRACTION
C
      XI=XX(I)
      YI=YY(I)-C0-C1*XI
      VI=VY(I)+V0+V1*XI**2+2*XI*COV
C
C     LORENTZ AND POLARIZATION CORRECTIONS
C
      TH=(XX(I)+TH0)*RAD
      VT=VX(I)*RAD**2
      CALL LZ (TH,VT,FL,VL)
      CALL PZ (TH,VT,FP,VP)
      CV=SQRT(VL*VP)
      YJ=YI
      IF (YI.EQ.0) YI=0.5*SQRT(VI)
      VI=(VI/YI**2+VL/FL**2+VP/FP**2+CV/(FL*FP))*(YI/(FL*FP))**2
      YI=YJ/(FL*FP)
C
C     INTEGRATED NET INTENSITY BY TRAPEZOIDAL INTEGRATION
C
      DX=0.50*(XX(I+1)-XX(I-1))
      VD=0.25*(VX(I+1)+VX(I-1))
      R=R+YI*DX
      VR=VR+YI**2*VD+DX**2*VI
    9 CONTINUE
      SIGR=SQRT(VR)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE BGLP2
C
C     STRUCTURED BACKGROUND UNDER LOW ANGLE REFLECTIONS DUE TO THE
C     ABSORPTION EDGE OF THE BETA FILTER AND THE PEAK IN THE WHITE
C     RADIATION BACKGROUND
C
C     R. J. NELMES (1975).  ACTA CRYST. A31, 273-279.
C
      PARAMETER (NMAX=100)
      DIMENSION X(NMAX),Y(NMAX)
      COMMON /BLOCK7/ NSTEP,XX(NMAX),YY(NMAX),VX(NMAX),VY(NMAX)
      COMMON /BLOCK8/ TH0,X0,SIGX0,XE,X1,X2,W1,W2,I1,I2,L1,L2,R,SIGR,EPS
      DATA RAD /0.017453293/
C
C     FIT STRAIGHT LINE SEGMENTS TO THE STRUCTURED BACKGROUND.
C
      A0=0
      A1=0
      VA0=0
      VA1=0
      COVA=0
      E0=0
      E1=0
      VE0=0
      VE1=0
      COVE=0
      B0=0
      B1=0
      VB0=0
      VB1=0
      COVB=0
      J1=NINDEX(XE-EPS)
      J2=NINDEX(XE+EPS)
      IF (J2.LE.I1) GO TO 99
C
C     FIT ONE STRAIGHT LINE TO THE LOW ANGLE BACKGROUND BELOW THE TAIL
C     OF THE BETA FILTER ABSORPTION EDGE.
C
      J1=MIN(J1,L1-1)
      J1=MAX(J1,I1)
      N=0
      DO 1 I=I1,J1
      N=N+1
      X(N)=XX(I)
      Y(N)=YY(I)
    1 CONTINUE
      IF (N.GE.3) THEN
        CALL FITLINE (N,X,Y,A0,A1,VA0,VA1,COVA)
      ELSE IF (N.EQ.2) THEN
        A0=0.5*(Y(1)+Y(2))
        VA0=0.5*A0
      ELSE
        A0=Y(1)
        VA0=A0
      END IF
C
C     FIT A SECOND STRAIGHT LINE TO THE ABSORPTION EDGE.
C
      J2=MIN(J2,L1-1)
      J2=MAX(J2,J1)
      N=0
      DO 2 I=J1,J2
      N=N+1
      X(N)=XX(I)
      Y(N)=YY(I)
    2 CONTINUE
      IF (N.GE.3) THEN
        CALL FITLINE (N,X,Y,E0,E1,VE0,VE1,COVE)
      ELSE IF (N.EQ.2) THEN
        E0=0.5*(Y(1)+Y(2))
        VE0=0.5*E0
      ELSE
        E0=Y(1)
        VE0=E0
      END IF
C
C     THE ABSORPTION EDGE SHOULD HAVE POSITIVE INTERCEPT AND SLOPE.
C
      IF (E0.LT.0.OR.E1.LT.0) GO TO 99
C
C     FIT A THIRD STRAIGHT LINE TO THE HIGH ANGLE, WHITE RADIATION
C     BACKGROUND ABOVE THE TAIL OF THE BRAGG PEAK.
C
      J2=MAX(J2,L2+1)
      J2=MIN(J2,I2)
      N=0
      DO 3 I=J2,I2
      N=N+1
      X(N)=XX(I)
      Y(N)=YY(I)
    3 CONTINUE
      IF (N.GE.3) THEN
        CALL FITLINE (N,X,Y,B0,B1,VB0,VB1,COVB)
      ELSE IF (N.EQ.2) THEN
        B0=0.5*(Y(1)+Y(2))
        VB0=0.5*B0
      ELSE
        B0=Y(1)
        VB0=B0
      END IF
C
C     THE HIGH ANGLE BACKGROUND SHOULD HAVE POSITIVE INTERCEPT AND
C     NEGATIVE SLOPE.
C
      IF (B0.LT.0.OR.B1.GT.0) GO TO 99
C
C     EVALUATE THE LOW ANGLE BACKGROUND AT XE - EPSILON AND THE HIGH
C     ANGLE BACKGROUND AT XE + EPSILON.
C
      J1=NINDEX(XE-EPS)
      J2=NINDEX(XE+EPS)
      J1=MAX(J1,I1)
      J2=MIN(J2,I2)
      X1=XX(J1)
      X2=XX(J2)
      Y1=A0+A1*X1
      Y2=B0+B1*X2
      V1=VA0+VA1*X1**2+X1*COVA
      V2=VB0+VB1*X2**2+X2*COVB
C
C     THE BACKGROUND AT XE - EPSILON SHOULD BE LOWER THAN THE
C     BACKGROUND AT XE + EPSILON.
C
      IF (Y1.GE.Y2) GO TO 99
C
C     CALCULATE THE COEFFICIENTS OF THE LINE JOINING THE LOW ANGLE
C     BACKGROUND AT X1 = XE - EPSILON TO THE HIGH ANGLE BACKGROUND AT
C     X2 = XE + EPSILON.
C
      IF (J1.LT.L1.AND.L1.LT.J2) THEN
C
C       THE LOW ANGLE PEAK LIMIT FALLS ON THE ABSORPTION EDGE.
C
        X1=XX(L1-1)
        Y1=E0+E1*X1
        V1=VE0+VE1*X1**2+X1*COVE
      END IF
      E1=(Y2-Y1)/(X2-X1)
      E0=Y1-E1*X1
      VE1=(V1+V2)/(X2-X1)**2
      VE0=(V1*X2**2+V2*X1**2)/(X2-X1)**2
      COVE=-(V1*X2+V2*X1)/(X2-X1)**2
C
C     SUBTRACT BACKGROUND POINT-BY-POINT, DIVIDE OUT LORENTZ AND 
C     POLARIZATION FACTORS POINT-BY-POINT, AND INTEGRATE THE PEAK
C     PROFILE BETWEEN PEAK LIMITS.
C
      R=0
      VR=0
      DO 9 I=MAX(I1+1,L1),MIN(L2,I2-1)
C
C     BACKGROUND SUBTRACTION
C
      XI=XX(I)
      IF (XI.LT.X1) THEN
        C0=A0
        C1=A1
        V0=VA0
        V1=VA1
        COV=COVA
      ELSE IF (X1.LE.XI.AND.XI.LE.X2) THEN
        C0=E0
        C1=E1
        V0=VE0
        V1=VE1
        COV=COVE
      ELSE
        C0=B0
        C1=B1
        V0=VB0
        V1=VB1
        COV=COVB
      END IF
      YI=YY(I)-C0-C1*XI
      VI=VY(I)+V0+V1*XI**2+2*XI*COV
C
C     PREVENT TRUNCATION OF THE LOW ANGLE TAIL OF THE PEAK BY THE
C     ESTIMATED ABSORPTION EDGE.
C
      IF (XI.LT.0.AND.YI.LT.0) YI=0
C
C     LORENTZ AND POLARIZATION CORRECTIONS
C
      TH=(XX(I)+TH0)*RAD
      VT=VX(I)*RAD**2
      CALL LZ (TH,VT,FL,VL)
      CALL PZ (TH,VT,FP,VP)
      CV=SQRT(VL*VP)
      YJ=YI
      IF (YI.EQ.0) YI=0.5*SQRT(VI)
      VI=(VI/YI**2+VL/FL**2+VP/FP**2+CV/(FL*FP))*(YI/(FL*FP))**2
      YI=YJ/(FL*FP)
C
C     INTEGRATED NET INTENSITY BY TRAPEZOIDAL INTEGRATION
C
      DX=0.50*(XX(I+1)-XX(I-1))
      VD=0.25*(VX(I+1)+VX(I-1))
      R=R+YI*DX
      VR=VR+YI**2*VD+DX**2*VI
    9 CONTINUE
      SIGR=SQRT(VR)
      IF (R.LT.CUTOFF*SIGR) GO TO 99
      RETURN
   99 CONTINUE
C
C     IF STRUCTURED BACKGROUND MODEL IS NOT APPROPRIATE, RESORT TO
C     LINEAR BACKGROUND.
C
      J2=MIN(NINDEX(XE+EPS),L1-1)
      I1=MAX(I1,J2)
      CALL BGLP1
      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 DPRINT (JJ,JH,JK,JL,JPRINT)
C
C     JPRINT = 1 FOR H00, 0K0, 00L, HH0, 0KK, H0H, AND HHH;
C     JPRINT = 0 OTHERWISE.
C
      IF (JJ.LT.0) GO TO 10
      J=ABS(JH)
      K=ABS(JK)
      L=ABS(JL)
      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 JPRINT=0
      RETURN
   11 JPRINT=1
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE EXTREM (IXMIN,IXMAX,IX)
      PARAMETER (N1=4,N2=12,N3=50)
      DIMENSION IXMIN(N1,N2,N3),IXMAX(N1,N2,N3),ITYPE(N1),IX(N2)
      DATA ITYPE /5, 6, 8, 12/
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(1),V(1)
      L=NINT(WINDOW*N)
      DO 1 I=L+1,N-L
      YB=0
      YI=0
      VB=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(1),Y(1)
      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
C
C     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       KAPPA IS MECHANICALLY LIMITED TO THE RANGE 0 TO 180 DEGREES.
C       WITH ALPHA = 50 DEGREES, THIS LIMITS CHI TO THE RANGE -100 TO
C       +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
        I=0
      ELSE
        H=A(I)
      END IF
      RETURN
      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,L+H=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         IABSENT = 1 FOR A FORBIDDEN REFLECTION.
C
      INTEGER H
      DIMENSION ICOND(38)
      IABSENT=0
      IF (NCOND.EQ.0) RETURN
      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.OR.MOD(L+H,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
      RETURN
   91 IABSENT=1
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE INPUT
C
C         READ CONTROL DATA AND EVALUATE VARIABLES.
C
      CHARACTER ATIME*8,ADATE*9,TITLE*80,FILE1*80,FILE2*80,HCOND*18,
     & ZTIME*8,ZDATE*9,TYPE1*4,TYPE2*4
      COMMON /BLOCK0/ IO1,IO2,IO3,ILP
      COMMON /BLOCKC/ NCOND,ICOND(38)
      COMMON /BLOCK1/ ATIME,ADATE,TITLE
      COMMON /BLOCK2/ GINV(3,3),WLE,WL1,WL0,WL2
      COMMON /BLOCK3/ WINDOW,CUTOFF,NSMOOTH,II1,II2,XT1,XT2,SC1,SC2,
     & XX0,XW1,XW2
      COMMON /BLOCK4/ C0KIN,C1KIN,C2KIN,C0DYN,C1DYN,C2DYN,FRACTD,VARFD,
     & INEUTRON
      COMMON /BLOCK5/ IDIFF,U(3,3),MODEL1,MODEL2,Q1(3,3),Q2(3,3),T1,T2,
     & F,SIGU(3,3),SIGQ1(3,3),SIGQ2(3,3),SIGT1,SIGT2,TANTHM
      COMMON /BLOCK6/ VARTHE,VARTIM,TAU,VARTAU,ATT,VARATT
      DIMENSION CELL(6)
      PI=ACOS(-1.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='bglp.dat',STATUS='OLD')
      READ (IO1,500) ZTIME
      READ (IO1,500) ZDATE
      READ (IO1,500) TITLE
      READ (IO1,500) FILE1
      READ (IO1,500) FILE2
      READ (IO1,501) CELL
      READ (IO1,502) IDIFF,INEUTRON
      READ (IO1,501) THM,RHOM,FRACTD,SIGFD
      READ (IO1,501) WLE,WL1,WL0,WL2
      READ (IO1,501) SIGTHE,SIGTIM
      READ (IO1,501) TAU,SIGTAU,ATT,SIGATT
      READ (IO1,501) WINDOW
      READ (IO1,501) CUTOFF
      READ (IO1,502) NSMOOTH
      READ (IO1,505) SC1,SC2
      READ (IO1,503) II1,II2,XT1,XT2
      READ (IO1,504)    ((U(I,J),J=1,3),I=1,3)
      READ (IO1,504) ((SIGU(I,J),J=1,3),I=1,3)
      READ (IO1,500) TYPE1
      READ (IO1,504)    ((Q1(I,J),J=1,3),I=1,3)
      READ (IO1,504) ((SIGQ1(I,J),J=1,3),I=1,3)
      READ (IO1,504)    T1
      READ (IO1,504) SIGT1
      READ (IO1,500) TYPE2
      READ (IO1,504)    ((Q2(I,J),J=1,3),I=1,3)
      READ (IO1,504) ((SIGQ2(I,J),J=1,3),I=1,3)
      READ (IO1,504)    T2
      READ (IO1,504) SIGT2
      F=0
      XX0=0.5
      XW1=0
      XW1=0
      READ (IO1,504,END=9) F
      READ (IO1,505,END=9) XX0,XW1,XW2
 500  FORMAT (1X,A)
 501  FORMAT (1X,6F10.5)
 502  FORMAT (1X,2I10)
 503  FORMAT (1X,2I10,2F10.3)
 504  FORMAT (3(1X,3E15.8/))
 505  FORMAT (1X,4F10.2)
 9    CLOSE (UNIT=IO1)
      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
C
C         PRINT CONTROL DATA.
C
      OPEN (UNIT=ILP,FILE='bglp.lp',STATUS='NEW',CARRIAGECONTROL=
     & 'FORTRAN')
      WRITE (ILP,610) ATIME,ADATE,TITLE
      WRITE (ILP,600) ZTIME,ZDATE,TITLE
      WRITE (ILP,611) FILE1
      WRITE (ILP,612) FILE2
      WRITE (ILP,629) 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)
      IF (INEUTRON.EQ.0) WRITE (ILP,651)
      IF (INEUTRON.EQ.1) WRITE (ILP,652)
      IF (THM.NE.0) WRITE (ILP,628) THM,RHOM,FRACTD,SIGFD
      WRITE (ILP,631) WLE,WL1,WL0,WL2
      WRITE (ILP,613) SIGTHE,SIGTIM
      WRITE (ILP,620) TAU,SIGTAU,ATT,SIGATT
      WRITE (ILP,622) II1,II2,XT1,XT2
      IF (SC1.GT.0.OR.SC2.LT.1) WRITE (ILP,621) SC1,SC2
      WRITE (ILP,634) WINDOW
      WRITE (ILP,626) CUTOFF
      IF (NSMOOTH.NE.2) WRITE (ILP,627) NSMOOTH
      WRITE (ILP,610) ATIME,ADATE,TITLE
      WRITE (ILP,623) ((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,617)
      IF (TYPE1.EQ.'LRNZ') WRITE (ILP,609)
      IF (TYPE1.EQ.'GAUS') WRITE (ILP,608)
      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,642)    T1
      WRITE (ILP,643) SIGT1
      WRITE (ILP,618)
      IF (TYPE2.EQ.'LRNZ') WRITE (ILP,609)
      IF (TYPE2.EQ.'GAUS') WRITE (ILP,608)
      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,642)    T2
      WRITE (ILP,643) SIGT2
C
C         CALCULATE RECIPROCAL SPACE METRIC TENSOR.
C
      CALL METRIC (CELL,GINV)
C
C         EVALUATE MONOCHROMATOR VARIABLES.
C
      THM=THM*PI/180
      RHOM=RHOM*PI/180
      COSSQ=COS(RHOM)**2
      SINSQ=SIN(RHOM)**2
      C0KIN=COS(2*THM)**2
      C1KIN=COSSQ+SINSQ*C0KIN
      C2KIN=SINSQ+COSSQ*C0KIN
      C0DYN=ABS(COS(2*THM))
      C1DYN=COSSQ+SINSQ*C0DYN
      C2DYN=SINSQ+COSSQ*C0DYN
      IF (FRACTD.NE.0.AND.SIGFD.EQ.0) SIGFD=0.05*FRACTD
      TANTHM=TAN(THM)
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
      VARFD=SIGFD**2
      VARTAU=SIGTAU**2
      VARATT=SIGATT**2
      VARTHE=SIGTHE**2
      VARTIM=SIGTIM**2
C
C         SET PEAK POSITION AND PEAK WIDTHS TO BE USED INSTEAD OF VALUES
C         CALCULATED FROM THE PEAK POSITION AND PEAK WIDTH PARAMETERS.
C
      IF (XX0.NE.0.OR.XW1.NE.0.OR.XW2.NE.0)
     & WRITE (ILP,610) ATIME,ADATE,TITLE
      IF (XX0.NE.0) WRITE (ILP,647) XX0
      IF (XW1.NE.0.OR.XW2.NE.0) WRITE (ILP,648) XW1,XW2
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,644)
      NCOND=0
 7    READ (IO1,506,END=8) HCOND
 506  FORMAT (A)
      WRITE (ILP,645) HCOND
      I=0
      CALL HKLCOND (HCOND,I)
      NCOND=NCOND+1
      ICOND(NCOND)=I
      GO TO 7
 8    IF (NCOND.EQ.0) THEN
        WRITE (ILP,646)
        CLOSE (UNIT=IO1,DISPOSE='DELETE')
      ELSE
        CLOSE (UNIT=IO1)
      END IF
C
C         OPEN INPUT AND OUTPUT REFLECTION DATA FILES.
C
      OPEN (UNIT=IO1,FILE=FILE1,STATUS='OLD',FORM='UNFORMATTED',
     & READONLY)
      OPEN (UNIT=IO2,FILE=FILE2,STATUS='NEW',FORM='UNFORMATTED')
      OPEN (UNIT=IO3,FILE='OUTPUT.DAT',STATUS='SCRATCH',FORM=
     & 'UNFORMATTED')
      RETURN
C
C         FORMATS FOR PRINTED OUTPUT
C
 600  FORMAT ('0INPUT DATA FROM PROGRAM REFPK JOB:'/
     &'0    TIME      DATE  TITLE'/' ',A,1X,A,2X,A/)
 610  FORMAT ('1'/'0PROGRAM BGLP.  ',A,1X,A,'.  'A/)
 611  FORMAT ('0INPUT FILE:   ',A)
 612  FORMAT ('0OUTPUT FILE:  ',A/)
 629  FORMAT ('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 ('0NICOLET (NEE SYNTEX) P3 DIFFRACTOMETER'/)
 616  FORMAT ('0ENRAF-NONIUS CAD4 DIFFRACTOMETER'/)
 651  FORMAT ('0X-RAY DATA'/)
 652  FORMAT ('0NEUTRON DATA'/)
 628  FORMAT ('0  THETA(M)    RHO(M) FRAC(DYN) SIG(FRAC)'/' ',2F10.2,
     &2E10.3/)
 631  FORMAT ('0FILT. EDGE ALPHA-ONE ALPHA-BAR ALPHA-TWO'/' ',4F10.5/)
 613  FORMAT ('0SIG(THETA) SIG(TIME) (THETA IN DEGREES AND TIME IN MICRO
     &SECONDS)'/' ',F10.4,F10.2/)
 620  FORMAT ('0       TAU  SIG(TAU)       ATT  SIG(ATT) (TAU IN MICROSE
     &CONDS)'/' ',2F10.2,2F10.2/)
 622  FORMAT ('0BEGINNING AND ENDING REFLECTION SERIAL NUMBERS AND X-RAY
     & EXPOSURE TIMES (HOURS):'/'0     IIMIN     IIMAX     XTMIN     XTM
     &AX'/' ',2I10,2F10.2/)
 621  FORMAT ('0SCAN LIMITS (DECIMAL FRACTION OF SCAN WIDTH):'/
     &'0        I1        I2'/' ',2F10.2/)
 634  FORMAT ('0WINDOW WIDTH FOR INITIAL SCAN FOR PEAK LIMITS (DECIMAL F
     &RACTION OF SCAN WIDTH)'/'0    WINDOW'/' ',F10.4)
 626  FORMAT ('0E.S.D. MULTIPLIER FOR STATISTICAL SIGNIFICANCE TESTS'/'0
     &    CUTOFF'/' ',F10.2)
 627  FORMAT ('0NUMBER OF PASSES OF DATA-SMOOTHING FOR PEAK LOCATION'/'0
     &   NSMOOTH'/' ',I10)
 623  FORMAT ('0PEAK POSITION MATRIX'/
     &'0ORIENTATION MATRIX RECALCULATED BY LEAST-SQUARES FIT TO OMEGA'/
     &' VALUES FOR THE INTENSITY-WEIGHTED PEAK CENTROIDS ALONG WITH'/
     &' THE DIFFRACTOMETER ANGLES CHI AND PHI.'/
     &'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)
 617  FORMAT ('0PEAK WIDTH TENSORS'/
     &'0W1 AND W2 ARE THE BASE WIDTHS OF THE HALF-PEAKS BELOW'/
     &' THETA(ALPHA-1) AND ABOVE THETA(ALPHA-2), RESPECTIVELY,'/
     &' IN UNITS OF DEGREES OF CRYSTAL ROTATION (THETA).'/
     &'0LOW ANGLE HALF-PEAK W1')
 618  FORMAT ('0HIGH ANGLE HALF-PEAK W2')
 609  FORMAT ('0LORENTZIAN MODEL'/
     &'0W = (Z(K)*Q(J,K)*Z(J))**(1/2) + T*TAN(THETA)')
 608  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 = ',F3.1)
 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)
 642  FORMAT ('0                    T         = ',F8.4)
 643  FORMAT ('0                    SIG       = ',F8.4)
 647  FORMAT ('0FOR WEAK REFLECTIONS WITH NO SIGNIFICANT PEAK ABOVE BACK
     &GROUND, THE'/' DEFAULT ESTIMATE OF PEAK POSITION WILL BE'/'0X0 =',
     &' ',F4.2,' (DECIMAL FRACTION OF SCAN WIDTH)'/'0INSTEAD OF THE VALU
     &E CALCULATED USING THE DIFFRACTOMETER ORIENTATION MATRIX.'/)
 648  FORMAT ('0FOR EACH REFLECTION, THE PEAK WIDTH WILL BE FIXED AS W =
     & L2 - L1, WHERE'/'0L1 = ',F4.2,' (DECIMAL FRACTION OF SCAN WIDTH)'
     &/' L2 = ',F4.2,', AND THUS'/'0W = ',F4.2,' CENTERED ABOUT THE INTE
     &NSITY-WEIGHTED PEAK CENTROID,'/'0INSTEAD OF THE VALUES CALCULATE',
     &'D USING THE ANISOTROPIC PEAK WIDTH TENSORS.'/)
 644  FORMAT ('0CONDITIONS LIMITING POSSIBLE REFLECTIONS'/)
 645  FORMAT (' ',2X,A)
 646  FORMAT (' NONE')
      END
C-----------------------------------------------------------------------
      SUBROUTINE LZ (TH,VT,FL,VL)
C
C     LORENTZ FACTOR.  EQUATORIAL DIFFRACTION GEOMETRY
C
      U=COS(2*TH)
      V=SIN(2*TH)
      FL=1/V
      VL=VT*(-2*U/V**2)**2
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE METRIC (A,GB)
C
C         CALCULATE RECIPROCAL LATTICE PARAMETERS AND DIRECT AND
C         RECIPROCAL LATTICE METRIC TENSORS FROM DIRECT LATTICE
C         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-----------------------------------------------------------------------
      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 /BLOCK7/ 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
      CHARACTER ATIME*8,ADATE*9,TITLE*80,HEADING*120
      COMMON /BLOCK0/ IO1,IO2,IO3,ILP
      COMMON /BLOCK1/ ATIME,ADATE,TITLE
      COMMON /BLOCK3/ WINDOW,CUTOFF,NSMOOTH,II1,II2,XT1,XT2,SC1,SC2,
     & XX0,XW1,XW2
      PARAMETER (M1=4,M2=12,M3=50,MI=26,MJ=13,MK=16)
      DIMENSION HEADING(M1),JX(M2),JXMIN(M1,M2,M3),JXMAX(M1,M2,M3)
      DIMENSION XI(MI),YI(MI),NI(MI),XJ(MJ),YJ(MJ),NJ(MJ),
     & XK(MK),YK(MK),NK(MK)
C
C         DATA FOR CATALOGING REFLECTIONS
C
      DATA XI /-1000, -300, -100, -50, -20, -10, -5, -2, -1, 0, 1, 2, 5,
     &         10, 20, 50, 100, 300, 1000, 3000, 10000, 30000, 100000,
     &         300000, 1000000, 0/
      DATA XJ /-10, -5, -2, -1, 0, 1, 2, 5, 10, 20, 50, 100, 0/
      DATA XK /0.5, 0.65, 0.8, 0.9, 1.0, 1.05, 1.1, 1.15, 1.175, 1.2,
     &         1.25, 1.3, 1.35, 1.4, 1.45, 1.5/
      DATA NTOT,NREF,NCEN,NSIG /4*0/
C
C         INITIALIZE ARRAYS FOR CATALOGING VARIOUS CLASSES OF DATA.
C
      DO 11 I=1,MI
      YI(I)=0
      NI(I)=0
 11   CONTINUE
      DO 12 J=1,MJ
      YJ(J)=0
      NJ(J)=0
 12   CONTINUE
      DO 13 K=1,MK
      YK(K)=0
      NK(K)=0
 13   CONTINUE
      DO 15 J1=1,M1
      DO 15 J2=1,M2
      DO 15 J3=1,M3
      JXMIN(J1,J2,J3)=+1E5
      JXMAX(J1,J2,J3)=-1E5
 15   CONTINUE
      WRITE (ILP,610) ATIME,ADATE,TITLE
      WRITE (ILP,635)
C
C         LOOP THROUGH REFLECTION DATA FILE.
C
      REWIND IO3
 1    READ (IO3,END=9) II,IH,IK,IL,S,R,SIGR,I0,IE,L1,L2,W,WEAK
      NTOT=NTOT+1
      IF (II.LT.0) GO TO 1
      NREF=NREF+1
      IF (WEAK.EQ.0) NCEN=NCEN+1
      IF (R.GE.CUTOFF*SIGR) NSIG=NSIG+1
C
C         PRINT SELECTED REFLECTION RECORDS.
C
      CALL DPRINT (II,IH,IK,IL,IPRINT)
      IF (IPRINT.EQ.1) WRITE (ILP,636) II,IH,IK,IL,I0,MAX(IE,0),L1,L2,
     & R,SIGR
C
C         CATALOG DATA IN CLASSES OF R = I/LP, R/SIG(R), AND S =
C         SIN(THETA)/LAMBDA.
C
      DO 31 I=1,MI-1
      IF (R.LT.XI(I)) GO TO 32
 31   CONTINUE
      I=MI
 32   NI(I)=NI(I)+1
      YI(I)=YI(I)+R/SIGR
      DO 33 J=1,MJ-1
      IF (R/SIGR.LT.XJ(J)) GO TO 34
 33   CONTINUE
      J=MJ
 34   NJ(J)=NJ(J)+1
      YJ(J)=YJ(J)+R/SIGR
      IF (R/SIGR.LT.CUTOFF) GO TO 1
      DO 35 K=1,MK-1
      IF (S.LT.XK(K)) GO TO 36
 35   CONTINUE
      K=MK
 36   NK(K)=NK(K)+1
      YK(K)=YK(K)+R/SIGR
C
C         CATALOG PROFILES WITH EXTREME FEATURES.
C
      JX(1)=II
      JX(2)=IH
      JX(3)=IK
      JX(4)=IL
      JX(5)=1000*S
      JX(6)=1000*R
      JX(7)=1000*SIGR
      JX(8)=I0
      JX(9)=MAX(IE,0)
      JX(10)=L1
      JX(11)=L2
      JX(12)=1000*W
      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)
      DO 95 J1=1,M1
      DO 95 J2=1,M2
      DO 95 J3=1,M3
      IF (JXMIN(J1,J2,J3).EQ.+1E5) JXMIN(J1,J2,J3)=0
      IF (JXMAX(J1,J2,J3).EQ.-1E5) JXMAX(J1,J2,J3)=0
 95   CONTINUE
C
C         PRINT CATALOG OF CLASSES OF DATA.
C
      WRITE (ILP,610) ATIME,ADATE,TITLE
      WRITE (ILP,641) NTOT,NREF,NCEN,NSIG,CUTOFF
      ISUM=0
      JSUM=0
      KSUM=0
      DO 37 I=1,MI
      IF (NI(I).EQ.0) GO TO 37
      ISUM=ISUM+NI(I)
      YI(I)=YI(I)/NI(I)
 37   CONTINUE
      DO 38 J=1,MJ
      IF (NJ(J).EQ.0) GO TO 38
      JSUM=JSUM+NJ(J)
      YJ(J)=YJ(J)/NJ(J)
 38   CONTINUE
      DO 39 K=1,MK
      IF (NK(K).EQ.0) GO TO 39
      KSUM=KSUM+NK(K)
      YK(K)=YK(K)/NK(K)
 39   CONTINUE
      WRITE (ILP,637) XI(1),NI(1),YI(1),(XI(I-1),XI(I),NI(I),
     &YI(I),I=2,(MI-1)),XI(MI-1),NI(MI),YI(MI)
      WRITE (ILP,640) ISUM
      WRITE (ILP,638) XJ(1),NJ(1),YJ(1),(XJ(J-1),XJ(J),NJ(J),
     &YJ(J),J=2,(MJ-1)),XJ(MJ-1),NJ(MJ),YJ(MJ)
      WRITE (ILP,640) JSUM
      WRITE (ILP,639) XK(1),NK(1),YK(1),(XK(K-1),XK(K),NK(K),
     &YK(K),K=2,(MK-1)),XK(MK-1),NK(MK),YK(MK)
      WRITE (ILP,640) KSUM
C
C         PRINT CATALOGS OF EXTREME CASES.
C
      WRITE (ILP,610) ATIME,ADATE,TITLE
      WRITE (ILP,644)
      DO 49 J1=1,M1
      WRITE (ILP,610) ATIME,ADATE,TITLE
      WRITE (ILP,645) HEADING(J1)
      DO 59 J3=1,MIN(M3,NSIG)
      F5=JXMIN(J1,5,J3)*0.001
      F6=JXMIN(J1,6,J3)*0.001
      F7=JXMIN(J1,7,J3)*0.001
      F12=JXMIN(J1,12,J3)*0.001
 59   WRITE (ILP,646) (JXMIN(J1,J2,J3),J2=1,4),F5,F6,F7,
     &(JXMIN(J1,J2,J3),J2=8,11),F12
      WRITE (ILP,610) ATIME,ADATE,TITLE
      WRITE (ILP,645) HEADING(J1)
      DO 69 J3=1,MIN(M3,NSIG)
      F5=JXMAX(J1,5,J3)*0.001
      F6=JXMAX(J1,6,J3)*0.001
      F7=JXMAX(J1,7,J3)*0.001
      F12=JXMAX(J1,12,J3)*0.001
 69   WRITE (ILP,646) (JXMAX(J1,J2,J3),J2=1,4),F5,F6,F7,
     &(JXMAX(J1,J2,J3),J2=8,11),F12
 49   CONTINUE
      RETURN
C
C         FORMAT STATEMENTS FOR PRINTED OUTPUT
C
 610  FORMAT ('1'/'0PROGRAM BGLP.  ',A,1X,A,'.  'A/)
 635  FORMAT ('0     I     H   K   L    X0    XE    L1  L2    Y=I/LP',
     &'        SIG(Y)'/)
 636  FORMAT (' ',I6,2X,3I4,2X,2(I4,2X),2I4,2F10.2)
 637  FORMAT ('0',12X'Y-INTERVAL         NUMBER  AVERAGE Y/SIG(Y), '
     &        'WHERE Y = I/LP'/
     &        '0',9X,  '   Y < ',F9.1,5X,I7,F9.1/
     &     24(' ',F9.1,' < Y < ',F9.1,5X,I7,F9.1/),
     &        ' ',F9.1,' < Y   ',9X,  5X,I7,F9.1)
 640  FORMAT (' ',32X,'TOTAL'/ ' ',30X,I7)
 638  FORMAT ('0',12X,'R-INTERVAL         NUMBER  AVERAGE R, WHERE '
     &        'R = Y/SIG(Y)'/
     &        '0',9X,  '   R < ',F9.1,5X,I7,F9.1/
     &     11(' ',F9.1,' < R < ',F9.1,5X,I7,F9.1/),
     &        ' ',F9.1,' < R   ',9X,  5X,I7,F9.1)
 639  FORMAT ('1',12X,'S-INTERVAL         NUMBER  AVERAGE R, WHERE '
     &        'S = SIN(THETA)/LAMBDA'/
     &        '0',9X,  '   S < ',F9.3,5X,I7,F9.1/
     &     14(' ',F9.3,' < S < ',F9.3,5X,I7,F9.1/),
     &        ' ',F9.3,' < S   ',9X,  5X,I7,F9.1)
 641  FORMAT ('0NTOT = ',I6,'  TOTAL REFLECTIONS (EXCLUDING SYMMETRY FOR
     &BIDDEN REFLECTIONS)'/
     &        ' NREF = ',I6,'  NOT STANDARD REFERENCE REFLECTIONS'/
     &        ' NCEN = ',I6,'  WITH A SIGNIFICANT INTENSITY-WEIGHTED PEA
     &K CENTROID'/
     &        ' NSIG = ',I6,'  WITH Y > CUTOFF*SIGMA(Y), WHERE CUTOFF =
     & ',F5.2)
 644  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 /
     &'S=SIN(THETA)/LAMBDA',
     &'Y=I/LP, LORENTZ AND POLARIZATION CORRECTED NET INTENSITY',
     &'X0, INTENSITY WEIGHTED PEAK CENTROID',
     &'W=(L2-L1+1)*DTH, BASE WIDTH OF REFLECTION (DEGREES THETA)'/
 645  FORMAT ('0',A/'0     I     H   K   L      S    Y=I/LP    SIG(Y)',
     &'    X0    XE    L1  L2         W'/)
 646  FORMAT (' ',I6,2X,3I4,2XF5.3,2F10.2,2XI4,2XI4,2X,2I4,F10.2)
      END
C-----------------------------------------------------------------------
      SUBROUTINE PROCESS
      PARAMETER (NMAX=100)
      CHARACTER ATIME*8,ADATE*9,TITLE*80
      COMMON /BLOCK0/ IO1,IO2,IO3,ILP
      COMMON /BLOCKC/ NCOND,ICOND(38)
      COMMON /BLOCK1/ ATIME,ADATE,TITLE
      COMMON /BLOCK2/ GINV(3,3),WLE,WL1,WL0,WL2
      COMMON /BLOCK3/ WINDOW,CUTOFF,NSMOOTH,II1,II2,XT1,XT2,SC1,SC2,
     & XX0,XW1,XW2
      COMMON /BLOCK5/ IDIFF,U(3,3),MODEL1,MODEL2,Q1(3,3),Q2(3,3),T1,T2,
     & F,SIGU(3,3),SIGQ1(3,3),SIGQ2(3,3),SIGT1,SIGT2,TANTHM
      COMMON /BLOCK6/ VARTHE,VARTIM,TAU,VARTAU,ATT,VARATT
      COMMON /BLOCK7/ NSTEP,XX(NMAX),YY(NMAX),VX(NMAX),VY(NMAX)
      COMMON /BLOCK8/ TH0,X0,SIGX0,XE,X1,X2,W1,W2,I1,I2,L1,L2,R,SIGR,EPS
      DIMENSION H(3),ANGLES(4),YYTEMP(NMAX),VYTEMP(NMAX)
      PI=ACOS(-1.0)
      DEG=180/PI
      RAD=PI/180
C
C     LOOP THROUGH PROCESSING 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     SKIP SYMMETRY FORBIDDEN REFLECTIONS.
C
      CALL HKLTEST (IH,IK,IL,NCOND,ICOND,IABSENT)
      IF (IABSENT.NE.0) GO TO 1
C
C     STORE INDICES IN 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     TRANSFORM TO SETTING ANGLES AS DEFINED BY WALTER HAMILTON,
C     INTERNATIONAL TABLES FOR X-RAY CRYSTALLOGRAPHY, VOL. IV, 1974, PP.
C     273-284.
C
      CALL GEOM (IDIFF,ANGLES,TWOTH,OMEGA,CHI,PHI)
C
C     SET 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     SMOOTH PROFILE FOR PEAK CENTROID LOCATION.
C
      DO 11 I=1,NSTEP
      YYTEMP(I)=YY(I)
      VYTEMP(I)=VY(I)
   11 CONTINUE
      IF (NSMOOTH.NE.0) THEN
        CALL SMOOTH (NSMOOTH,NSTEP,YYTEMP)
        CALL SMOOTH (NSMOOTH,NSTEP,VYTEMP)
      END IF
C
C     FIND INTENSITY-WEIGHTED CENTROID OF PEAK ABOVE BACKGROUND.
C
      CALL CENTROID (WINDOW,CUTOFF,XX,YYTEMP,VX,VYTEMP,I1,I2,X0,SIGX0)
C
C     TREAT SCANS THAT HAVE NO SIGNIFICANT PEAK ABOVE BACKGROUND.
C
      IF (X0.EQ.90) THEN
        WEAK=1
        CALL XCALC (XX0,XX,I1,I2,H,U,SIGU,GINV,TWOTH,OMEGA,CHI,PHI,
     &   X0,SIGX0)
      ELSE
        WEAK=0
      END IF
C
C     CALCULATE POSITIONS OF SPECTRAL FEATURES.
C
      XE=X0+(THE-TH0)
      X1=X0+(TH1-TH0)
      X2=X0+(TH2-TH0)
      IF (WEAK.EQ.0) CALL XCHECK (XE,X1,X0,X2,XX,YYTEMP,I1,I2)
C
C     CALCULATE PEAK WIDTHS.
C
      TANTH=TAN(TH0*RAD)
      CALL WCALC (XW1,XX,I1,I2,CHI,PHI,MODEL1,Q1,T1,F,TANTH,TANTHM,
     & SIGQ1,SIGT1,W1,SIGW1)
      CALL WCALC (XW2,XX,I1,I2,CHI,PHI,MODEL2,Q2,T2,F,TANTH,TANTHM,
     & SIGQ2,SIGT2,W2,SIGW2)
C
C     LIMIT THE ESTIMATED UNCERTAINTIES IN PEAK POSITION AND PEAK WIDTHS
C     TO ONE STEP.
C
      DX=(XX(I2)-XX(I1))/(I2-I1+1)
      SIGX0=MIN(SIGX0,DX)
      SIGW1=MIN(SIGW1,DX)
      SIGW2=MIN(SIGW2,DX)
C
C     CALCULATE PEAK LIMITS.
C
      L1=NINDEX(X1-SIGX0-W1-SIGW1)
      L2=NINDEX(X2+SIGX0+W2+SIGW2)
C
C     PERFORM BACKGROUND SUBTRACTION AND LORENTZ AND POLARIZATION
C     CORRECTIONS.
C
      IF (WLE.EQ.0.OR.WEAK.EQ.1) THEN
        CALL BGLP1
      ELSE
C
C       EPSILON IS AN ESTIMATE OF THE HALF-WIDTH OF THE BETA-FILTER
C       ABSORPTION EDGE BROADENED DUE TO FINITE INSTRUMENTAL RESOLUTION.
C
        EPS=0.5*(W1+W2)
        J1=NINDEX(XE-EPS)
        J2=NINDEX(XE+EPS)
        IF (J2.LT.L1-1) THEN
          I1=MAX(I1,J2)
          CALL BGLP1
        ELSE
          CALL BGLP2
        END IF
      END IF
C
C     WRITE REFLECTION RECORD TO OUTPUT FILES.
C
      CALL WRITE1 (IO2,II,IH,IK,IL,ANGLES,R,SIGR,XTIME)
      I0=NINDEX(X0)
      IE=NINDEX(XE)
      IF (WLE.EQ.0) IE=0
      W=XX(L2)-XX(L1)
      WRITE (IO3) II,IH,IK,IL,S,R,SIGR,I0,IE,L1,L2,W,WEAK
C
C     LOOP BACK TO PROCESS NEXT REFLECTION.
C
      GO TO 1
C
C     END LOOP THROUGH REFLECTION DATA.
C
    9 CLOSE (UNIT=IO2)
      ENDFILE IO3
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE PZ (TH,VT,FP,VP)
C
C     POLARIZATION FACTOR.  WITH OR WITHOUT INCIDENT BEAM CRYSTAL
C     MONOCHROMATOR
C
      COMMON /BLOCK4/ C0K,C1K,C2K,C0D,C1D,C2D,FD,VF,INEUTRON
C
C     C0K = COS(2*THM)**2   KINEMATIC (MOSAIC) MONOCHROMATOR CRYSTAL
C     C0D = ABS(COS(2*THM)) DYNAMIC (PERFECT) MONOCHROMATOR CRYSTAL
C
C     C1 = COS(RHO)**2 + SIN(RHO)**2*C0
C     C2 = SIN(RHO)**2 + COS(RHO)**2*C0
C
C     RHO =  0      PARALLEL MONOCHROMATOR-DIFFRACTOMETER GEOMETRY
C     RHO = 90 DEG  PERPENDICULAR GEOMETRY
C
C     FD = 0  ALL KINEMATIC, NO DYNAMIC MONOCHROMATOR DIFFRACTION
C     FD = 1  ALL DYNAMIC, NO KINEMATIC MONOCHROMATOR DIFFRACTION
C
C     VF = VAR(FD), THE VARIANCE OF THE DYNAMIC DIFFRACTION FRACTION
C
C     FOR AN UNPOLARIZED INCIDENT BEAM (I.E., NO MONOCHROMATOR),
C     THM = 0, C0 = C1 = C2 = 1, AND FD = VF = 0.
C
C     FOR NEUTRON DIFFRACTION DATA, INEUTRON = 1.
C
      IF (INEUTRON.EQ.1) THEN
        FP=1
        VP=0
      ELSE
        U=COS(2*TH)
        V=SIN(2*TH)
        PK=(C1K+C2K*U**2)/(1+C0K)
        VK=VT*(-4*C2K*U*V/(1+C0K))**2
        IF (FD.EQ.0) THEN
          FP=PK
          VP=VK
        ELSE
          PD=(C1D+C2D*U**2)/(1+C0D)
          VD=VT*(-4*C2D*U*V/(1+C0D))**2
          FP=FD*PD+(1-FD)*PK
          VP=FD**2*VD+(1-FD)**2*VK+(PD**2+PK**2)*VF
        END IF
      END IF
      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 WCALC (FW,X,I1,I2,CHI,PHI,MODEL,Q,T,F,TANTH,TANTHM,
     & SIGQ,SIGT,W,SIGW)
C
C     CALCULATED PEAK WIDTHS
C
      DIMENSION X(1),Q(3,3),SIGQ(3,3),Z(3)
      DATA RAD /0.017453293/
      IF (FW.NE.0) THEN
        W=FW*(X(I2)-X(I1))
        SIGW=0
        RETURN
      END IF
C
C     GET THE COMPONENTS ALONG CARTESIAN AXES FIXED IN THE CRYSTAL OF A
C     UNIT VECTOR, Z, NORMAL TO THE EQUATORIAL PLANE OF THE INCIDENT AND
C     DIFFRACTED BEAM VECTORS.
C
      COSCH=COS(CHI*RAD)
      SINCH=SIN(CHI*RAD)
      COSPH=COS(PHI*RAD)
      SINPH=SIN(PHI*RAD)
      Z(1)=SINPH*SINCH
      Z(2)=COSPH*SINCH
      Z(3)=COSCH
C
C     CALCULATE PEAK WIDTHS.
C
      U=0
      V=0
      DO 10 I=1,3
      DO 10 J=1,3
      U=U+Z(I)*Z(J)*Q(I,J)
      V=V+(Z(I)*Z(J)*SIGQ(I,J))**2
   10 CONTINUE
      IF (MODEL.EQ.1) THEN
C
C       LORENTZIAN CONVOLUTION MODEL
C
        W=SQRT(U)+T*TANTH
        SIGW=SQRT(V/(2*U)**2+(SIGT*TANTH)**2)
      ELSE IF (MODEL.EQ.2) THEN
C
C       GAUSSIAN CONVOLUTION MODEL
C
        W=SQRT(U+(T*TANTH)**2)
        SIGW=SQRT(V+(2*T*TANTH**2*SIGT)**2)/(2*W)
      ELSE IF (MODEL.EQ.3) THEN
C
C       GAUSSIAN CONVOLUTION MODEL FOR PARALLEL MONOCHROMATOR GEOMETRY
C
        W=SQRT(ABS(U-2*F*T*TANTH/TANTHM+T*(TANTH/TANTHM)**2))
        SIGW=SQRT(V/(2*U)**2+(2*F*SIGT*TANTH/TANTHM)**2
     &         +(SIGT*(TANTH/TANTHM)**2)**2)
      END IF
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE WRITE1 (IFILE,II,IH,IK,IL,ANGLES,R,SIGR,XTIME)
      DIMENSION ANGLES(4)
      WRITE (IFILE) II,IH,IK,IL,ANGLES,R,SIGR,XTIME
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE XCALC (FX,X,I1,I2,H,U,SIGU,GINV,TWOTH,OMEGA,CHI,PHI,
     & X0,SIGX0)
C
C     EXPECTED PEAK POSITION
C
      DIMENSION X(1),H(3),U(3,3),SIGU(3,3),GINV(3,3),Y(3),V(3)
      DATA RAD,DEG /0.017453293, 57.2957795/
      IF (FX.NE.0) THEN
        X0=X(I1)+FX*(X(I2)-X(I1))
        SIGX0=(X(I2)-X(I1))/(I2-I1+1)
        RETURN
      END IF
C
C     GET THE COMPONENTS ALONG CARTESIAN AXES FIXED IN THE CRYSTAL OF A
C     UNIT VECTOR, Y, PARALLEL TO THE RECIPROCAL LATTICE VECTOR, DSTAR,
C     IN THE DIFFRACTING CONDITION.
C
      DO 1 I=1,3
      Y(I)=0
      V(I)=0
      DO 1 J=1,3
      Y(I)=Y(I)+H(J)*U(J,I)
      V(I)=V(I)+(H(J)*SIGU(J,I))**2
    1 CONTINUE
      DSTAR=2*SINTHL(H,GINV)
      DO 2 I=1,3
      Y(I)=Y(I)/DSTAR
      V(I)=V(I)/DSTAR**2
    2 CONTINUE
      IF (TWOTH.LT.0) THEN
        DO 3 I=1,3
        Y(I)=-Y(I)
    3   CONTINUE
      END IF
C
C     CALCULATE EXPECTED PEAK DISPLACEMENT, DELTA(OMEGA).
C
      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
C
C     THE POSITIVE SENSE OF THE OMEGA ROTATION HAS BEEN DEFINED SUCH
C     THAT DELTA(X0) = -DELTA(OMEGA).
C
      X0=0.5*(X(I1)+X(I2))-DELTA
      VARSIN=V(1)*COSPH**2+V(2)*SINPH**2
      VARCOS=(V(1)*SINPH**2+V(2)*COSPH**2)*COSCH**2+V(3)*SINCH**2
      OMEGA=ARCTAN(SINOM,COSOM)
      VAROM=(COS(OMEGA)**2/COSOM)**2*(TAN(OMEGA)**2*VARCOS+VARSIN)
      SIGX0=SQRT(VAROM)*DEG
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE XCHECK (XE,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
        XE=XE+D
        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 /BLOCK6/ 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-----------------------------------------------------------------------
