	SUBROUTINE PFIT(X, Y, NPTS, NORDER, A, B, NS, S)
C
C
C	PFIT PERFORMS A LEAST SQUARES POLYNOMIAL REGRESSION
C	ON THE GIVEN DATA TO DETERMINE THE BEST FIT
C	EQUATION TO THE DATA.
C
C
C	PROGRAMMER:
C
C
C	INPUT:
C
C	X - ABSCISSA VALUES FOR THE POLYNOMIAL FIT.
C
C	Y - ORIGINAL ORDINATE VALUES FOR EACH RESPECTIVE X.
C
C	NPTS - THE NUMBER OF POINTS TO BE INCLUDED IN THE
C	       INTERVAL BEING SMOOTHED. (IE. N-POINT SMOOTH)
C
C	NORDER - ORDER OF THE POLYNOMIAL TO BE GENERATED.
C	    NEGITAVE - EXAMINE POLYNOMIALS WITH ORDERS
C	              FROM 1 TO IABS(NORDER) AND STOP
C	              WHEN THE MINIMUM STANDARD DEVIATION
C	              IS FOUND, EVEN IF LESS THAN NORDER.
C	              (IE. 1,2,3,....,NORDER)
C
C	    ZERO - IGNORE THE CALL AND PROCEDE WITH PROGRAM.
C
C	    POSITIVE - EXAMINE THE POLYNOMIAL WITH ORDER
C                 NORDER ONLY. (IE. 2ND DEGREE ONLY)
C
C	OUTPUT:
C
C	A - Y INTERCEPT OF THE POLYNOMIAL.
C
C	B - ARRAY CONTAINING THE COEFFICIENTS OF THE
C	    POLYNOMIAL.
C
C	NS - ORDER OF THE POLYNOMIAL WITH THE MINIMUM
C	     STANDARD DEVIATION FROM THE DATA.
C
C	S - THE STANDARD DEVIATION OF THE POLYNOMIAL.
C
C
C	*** NOTE ***
C
C	THE POLYNOMIAL HAS THE FORM OF:
C
C	ORDINATE = A + B(I) * (ABSCISSA ** I)
C	    FOR I = 1 TO NS
C
C
	IMPLICIT REAL*8 (A-H,O-Z)
	REAL*8 A, B, C, D, EPS, REGR, X, Y
	DIMENSION B(10),D(10),X(100),Y(100)
	DATA EPS/1.0D-20/
	DATA D/10*0.D0/
C
C	... WHICH OPTION IS IN EFFECT ...
C
	IF(NORDER) 20,50,10
C
C	... FIT WITH POLYNOMIAL OF ORDER NORDER ONLY ...
C
10	S = REGR(NPTS,X,Y,NORDER,A,B)
	GO TO 50
C
C	... CHECK ALL POLYNOMIALS FROM 1 TO NORDER
C	    AND RETAIN THE ONE WITH LOWEST S VALUE ...
C
20	N = IABS(NORDER)
	NS = 1
C
C	... GENERATE FIRST ORDER AS A BASIS ...
C
	SMALL = REGR(NPTS,X,Y,NS,A,B)
	IF(N.EQ.1) GO TO 40
C
C	... SCAN OTHER POLYNOMIALS FOR MINIMUM S ...
C
	DO 35 I1 = 2,N
	I=I1
	IF(SMALL.LE.EPS) GO TO 40
	S = REGR(NPTS,X,Y,I,C,D)
	IF(S.GE.SMALL) GO TO 35
C
C	... SAVE DETAILS ABOUT THE POLYNOMIAL ...
C
	SMALL = S
	NS = I
C
C	... SAVE COEFFICIENTS ...
C
	A = C
	DO 30 J =1,N
	B(J) = D(J)
30	CONTINUE
35	CONTINUE
40	S = SMALL
50	RETURN
	END
      FUNCTION REGR(M, X, Y, N, A, B)
