      PROGRAM TSCALE
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 (MSTD=20,MCUT=10,MDAT=500,MOMIT=20,MREJ=50)
      CHARACTER TITLE*80,ATIME*8,ADATE*9,FILE1*40,FILE2*40
      COMMON /BLOCK1/ TITLE,ATIME,ADATE
      COMMON /BLOCK2/ IO1,IO2,ILP,ITABL,IPLOT,IQSUM,IGRAPH,
     &                IAPPLY,IANISO,IYDIFF,ISDIFF,A(6),GINV(3,3)
      COMMON /BLOCK3/ NSTD,IHKL(MSTD,3),NCUT(MSTD),XCUT(MSTD,MCUT),XZERO
      COMMON /BLOCK4/ NIOMIT(MSTD),IOMIT(MSTD,MOMIT,2),
     &                NXOMIT(MSTD),XOMIT(MSTD,MOMIT,2)
      COMMON /BLOCK5/ NIREJ,IREJ(MREJ,2),
     &                NXREJ,XREJ(MREJ,2)
      COMMON /BLOCK6/ NDAT(MSTD,MCUT),XX(MSTD,MCUT,MDAT),
     &                YY(MSTD,MCUT,MDAT),SS(MSTD,MCUT,MDAT)
      COMMON /BLOCK7/ PAR(MSTD,MCUT,4),COV(MSTD,MCUT,4,4),
     &                XLIM(MSTD,MCUT,2),YZERO(MSTD),P
      CALL INPUT
      CALL CULL
      CALL ANALYZE
      CALL APPLY
      STOP 'PROGRAM TSCALE FINIS!'
      END
C-----------------------------------------------------------------------
      SUBROUTINE ANALYZE
      PARAMETER (MSTD=20,MCUT=10,MDAT=500,MOMIT=20,MREJ=50)
      PARAMETER (NX=101,NY=51)
      CHARACTER TITLE*80,ATIME*8,ADATE*9
      COMMON /BLOCK1/ TITLE,ATIME,ADATE
      COMMON /BLOCK2/ IO1,IO2,ILP,ITABL,IPLOT,IQSUM,IGRAPH,
     &                IAPPLY,IANISO,IYDIFF,ISDIFF,A(6),GINV(3,3)
      COMMON /BLOCK3/ NSTD,IHKL(MSTD,3),NCUT(MSTD),XCUT(MSTD,MCUT),XZERO
      COMMON /BLOCK6/ NDAT(MSTD,MCUT),XX(MSTD,MCUT,MDAT),
     &                YY(MSTD,MCUT,MDAT),SS(MSTD,MCUT,MDAT)
      COMMON /BLOCK7/ PAR(MSTD,MCUT,4),COV(MSTD,MCUT,4,4),
     &                XLIM(MSTD,MCUT,2),YZERO(MSTD),P
      DIMENSION X(MDAT),Y(MDAT),S(MDAT),C(4),CC(4,4)
      DIMENSION SZERO(MSTD),PP(MSTD,MCUT)
C
C         ANALYZE STANDARD REFLECTION INTENSITIES.
C
      DO 10 I=1,NSTD
      DO 10 J=1,NCUT(I)+1
      N=NDAT(I,J)
      IF (N.EQ.0) GO TO 10
      DO 50 K=1,4
      PAR(I,J,K)=0
      DO 50 L=1,4
      COV(I,J,K,L)=0
 50   CONTINUE
      XMIN=+1E6
      XMAX=-1E6
      YMIN=+1E6
      YMAX=-1E6
      SUMN=0
      SUMW=0
      SUMWY=0
      SUMWY2=0
C
C         PICK OUT DATA AND EVALUATE STATISTICS FOR JTH SEGMENT OF ITH
C         STANDARD.
C
      DO 11 K=1,N
      X(K)=XX(I,J,K)
      Y(K)=YY(I,J,K)
      S(K)=SS(I,J,K)
      XI=X(K)
      YI=Y(K)
      SIGYI=S(K)
      WI=1/SIGYI**2
      SUMN=SUMN+1
      SUMW=SUMW+WI
      SUMWY=SUMWY+WI*YI
      SUMWY2=SUMWY2+WI*YI**2
      XMIN=MIN(XMIN,XI)
      XMAX=MAX(XMAX,XI)
 11   CONTINUE
      XLIM(I,J,1)=XMIN
      XLIM(I,J,2)=XMAX
      YMEAN=SUMWY/SUMW
      RMSD=SQRT(ABS(SUMWY2/SUMW-YMEAN**2))
      ESD=SQRT(1/SUMW)
      IF (SUMN.LT.2) GO TO 12
      RMSD=RMSD*SQRT(SUMN/(SUMN-1))
      ESD=ESD*SQRT(SUMN)
 12   IF (ITABL.EQ.0) GO TO 13
      WRITE (ILP,102) ATIME,ADATE,TITLE
      WRITE (ILP,101) (IHKL(I,K),K=1,3),YMEAN,ESD,RMSD,
     & (X(K),Y(K),K=1,N)
 13   CONTINUE
C
C         FIT A CURVE TO THE DATA (X,Y) FOR THE SEGMENT.
C
      IF (N.GT.2) GO TO 14
C
C         FOR SEGMENTS WITH ONLY ONE OR TWO DATA POINTS, ASSUME Y(X) =
C         CONSTANT.
C
      C(1)=YMEAN
      CC(1,1)=MAX(RMSD,ESD)
      DO 130 K=2,4
      C(K)=0
      DO 130 L=2,4
      CC(K,L)=0
 130  CONTINUE
      GO TO 140
 14   CONTINUE
C
C         FOR SEGMENTS WITH THREE OR MORE DATA POINTS, FIT A POWER
C         SERIES POLYNOMIAL, Y(X) = A0 + A1*X + A2*X**2 + A3*X**3, TO
C         THE SEGMENT.
C
      CALL REGRESS (N,X,Y,S,C,CC)
 140  DO 141 K=1,4
      PAR(I,J,K)=C(K)
      DO 141 L=1,4
      COV(I,J,K,L)=CC(K,L)
 141  CONTINUE
C
C         CALCULATE PROPORTIONALITY FACTOR, P, FOR INSTRUMENTAL
C         INSTABILITIES FOR THE SEGMENT.
C
      CALL PCALC (N,X,Y,S,C,PSQ)
      PP(I,J)=PSQ
C
C         PLOT DATA AND FITTED CURVE FOR THE SEGMENTS WITH TWO OR MORE
C         DATA POINTS.
C
      IF (N.LT.2) GO TO 10
      IF (IPLOT.EQ.0) GO TO 17
      DO 15 K=1,N
      YI=Y(K)
      SIGYI=S(K)
      YMIN=MIN(YMIN,YI-2*SIGYI)
      YMAX=MAX(YMAX,YI+2*SIGYI)
 15   CONTINUE
      YMID=(YMIN+YMAX)/2
      DELY=(NY-1)*ESD
      Y1=YMID-DELY/2
      Y2=YMID+DELY/2
      YMIN=MIN(YMIN,Y1)
      YMAX=MAX(YMAX,Y2)
      DELX=(XMAX-XMIN)/(NX-1)
      DELY=(YMAX-YMIN)/(NY-1)
      WRITE (ILP,110) ATIME,ADATE,TITLE
      CALL PLOT1 (N,X,Y,XMIN,XMAX,YMIN,YMAX,C)
      WRITE (ILP,100) (IHKL(I,K),K=1,3),(PAR(I,J,K),K=1,4),
     & (SQRT(COV(I,J,K,K)),K=1,4),SQRT(MAX(PP(I,J),0.0)),DELX,DELY,ESD
 17   CONTINUE
C
C         EVALUATE AND PLOT Q-SUMS FOR THE SEGMENT.
C
      IF (IQSUM.EQ.0) GO TO 10
      QSUM=0
      YMIN=+1E6
      YMAX=-1E6
      DO 16 K=1,N
      QSUM=QSUM+(Y(K)-YMEAN)
      Y(K)=QSUM
      YMIN=MIN(YMIN,QSUM)
      YMAX=MAX(YMAX,QSUM)
 16   CONTINUE
      YMID=(YMIN+YMAX)/2
      DELY=(NY-1)*ESD
      Y1=YMID-DELY/2
      Y2=YMID+DELY/2
      YMIN=MIN(YMIN,Y1)
      YMAX=MAX(YMAX,Y2)
      DELY=(YMAX-YMIN)/(NY-1)
      WRITE (ILP,110) ATIME,ADATE,TITLE
      CALL PLOT2 (N,X,Y,XMIN,XMAX,YMIN,YMAX)
      WRITE (ILP,103) (IHKL(I,K),K=1,3),YMEAN,ESD,DELX,DELY
 10   CONTINUE
C
C         CALCULATE Y0 = Y(X0) AND SIGMA(Y0) FOR EACH STANDARD.
C
      WRITE (ILP,102) ATIME,ADATE,TITLE
      WRITE (ILP,104) XZERO
      DO 30 I=1,NSTD
      IF (NCUT(I).EQ.0.AND.NDAT(I,1).EQ.0) GO TO 30
C
C         FIND THE SEGMENT FOR THE ITH STANDARD THAT INCLUDES X = X0 IN
C         ITS DOMAIN.
C
      N=0
      DO 31 J=1,NCUT(I)+1
      IF (NDAT(I,J).EQ.0) GO TO 31
      N=N+1
      X1=XLIM(I,J,1)
      X2=XLIM(I,J,2)
      IF (X1.LE.XZERO.AND.XZERO.LE.X2) GO TO 32
 31   CONTINUE
      IF (N.EQ.0) GO TO 30
C
C         IF NO SEGMENT OF THE ITH STANDARD INCLUDES X0, THEN RESORT TO
C         THE SEGMENT NEAREST X0, AND EXTRAPOLATE TO X0.
C
      DMIN=1E6
      DO 33 J=1,NCUT(I)+1
      X1=XLIM(I,J,1)
      X2=XLIM(I,J,2)
      D1=ABS(X1-XZERO)
      D2=ABS(X2-XZERO)
      D=MIN(D1,D2)
      IF (D.GE.DMIN) GO TO 33
      DMIN=D
      JMIN=J
 33   CONTINUE
      J=JMIN
 32   CONTINUE
C
C         CALCULATE Y0 AND SIGMA(Y0).
C
      YZERO(I)=0
      DO 36 K=1,4
      YZERO(I)=YZERO(I)+PAR(I,J,K)*XZERO**(K-1)
 36   CONTINUE
      SZERO(I)=0
      DO 37 K=1,4
      DO 37 L=1,4
      SZERO(I)=SZERO(I)+XZERO**(K+L-2)*COV(I,J,K,L)
 37   CONTINUE
      SZERO(I)=SQRT(ABS(SZERO(I)))
      WRITE (ILP,105) I,(IHKL(I,K),K=1,3),YZERO(I),SZERO(I)
 30   CONTINUE
