10 REM
20 REM   POLYNOMIAL  CURVFIT
30 REM   By Ralph Johnson 1-1982
40 REM
50 REM  MODIFIED BY BOB SCOTT MAR 1983
60 REM  TO DRIVE THE RADIO SHACK CGP-115 
70 REM  PRINTER/PLOTTER.  PROGRAM NOW
80 REM  PRIVIDES GRAPHIC CURVE OUTPUT
90 REM  AND PLOTS INPUT DATA points.
100 PRINT "PRESS RETURN WHEN PLOTTER READY":INPUT AZ$
101 POKE 3,(PEEK(3) AND 63) OR 64
110       REM  MAIN PROGRAM
120 GOSUB 260 	REM SET UP ARRAYS
130 GOSUB 390  REM INPUT DATA
140 GOSUB 680  REM CLEAR AND LOAD ARRAYS
150 GOSUB 870  REM SIMUL. EQS. SOLUTION
160 GOSUB 1340  REM SUM OF ERRORS SQRD
170 GOSUB 1470  REM OUTPUT
180 GOTO 1750   REM DECISIONS
190 GOSUB 1890  REM EVALUATE
200 GOSUB 2220  REM PLOT
210 GOTO 1750
220	REM END MAIN
230 REM
240	REM SET UP ARRAY MEMORY
250 REM
260 DIM A(20,20),B(20),C(20),X(20),Y(20)
270 DIM XX(120),YY(120),Q(45),YZ(20)
280 REM	GREETING
290 PRINT CHR$(&H1A)
300 PRINT "******* POLYNOMIAL  CURVE  FIT *******"
310 PRINT:PRINT
320 PRINT"    THIS PROGRAM PERFORMS A LEAST SQUARED"
330 PRINT"CURVE FIT FOR A GENERAL POLYNOMIAL IN X AND Y"
340 PRINT"IT WILL ALSO EVALUATE AND PLOT THE FUNCTION"
350 PRINT"IF YOU WISH"
360 PRINT:PRINT:PRINT
370 RETURN
380	REM END SETUP MEMORY
390    REM INPUT SUBROUTINE 
400	REM RESET ARRAYS
410 FOR I%=0 TO 20 
420	B(I%)=0
430	Q(I%)=0
440 NEXT I%
450 IF Z% =2 THEN GOTO 570
460 INPUT"HOW MANY DATA PTS. ARE THERE";DATAPTS%
470 PRINT CHR$(&H1A)
480 FOR I%=1 TO 5
490    PRINT
500 NEXT I%
510 PRINT"NOW INPUT THE DATA.":PRINT:PRINT
520 FOR I%=1 TO DATAPTS%
530    PRINT"X(";I%;"), Y(";I%;")";
540	INPUT X(I%),Y(I%)
550 NEXT I%
560 REM
570 INPUT"WHAT ORDER FIT WOULD YOU LIKE";ORDER%
580     N%=ORDER%+1
590 REM 
600	REM VALIDITY CHECK
610 REM
620 IF ORDER%>=DATAPTS% THEN PRINT"NOT ENOUGH DATA FOR THAT ORDER FIT": GOTO 570
630 PRINT CHR$(&H1A)
640 RETURN
650	REM END INPUT
660    REM MATRIX LOAD
670  REM GENERATE VALUES FOR ARRAY
680 FOR I%=0 TO ORDER%*2
690	FOR J%=1 TO DATAPTS%
700	    Q(I%)=X(J%)^I%+Q(I%)
710	NEXT J%
720    Q(I%)=2*Q(I%)
730 NEXT I%
740 REM LOAD
750 FOR I%=0 TO ORDER%
760	FOR J%=0 TO ORDER%
770	    A(I%+1,J%+1)=Q(I%+J%)
780	NEXT J%
790 NEXT I%
800 FOR I%=0 TO ORDER%
810    FOR J%=1 TO DATAPTS%
820	    B(I%+1)=2*X(J%)^I%*Y(J%)+B(I%+1)
830	NEXT J%
840 NEXT I%
850 RETURN
860    REM END MATRIX LOAD
870    REM SIMULTANEOUS EQU SOLN
880 FLAG%=0
890  IF N%<>0 THEN GOTO 930
900  IF A(1,1)=0 THEN FLAG%=1: RETURN
910  X(1)=B(1)/A(1,1)
920  RETURN
930  M%=N%-1
940  FOR I%=1 TO M%
950    H=ABS(A(I%,I%))
960	L%=I%
970	I1%=I%+1
980	FOR J%=I1% TO N%
990	   IF ABS(A(J%,I%)) < H THEN 1020
1000	   H=ABS(A(J%,I%))
1010	   L%=J%
1020    NEXT J%
1030	IF H=0 THEN FLAG%=1: RETURN
1040	IF L%=I% THEN 1130
1050	FOR J%=1 TO N%
1060	   G=A(L%,J%)
1070	   A(L%,J%)=A(I%,J%)
1080	   A(I%,J%)=G
1090	NEXT J%
1100	G=B(L%)
1110	B(L%)=B(I%)
1120	B(I%)=G
1130	FOR J%=I1% TO N%
1140       T=A(J%,I%)/A(I%,I%)
1150	   FOR K%=I1% TO N%
1160		A(J%,K%)=A(J%,K%)-T*A(I%,K%)
1170	   NEXT K%
1180	   B(J%)=B(J%)-T*B(I%)
1190	NEXT J%
1200  NEXT I%
1210 IF A(N%,N%)=0 THEN FLAG%=1: RETURN
1220 C(N%)=B(N%)/A(N%,N%)
1230 I%=N%-1
1240 S=0
1250 I1%=I%+1
1260 FOR J%=I1% TO N%
1270    S=S+A(I%,J%)*C(J%)
1280 NEXT J%
1290 C(I%)=(B(I%)-S)/A(I%,I%)
1300 I%=I%-1
1310 IF I%>0 THEN 1240
1320 RETURN
1330     REM END SIMULTANEOUS EQU
1340     REM SUM OF ERRORS SQRD
1350 FOR I%=1 TO DATAPTS%
1360	 YZ(I%)=0
1370 NEXT I%
1380  E=0
1390  FOR I%=1 TO DATAPTS%
1400     FOR J%=1 TO N%
1410	    YZ(I%)=C(J%)*X(I%)^(J%-1)+YZ(I%)
1420	 NEXT J%
1430	 E=E+(Y(I%)-YZ(I%))^2
1440 NEXT I%
1450 RETURN
1460 REM END SUM ERRORS SQRD
1470     REM OUTPUT STATEMENTS
1480 IF FLAG%=1 THEN PRINT "AMBIGUOUS DATA. NO SOLUTION FOUND":RETURN
1490 PRINT"LIST OF COEFFICIENTS":PRINT
1500 PRINT"POWER OF X","COEFFICIENT":PRINT
1510 FOR I%=N% TO 1 STEP -1
1520	 PRINT "X^";I%-1,C(I%)
1530 NEXT I%
1540 PRINT "SUM OF ERRORS SQRD=";E
1550 PRINT
1560 FOR I%=1 TO DATAPTS%
1570	PRINT"X=";X(I%),"Y=";YZ(I%),"DESIRED Y:";Y(I%)
1580 NEXT I%
1590	REM Printer
1600 PRINT:PRINT
1610 INPUT "Do you want a print out of this (y/n)";PR$
1620 IF PR$<>"y" AND PR$<>"Y" THEN GOTO 1740
1630 LPRINT "LIST OF COEFFICIENTS": LPRINT " "
1640 LPRINT " ": LPRINT "POWER OF X", "COEFFICIENT": LPRINT" "
1650 FOR I%=N% TO 1 STEP -1
1660	LPRINT "X^";I%-1,C(I%)
1670 NEXT I%
1680 LPRINT : LPRINT "SUM OF ERRORS SQRD=";E
1690 LPRINT" "
1700 FOR I%=1 TO DATAPTS%
1710	LPRINT "x=";X(I%),"y=";YZ(I%):LPRINT " ","Desired y:";Y(I%)
1720 NEXT I%
1730 RETURN
1740 REM end output
1750 REM OPTIONS
1760 REM
1770 INPUT"HIT ENTER TO CONTINUE";A$
1780 PRINT CHR$(&H1A): PRINT: PRINT
1790 PRINT"Enter the desired option": PRINT: PRINT
1800 PRINT"	1)EVALUATE THE FUNCTION"
1810 PRINT"	2)TRY A DIFFERENT ORDER FIT"
1820 PRINT"	3)TRY A NEW DATA SET"
1830 PRINT"	4)EXIT THE PROGRAM"
1840 PRINT"        5)PLOT THE FUNCTION VS DATA"
1850 PRINT:INPUT"Enter your choice";Z%
1860 ON Z% GOTO 190,130,130,1870,2220
1870 POKE 3,(PEEK(3) AND 63) OR 128
1871 STOP
1880 REM    END OPTIONS
1890      REM EVALUATE
1900 INPUT"SELECT THE RANGE YOU WANT TO EVALUATE min,max";XINIT,XFIN
1910 INCR=(XFIN-XINIT)/120
1920 ZO=0
1930 FOR I%=0 TO 120
1940     XX(I%)=XINIT+INCR*I%
1950     YY(I%)=0
1960	FOR J%=1 TO N%
1970	    YY(I%)=C(J%)*XX(I%)^(J%-1)+YY(I%)
1980	NEXT J%
1990    IF I%=0 THEN YMIN=YY(0):YMAX=YMIN
2000	IF YY(I%)>YMAX THEN YMAX=YY(I%)
2010	IF YY(I%)<YMIN THEN YMIN=YY(I%)
2020 NEXT I%
2030	REM output evaluate
2040 PRINT CHR$(&H1A)
2050 PRINT"X AND Y VALUES":PRINT
2060 PRINT"X","Y","X","Y"
2070 FOR I%=0 TO 60 STEP 6
2080	PRINT XX(I%),YY(I%),XX(I%+60),YY(I%+60)
2090 NEXT I%
2100 PRINT
2110 PRINT"Y MAX.:";YMAX,"Y MIN.:";YMIN
2120 INPUT "DO YOU WANT THIS PRINTED OUT (y/n)";PR$
2130 IF PR$<>"Y" AND PR$<>"y" THEN RETURN
2140 LPRINT"X and Y values" : LPRINT " "
2150 LPRINT"X","Y","X","Y"
2160 FOR I%=0 TO 60 STEP 6
2170	LPRINT XX(I%),YY(I%),XX(I%+60),YY(I%+60)
2180 NEXT I%
2190 LPRINT" "
2200 LPRINT "Y MIN=";YMIN,"Y MAX=";YMAX
2210 RETURN
2220	REM PLOT
2230 INPUT "PLOTTER READY ";AN$
2240 IF AN$="N" OR AN$="n" THEN RETURN
2250 PRINT CHR$(26);"****** PLOTTING ******"
2260 XMAX=-9E+25
2270 YMAX=-9E+25
2280 XMIN=9E+25
2290 YMIN=9E+25
2300 FOR I%=1 TO DATAPTS%
2310 IF X(I%)>XMAX THEN XMAX=X(I%)
2320 IF Y(I%)>YMAX THEN YMAX=Y(I%)
2330 IF X(I%)<XMIN THEN XMIN=X(I%)
2340 IF Y(I%)<YMIN THEN YMIN=Y(I%)
2350 NEXT I%
2360 LPRINT:LPRINT
2370  LPRINT CHR$(18)
2380 LPRINT "M0,-500"
2390 LPRINT "I"
2400 LPRINT "X0,20,24"
2410 LPRINT "H"
2420 LPRINT "X1,20,24"
2430 LPRINT "H"
2440 REM PLOT DATA POINTS...DETERMINE OFFSET/SCALE
2450 YOFF=YMIN*-1
2460 XOFF=XMIN*-1
2470 YSCALE=480/(YMAX+YOFF)*.99
2480 XSCALE=480/(XMAX+XOFF)*.99
2490 SCALE=YSCALE
2500 IF XSCALE<YSCALE THEN SCALE=XSCALE
2510 LPRINT "C2"
2520 FOR I%=1 TO DATAPTS%
2530 XPT=(X(I%)+XOFF)*XSCALE
2540 YPT=(Y(I%)+YOFF)*YSCALE
2550 LPRINT "M";XPT;",";YPT
2560 LPRINT "P";"*"
2570 NEXT I%
2580 LPRINT "H"
2590 LPRINT "C3"
2600 REM PLOT COMPUTED FUNCTION
2610 FOR I=XMIN TO XMAX STEP (XMAX-XMIN)/120
2620 YPT=0
2630 FOR J%=1 TO N%
2640 YPT=C(J%)*I^(J%-1)+YPT
2650 NEXT J%
2660 YPT=YPT+YOFF
2670 XPT=I+XOFF
2680 IF I=XMIN THEN LPRINT "M";XPT*XSCALE;",";YPT*YSCALE
2690 LPRINT "D";XPT*XSCALE;",";YPT*YSCALE
2700 NEXT I
2710 NTIC=ABS(XMAX)
2720 IF ABS(YMAX)>NTIC THEN NTIC=ABS(YMAX)
2730 LPRINT "C1"
2740 LPRINT "M";480-(LEN(STR$(XMAX))*15);",";10
2750 LPRINT "P";XMAX
2760 LPRINT "M10,480"
2770 LPRINT "P";YMAX
2780 LPRINT "H"
2785 LPRINT "C0"
2790 LPRINT "A"
2800 LPRINT
2810 LPRINT "X TIC INTERVAL=";XMAX/24
2820 LPRINT "Y TIC INTERVAL=";YMAX/24
2830 LPRINT "ORIGIN=";XMIN;",";YMIN
2840 GOTO 1780
 LPRINT "X TIC INTERV