C
C      THIS PROGRAM ORIGINALLY CAME FROM THE BOOK:
C
C      B. CARNAHAN, H. LUTHER, AND J. WILKES, APPLIED NUMERICAL
C      METHODS, JOHN WILEY AND SONS, NEW YORK, 1969, P 580
C
C
C              FUNCTION REGR, WHICH CARRIES OUT AN NTH-ORDER
C              POLYNOMIAL REGRESSION ON M DATA POINTS CONTAINED
C              IN THE X AND Y ARRAYS.  THE FUNCTION NORMALLY
C              RETURNS THE STANDARD DEVIATION S OF THE POINTS
C              ABOUT THE REGRESSION CURVE.  THE REGRESSION
C              COEFFICIENTS ARE PLACED IN A AND B(1)...B(N),.
C              HOWEVER, IF THE SIMULTANEOUS EQUATION SOLVING 
C              ROUTINE SIMUL ENCOUNTERS A NEAR-SINGULAR MATRIX,
C             THE FUNCTION RETURNS THE VALUE 0.0 AS REGR.
C
C
C      PARAMETER DEFINITIONS:
C
C      M - THE NUMBER OF DATA POINTS IN X AND Y.
C
C      X,Y - THE ARRAYS CONTAINING THE ABSCISSA AND ORDINATE
C            DATA VALUES RESPECTIVELY.
C
C      N - THE ORDER OF THE POLYNOMIAL DESIRED FOR THE FIT.
C
C      A - Y INTERCEPT OF THE POLYNOMIAL.
C
C      B - ARRAY CONTAINING THE COEFFICIENTS OF THE POLYNOMIAL.
C
C      REGR - RETURNS AS THE VALUE OF THE STANDARD
C             DEVIATION OF THE FIT TO THE DATA.  IF
C             A NEAR SINGULAR OR IF A SINGULAR
C             MATRIX IS ENCOUNTERED THE VALUE OF
C             0 (ZERO) IS RETURNED AS REGR.
C
C
C      **** NOTE ****
C
C      THE POLYNOMIAL HAS THE FORM OF:
C
C      ORDINATE = A + B(I) * (ABSCISSA**I)
C          FOR I=1 TO N
C
C
      IMPLICIT REAL*8 (A - H, O-Z)
      REAL*8 A, B, REGR, X, Y
      DIMENSION C(11,11), SX(20), SYX(10), CYX(10), X(100), Y(100),
     1B(10)
      DATA EPS/ 1.0D-20/
C
C      ..... COMPUTE SUMS OF POWERS AND PRODUCTS .....
      NTWO = 2*N
      NP1 =N + 1
      SY = 0.0
      SYY = 0.0
      DO 1 I=1,N
      NPI = N + I
      SX(I) = 0.0
      SX( NPI) = 0.0
      SYX(I) = 0.0
1      CONTINUE
      DO 3 I = 1,M
      SY = SY + Y(I)
      SYY = SYY + Y( I)**2
      DUM=1.0
      DO 2 J=1,N
      DUM = DUM*X(I)
      SX(J) = SX(J) + DUM
      SYX(J) = SYX(J) + Y(I) *DUM
2      CONTINUE
      DO 3 J=NP1,NTWO
      DUM = DUM*X(I)
      SX(J) = SX(J) + DUM
3      CONTINUE
C
C      ..... COMPUTE COEFFICIENTS C(I,J) .....
C
      FM = M
      CYY = SYY - SY*SY/FM
      DO 4 I =1,N
      CYX(I) =SYX(I) - SY*SX(I)/FM
      C(I,NP1) = CYX(I)
      DO 4 J=1,N
      IPJ = I + J
      C(I,J) = SX(IPJ) - SX(I)*SX(J)/FM
4      CONTINUE
C
C      ..... CALL ON SIMUL TO SOLVE SIMULTANEOUS EQUATIONS ..
C
      DET = SIMUL (N, C, B, EPS, 1, 11)
C
C      I INSERTED THE COMMENTS INFRONT OF THE WRITE
C      STATEMENT BECAUSE I DIDN'T WANT TO PRINT IT OUT.
C
C      WRITE(5,200) DET
C200      FORMAT(13H0       DET = ,E14.8)
      IF(DET.NE.0.0) GO TO 6
      REGR = 0.0
      RETURN
C
C      ..... COMPUTE INTERCEPT A AND STANDARD DEVIATION S .....
C
6      DUM = SY
      TEMP = CYY
      DO 7 I = 1,N
      DUM = DUM - B(I)*SX(I)
      TEMP = TEMP - B(I)*CYX(I)
7      CONTINUE
      A = DUM/FM
      DENOM = M - N - 1