C
C         CONVERT TO NORMALIZED INVERSE SCALING FACTORS Y/Y0 FOR EACH
C         SEGMENT.
C
      WRITE (ILP,107)
      DO 20 I=1,NSTD
      DO 20 J=1,NCUT(I)+1
      IF (NDAT(I,J).EQ.0) GO TO 20
      DO 21 K=1,4
      PAR(I,J,K)=PAR(I,J,K)/YZERO(I)
      DO 21 L=1,K
      T=(COV(I,J,K,L)+PAR(I,J,K)*PAR(I,J,L)*SZERO(I)**2)/YZERO(I)**2
      COV(I,J,K,L)=T
      COV(I,J,L,K)=T
 21   CONTINUE
      DO 22 K=1,NDAT(I,J)
      YI=YY(I,J,K)/YZERO(I)
      SI=YI**2*((SS(I,J,K)/YY(I,J,K))**2+(SZERO(I)/YZERO(I))**2)
      SI=SQRT(SI)
      YY(I,J,K)=YI
      SS(I,J,K)=SI
 22   CONTINUE
      WRITE (ILP,108) I,J,(IHKL(I,K),K=1,3),NDAT(I,J),(XLIM(I,J,K),K=1,
     & 2),(PAR(I,J,K),K=1,4),(SQRT(COV(I,J,K,K)),K=1,4)
 20   CONTINUE
C
C         EVALUATE N-WEIGHTED, SQUARED-INTENSITY-WEIGHTED ROOT-MEAN-
C         SQUARE P-VALUE
C
      WRITE (ILP,102) ATIME,ADATE,TITLE
      WRITE (ILP,109)
      PSQMIN=+1
      PSQMAX=-1
      SUMW=0
      SUMWPSQ=0
      DO 25 I=1,NSTD
      DO 25 J=1,NCUT(I)+1
      N=NDAT(I,J)
      IF (N.EQ.0) GO TO 25
      Y0=YZERO(I)
      PSQ=PP(I,J)
      WRITE (ILP,120) I,J,(IHKL(I,K),K=1,3),N,Y0,PSQ,SQRT(MAX(PSQ,0.0))
      PSQMIN=MIN(PSQMIN,PSQ)
      PSQMAX=MAX(PSQMAX,PSQ)
      W=N*Y0**2
      SUMW=SUMW+W
      SUMWPSQ=SUMWPSQ+W*PSQ
 25   CONTINUE
      PSQ=SUMWPSQ/SUMW
      P=SQRT(MAX(PSQ,0.0))
      WRITE (ILP,106) P,PSQ,PSQMIN,PSQMAX
C
C         PLOT NORMALIZED DATA POINTS AND INVERSE SCALING FUNCTIONS FOR
C         EACH STANDARD AND COMPOSITE INVERSE SCALING FUNCTION AND ITS
C         FRACTIONAL ERROR.
C
      IF (IGRAPH.NE.0) CALL GRAPH
      RETURN
 110  FORMAT('1'/'0',10X,'PROGRAM TSCALE.  ',A,1X,A,'.  ',A)
 100  FORMAT('0',10X,'REFLECTION INTENSITIES.  Y = A0 + A1*X + A2*X**2',
     &' + A3*X**3, WHERE Y = I/LP AND X = T(HR).'/
     &'0',10X,'  H  K  K        A0        A1        A2        A3  SIG(A0
     &)  SIG(A1)  SIG(A2)  SIG(A3)           P'/
     &' ',10X,3I3,4E10.3,4E9.2,2x,E10.3/
     &'0',10X,'DELTA(X) = ',F10.3,' HR AND DELTA(Y) = ',E10.3,' PER GRID
     & POINT.  SIGMA(Y) = ',E10.3)
 101  FORMAT('0','STANDARD REFERENCE REFLECTION INTENSITY DATA Y = ',
     &           'I/LP AS A FUNCTION OF X = T(HR):'/
     &       '0','YMEAN = SUM(W*Y)/SUM(W)'/
     &       ' ','ESD   = SQRT(N/SUM(W))'/
     &       ' ','RMSD  = SQRT((N/(N-1))*(SUM(W*Y**2)/SUM(W) - ',
     &           'YMEAN**2))'/
     &       '0','  H  K  L     YMEAN       ESD      RMSD'/
     &       ' ',3I3,3F10.2/'0',5(9X,'X',9X,'Y',5X)/
     &      ('0',5(F10.3,F10.2,5X)))
 102  FORMAT('1'/'0PROGRAM TSCALE.  ',A,1X,A,'.  ',A)
 103  FORMAT('0',10X,'Q-SUMS.  Q(X(J)) = SUM(I=1,J)(Y(X(I)) - YMEAN) ',
     &               'AS A FUNCTION OF X = T(HR).'/
     &       '0',10X,'  H  K  L     YMEAN       ESD'/
     &       ' ',10X,3I3,2F10.2/
     &       '0',10X,'DELTA(X) = ',F10.3,' HR AND DELTA(Q) = ',E10.3,
     &               ' PER GRID POINT.')
 109  FORMAT('0','PROPORTIONALITY FACTORS FOR INSTRUMENTAL INSTABILITIES
     & FOR EACH SEGMENT:'/
     &       '0','PSQ = (SUM(DELTA(Y)**2) - SUM(SIGMA(Y)**2))/SUM(Y**2)'
     &/
     &       '0','PSQ = (SUM((YI - Y(XI))**2) - SUM(SIGMA(YI)**2)/SUM(YI
     &**2))'/
     &       '0','  I  J    H  K  L      N        Y0         PSQ',
     &'         P = SQRT(MAX(PSQ,0.0))'/)
 120  FORMAT(' ',2I3,2X,3I3,2X,I5,F10.2,E12.3,F10.5)
 106  FORMAT('0','N-WEIGHTED, SQUARED-INTENSITY-WEIGHTED ROOT-MEAN-SQUAR
     &E P-VALUE:'/
     &       '0','P      = ',F8.5/
     &       '0','PSQ    = ',E10.3/
     &       '0','PSQMIN = ',E10.3/
     &       ' ','PSQMAX = ',E10.3)
 104  FORMAT('0','REFERENCE INTENSITY VALUES Y=I/LP CALCULATED AT THE ',
     &           'REFERENCE TIME OF THE EXPERIMENT X0 = ',F10.3,' HR'/
     &       '0','  I       H  K  L  Y0=Y(X0) SIGMA(Y0)'/)
 105  FORMAT(' ',I3,5X,3I3,2F10.2)
 107  FORMAT('0','NORMALIZED INVERSE SCALING FACTORS.  Y/Y0 = A0 + A1*',
     &'X + A2*X**2 + A3*X**3, WHERE X = T(HR).'/
     &'0','  I  J    H  K  L       N    XMIN    XMAX         A0       ',
     &'  A1         A2         A3  SIG(A0)  SIG(A1)  SIG(A2)  SIG(A3)'/)
 108  FORMAT(' ',2I3,2X,3I3,2X,I6,2F8.2,4E11.3,4E9.2)
      END
C-----------------------------------------------------------------------
      SUBROUTINE APPLY
C
C         SCALES REFLECTION DATA ACCORDING TO THE SCALING FUNCTION
C         FROM SUBROUTINE ANALYZE.
C
      PARAMETER (MSTD=20,MCUT=10,MDAT=500,MOMIT=20,MREJ=50)
      CHARACTER TITLE*80,ATIME*8,ADATE*9
      COMMON /BLOCK1/ TITLE,ATIME,ADATE
      COMMON /BLOCK2/ IO1,IO2,ILP,ITABL,IPLOT,IQSUM,IGRAPH,
     &                IAPPLY,IANISO,IYDIFF,ISDIFF,A(6),GINV(3,3)
      COMMON /BLOCK3/ NSTD,IHKL(MSTD,3),NCUT(MSTD),XCUT(MSTD,MCUT),XZERO
      COMMON /BLOCK5/ NIREJ,IREJ(MREJ,2),
     &                NXREJ,XREJ(MREJ,2)
      COMMON /BLOCK6/ NDAT(MSTD,MCUT),XX(MSTD,MCUT,MDAT),
     &                YY(MSTD,MCUT,MDAT),SS(MSTD,MCUT,MDAT)
      COMMON /BLOCK7/ PAR(MSTD,MCUT,4),COV(MSTD,MCUT,4,4),
     &                XLIM(MSTD,MCUT,2),YZERO(MSTD),P
      DIMENSION U(3),V(3)
      IF (IAPPLY.EQ.0) RETURN
C
C         TEST WHETHER SCALING IS TO BE WEIGHTED ANISOTROPICALLY OR BY
C         INTENSITY OR SIN(THETA)/LAMBDA DIFFERENCES.
C
      IF (IANISO.EQ.0.AND.IYDIFF.EQ.0.AND.ISDIFF.EQ.0) GO TO 40
C
C         PRINT MESSAGES.
C
      WRITE (ILP,300) ATIME,ADATE,TITLE
      IF (IANISO.NE.0.OR.ISDIFF.NE.0) WRITE (ILP,400) (A(I),I=1,6),
     & ((GINV(I,J),J=1,3),I=1,3)
      IF (IANISO.NE.0) C1=IANISO
      IF (IYDIFF.NE.0) C2=IYDIFF
      IF (ISDIFF.NE.0) C3=ISDIFF
      IF (C1.LT.0) C1=-1/C1
      IF (C2.LT.0) C2=-1/C2
      IF (C3.LT.0) C3=-1/C3
      IF (IANISO.NE.0) WRITE (ILP,401) C1
      IF (IYDIFF.NE.0) WRITE (ILP,402) C2
      IF (ISDIFF.NE.0) WRITE (ILP,403) C3
 40   CONTINUE
C
C         INITIALIZE SOME VARIABLES.
C
      WRITE (ILP,309) ATIME,ADATE,TITLE
      NSUM=0
      SUMF=0
      SUME=0
      FMIN=+1E6
      EMIN=+1E6
      FMAX=-1E6
      EMAX=-1E6
      NREAD=0
      NWRIT=0
      NPRINT=0
C
C         LOOP THROUGH REFLECTION DATA FILE.
C
      REWIND IO1
      IEND=0
 1    CALL READ1 (IO1,IEND,II,IH,IK,IL,A1,A2,A3,A4,Y,SIGY,XTIME)
      IF (IEND.NE.0) GO TO 9
      NREAD=NREAD+1
C
C         X**N IS UNDEFINED FOR X = N = 0.
C
      IF (XTIME.EQ.0) XTIME=0.0001
      X=XTIME
C
C         IS THIS REFLECTION TO BE REJECTED FROM THE DATA SET?
C
      IF (NIREJ.EQ.0) GO TO 25
      DO 20 I=1,NIREJ
      I1=IREJ(I,1)
      I2=IREJ(I,2)
      IF (I1.LE.ABS(II).AND.ABS(II).LE.I2) GO TO 1
 20   CONTINUE
 25   IF (NXREJ.EQ.0) GO TO 29
      DO 21 I=1,NXREJ
      X1=XREJ(I,1)
      X2=XREJ(I,2)
      IF (X1.LT.X.AND.X.LT.X2) GO TO 1
 21   CONTINUE
 29   CONTINUE
