SUBROUTINE POLYFIT ( X, Y, ERROR, YFIT, N,  XLO, XHI, YLO, YHI,
     $                     ORDER, XLINE, YLINE, NLINE, RMS )
C
C  Do a weighted fit of a polynomial to a data set.
C
C  Input data                 X(i), Y(i), i=1, N
C  Values fitted to the data  Y(i), i=1, N
C  Range for valid x data     XLO, XHI
C  Range for valid y data     YLO, YHI
C  Order of fit       ORDER
C
C  This subroutine does not output the coefficients of the fit.
C  However, to facilitate plotting the function, it does output
C  the value of the polynomial, YLINE(i), on a grid of x values
C  XLINE(i), i=1, NLINE input by the user.
C
C  The quality of the fit is evaluated by RMS, the rms of the difference
C  between the data and the fit.
C  
      IMPLICIT UNDEFINED (A-Z)
      INTEGER *4 N, ORDER, NLINE
      REAL *4    X(*), Y(*), ERROR(*), XLINE(*), YLINE(*), YFIT(*)
      REAL *4 YLO, YHI, RMS, MEAN, XLO, XHI, DELT
      INTEGER *4 MAXPAR
      PARAMETER (MAXPAR=20)
      REAL *8 A(MAXPAR,MAXPAR), AA(MAXPAR,MAXPAR), WEIGHT
      REAL *8 CC(MAXPAR), DD(MAXPAR), EE(MAXPAR), BB, VV
      INTEGER *4 I, J, NUNK, COUNT
      SAVE
C.......................................................................
      NUNK = ORDER + 1
      IF ( NUNK+2 .GT. MAXPAR ) THEN
          WRITE(0,*) ' TOO many parameters.  Limit is ',MAXPAR
          NUNK = MAXPAR-2
      END IF
C               LS mode=1, initialize for least squares fit
      WEIGHT = 1.
      CALL LS(1,NUNK,BB,CC,DD,EE,VV, WEIGHT, A, AA, MAXPAR)
      COUNT = 0
      MEAN  = 0.
      RMS   = 0.
      DO 200 I = 1, N
          IF ( (Y(I) .LT. YLO).OR.(Y(I).GT.YHI) ) GO TO 200
          IF ( (X(I) .LT. XLO).OR.(X(I).GT.XHI) ) GO TO 200
          COUNT = COUNT + 1
          MEAN  = MEAN + Y(I)
          RMS   = RMS  + Y(I)*Y(I)
          BB    = Y(I)
          CC(1) = 1.D0
          WEIGHT = 1.D0 / ( ERROR(I) * ERROR(I) )
          DO 150 J = 1, NUNK - 1
              CC(J+1) = CC(J) * X(I)
  150     CONTINUE
C
C              LS mode = 2, increment coefficient matrix for least sq fit.
C                 mode = 3, increment solution vector    for least sq fit.
C                           The weight is included as real*8
C
          CALL LS(2,NUNK,BB,CC,DD,EE,VV, WEIGHT, A, AA, MAXPAR )
          CALL LS(3,NUNK,BB,CC,DD,EE,VV, WEIGHT, A, AA, MAXPAR )
  200 CONTINUE
C
C          LS mode = 4, do the least squares fit
C
      DO 220 J = 2, NUNK
          DD(J) = 0.
  220 CONTINUE
      IF ( COUNT .EQ. 0 ) THEN
          DD(1) = 0.
      ELSE IF ( COUNT .LT. NUNK ) THEN
          DD(1) = MEAN / FLOAT(COUNT)
      ELSE 
          CALL LS(4,NUNK,BB,CC,DD,EE,VV, WEIGHT, A, AA, MAXPAR )
      END IF
      DO 300 I = 1, NLINE
          YLINE(I) = DD(NUNK)
          DO 250 J = 1, NUNK - 1
              YLINE(I) = YLINE(I)*XLINE(I) + DD(NUNK-J)
  250     CONTINUE
  300 CONTINUE
      MEAN = 0.
      RMS  = 0.
      COUNT  = 0
      DO 330 I = 1, N
          YFIT(I) = DD(NUNK)
          DO 320 J = 1, NUNK - 1
              YFIT(I) = YFIT(I)*X(I) + DD(NUNK-J)
  320     CONTINUE
          IF ( Y(I) .LT. YLO ) GO TO 330
          IF ( Y(I) .GT. YHI ) GO TO 330
          IF ( X(I) .LT. XLO ) GO TO 330
          IF ( X(I) .GT. XHI ) GO TO 330
          DELT = (Y(I)-YFIT(I) )
          MEAN = MEAN + DELT
          RMS  = RMS  + DELT*DELT
          COUNT  = COUNT  + 1
  330 CONTINUE
      IF ( COUNT .LE. 1 ) THEN
          RMS  = 0.
      ELSE
          MEAN = MEAN / FLOAT(COUNT)
          RMS  = SQRT ( RMS / FLOAT(COUNT) - MEAN*MEAN )
      END IF
      RETURN
      END
      SUBROUTINE EVENPOLY ( X, Y, ERROR, YFIT, N, ORDER, COEF,
     $                      XLINE, YLINE, NLINE, RMS, CHISQ )