C
C      ***** NOTE *****
C
C      I INSERTED THE FOLLOWING IF STATEMENT TO ELIMINATE
C      TEMP BEING NEAR ZERO BUT SLIGHTLY NEGATIVE.
C      AN EXAMPLE IS IF TEMP EQUALS -.12345E-15.
C
C      THIS STATEMENT IS NOT SUPPLIED IN
C      THE BOOK VERSION, BUT RATHER IS A HOMEBREW PATCH.
C
C      ****************
      IF(TEMP.LE.EPS) TEMP=EPS
      S = DSQRT(TEMP/DENOM)
      REGR = S
C:BUG      WRITE( 6, 91) N, REGR, A, B
   91 FORMAT(' REGR: N, REGR, A: ', I4, 2E15.5,/,
     &        ' B: ',10E10.3)
      RETURN
      END

      FUNCTION SIMUL( N, A, X, EPS, INDIC, NRC)
C
C      THIS PROGRAM ORIGINALLY CAME FROM THE BOOK:
C
C      B. CARNAHAN, H. LUTHER, AND J. WILKES, APPLIED NUMERICAL
C      METHODS, JOHN WILEY AND SONS, NEW YORK, 1969, P290
C
C
C          WHEN INDIC IS NEGATIVE, SIMUL COMPUTES THE INVERSE
C          OF THE N BY N MATRIX A IN PLACE.  WHEN INDIC IS ZERO,
C          SIMUL COMPUTES THE N SOLUTIONS X(1)...X(N) CORRESPONDING
C          TO THE SET OF LINEAR EQUATIONS WITH AUGEMENTED MATRIX
C          OF COEFFICIENTS IN THE N BY N+1 ARRAY A AND IN
C          ADDITION COMPUTES TH INVERSE OF THE COEFFICIENT
C          MATRIX IN PLACE AS ABOVE.  IF INDIC IS POSITIVE,
C          THE SET OF LINERAR EQUATIONS IS SOLVED BUT THE
C          INVERSE IS NOT COMPUTED IN PLACE.  THE GAUSS-
C          JORDAN COMPLETE ELIMINATION METHOD IS EMPLOYED
C          WITH THE MAXIMUM PIVOT STRATEGY.  ROW AND
C          COLUMN SUBSCRIPTS OF SUCCESIVE PIVOT ELEMENTS ARE
C          SAVED IN ORDER IN THE IROW AND JCOL ARRAYS
C          RESPECTIVELY.  K IS THE PIVOT COUNTER, PIVOT IS
C          THE ALGEBRAIC VALUE OF THE PIVOT ELEMENT.  MAX
C          THE NUMBER OF COLUMNS IN A AND DETER THE DETERMINANT
C          OF THE COEFFICIENT MATRIX.  THE SOLUTIONS ARE
C          DETERMINED IN THE N+1TH COLUMN OF A AND THEN
C          UNSCRAMBLED AND PUT IN ORDER IN X(1)...X(N)
C          USING THE PIVOT SUBSCRIPT INFORMATION AVAILABLE
C          IN THE IROW AND JCOL ARRAYS.  THE SIGN OF THE
C          DETERMINANT IS ADJUSTED, IF NECESSARY, BY DETERMINING
C          IF AN EVEN OR ODD NUMBER OF PAIRWISE INTERCHANGES IS
C      REQUIRED TO PUT THE ELEMENTS OF THE JORD ARRAY IN
C          ASCENDING SEQUENCE WHERE JORD(IROW(I)) = JCOL(I).
C          IF THE UNIVERSE IS REQUIRED, IT IS UNSCRAMBLED IN PLACE
C          USING Y(1)...Y(N) AS TEMPORARY STORAGE.  THE VALUE OF
C          THE DETERMINANT IS GIVEN AS THE VALUE OF THE FUNCTION.
C          SHOULD THE POTENTIAL PIVOT OF LARGEST MAGNITUDE BE
C          SMALLER IN MAGNITUDE THAN EPS, THE MATRIX IS
C          CONSIDERED TO BE SINGULAR AND A TRUE ZERO IS RETURNED
C          AS THE VALUE OF THE FUNCTION.
C
      IMPLICIT REAL*8(A-H, O-Z)
      REAL*8 A, X, EPS, SIMUL
      DIMENSION IROW(50),JCOL(50),JORD(50),Y(50),A(NRC,NRC),X(N)
C
      MAX = N
      IF(INDIC.GE.0) MAX = N+1