C
C         FIND THE INVERSE SCALING FACTOR APPLICABLE AT XTIME = X BY
C         AVERAGING OVER THE POLYNOMIAL CURVE SEGMENTS WHOSE DOMAIN
C         INCLUDES X = XTIME.
C
      SUMN=0
      SUMW=0
      SUMWF=0
      SUMWF2=0
      SUMW2V=0
C
C         Q IS AN EXPANSION FACTOR USED TO DEAL WITH THE (RARE) CASE OF
C         A MEASUREMENT TIME NOT IN THE DOMAIN OF ANY SEGMENT.
C         EXTRAPOLATE NO MORE THAN HALF THE AVERAGE TIME BETWEEN
C         MEASUREMENTS OF STANDARDS.
C
      Q=0
 10   IF (Q.GT.0.5) GO TO 1
      DO 11 I=1,NSTD
      DO 11 J=1,NCUT(I)+1
      IF (NDAT(I,J).EQ.0) GO TO 11
      N=NDAT(I,J)
      X1=XLIM(I,J,1)
      X2=XLIM(I,J,2)
      DX=Q*(X2-X1)/N
      X1=X1-DX
      X2=X2+DX
      IF (X.LT.X1) GO TO 11
      IF (X.GT.X2) GO TO 11
      F=0
      VAR=0
      DO 12 K=1,4
      F=F+PAR(I,J,K)*X**(K-1)
      DO 12 L=1,4
      VAR=VAR+X**(K+L-2)*COV(I,J,K,L)
 12   CONTINUE
C
C         VARIANCE CAN BE ZERO OR, DUE TO NUMERICAL IMPRECISION IN
C         COV(L,M), HAVE A SMALL NEGATIVE VALUE WHEN X IS NEAR X = XZERO
C
      VAR=MAX(VAR,1E-9)
C
C         RELATIVE WEIGHT FOR AVERAGING W = W0*W1*W2*W3.
C
C         W0 FROM VARIANCE DUE TO POLYNOMIAL COEFFICIENTS.
C
C         W1 FROM ANGLE BETWEEN RECIPROCAL LATTICE VECTORS FOR GIVEN
C            REFLECTION AND REFERENCE REFLECTION.
C
C         W2 FROM DIFFERENCE BETWEEN INTENSITIES OF GIVEN REFLECTION
C            AND REFERENCE REFLECTION.
C
C         W3 FROM DIFFERENCE BETWEEN SIN(THETA)/LAMBDA VALUES OF GIVEN
C            REFLECTION AND REFERENCE REFLECTION.
C
      W0=1/VAR
      W1=1
      W2=1
      W3=1
      IF (IANISO.EQ.0.AND.ISDIFF.EQ.0) GO TO 50
      U(1)=IHKL(I,1)
      U(2)=IHKL(I,2)
      U(3)=IHKL(I,3)
      V(1)=IH
      V(2)=IK
      V(3)=IL
 50   CONTINUE
C
C         FUNCTION VTMV CALCULATES THE BILINEAR FORM,
C         SCALAR = (VECTOR TRANSPOSE)*(MATRIX)*(VECTOR).
C
      IF (IANISO.EQ.0) GO TO 51
      COSPH=VTMV(U,GINV,V)/SQRT(VTMV(U,GINV,U)*VTMV(V,GINV,V))
      W1=1/(1+C1**2*(1-COSPH**2))
 51   IF (IYDIFF.EQ.0) GO TO 52
      DELTA=Y-YZERO(I)*F
      W2=1/(1+(C2*DELTA)**2)
 52   IF (ISDIFF.EQ.0) GO TO 53
      DELTA=SQRT(VTMV(U,GINV,U))-SQRT(VTMV(V,GINV,V))
      W3=1/(1+(C3*DELTA)**2)
 53   W=W0*W1*W2*W3
      SUMN=SUMN+1
      SUMW=SUMW+W
      SUMWF=SUMWF+W*F
      SUMWF2=SUMWF2+W*F**2
      SUMW2V=SUMW2V+W**2*VAR
 11   CONTINUE
      Q=Q+0.5
      IF (SUMN.EQ.0) GO TO 10
C
C         CALCULATE MEAN INVERSE SCALING FACTOR AND ITS STANDARD
C         DEVIATION.
C
      F=SUMWF/SUMW
      WMSD=SUMWF2/SUMW-F**2
      IF (SUMN.GT.1) WMSD=WMSD/(SUMN-1)
      VAR=SUMW2V/SUMW**2
      VAR=MAX(VAR,WMSD)
      SIGF=SQRT(VAR)
C
C         INVERT INVERSE SCALING FACTOR.
C
      SIGF=SIGF/F**2
      F=1/F
C
C         ACCUMULATE SUMS FOR AVERAGE SCALING FACTOR AND AVERAGE
C         FRACTIONAL ERROR, AND TABULATE EXTREME VALUES.
C
      NSUM=NSUM+1
      SUMF=SUMF+F
      SUME=SUME+SIGF/F
      FMIN=MIN(FMIN,F)
      FMAX=MAX(FMAX,F)
      EMIN=MIN(EMIN,(SIGF/F))
      EMAX=MAX(EMAX,(SIGF/F))
C
C         APPLY SCALING FACTOR TO Y = I/LP VALUE.
C
      VAR=(Y*SIGF)**2+F**2*(SIGY**2+(P*Y)**2)
      SIGFY=SQRT(VAR)
      FY=F*Y
