' "POLYNOMIAL CURVE FITTER"
' "FROM WILLIAM G. HOOD, BYTE, JUNE 1987"
' PRINT "MODIFIED BY DALE HOLT"


 LN=1000 : LD=11 ' LN=MAX DATA POINTS; LD=HIGHEST DEGREE + 1
 
 N=0: M=0: R2=0: MF=0
 S1=0: S2=0: S3=0: S4=0: P1=0: P2=0
 P3=0: I=0: J=0: J1=0: K=0
 VR=0: MM=0: WT=0: P=0
 DIM X(LN), Y(LN), W(LN), C(LD)
 DIM D1(LD), D2(LD), D3(LD), D4(LD), D5(LD), D6(LD)

' GOSUB MATH 'CALLS THE SUBROUTINE
' INPUT:
' N=# OF DATA POINTS
' X( )=X-COORDINATES OF THE DATA POINTS
' Y( )=Y-COORDINATES OF THE DATA POINTS
' W( )=WEIGHTING FACTORS OF THE DATA POINTS
' M=DEGREE OF POLYNOMIAL
' MF=0 IF NEW DATA, MF=1 IF OLD DATA BUT HIGHER DEGREE
 
' OUTPUT:
' C=ARRAY OF M+1 COEFFICIENTS
' S2=RESIDUAL VARIANCE
' R2=COEFFICIENT OF DETERMINATION

START:
  CLS
  PRINT "POLYFIT. COPYRIGHT (C) 1986 BY WILLIAM G. HOOD"
  PRINT "MODIFIED IN 1987 BY DALE HOLT"
  PRINT
  PRINT "THIS PROGRAM FINDS THE COEFFICIENTS OF THE NTH DEGREE POLY";
  PRINT "NOMIAL"
  PRINT
  PRINT "Y = C(1)*X^N + C(2)*X^(N-1) + ... + C(N)*X + C(N+1)
  PRINT
  PRINT "THAT FITS A SET OF DATA POINTS IN A LEAST-SQUARES SENSE."
  PRINT"EACH DATA POINT MUST CONSIST OF AN X-VALUE, A Y-VALUE, AND AN
  PRINT"OPTIONAL WEIGHT, SEPARATED BY COMMAS AND TERMINATED BY A";
  PRINT " RETURN."
  PRINT "DATA IS READ FROM A DISK FILE UNTIL THE EOF IS REACHED, ";
  PRINT "OR A"
  PRINT "SPECIFIED NUMBER OF DATA POINTS ARE READ IN FROM THE ";
  PRINT "KEYBOARD."
  PRINT: PRINT "MODIFICATIONS BY DALE HOLT PERMIT SAVING A FILE";
  PRINT " OF POINTS"
  PRINT "GENERATED BY THE POLYNOMIAL WHICH HAS BEEN CALCULATED."
  PRINT "THE ORIGINAL DATA POINTS CAN THEN BE PLOTTED ALONG WITH";
  PRINT " THE CURVES"
  PRINT "GENERATED BY THE LEAST SQUARES FIT FOR A VISUAL CHECK OF FIT
  PRINT
  
  INPUT "NAME OF THE DATA FILE (KYBD: = KEYBOARD)? ",FI$
  INPUT "IS THE DATA WEIGHTED (Y/N)? ",W$
  IF W$<>"Y" AND W$<>"y" THEN W$="N"
  IF FI$<>"KYBD:" THEN L1190
L810:
  PRINT "Keyboard data entry"
  PRINT "How many data points (2 -";LN;")";: INPUT N
  IF N=0 THEN N=LN
  IF N<2 THEN PRINT "INVALID! ";: GOTO L810
  FOR K=1 TO N
    IF W$<>"N" THEN INPUT "X,Y,W? ",X(K),Y(K),W(K): W(K)=ABS(W(K))
    IF W$="N" THEN INPUT "X,Y? ",X(K),Y(K): W(K)=1
  NEXT K
L890:
 PM=LD-1: IF N-1<PM THEN PM=N-1
L900:
  PRINT "Degree of polynomial fit (1-";PM;")? ";
  INPUT "",M : M=ABS(M)
  IF M<1 OR M>PM THEN PRINT "INVALID! ";: GOTO L900
  GOSUB MATH
  PRINT: PRINT "Power Series Coefficients:" 
  FOR J=M+1 TO 1 STEP -1
  PRINT "    ";C(J);TAB(22);"X^";M+1-J
  NEXT J
  PRINT: PRINT "Residual variance: ";S2
   
  PRINT: PRINT "Coefficient of Determination:";R2