C
C      .....IS N LARGER THAN 50 .....
      IF(N.LE.50) GO TO 5
      WRITE(5,200)
      SIMUL = 0.0
      RETURN
C
C      ..... BEGIN ELIMINATION PROCEDURE .....
5      DETER = 1.
      DO 18 K = 1,N
      KM1 = K -1
C
C      ..... SEARCH FOR THE PIVOT ELEMENT .....
      PIVOT = 0.0
      DO 11I = 1,N
      DO 11 J = 1,N
C      ..... SCAN IROW AND JCOL ARRAYS FOR INVALID PIVOT SUBSCRIPTS....
      IF(K.EQ.1) GO TO 9
      DO 8 ISCAN = 1,KM1

      DO 8 JSCAN = 1,KM1
      IF(I.EQ.IROW(ISCAN)) GO TO 11
      IF(J.EQ.JCOL(JSCAN)) GO TO 11
8      CONTINUE
9      IF(DABS(A(I,J)).LE.DABS(PIVOT)) GO TO 11
      PIVOT = A(I,J)
      IROW(K) = I
      JCOL(K) = J
11      CONTINUE
C
C      ..... INSURE THAT SELECTED PIVOT IS LARGER THAN EPS ....
      IF(DABS(PIVOT).GT.EPS) GO TO 13
      SIMUL = 0.0
      RETURN
C
C      ..... UPDATE THE DETERMINANT VALUE .....
13      IROWK = IROW(K)
      JCOLK = JCOL(K)
      DETER = DETER*PIVOT
C
C      ..... NORMALIZE PIVOT ROW ELEMENTS .....
      DO 14 J = 1,MAX
      A(IROWK,J) = A (IROWK,J)/PIVOT
14      CONTINUE
C
C      ..... CARRY OUT ELIMINATION AND DEVELOP INVERSE .....
      A(IROWK,JCOLK) = 1./PIVOT
      DO 18 I = 1,N
      AIJCK = A(I,JCOLK)
      IF(I.EQ.IROWK) GO TO 18
      A(I,JCOLK) = - AIJCK/PIVOT
      DO 17 J = 1,MAX
      IF(J.NE.JCOLK) A(I,J) = A(I,J) -AIJCK*A(IROWK,J)
17      CONTINUE
18      CONTINUE
C
C      ..... ORDER SOLUTION VALUES (IF ANY) AND CREATE JORD ARRAY .....
      DO 20 I = 1,N
      IROWI = IROW(I)
      JCOLI = JCOL(I)
      JORD(IROWI) = JCOLI
      IF(INDIC.GE.0) X(JCOLI) = A(IROWI,MAX)
20      CONTINUE
C
C      ..... ADJUST SIGN OF DETERMINANT .....
      INTCH = 0
      NM1 = N -1
      DO 22 I = 1,NM1
      IP1 = I + 1
      DO 22 J = IP1,N
      IF(JORD(J).GE.JORD(I)) GO TO 22
      JTEMP = JORD(J)
      JORD(J) = JORD(I)
      JORD(I) = JTEMP
22      CONTINUE
      IF(INTCH.NE.INTCH/2*2) DETER = -DETER
C
C      ..... IF INDIC IS POSITIVE RETURN WITH RESULTS.....
24      IF(INDIC.LE.0) GO TO 26
      SIMUL = DETER
      RETURN
C
C      .... IF INDIC IS NEGATIVE OR ZERO, UNSCRAMBLE THE INVERSE
C            FIRST BY ROWS .....
26      DO 28 J = 1,N
      DO 27 I = 1,N
      IROWI = IROW(I)
      JCOLI = JCOL(I)
      Y(JCOLI) = A(IROWI,J)
27      CONTINUE
      DO 28 I = 1,N
      A(I,J) = Y(I)
28      CONTINUE
C      ..... THEN BY COLUMNS ....
      DO 30 I = 1,N
      DO 29 J = 1,N
      IROWJ = IROW(J)
      JCOLJ = JCOL(J)
      Y(IROWJ) = A(I,JCOLJ)
29      CONTINUE
      DO 30 J = 1,N
      A(I,J) = Y(J)
30      CONTINUE
C
C      ..... RETURN FOR INDIC NEGATIVE OR ZERO ...
      SIMUL = DETER
      RETURN
C
C      ..... FORMAT FOR OUTPUT STATEMENT .....
200      FORMAT(10HON TOO BIG)
      END