C
C         OUTPUT
C
      CALL WRITE1 (IO2,II,IH,IK,IL,A1,A2,A3,A4,FY,SIGFY,XTIME)
      NWRIT=NWRIT+1
      IF (II.LT.0) GO TO 1
      J=ABS(IH)
      K=ABS(IK)
      L=ABS(IL)
      IF (J.EQ.0.AND.K.EQ.0) GO TO 39
      IF (J.EQ.0.AND.L.EQ.0) GO TO 39
      IF (K.EQ.0.AND.L.EQ.0) GO TO 39
      IF (J.EQ.K.AND.L.EQ.0) GO TO 39
      IF (J.EQ.L.AND.K.EQ.0) GO TO 39
      IF (K.EQ.L.AND.J.EQ.0) GO TO 39
      IF (J.EQ.K.AND.K.EQ.L) GO TO 39
      GO TO 1
 39   WRITE (ILP,310) II,IH,IK,IL,Y,SIGY,FY,SIGFY,XTIME,F,SIGF
      NPRINT=NPRINT+1
      IF (MOD(NPRINT,50).NE.0) GO TO 1
      WRITE (ILP,309) ATIME,ADATE,TITLE
      GO TO 1
 9    CONTINUE
      FAV=SUMF/NSUM
      EAV=SUME/NSUM
      WRITE (ILP,311) FMIN,FMAX,FAV,EMIN,EMAX,EAV
      IF (NWRIT.NE.NREAD) WRITE (ILP,312) NREAD,NWRIT,NREAD-NWRIT
      RETURN
 300  FORMAT('1'/'0PROGRAM TSCALE.  ',A,1X,A,'.  ',A)
 400  FORMAT('0LATTICE PARAMETERS:'/'0        A         B         C',
     &'     ALPHA      BETA     GAMMA'/3F10.4,3F10.3/'0RECIPROCAL SPAC',
     &'E METRIC TENSOR:'/'0G11 G12 G13   ',3F10.6/' G21 G22 G23 = ',
     &3F10.6/' G31 G23 G33   ',3F10.6)
 401  FORMAT('0FOR EACH REFLECTION, AVERAGED SCALE FACTOR WILL BE WEIGHT
     &ED ANISOTROPICALLY.'/'0W1 = 1/(1 + (C1*(1 - COS(PHI)**2)**2), WHER
     &E PHI = ANGLE(H,H0) AND'/'0C1 = ',E10.3)
 402  FORMAT('0FOR EACH REFLECTION, AVERAGED SCALE FACTOR WILL BE WEIGHT
     &ED BY INTENSITY DIFFERENCES.'/'0W2 = 1/(1 + (C2*(I - I0)**2)'/'0C2
     & = ',E10.3)
 403  FORMAT('0FOR EACH REFLECTION, AVERAGED SCALE FACTOR WILL BE WEIGHT
     &ED BY SIN(THETA)/LAMBDA DIFFERENCES.'/'0W3 = 1/(1 + (C3*(S - S0)**
     &)'/'0C3 = ',E10.3)
 309  FORMAT('1'/'0PROGRAM TSCALE.  ',A,1X,A,'.  ',A/'0     I   H   K ',
     &'  L  Y = I/LP  SIGMA(Y) Z = K(T)*Y SIGMA(Z) XTIME = T  K(T) SIGMA
     &(K)'/'      -   -   -   -  --------  -------- ---------- -------',
     &'- ---------  ---- --------')
 310  FORMAT(' ',I6,3I4,4F10.2,F10.3,2F6.3)
 319  FORMAT('+',77X,'  WARNING:  GAP IN ONE OR MORE SCALING CURVES.')
 311  FORMAT('0MINIMUM SCALING FACTOR              KMIN = ',F7.3/
     &       ' MAXIMUM SCALING FACTOR              KMAX = ',F7.3/
     &       '0AVERAGE SCALING FACTOR                 K = ',F7.3/
     &       '0MINIMUM FRACTIONAL ERROR            EMIN = ',F7.3/
     &       ' MAXIMUM FRACTIONAL ERROR            EMAX = ',F7.3/
     &       '0AVERAGE FRACTIONAL ERROR  SIGMA(K)/K = E = ',F7.3)
 312  FORMAT('0NIN  = ',I5,' REFLECTIONS READ FROM THE INPUT FILE.'/' NO
     &UT = ',I5,' REFLECTIONS WRITTEN TO THE OUTPUT FILE.'/'0NREJ = ',
     &I5,' REFLECTIONS REJECTED BASED ON INPUT SERIAL NUMBERS OR EXPOSUR
     &E TIMES,'/'             OR DUE TO A GAP BETWEEN SEGMENTS OF THE SC
     &ALING FUNCTION.')
      END
C-----------------------------------------------------------------------
      SUBROUTINE CULL
C
C         CULLS THE STANDARD REFERENCE REFLECTIONS FROM THE INPUT
C         REFLECTION DATA FILE AND COMPILES SUMMARY STATISTICS ON THE
C         DISTRIBUTIONS OF THE REFERENCE INTENSITIES.
C
      PARAMETER (MSTD=20,MCUT=10,MDAT=500,MOMIT=20,MREJ=50)
      CHARACTER TITLE*80,ATIME*8,ADATE*9
      COMMON /BLOCK1/ TITLE,ATIME,ADATE
      COMMON /BLOCK2/ IO1,IO2,ILP,ITABL,IPLOT,IQSUM,IGRAPH,
     &                IAPPLY,IANISO,IYDIFF,ISDIFF,A(6),GINV(3,3)
      COMMON /BLOCK3/ NSTD,IHKL(MSTD,3),NCUT(MSTD),XCUT(MSTD,MCUT),XZERO
      COMMON /BLOCK4/ NIOMIT(MSTD),IOMIT(MSTD,MOMIT,2),
     &                NXOMIT(MSTD),XOMIT(MSTD,MOMIT,2)
      COMMON /BLOCK6/ NDAT(MSTD,MCUT),XX(MSTD,MCUT,MDAT),
     &                YY(MSTD,MCUT,MDAT),SS(MSTD,MCUT,MDAT)
      DIMENSION SUMN(MSTD),SUMW(MSTD),
     &          SUMY(MSTD),SUMWY(MSTD),SUMY2(MSTD),SUMWY2(MSTD),
     &          IMIN(MSTD),TMIN(MSTD),XMIN(MSTD),YMIN(MSTD),SMIN(MSTD),
     &          IMAX(MSTD),TMAX(MSTD),XMAX(MSTD),YMAX(MSTD),SMAX(MSTD),
     &          H(3)
      DO 10 I=1,MSTD
      SUMN(I)=0
      SUMW(I)=0
      SUMY(I)=0
      SUMWY(I)=0
      SUMY2(I)=0
      SUMWY2(I)=0
      XMIN(I)=+1E6
      XMAX(I)=-1E6
      YMIN(I)=+1E6
      YMAX(I)=-1E6
      DO 10 J=1,MCUT
      NDAT(I,J)=0
 10   CONTINUE
      XTIME1=+1E5
      XTIME2=-1E5
      IF (ITABL.EQ.0) GO TO 11
      WRITE (ILP,102) ATIME,ADATE,TITLE
      WRITE (ILP,101)
      JJ=0
 11   CONTINUE
C
C         LOOP THROUGH REFLECTION DATA FILE.
C
      REWIND IO1
      IEND=0
 1    CALL READ1 (IO1,IEND,II,IH,IK,IL,A1,A2,A3,A4,Y,SIGY,XTIME)
      IF (IEND.NE.0) GO TO 9
C
C         X**N IS UNDEFINED FOR X = N = 0.
C
      IF (XTIME.EQ.0) XTIME=0.0001
      XTIME1=MIN(XTIME1,XTIME)
      XTIME2=MAX(XTIME2,XTIME)
C
C         IS THIS A STANDARD REFERENCE REFLECTION?
C
      IF (II.GT.0) GO TO 1
C
C         LIST MEASUREMENTS OF SETS OF STANDARDS.
C
      IF (ITABL.EQ.0) GO TO 15
      IF (II.NE.JJ-1) WRITE (ILP,109)
      WRITE (ILP,110) II,IH,IK,IL,Y,SIGY,XTIME
      JJ=II
 15   CONTINUE
C
C         WHAT TYPE OF STANDARD IS IT?
C
      IF (NSTD.EQ.0) GO TO 50
      DO 20 I=1,NSTD
      JH=IHKL(I,1)
      JK=IHKL(I,2)
      JL=IHKL(I,3)
      IF (IH.NE.JH) GO TO 20
      IF (IK.NE.JK) GO TO 20
      IF (IL.NE.JL) GO TO 20
C
C         IS IT TO BE OMITTED FROM THE CALCULATION OF SCALING FACTORS?
C
      IF (NIOMIT(I).EQ.0) GO TO 21
      DO 22 J=1,NIOMIT(I)
      I1=IOMIT(I,J,1)
      I2=IOMIT(I,J,2)
      IF (I1.LE.ABS(II).AND.ABS(II).LE.I2) GO TO 1
 22   CONTINUE
 21   IF (NXOMIT(I).EQ.0) GO TO 30
      DO 23 J=1,NXOMIT(I)
      X1=XOMIT(I,J,1)
      X2=XOMIT(I,J,2)
      IF (X1.LT.XTIME.AND.XTIME.LT.X2) GO TO 1
 23   CONTINUE
C
C         IN WHICH TIME SEGMENT IS IT?
C
 30   J=1
      IF (NCUT(I).EQ.0) GO TO 40
      DO 31 J=1,NCUT(I)
      IF (XTIME.LE.XCUT(I,J)) GO TO 40
 31   CONTINUE
      J=NCUT(I)+1
      GO TO 40
 20   CONTINUE
 50   CONTINUE
C
C         IF IT IS A NEW TYPE, EXTEND THE TABLE OF TYPES.
C
      NSTD=NSTD+1
      I=NSTD
      IHKL(I,1)=IH
      IHKL(I,2)=IK
      IHKL(I,3)=IL
      J=1
 40   CONTINUE
C
C         TABULATE DATA FOR SUBSEQUENT ANALYSIS.
C
      NDAT(I,J)=NDAT(I,J)+1
      IF (NDAT(I,J).GT.MDAT) STOP 'TOO MANY MEASUREMENTS OF STANDARDS'
      K=NDAT(I,J)
      XX(I,J,K)=XTIME
      YY(I,J,K)=Y
      SS(I,J,K)=SIGY
C
C         COMPILE SUMMARY STATISTICS FOR EACH TYPE OF STANDARD.
C
      W=1/SIGY**2
      SUMN(I)=SUMN(I)+1
      SUMW(I)=SUMW(I)+W
      SUMY(I)=SUMY(I)+Y
      SUMWY(I)=SUMWY(I)+W*Y
      SUMY2(I)=SUMY2(I)+Y**2
      SUMWY2(I)=SUMWY2(I)+W*Y**2
      XMIN(I)=MIN(XMIN(I),XTIME)
      XMAX(I)=MAX(XMAX(I),XTIME)
      IF (YMIN(I).LE.Y) GO TO 41
      IMIN(I)=II
      TMIN(I)=XTIME
      YMIN(I)=Y
      SMIN(I)=SIGY
 41   CONTINUE
      IF (YMAX(I).GE.Y) GO TO 42
      IMAX(I)=II
      TMAX(I)=XTIME
      YMAX(I)=Y
      SMAX(I)=SIGY
 42   CONTINUE
C
C         LOOP BACK FOR NEXT MEASUREMENT.
C
      GO TO 1
 9    CONTINUE
C
C         CALCULATE AND PRINT REFERENCE TIME OF THE EXPERIMENT.
C
      WRITE (ILP,213) XTIME1,XTIME2,0.5*(XTIME1+XTIME2)
      IF (XZERO.EQ.0) XZERO=0.5*(XTIME1+XTIME2)
C
C         CHECK THAT THE DOMAIN OF EACH STANDARD INCLUDES THE REFERENCE
C         TIME OF THE EXPERIMENT.
C
      DO 19 I=1,NSTD
      IF (SUMN(I).EQ.0) GO TO 19
      IF (XMIN(I).LE.XZERO.AND.XZERO.LE.XMAX(I)) GO TO 19
      WRITE (ILP,214) I,(IHKL(I,K),K=1,3),XMIN(I),XMAX(I)
 19   CONTINUE
C
C         EVALUATE AND PRINT SUMMARY STATISTICS.
C
      WRITE (ILP,102) ATIME,ADATE,TITLE
      WRITE (ILP,210)
      DO 39 I=1,NSTD
      H(1)=IHKL(I,1)
      H(2)=IHKL(I,2)
      H(3)=IHKL(I,3)
      SINTHL=0.5*SQRT(VTMV(H,GINV,H))
      N=SUMN(I)
      IF (N.EQ.0) GO TO 38
      YMEAN=SUMY(I)/SUMN(I)
      WYMEAN=SUMWY(I)/SUMW(I)
      ESD=SQRT(SUMN(I)/SUMW(I))
      IF (N.GT.1) GO TO 37
      RMSD=0
      RWMSD=0
      GO TO 36
 37   RMSD=SQRT(ABS(SUMY2(I)/SUMN(I)-YMEAN**2)*SUMN(I)/(SUMN(I)-1))
      RWMSD=SQRT(ABS(SUMWY2(I)/SUMW(I)-WYMEAN**2)*SUMN(I)/(SUMN(I)-1))
 36   WRITE (ILP,211) I,(IHKL(I,K),K=1,3),SINTHL,
     & IMIN(I),TMIN(I),YMIN(I),SMIN(I),
     & N,YMEAN,RMSD,WYMEAN,RWMSD,ESD,
     & IMAX(I),TMAX(I),YMAX(I),SMAX(I)
      GO TO 39
 38   WRITE (ILP,212) I,(IHKL(I,K),K=1,3),SINTHL
 39   CONTINUE
      RETURN
 102  FORMAT('1'/'0PROGRAM TSCALE.  ',A,1X,A,'.  ',A)
 101  FORMAT ('0STANDARD REFERENCE REFLECTIONS:'/'0         J     H  K',
     &'  L           Y      SIGY     XTIME'/)
 109  FORMAT (' ')
 110  FORMAT (' ',I10,3X,3I3,2X,2F10.2,F10.3)
 210  FORMAT('0SUMMARY STATISTICS ON STANDARD REFERENCE REFLECTIONS:'/
     &'0Y = I/LP'/
     &' W = 1/SIGMA(Y)**2'/
     &'0MEAN  = SUM(Y)/N'/
     &' WMEAN = SUM(W*Y)/SUM(W)'/
     &'0RMSD  = SQRT((N/(N-1))*(SUM(Y**2)/N - (SUM(Y)/N)**2))'/
     &' RWMSD = SQRT((N/(N-1))*(SUM(W*Y**2)/SUM(W) - (SUM(W*Y)/SUM(W))**
     &2))'/
     &' ESD   = SQRT(N/SUM(W))'/
     &'0  I     H  K  L SIN(TH)/L         J     T(HR)         Y  SIGMA(Y
     &)         N      MEAN      RMSD     WMEAN     RWMSD       ESD'/
     &' ',48X,'MIN/MAX'/)
 211  FORMAT(' ',I3,3X,3I3,F10.4,I10,F10.3,2F10.2,I10,5F10.2/' ',25X,
     &I10,F10.3,2F10.2/)
 212  FORMAT(' ',I3,3X,3I3,F10.4,49X,'0'/)
 213  FORMAT('0TMIN = ',F10.3,' HR'/
     &       ' TMAX = ',F10.3,' HR'/
     &       '0MID-TIME OF THE EXPERIMENT:'/
     &       '0TZERO = (TMIN + TMAX)/2 = ',F10.3,' HR')
 214  FORMAT('0WARNING:  THE DOMAIN OF THE FOLLOWING STANDARD REFLECTION
     & DOES NOT INCLUDE TZERO.'/' PERHAPS THE STANDARD SHOULD BE OMITTED
     & FROM THE CALCULATION OF SCALING FACTORS,'/' OR A DIFFERENT REFERE
     &NCE TIME SHOULD BE CHOSEN.'/'   I     H  K  L      TMIN      TMAX'
     &/' ',I3,3X,3I3,2F10.3)
      END
C-----------------------------------------------------------------------
      FUNCTION FCRIT (NFREE)
C
C     TABLE OF F-VALUES FOR ONE VS. NFREE DEGREES OF FREEDOM.
C
C     PROBABILITY, P = P(F, 1, NFREE), THAT A RANDOM SAMPLE OF DATA
C     WOULD GIVE A RATIO OF CHI-SQUARED VALUES LARGER THAN F.
C
C     P. R. BEVINGTON, DATA REDUCTION AND ERROR ANALYSIS FOR THE
C     PHYSICAL SCIENCES, MC GRAW-HILL BOOK CO., NEW YORK, 1969,
C     P. 318.
C
      DIMENSION N(20),F(20)
      DATA N /1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
     &        15, 20, 24, 30, 40, 60, 120, 1000000/
C
C     P = 0.1
C
      DATA F /39.9, 8.53, 5.54, 4.54, 4.06, 3.78, 3.59, 3.46, 3.36,
     &        3.28, 3.23, 3.18, 3.07, 2.97, 2.93, 2.88, 2.84, 2.79,
     &        2.75, 2.71/
C
C     P = 0.05
C
C     DATA F /161,  18.5, 10.1, 7.71, 6.61, 5.99, 5.59, 5.32, 5.12,
C    &        4.96, 4.84, 4.75, 4.54, 4.35, 4.26, 4.17, 4.08, 4.00,
C    &        3.92, 3.84/
C
C     P = 0.01
C
C     DATA F /4050, 98.5, 34.1, 21.2, 16.3, 13.7, 12.2, 11.3, 10.6,
C    &        10.0, 9.65, 9.33, 8.68, 8.10, 7.82, 7.56, 7.31, 7.08,
C    &        6.85, 6.63/
      IF (NFREE.GT.12) GO TO 5
      FCRIT=F(NFREE)
      RETURN
    5 DO 1 I=12,19
      IF (NFREE.GT.N(I+1)) GO TO 1
      DF=F(I+1)-F(I)
      DN=N(I+1)-N(I)
      DFDN=DF/DN
      DN=NFREE-N(I)
      FCRIT=F(I)+DN*DFDN
      RETURN
    1 CONTINUE
      END
C-----------------------------------------------------------------------
      SUBROUTINE GRAPH
C
C         CALCULATES AND PLOTS POINTS OF INVERSE SCALING FUNCTION
C         R = Y/Y0 AND ITS FRACTIONAL ERROR SIGMA(R)/R.
C
      PARAMETER (MSTD=20,MCUT=10,MDAT=500,MOMIT=20,MREJ=50)
      PARAMETER (NX=101,NY=51)
      CHARACTER TITLE*80,ATIME*8,ADATE*9
      COMMON /BLOCK1/ TITLE,ATIME,ADATE
      COMMON /BLOCK2/ IO1,IO2,ILP,ITABL,IPLOT,IQSUM,IGRAPH,
     &                IAPPLY,IANISO,IYDIFF,ISDIFF,A(6),GINV(3,3)
      COMMON /BLOCK3/ NSTD,IHKL(MSTD,3),NCUT(MSTD),XCUT(MSTD,MCUT),XZERO
      COMMON /BLOCK6/ NDAT(MSTD,MCUT),XX(MSTD,MCUT,MDAT),
     &                YY(MSTD,MCUT,MDAT),SS(MSTD,MCUT,MDAT)
      COMMON /BLOCK7/ PAR(MSTD,MCUT,4),COV(MSTD,MCUT,4,4),
     &                XLIM(MSTD,MCUT,2),YZERO(MSTD),P
      DIMENSION X(MDAT),Y(MDAT),S(MDAT),C(4)
C
C         FIND PLOT LIMITS.
C
      XMIN=+1E6
      XMAX=-1E6
      YMIN=+1E6
      YMAX=-1E6
      SUMSSQ=0
      N=0
      DO 20 I=1,NSTD
      DO 20 J=1,NCUT(I)+1
      IF (NDAT(I,J).EQ.0) GO TO 20
      DO 21 K=1,NDAT(I,J)
         XI=XX(I,J,K)
         YI=YY(I,J,K)
      SIGYI=SS(I,J,K)
      XMIN=MIN(XMIN,XI)
      XMAX=MAX(XMAX,XI)
      YMIN=MIN(YMIN,YI-2*SIGYI)
      YMAX=MAX(YMAX,YI+2*SIGYI)
      SUMSSQ=SUMSSQ+SIGYI**2
      N=N+1
 21   CONTINUE
 20   CONTINUE
      RMSSIG=SQRT(SUMSSQ/N)
      YMID=(YMIN+YMAX)/2
      DELY=(NY-1)*RMSSIG
      Y1=YMID-DELY/2
      Y2=YMID+DELY/2
      YMIN=MIN(YMIN,Y1)
      YMAX=MAX(YMAX,Y2)
      DELX=(XMAX-XMIN)/(NX-1)
      DELY=(YMAX-YMIN)/(NY-1)
C
C         PLOT DATA AND SCALING FUNCTION FOR FOR EACH STANDARD.
C
      DO 31 I=1,NSTD
      N=0
      DO 32 J=1,NCUT(I)+1
      IF (NDAT(I,J).EQ.0) GO TO 32
      DO 33 K=1,NDAT(I,J)
      N=N+1
      X(N)=XX(I,J,K)
      Y(N)=YY(I,J,K)
 33   CONTINUE
 32   CONTINUE
      IF (N.EQ.0) GO TO 31
      WRITE (ILP,110) ATIME,ADATE,TITLE
      CALL PLOT3 (I,N,X,Y,XMIN,XMAX,YMIN,YMAX)
      WRITE (ILP,101) (IHKL(I,K),K=1,3),DELX,DELY,RMSSIG
 31   CONTINUE
C
C         PLOT COMPOSITE SCALING FUNCTION AND ITS FRACTIONAL ERROR.
C
      FMIN=+1E6
      FMAX=-1E6
      SUMF=0
      DO 11 I=1,NX
      SUMN=0
      SUMW=0
      SUMWY=0
      SUMWY2=0
      X(I)=0
      Y(I)=0
      S(I)=0
      XI=XMIN+(I-1)*DELX
C
C         Q IS AN EXPANSION FACTOR USED TO DEAL WITH THE (RARE) CASE OF
C         A TIME X NOT IN THE DOMAIN OF ANY SEGMENT.  EXTRAPOLATE NO
C         MORE THAN HALF THE AVERAGE TIME BETWEEN MEASUREMENTS OF
C         STANDARDS.
C
      Q=0
 10   IF (Q.GT.0.5) GO TO 11
      DO 12 J=1,NSTD
      DO 12 K=1,NCUT(J)+1
      IF (NDAT(J,K).EQ.0) GO TO 12
      N=NDAT(J,K)
      X1=XLIM(J,K,1)
      X2=XLIM(J,K,2)
      DX=Q*(X2-X1)/N
      X1=X1-DX
      X2=X2+DX
      IF (XI.LT.X1) GO TO 12
      IF (XI.GT.X2) GO TO 12
      YI=0
      VAR=0
      DO 13 L=1,4
      YI=YI+PAR(J,K,L)*XI**(L-1)
      DO 13 M=1,4
      VAR=VAR+XI**(L+M-2)*COV(J,K,L,M)
 13   CONTINUE
C
C         VARIANCE CAN BE ZERO OR, DUE TO NUMERICAL IMPRECISION IN
C         COV(L,M), HAVE A SMALL NEGATIVE VALUE WHEN X IS NEAR X = XZERO
C
      VAR=MAX(VAR,1E-9)
      WI=1/VAR
      SUMN=SUMN+1
      SUMW=SUMW+WI
      SUMWY=SUMWY+WI*YI
      SUMWY2=SUMWY2+WI*YI**2
 12   CONTINUE
      Q=Q+0.1
      IF (SUMN.EQ.0) GO TO 10
C
C         CALCULATE MEAN INVERSE SCALING FACTOR AND ITS STANDARD
C         DEVIATION.
C
      YI=SUMWY/SUMW
      VAR=1/SUMW
      WMSD=SUMWY2/SUMW-YI**2
      IF (SUMN.GT.1) WMSD=WMSD/(SUMN-1)
      VAR=MAX(VAR,WMSD)
      SIG=SQRT(VAR)
      X(I)=XI
      Y(I)=YI
      S(I)=SIG/YI
      FMIN=MIN(FMIN,YI)
      FMAX=MAX(FMAX,YI)
      SUMF=SUMF+YI
 11   CONTINUE
      FAVG=SUMF/NX
C
C         PLOT COMPOSITE SCALING FUNCTION.
C
      WRITE (ILP,110) ATIME,ADATE,TITLE
      CALL PLOT2 (NX,X,Y,XMIN,XMAX,YMIN,YMAX)
      WRITE (ILP,102) FAVG,FMIN,FMAX,DELX,DELY
C
C         PLOT FRACTIONAL ERROR IN COMPOSITE SCALING FUNCTION.
C
      SUMY=0
      YMIN=0
      YMAX=0
      DO 14 I=1,NX
      YI=S(I)
      SUMY=SUMY+YI
      YMAX=MAX(YMAX,YI)
 14   CONTINUE
      YAVG=SUMY/NX
      DELY=(YMAX-YMIN)/(NY-1)
      WRITE (ILP,110) ATIME,ADATE,TITLE
      CALL PLOT2 (NX,X,S,XMIN,XMAX,YMIN,YMAX)
      WRITE (ILP,103) YAVG,YMAX,DELX,DELY
      RETURN
 110  FORMAT('1'/'0',10X,'PROGRAM TSCALE.  ',A,1X,A,'.  ',A)
 101  FORMAT('0',10X,'Y = INVERSE SCALING FACTOR VERSUS X = TIME (HR).'/
     &       '0',10X,'  H  K  L'/' ',10X,3I3/
     &       '0',10X,'DELTA(X) = ',F7.3,' HR AND DELTA(Y) = ',E10.3,
     &               ' PER GRID POINT.  R.M.S. SIGMA(Y) = 'E10.3'.')
 102  FORMAT('0',10X,'Y = AVERAGE OF INVERSE SCALING FACTORS VERSUS X ='
     &               ' TIME (HR).'/
     &       '0',10X,'AVERAGE SCALING FACTOR = ',F7.3,'; MINIMUM'
     &               ' SCALING FACTOR = ',F7.3,'; MAXIMUM SCALING'
     &               ' FACTOR = ',F7.3,'.'/
     &       '0',10X,'DELTA(X) = ',F7.3,' HR AND DELTA(Y) = ',E10.3,
     &               ' PER GRID POINT.')
 103  FORMAT('0',10X,'Y = FRACTIONAL ERROR IN AVERAGED INVERSE SCALING',
     &               ' FACTOR VERSUS X = TIME (HR).'/
     &       '0',10X,'AVERAGE FRACTIONAL ERROR = ',F7.3,'; MAXIMUM'
     &               ' FRACTIONAL ERROR = ',F7.3,'.'/
     &       '0',10X,'DELTA(X) = ',F7.3,' HR AND DELTA(Y) = ',E10.3,
     &               ' PER GRID POINT.')
      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 METRIC (A,GINV)
C
C        CALCULATE RECIPROCAL METRIC TENSOR.
C
      DIMENSION A(6),GINV(3,3)
      PI=ACOS(-1.0)
      CA=COS(A(4)*PI/180)
      CB=COS(A(5)*PI/180)
      CG=COS(A(6)*PI/180)
      SA=SIN(A(4)*PI/180)
      SB=SIN(A(5)*PI/180)
      SG=SIN(A(6)*PI/180)
      V=A(1)*A(2)*A(3)*SQRT(1-CA**2-CB**2-CG**2+2*CA*CB*CG)
      B1=A(2)*A(3)*SA/V
      B2=A(1)*A(3)*SB/V
      B3=A(1)*A(2)*SG/V
      B4=(CB*CG-CA)/(SB*SG)
      B5=(CA*CG-CB)/(SA*SG)
      B6=(CA*CB-CG)/(SA*SB)
      GINV(1,1)=B1*B1
      GINV(1,2)=B1*B2*B6
      GINV(1,3)=B1*B3*B5
      GINV(2,1)=GINV(1,2)
      GINV(2,2)=B2*B2
      GINV(2,3)=B2*B3*B4
      GINV(3,1)=GINV(1,3)
      GINV(3,2)=GINV(2,3)
      GINV(3,3)=B3*B3
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE PCALC (N,X,Y,S,C,PSQ)
C
C         CALCULATES 'INSTABILITY CONSTANT' P AS DESCRIBED BY
C         MC CANDLISH, STOUT, AND ANDREWS, ACTA CRYST., A31,
C         245-249 (1975).
C
      DIMENSION X(1),Y(1),S(1),C(1)
      SUMDEL=0
      SUMSIG=0
      SUMOBS=0
      DO 1 I=1,N
      SIGY=S(I)
      YOBS=Y(I)
      YCALC=0
      DO 2 J=1,4
      YCALC=YCALC+C(J)*X(I)**(J-1)
    2 CONTINUE
      SUMDEL=SUMDEL+(YOBS-YCALC)**2
      SUMSIG=SUMSIG+SIGY**2
      SUMOBS=SUMOBS+YOBS**2
    1 CONTINUE
      PSQ=(SUMDEL-SUMSIG)/SUMOBS
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE PLOTA (AA)
      PARAMETER (NX=101,NY=51)
      CHARACTER AA*1
      DIMENSION AA(NX,NY)
      DO 10 IX=1,NX
      DO 10 IY=1,NY
      AA(IX,IY)=' '
 10   CONTINUE
      DO 21 IX=1,NX
      AA(IX, 1)='.'
      AA(IX,NY)='.'
 21   CONTINUE
      DO 22 IX=1,NX,10
      AA(IX, 1)='+'
      AA(IX,NY)='+'
 22   CONTINUE
      DO 23 IY=1,NY
      AA( 1,IY)='.'
      AA(NX,IY)='.'
 23   CONTINUE
      DO 24 IY=1,NY,10
      AA( 1,IY)='+'
      AA(NX,IY)='+'
 24   CONTINUE
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE PLOTB (XMIN,XMAX,YMIN,YMAX,AA)
      PARAMETER (NX=101,NY=51)
      CHARACTER AA*1
      DIMENSION AA(NX,NY),XGRID(11)
      COMMON /BLOCK2/ IO1,IO2,ILP,ITABL,IPLOT,IQSUM,IGRAPH,
     &                IAPPLY,IANISO,IYDIFF,ISDIFF,A(6),GINV(3,3)
      DELX=XMAX-XMIN
      DELY=YMAX-YMIN
      F=0
      DO 90 IY=1,NY
      IF (MOD(IY,10).EQ.1) THEN
        YGRID=YMAX-F*DELY
        F=F+0.2
        WRITE (ILP,101) YGRID,(AA(IX,IY),IX=1,NX)
      ELSE
        WRITE (ILP,100)       (AA(IX,IY),IX=1,NX)
      END IF
 90   CONTINUE
 101  FORMAT(' ',E10.3,101A1)
 100  FORMAT(' ',10X,  101A1)
      F=0
      DO 91 I=1,11
      XGRID(I)=XMIN+F*DELX
      F=F+0.1
 91   CONTINUE
      WRITE (ILP,102) XGRID
 102  FORMAT(' ',1X,11F10.2)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE PLOT1 (N,X,Y,XMIN,XMAX,YMIN,YMAX,A)
C
C         PLOT DATA POINTS AND FITTED CURVE FOR THE JTH SEGMENT OF THE
C         ITH STANDARD.
C
      PARAMETER (NX=101,NY=51)
      CHARACTER AA*1
      DIMENSION X(1),Y(1),A(1),AA(NX,NY)
C
C         CLEAR THE PLOT ARRAY AND FILL IN THE GRID OUTLINE.
C
      CALL PLOTA (AA)
C
C         FILL IN THE CALCULATED POINTS.
C
      DELX=XMAX-XMIN
      DELY=YMAX-YMIN
      IF (A(1).EQ.0.AND.A(2).EQ.0.AND.A(3).EQ.0.AND.A(4).EQ.0) GO TO 49
      DX=DELX/(NX-1)
      DO 41 I=1,NX
      XI=XMIN+(I-1)*DX
      YI=0
      DO 42 J=1,4
      YI=YI+A(J)*XI**(J-1)
 42   CONTINUE
      IX=1+((XI-XMIN)/DELX)*(NX-1)+0.49
      IY=1+((YMAX-YI)/DELY)*(NY-1)+0.49
      IF (IX.LT.1.OR.IX.GT.NX) GO TO 41
      IF (IY.LT.1.OR.IY.GT.NY) GO TO 41
      AA(IX,IY)='*'
 41   CONTINUE
 49   CONTINUE
C
C         FILL IN THE EXPERIMENTAL DATA POINTS.
C
      DO 30 I=1,N
      IX=1+((X(I)-XMIN)/DELX)*(NX-1)+0.49
      IY=1+((YMAX-Y(I))/DELY)*(NY-1)+0.49
      IX=MAX(IX,1)
      IX=MIN(IX,NX)
      IY=MAX(IY,1)
      IY=MIN(IY,NY)
      AA(IX,IY)='O'
 30   CONTINUE
C
C         PRINT THE PLOT ARRAY.
C
      CALL PLOTB (XMIN,XMAX,YMIN,YMAX,AA)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE PLOT2 (N,X,Y,XMIN,XMAX,YMIN,YMAX)
C
C         PLOT THE Q-SUM FUNCTION, THE COMPOSITE SCALING FUNCTION, OR
C         THE FRACTIONAL ERROR IN THE COMPOSITE SCALING FUNCTION.
C
      PARAMETER (NX=101,NY=51)
      CHARACTER AA*1
      DIMENSION X(1),Y(1),AA(NX,NY)
C
C         CLEAR THE PLOT ARRAY AND FILL IN THE GRID OUTLINE.
C
      CALL PLOTA (AA)
C
C         FILL IN THE CALCULATED POINTS.
C
      DELX=XMAX-XMIN
      DELY=YMAX-YMIN
      DO 30 I=1,N
      IX=1+((X(I)-XMIN)/DELX)*(NX-1)+0.49
      IY=1+((YMAX-Y(I))/DELY)*(NY-1)+0.49
      IF (IX.LT.1.OR.IX.GT.NX) GO TO 30
      IF (IY.LT.1.OR.IY.GT.NY) GO TO 30
      AA(IX,IY)='*'
 30   CONTINUE
C
C         PRINT THE PLOT ARRAY.
C
      CALL PLOTB (XMIN,XMAX,YMIN,YMAX,AA)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE PLOT3 (ISTD,NDATA,X,Y,XMIN,XMAX,YMIN,YMAX)
C
C         PLOT DATA POINTS AND THE (SEGMENTS OF THE) SCALING FUNCTION
C         FOR THE ITH STANDARD.
C
      PARAMETER (MSTD=20,MCUT=10,MDAT=500,MOMIT=20,MREJ=50)
      PARAMETER (NX=101,NY=51)
      CHARACTER AA*1
      DIMENSION X(1),Y(1),C(4),AA(NX,NY)
      COMMON /BLOCK2/ IO1,IO2,ILP,ITABL,IPLOT,IQSUM,IGRAPH,
     &                IAPPLY,IANISO,IYDIFF,ISDIFF,A(6),GINV(3,3)
      COMMON /BLOCK3/ NSTD,IHKL(MSTD,3),NCUT(MSTD),XCUT(MSTD,MCUT),XZERO
      COMMON /BLOCK6/ NDAT(MSTD,MCUT),XX(MSTD,MCUT,MDAT),
     &                YY(MSTD,MCUT,MDAT),SS(MSTD,MCUT,MDAT)
      COMMON /BLOCK7/ PAR(MSTD,MCUT,4),COV(MSTD,MCUT,4,4),
     &                XLIM(MSTD,MCUT,2),YZERO(MSTD),P
C
C         CLEAR THE PLOT ARRAY AND FILL IN THE GRID OUTLINE.
C
      CALL PLOTA (AA)
C
C         FILL IN THE CALCULATED POINTS.
C
      DELX=XMAX-XMIN
      DELY=YMAX-YMIN
      DX=DELX/(NX-1)
      DO 40 I=1,NX
      XI=XMIN+(I-1)*DX
C
C         FIND THE SEGMENT THAT INCLUDES X = XI IN ITS DOMAIN.
C
      DO 41 JSEG=1,NCUT(ISTD)+1
      X1=XLIM(ISTD,JSEG,1)
      X2=XLIM(ISTD,JSEG,2)
      IF (X1.LE.XI.AND.XI.LE.X2) GO TO 42
 41   CONTINUE
      GO TO 40
 42   CONTINUE
      DO 43 K=1,4
      C(K)=PAR(ISTD,JSEG,K)
 43   CONTINUE
C
C         CALCULATE YI = YI(XI).
C
      YI=0
      DO 44 K=1,4
      YI=YI+C(K)*XI**(K-1)
 44   CONTINUE
      IX=1+((XI-XMIN)/DELX)*(NX-1)+0.49
      IY=1+((YMAX-YI)/DELY)*(NY-1)+0.49
      IF (IX.LT.1.OR.IX.GT.NX) GO TO 40
      IF (IY.LT.1.OR.IY.GT.NY) GO TO 40
      AA(IX,IY)='*'
 40   CONTINUE
C
C         FILL IN THE EXPERIMENTAL DATA POINTS.
C
      DO 30 I=1,NDATA
      IX=1+((X(I)-XMIN)/DELX)*(NX-1)+0.49
      IY=1+((YMAX-Y(I))/DELY)*(NY-1)+0.49
      IX=MAX(IX,1)
      IX=MIN(IX,NX)
      IY=MAX(IY,1)
      IY=MIN(IY,NY)
      AA(IX,IY)='O'
 30   CONTINUE
C
C         PRINT THE PLOT ARRAY.
C
      CALL PLOTB (XMIN,XMAX,YMIN,YMAX,AA)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE INPUT
C
C         READ CONTROL DATA.
C
      PARAMETER (MSTD=20,MCUT=10,MDAT=500,MOMIT=20,MREJ=50)
      CHARACTER TITLE*80,ATIME*8,ADATE*9,FILE1*40,FILE2*40
      COMMON /BLOCK1/ TITLE,ATIME,ADATE
      COMMON /BLOCK2/ IO1,IO2,ILP,ITABL,IPLOT,IQSUM,IGRAPH,
     &                IAPPLY,IANISO,IYDIFF,ISDIFF,A(6),GINV(3,3)
      COMMON /BLOCK3/ NSTD,IHKL(MSTD,3),NCUT(MSTD),XCUT(MSTD,MCUT),XZERO
      COMMON /BLOCK4/ NIOMIT(MSTD),IOMIT(MSTD,MOMIT,2),
     &                NXOMIT(MSTD),XOMIT(MSTD,MOMIT,2)
      COMMON /BLOCK5/ NIREJ,IREJ(MREJ,2),
     &                NXREJ,XREJ(MREJ,2)
      IO1=1
      IO2=2
      ILP=3
      OPEN (UNIT=IO1,NAME='tscale.dat',STATUS='OLD')
      READ (IO1,101) TITLE
      READ (IO1,101) FILE1
      READ (IO1,101) FILE2
      ITABL=1
      IPLOT=1
      IQSUM=0
      IGRAPH=1
      IAPPLY=1
      IANIS0=0
      IYDIFF=0
      ISDIFF=0
      A(1)=0
      NSTD=0
      DO 1000 I=1,MSTD
      NIOMIT(I)=0
      NXOMIT(I)=0
      NCUT(I)=0
 1000 CONTINUE
      NIREJ=0
      NXREJ=0
      ICUT=0
      IIOMIT=0
      IXOMIT=0
      IIREJ=0
      IXREJ=0
      XZERO=0
      READ (IO1,104,END=9) 
     & ITABL,IPLOT,IQSUM,IGRAPH,IAPPLY,IANISO,IYDIFF,ISDIFF
C
C         FOR ANISOTROPIC SCALING OR SIN(THETA)/LAMBDA DEPENDENT
C         SCALING, READ LATTICE PARAMETERS.
C
      READ (IO1,105,END=9) A
C
C         TABULATE (1) ANY STANDARDS THAT ARE TO BE ANALYZED IN SEG-
C         MENTS, (2) AND (3) ANY STANDARDS THAT ARE TO BE OMITTED FROM
C         THE ANALYSIS, AND (4) AND (5) ANY DATA THAT ARE TO BE REJECTED
C         RATHER THAN SCALED.
C
C         EACH OF THE FIVE INPUT LISTS IS ENDED BY A BLANK RECORD.
C
 70   READ (IO1,103,END=9) JH,JK,JL,XTCUT
      IF (JH.EQ.0.AND.JK.EQ.0.AND.JL.EQ.0) GO TO 79
      ICUT=1
      CALL TABULATE (JH,JK,JL,I)
      NCUT(I)=NCUT(I)+1
      IF (NCUT(I).EQ.MCUT) STOP 'TOO MANY EXPOSURE TIME CUTS'
      J=NCUT(I)
C
C          X-RAY EXPOSURE TIMES TO THE NEAREST 0.001 HR = 3.6 SEC
C
      XCUT(I,J)=XTCUT+0.00049
      GO TO 70
 79   CONTINUE
 10   READ (IO1,102,END=9) JH,JK,JL,I1,I2
      IF (JH.EQ.0.AND.JK.EQ.0.AND.JL.EQ.0) GO TO 19
      IF (I2.EQ.0) I2=I1
      IIOMIT=1
      CALL TABULATE (JH,JK,JL,I)
      NIOMIT(I)=NIOMIT(I)+1
      IF (NIOMIT(I).GT.MOMIT) STOP 'TOO MANY IOMIT RANGES'
      J=NIOMIT(I)
      IOMIT(I,J,1)=ABS(I1)
      IOMIT(I,J,2)=ABS(I2)
      GO TO 10
 19   CONTINUE
 20   READ (IO1,103,END=9) JH,JK,JL,X1,X2
      IF (JH.EQ.0.AND.JK.EQ.0.AND.JL.EQ.0) GO TO 29
      IF (X2.EQ.0) X2=X1
      X1=X1-0.00051
      X2=X2+0.00051
      IXOMIT=1
      CALL TABULATE (JH,JK,JL,I)
      NXOMIT(I)=NXOMIT(I)+1
      IF (NXOMIT(I).GT.MOMIT) STOP 'TOO MANY XOMIT RANGES'
      J=NXOMIT(I)
      XOMIT(I,J,1)=X1
      XOMIT(I,J,2)=X2
      GO TO 20
 29   CONTINUE
 30   READ (IO1,104,END=9) I1,I2
      IF (I1.EQ.0.AND.I2.EQ.0) GO TO 39
      IF (I2.EQ.0) I2=I1
      IIREJ=1
      NIREJ=NIREJ+1
      IF (NIREJ.GT.MREJ) STOP 'TOO MANY IREJECT RANGES'
      J=NIREJ
      IREJ(J,1)=ABS(I1)
      IREJ(J,2)=ABS(I2)
      GO TO 30
 39   CONTINUE
 40   READ (IO1,105,END=9) X1,X2
      IF (X1.EQ.0.AND.X2.EQ.0) GO TO 49
      IF (X2.EQ.0) X2=X1
C
C          X-RAY EXPOSURE TIMES TO THE NEAREST 0.001 HR = 3.6 SEC
C
      X1=X1-0.00051
      X2=X2+0.00051
      IXREJ=1
      NXREJ=NXREJ+1
      IF (NXREJ.GT.MREJ) STOP 'TOO MANY XREJECT RANGES'
      J=NXREJ
      XREJ(J,1)=X1
      XREJ(J,2)=X2
      GO TO 40
 49   CONTINUE
C
C         READ IN A VALUE OF XZERO, THE REFERENCE TIME FOR SCALING.
C         IF THE VALUE IS XZERO = 0, THE PROGRAM WILL DEFAULT TO THE
C         MID-TIME OF THE EXPERIMENT, XZERO = 0.5*(XMIN + XMAX).
C
      READ (IO1,105,END=9) XZERO
 9    CLOSE (UNIT=IO1)
      OPEN (UNIT=IO1,NAME=FILE1,STATUS='OLD',FORM='UNFORMATTED',
     & READONLY)
      IF (IAPPLY.NE.0)
     & OPEN (UNIT=IO2,NAME=FILE2,STATUS='NEW',FORM='UNFORMATTED')
      OPEN (UNIT=ILP,NAME='tscale.lp',STATUS='NEW',CARRIAGECONTROL=
     &'FORTRAN')
C
C         PRINT CONTROL DATA.
C
      CALL TIME (ATIME)
      CALL DATE (ADATE)
      WRITE (ILP,201) ATIME,ADATE,TITLE
      WRITE (ILP,202) FILE1
      WRITE (ILP,203) FILE2
      WRITE (ILP,205) 
     & ITABL,IPLOT,IQSUM,IGRAPH,IAPPLY,IANISO,IYDIFF,ISDIFF
      IF (A(1).NE.0) THEN
        WRITE (ILP,212) A
        CALL METRIC (A,GINV)
      END IF
      IF (ICUT.EQ.0) GO TO 61
      WRITE (ILP,220)
      DO 62 I=1,NSTD
      IF (NCUT(I).EQ.0) GO TO 62
      WRITE (ILP,221) I,(IHKL(I,J),J=1,3),(XCUT(I,J),J=1,NCUT(I))
 62   CONTINUE
 61   IF (IIOMIT.EQ.0.AND.IXOMIT.EQ.0) GO TO 51
      WRITE (ILP,204)
      IF (IIOMIT.EQ.0) GO TO 52
      DO 53 I=1,NSTD
      IF (NIOMIT(I).EQ.0) GO TO 53
      DO 43 J=1,NIOMIT(I)
      WRITE (ILP,206) I,(IHKL(I,K),K=1,3),(IOMIT(I,J,K),K=1,2)
 43   CONTINUE
 53   CONTINUE
 52   IF (IXOMIT.EQ.0) GO TO 51
      DO 54 I=1,NSTD
      IF (NXOMIT(I).EQ.0) GO TO 54
      DO 44 J=1,NXOMIT(I)
      WRITE (ILP,207) I,(IHKL(I,K),K=1,3),(XOMIT(I,J,K),K=1,2)
 44   CONTINUE
 54   CONTINUE
 51   IF (IIREJ.EQ.0.AND.IXREJ.EQ.0) GO TO 55
      WRITE (ILP,208)
      IF (IIREJ.EQ.0) GO TO 56
      WRITE (ILP,209) ((IREJ(I,J),J=1,2),I=1,NIREJ)
 56   IF (IXREJ.EQ.0) GO TO 55
      WRITE (ILP,210) ((XREJ(I,J),J=1,2),I=1,NXREJ)
 55   IF (XZERO.EQ.0) GO TO 59
      WRITE (ILP,211) XZERO
 59   RETURN
 101  FORMAT(A)
 102  FORMAT(3I4,2I10)
 103  FORMAT(3I4,2F10.0)
 104  FORMAT(8I10)
 105  FORMAT(8F10.0)
 201  FORMAT('1'/'0PROGRAM TSCALE.  ',A,1X,A,'.  ',A)
 202  FORMAT('0INPUT REFLECTION DATA FILE TO BE SCALED:  ',A)
 203  FORMAT('0OUTPUT SCALED REFLECTION DATA FILE:       ',A)
 205  FORMAT('0JOB CONTROL VARIABLES:'/'0   ITABL   IPLOT   IQSUM  IGRAP
     &H  IAPPLY  IANISO  IYDIFF  ISDIFF'/' ',8I8)
 212  FORMAT('0LATTICE PARAMETERS:'/'0         A         B         C',
     &'     ALPHA      BETA     GAMMA'/' ',3F10.4,3F10.3)
 204  FORMAT('0STANDARD REFERENCE REFLECTION MEASUREMENTS TO BE OMITTED
     & FROM THE ANALYSIS:'/
     &'0  I     H  K  L      FROM        TO'/)
 206  FORMAT(' ',I3,3X,3I3,2I10)
 207  FORMAT(' ',I3,3X,3I3,2F10.3,' HR')
 208  FORMAT('0REFLECTION MEASUREMENTS TO BE REJECTED FROM THE DATA SET:
     &'/
     &'0                     FROM        TO'/)
 209  FORMAT(' ',15X,2I10)
 210  FORMAT(' ',15X,2F10.3,' HR')
 220  FORMAT('0STANDARD REFERENCE REFLECTION DATA TO BE ANALYZED IN TIME
     & SEGMENTS:'/
     &'0  I     H  K  L XTCUT(HR)'/)
 221  FORMAT(' ',I3,3X,3I3,10F10.3)
 211  FORMAT('0INPUT VALUE OF XZERO:'/'0XZERO = ',F10.3,' HR'/'0WILL B',
     &'E USED INSTEAD OF THE MID-TIME OF THE EXPERIMENT AS THE REFERENCE
     & TIME FOR THE SCALING.')
      END
C-----------------------------------------------------------------------
      SUBROUTINE READ1 (IFILE,IEND,II,IH,IK,IL,A1,A2,A3,A4,Y,SIGY,XTIME)
      IEND=0
      READ (IFILE,END=9,ERR=9) II,IH,IK,IL,A1,A2,A3,A4,Y,SIGY,XTIME
C
C     NICOLET (SYNTEX) TIME REGISTER OVERFLOWS AFTER (2**16 - 1) MIN =
C     1092.25 HR.
C
      IF (XTIME.LT.0) XTIME=XTIME+1092.25
      RETURN
    9 IEND=1
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE REGRESS (NDAT,X,Y,SIG,PAR,COV)
C
C         LEAST-SQUARES FIT OF A POWER SERIES POLYNOMIAL BY MULTIPLE
C         LINEAR REGRESSION.
C
C         Y = A0 + A1*X + A2*X**2 + ... + AN*X**N
C
C         RETURNS SCALED INVERSE MATRIX AS VARIANCE-COVARIANCE MATRIX.
C
      DIMENSION X(1),Y(1),SIG(1)
      DIMENSION XMEAN(3),SUMWX(3),SUMWX2(3),SIGX(3),RXY(3),RXX(3,3)
      DIMENSION AXX(3,3),A(3),COVA0(3),PAR(4),COV(4,4)
      DOUBLE PRECISION AXX
      NMAX=3
      DO 1 J=1,NMAX+1
      PAR(J)=0
      DO 1 K=1,NMAX+1
      COV(J,K)=0
    1 CONTINUE
      NMAX=MIN(NMAX,NDAT-2)
C
C         CALCULATE WEIGHTED MEANS AND STANDARD DEVIATIONS OF YOBS AND
C         XOBS**J IN ORDER TO SEPARATE THE A0 PARAMETER AND SCALE THE
C         LEAST-SQUARES NORMAL MATRIX AND VECTOR.
C
      SUMW=0
      SUMWY=0
      SUMWY2=0
      DO 10 J=1,NMAX
      SUMWX(J)=0
      SUMWX2(J)=0
   10 CONTINUE
      DO 11 I=1,NDAT
      WI=1/SIG(I)**2
      YI=Y(I)
      SUMW=SUMW+WI
      SUMWY=SUMWY+WI*YI
      SUMWY2=SUMWY2+WI*YI**2
      DO 11 J=1,NMAX
      XIJ=X(I)**J
      SUMWX(J)=SUMWX(J)+WI*XIJ
      SUMWX2(J)=SUMWX2(J)+WI*XIJ**2
   11 CONTINUE
      YMEAN=SUMWY/SUMW
      SIGY=SQRT(ABS((SUMWY2/SUMW)-YMEAN**2)*NDAT/(NDAT-1))
      DO 19 J=1,NMAX
      XMEAN(J)=SUMWX(J)/SUMW
      SIGX(J)=SQRT(ABS((SUMWX2(J)/SUMW)-XMEAN(J)**2)*NDAT/(NDAT-1))
   19 CONTINUE
C
C         ACCUMULATE NORMAL MATRIX RXX(NMAX,NMAX) AND VECTOR RXY(NMAX).
C
      DO 30 J=1,NMAX
      RXY(J)=0
      DO 30 K=1,J
      RXX(J,K)=0
   30 CONTINUE
      DO 35 I=1,NDAT
      WI=1/SIG(I)**2
      YI=Y(I)
      DO 35 J=1,NMAX
      XIJ=X(I)**J
      RXY(J)=RXY(J)+WI*(XIJ-XMEAN(J))*(YI-YMEAN)
      DO 35 K=1,J
      XIK=X(I)**K
      RXX(J,K)=RXX(J,K)+WI*(XIJ-XMEAN(J))*(XIK-XMEAN(K))
   35 CONTINUE
      DO 39 J=1,NMAX
      RXY(J)=(RXY(J)/SUMW)*NDAT/(NDAT-1)
      RXY(J)=RXY(J)/(SIGX(J)*SIGY)
      DO 39 K=1,J
      RXX(J,K)=(RXX(J,K)/SUMW)*NDAT/(NDAT-1)
      RXX(J,K)=RXX(J,K)/(SIGX(J)*SIGX(K))
      RXX(K,J)=RXX(J,K)
   39 CONTINUE
C
C         TRY POLYNOMIALS OF SUCESSIVELY LOWER DEGREE AND PERFORM THE F-
C         TEST ON CHI-SQUARED VALUES FOR THE SIGNIFICANCE OF EACH HIGHER
C         DEGREE TERM.
C
      CHISQ=9E9
      DO 5 NPAR=NMAX,1,-1
      CHISQ1=CHISQ
C
C         INVERT THE NORMAL MATRIX RXX(NPAR,NPAR).
C
      DO 21 J=1,NPAR
      DO 21 K=1,NPAR
      AXX(J,K)=RXX(J,K)
   21 CONTINUE
      CALL MATINV (NPAR,NMAX,AXX,DET)
      IF (DET.EQ.0) RETURN
C
C         CALCULATE PARAMETERS.
C
      A0=YMEAN
      DO 40 J=1,NPAR
      A(J)=0
      DO 41 K=1,NPAR
      A(J)=A(J)+RXY(K)*AXX(J,K)
   41 CONTINUE
      A(J)=A(J)*SIGY/SIGX(J)
      A0=A0-A(J)*XMEAN(J)
   40 CONTINUE
C
C         CALCULATE CHI-SQUARED.
C
      CHISQ=0
      DO 50 I=1,NDAT
      W=1/SIG(I)**2
      YOBS=Y(I)
      YCALC=A0
      DO 51 J=1,NPAR
      YCALC=YCALC+A(J)*X(I)**J
   51 CONTINUE
      CHISQ=CHISQ+W*(YOBS-YCALC)**2
   50 CONTINUE
C
C         TEST SIGNIFICANCE OF HIGHER DEGREE TERM BY F-TEST USING
C         TABLE FUNCTION FCRIT.
C
      NFREE=NDAT-NPAR-1
      FTEST=(CHISQ-CHISQ1)/(CHISQ1/NFREE)
C
C         IF F-TEST PASSED, ACCEPT HIGHER DEGREE POLYNOMIAL.
C
      IF (FTEST.GE.FCRIT(NFREE)) GO TO 6
    5 CONTINUE
C
C         TRY A CONSTANT, ZEROTH DEGREE POLYNOMIAL, Y = A0 = YMEAN.
C
      NPAR=0
      CHISQ1=CHISQ
      CHISQ=0
      DO 20 I=1,NDAT
      CHISQ=CHISQ+((Y(I)-YMEAN)/SIG(I))**2
   20 CONTINUE
      NFREE=NDAT-1
      FTEST=(CHISQ-CHISQ1)/(CHISQ1/NFREE)
      IF (FTEST.LT.FCRIT(NFREE)) GO TO 99
    6 CONTINUE
C
C         EVALUATE FINAL POLYNOMIAL.
C
      NPAR=NPAR+1
      DO 121 J=1,NPAR
      DO 121 K=1,NPAR
      AXX(J,K)=RXX(J,K)
  121 CONTINUE
      CALL MATINV (NPAR,NMAX,AXX,DET)
      IF (DET.EQ.0) RETURN
C
C         CALCULATE PARAMETERS.
C
      A0=YMEAN
      DO 140 J=1,NPAR
      A(J)=0
      DO 141 K=1,NPAR
      A(J)=A(J)+RXY(K)*AXX(J,K)
  141 CONTINUE
      A(J)=A(J)*SIGY/SIGX(J)
      A0=A0-A(J)*XMEAN(J)
  140 CONTINUE
C
C         CALCULATE CHI-SQUARED.
C
      CHISQ=0
      DO 150 I=1,NDAT
      W=1/SIG(I)**2
      YOBS=Y(I)
      YCALC=A0
      DO 151 J=1,NPAR
      YCALC=YCALC+A(J)*X(I)**J
  151 CONTINUE
      CHISQ=CHISQ+W*(YOBS-YCALC)**2
  150 CONTINUE
C
C         CALCULATE PARAMETER ERROR ESTIMATES.
C
      NFREE=NDAT-NPAR-1
      CHISQ=CHISQ/NFREE
      DO 89 J=1,NPAR
      DO 89 K=1,J
      COV(J,K)=CHISQ*(AXX(J,K)/SUMW)*NDAT/(NDAT-1)
      COV(J,K)=COV(J,K)/(SIGX(J)*SIGX(K))
      COV(K,J)=COV(J,K)
   89 CONTINUE
C
C         SET UP PARAMETER AND VARIANCE-COVARIANCE ARRAYS TO INCLUDE THE
C         A0 ELEMENTS FOR RETURN TO THE CALLING PROGRAM.
C
      DO 90 I=(NPAR+1),2,-1
      PAR(I)=A(I-1)
      DO 90 J=I,2,-1
      COV(I,J)=COV(I-1,J-1)
      COV(J,I)=COV(I,J)
   90 CONTINUE
      PAR(1)=A0
      VARA0=0
      DO 91 J=1,NPAR
      COVA0(J)=0
      DO 91 K=1,NPAR
      COVA0(J)=COVA0(J)-XMEAN(K)*COV(K+1,J+1)
      VARA0=VARA0+XMEAN(J)*XMEAN(K)*COV(J+1,K+1)
   91 CONTINUE
      COV(1,1)=VARA0
      DO 92 I=1,NPAR
      J=1+I
      COV(1,J)=COVA0(I)
      COV(J,1)=COV(1,J)
   92 CONTINUE
      RETURN
   99 CONTINUE
C
C         RETURN A CONSTANT, ZEROTH DEGREE POLYNOMIAL, Y = A0.
C
      PAR(1)=YMEAN
      COV(1,1)=MAX(1/SUMW,((SUMWY2/SUMW)-YMEAN**2)/(NDAT-1))
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE TABULATE (JH,JK,JL,I)
      PARAMETER (MSTD=20,MCUT=10,MDAT=500,MOMIT=20,MREJ=50)
      COMMON /BLOCK3/ NSTD,IHKL(MSTD,3),NCUT(MSTD),XCUT(MSTD,MCUT),XZERO
      IF (NSTD.EQ.0) GO TO 1
      DO 2 I=1,NSTD
      IH=IHKL(I,1)
      IK=IHKL(I,2)
      IL=IHKL(I,3)
      IF (JH.NE.IH) GO TO 2
      IF (JK.NE.IK) GO TO 2
      IF (JL.NE.IL) GO TO 2
      GO TO 3
    2 CONTINUE
    1 CONTINUE
      NSTD=NSTD+1
      IF (NSTD.GT.MSTD) STOP 'TOO MANY TYPES OF STANDARDS'
      I=NSTD
      IHKL(I,1)=JH
      IHKL(I,2)=JK
      IHKL(I,3)=JL
    3 RETURN
      END
C-----------------------------------------------------------------------
      FUNCTION VTMV(U,G,V)
C
C     BILINEAR FORM = (VECTOR TRANSPOSE)*(MATRIX)*(VECTOR) = SCALAR
C     BILINEAR FORM = (ROW MATRIX)*(SQUARE MATRIX)*(COLUMN MATRIX)
C
      DIMENSION U(3),G(3,3),V(3)
      VTMV=0
      DO 1 I=1,3
      DO 1 J=1,3
      VTMV=VTMV+U(J)*G(I,J)*V(I)
    1 CONTINUE
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE WRITE1 (IFILE,II,IH,IK,IL,A1,A2,A3,A4,Y,SIGY,XTIME)
      WRITE (IFILE) II,IH,IK,IL,A1,A2,A3,A4,Y,SIGY,XTIME
      RETURN
      END