L1010:
 INPUT "Make a file of calculated points (Y/N)? ",A$
  IF A$="" THEN A$="N"
  A$=UCASE$(A$)
  IF A$="N" THEN L1130 
  IF A$<>"Y" THEN PRINT "INVALID! ";: GOTO L1010
L1040:
 INPUT "File name? ",F$
 IF F$="" THEN PRINT "INVALID! ";: GOTO L1040
 OPEN F$ FOR OUTPUT AS #5
 FOR I=1 TO N: P=C(1)
   FOR K=1 TO M
     P=P*X(I)+C(K+1)
   NEXT K
   ' PRINT X(I);",";P       ' PRINT TO SCREEN ALSO
   PRINT #5, X(I);",";P
 NEXT I
 CLOSE 5
L1130:
  INPUT "Try another degree? ",A$
  IF A$="" THEN A$="N"
  A$=UCASE$(A$)
  IF A$="N" THEN END
  IF A$="Y" THEN GOTO L900
  PRINT "INVALID! ";: GOTO L1130

L1190:
  PRINT "Reading from file: ";FI$
  OPEN FI$ FOR INPUT AS #1
  N=0
  PRINT "X";TAB(15);"Y";TAB(28);"W";:PRINT
  WHILE NOT EOF(1) 
    N=N+1
    IF W$<>"N" THEN INPUT #1,X(N),Y(N),W(N): W(N)=ABS(W(N))
    IF W$="N" THEN INPUT #1,X(N),Y(N): W(N)=1
    PRINT X(N);TAB(14);Y(N);TAB(28);W(N)
  WEND
  CLOSE 1
  PRINT: PRINT N;" points were read from the file."
  GOTO L890
  
  GOTO START
 
 
' CALCULATION SUBROUTINE
MATH:
 IF MF>0 AND M>MM THEN J1=MM+1: MM=M: GOTO L200
 J1=1: MM=M: S1=0: S2=0: S3=0: S4=0
 FOR I=1 TO N
   WT=W(I)
   S1=S1+WT*X(I): S2=S2+WT: S3=S3+WT*Y(I): S4=S4+WT*Y(I)*Y(I)
 NEXT I
 D4(1)=S1/S2: D5(1)=0: D6(1)=S3/S2: D1(1)=0: D2(1)=1
 VR=S4-S3*D6(1)
L200:
 FOR J=J1 TO MM
   S1=0: S2=0: S3=0: S4=0
   FOR I=1 TO N
     P1=0: P2=1
     FOR K=1 TO J
       P=P2: P2=(X(I)-D4(K))*P2-D5(K)*P1: P1=P
     NEXT K
     WT=W(I): P=WT*P2*P2
     S1=S1+P*X(I): S2=S2+P: S3=S3+WT*P1*P1: S4=S4+WT*Y(I)*P2
   NEXT I
   D4(J+1)=S1/S2: D5(J+1)=S2/S3: D6(J+1)=S4/S2
   D3(1)=-D4(J)*D2(1)-D5(J)*D1(1)
   IF J=>4 THEN 
     FOR K=2 TO J-2
       D3(K)=D2(K-1)-D4(J)*D2(K)-D5(J)*D1(K)
     NEXT K
   END IF
   IF J>2 THEN D3(J-1)=D2(J-2)-D4(J)*D2(J-1)-D5(J)
   IF J>1 THEN D3(J)=D2(J-1)-D4(J)
   FOR K=1 TO J
     D1(K)=D2(K): D2(K)=D3(K): D6(K)=D6(K)+D3(K)*D6(J+1)
   NEXT K
 NEXT J
 FOR J=1 TO M+1
   C(J)=D6(M+2-J)
 NEXT J
 P2=0
 FOR I=1 TO N
   P=C(1)
   FOR J=1 TO M
     P=P*X(I)+C(J+1)
   NEXT J
   P=P-Y(I): P(2)=P2+W(I)*P*P
 NEXT I
 S2=0: IF N>M+1 THEN S2=P2/(N-M-1)
 R2=1: IF VR<>0 THEN R2=1-P2/VR: IF R2<0 THEN R2=0
 RETURN

 END