C
C  Do a weighted fit of an even polynomial to a data set.
C
C  Input data                 X(i), Y(i), i=1, N
C  Values fitted to the data  Y(i), i=1, N
C  Errors in y data           ERROR(i),i=1,N
C  Order of fit               ORDER
C  Coefficients of fit        COEF(i), i=1, ORDER
C
C  However, to facilitate plotting the function, it does output
C  the value of the polynomial, YLINE(i), on a grid of x values
C  XLINE(i), i=1, NLINE input by the user.
C
C  The weights are determined as 1/(error) squared
C
C  The quality of the fit is evaluated by RMS, the rms of the difference
C  between the data and the fit and by the chi square.
C  
      IMPLICIT UNDEFINED (A-Z)
      INTEGER *4 N, ORDER, NLINE, NDF
      REAL *4    X(*), Y(*), XLINE(*), YLINE(*), YFIT(*), ERROR(*)
      REAL *4 RMS, MEAN, DELT, CHISQ, RELERR, MINERR, COEF(*)
      INTEGER *4 MAXPAR
      PARAMETER (MAXPAR=20)
      REAL *8 A(MAXPAR,MAXPAR), AA(MAXPAR,MAXPAR), WEIGHT, NEWX
      REAL *8 CC(MAXPAR), DD(MAXPAR), EE(MAXPAR), BB, VV, XMAX
      INTEGER *4 I, J, NUNK, COUNT
      SAVE
      DATA RELERR / .01 /
C.......................................................................
      NUNK = ORDER/2 + 1
      IF ( NUNK+2 .GT. MAXPAR ) THEN
          WRITE(0,*) ' TOO many parameters.  Limit is ',MAXPAR
          NUNK = MAXPAR-2
      END IF
C               LS mode=1, initialize for least squares fit
      WEIGHT = 1.
      CALL LS(1,NUNK,BB,CC,DD,EE,VV, WEIGHT, A, AA, MAXPAR)
C
C  Find maximum value of x and rescale the x axis
C
      XMAX = X(1)
      DO I = 2, N
          XMAX = MAX ( X(I), XMAX )
      END DO
      COUNT = 0
      MEAN  = 0.
      RMS   = 0.
      DO I = 1, N
          COUNT = COUNT + 1
          MEAN  = MEAN + Y(I)
          RMS   = RMS  + Y(I)*Y(I)
          BB    = Y(I)
          MINERR = RELERR * Y(I)
          WEIGHT = 1./(MAX(ERROR(I), MINERR))**2
          CC(1) = 1.D0
          NEWX = X(I) / XMAX
          DO J = 1, NUNK - 1
              CC(J+1) = CC(J) * NEWX * NEWX
          END DO
C
C              LS mode = 2, increment coefficient matrix for least sq fit.
C                 mode = 3, increment solution vector    for least sq fit.
C                           The weight is included as real*8
C
          CALL LS(2,NUNK,BB,CC,DD,EE,VV, WEIGHT, A, AA, MAXPAR )
          CALL LS(3,NUNK,BB,CC,DD,EE,VV, WEIGHT, A, AA, MAXPAR )
      END DO
C
C          LS mode = 4, do the least squares fit
C
      DO J = 2, NUNK
          DD(J) = 0.
      END DO
      IF ( COUNT .EQ. 0 ) THEN
          DD(1) = 0.
      ELSE IF ( COUNT .LT. NUNK ) THEN
          DD(1) = MEAN / FLOAT(COUNT)
      ELSE 
          CALL LS(4,NUNK,BB,CC,DD,EE,VV, WEIGHT, A, AA, MAXPAR )
      END IF
      COEF(1) = DD(1)
      DO I = 1, NUNK-1
          COEF(I+1) = DD(I+1) / (XMAX**(2*I))
      END DO
C      DO I = 1, NLINE
C          YLINE(I) = DD(1)
C          DO J = 1, NUNK-1
C              YLINE(I) = YLINE(I) + DD(J+1) * XLINE(I)**(2*J)
C          END DO
C      END DO
      DO I = 1, NLINE
          YLINE(I) = DD(NUNK)
          NEWX = XLINE(I) / XMAX
          DO J = 1, NUNK-1
              YLINE(I) = DD(NUNK-J) + YLINE(I) * NEWX * NEWX
          END DO
      END DO
      MEAN   = 0.
      RMS    = 0.
      CHISQ  = 0.
      COUNT  = 0
      DO I = 1, N
          YFIT(I) = DD(NUNK)
          NEWX = X(I) / XMAX
          DO J = 1, NUNK - 1
              YFIT(I) = DD(NUNK-J) + YFIT(I) * NEWX * NEWX
          END DO
          DELT = (Y(I)-YFIT(I) )
          MEAN = MEAN + DELT
          RMS  = RMS  + DELT*DELT
          MINERR = RELERR * Y(I)
          DELT = DELT / MAX( ERROR(I), MINERR )
          CHISQ = CHISQ + DELT*DELT
          COUNT  = COUNT  + 1
      END DO
      IF ( COUNT .LE. 1 ) THEN
          RMS  = 0.
      ELSE
          MEAN = MEAN / FLOAT(COUNT)
          RMS  = SQRT ( RMS / FLOAT(COUNT) - MEAN*MEAN )
      END IF
      NDF = COUNT - ORDER - 1
      IF ( NDF .LE. 0 ) THEN
          CHISQ = 99.999
      ELSE
          CHISQ = CHISQ / FLOAT(NDF)
      END IF
      RETURN
      END