home *** CD-ROM | disk | FTP | other *** search
/ Oakland CPM Archive / oakcpm.iso / cpm / graphics / basplot.lbr / PLOT.BZS / PLOT.BAS
Encoding:
BASIC Source File  |  1991-06-25  |  21.6 KB  |  285 lines
  1. '                            <<<   PLOT.BAS   >>>
  2.  
  3. ' This program was written by J. V. Brancik of Univ. N. S. W. with intention to provide students (or anyone)
  4. ' with a simple program to plot experimantal data and do some basic manipulation (like linear regression analysis).
  5. ' The whole program consists of three sub-programs which are chained (called) together.  The sub-programs are:
  6. ' PLOT.bas, PLOT01.bas, and PLOT02.bas.  The source code is for Microsoft Basic Compiler BASCOM, version 5.3.
  7. ' The spline functions and regression analysis subroutines are based on various Hewlett-Packard publications.
  8.  
  9. ' PLOT.bas:
  10. ' --------
  11. ' Program to enter experimental data and, when required, calculate spline functions and linear regression coefficients,
  12. ' and then to chain complimentary programs PLOT01.BAS and PLOT02.BAS, which will plot the data and required functions
  13. ' on the graphics screen or on the HP or FACIT plotter.  The data file created by this sub-program are as follows:
  14. ' PLOT.$$$ - 
  15. ' is a data file which contains data entered (experimental poins), and which is random access, rec. length 4 bytes
  16. ' (real numbers), max. 40 x and y data pairs.
  17. ' PLOT.FLG - 
  18. ' is a data file which contains information shared between the chaining programs (flags).  It is random access, rec.
  19. ' length 2 byte there are 12 records (integers).
  20. ' PLOT02.DAT - 
  21. ' contains function points for the spline functions, rec. length is 4 bytes (real numbers), there are four sections;
  22. ' each section has 400 recs., first conatins x values for complimentary second and third and fourth sections, which
  23. ' contain y values of f(x), f'(x), and f''(x).
  24. ' PLOT03.DAT - 
  25. ' is similar to the preceeding one, but contains nine sections; first contains x values common to all other sections,
  26. ' second to nineth contain f(x) values of corresponding eight regression models.  File PLOT03.DOC is a sequential file
  27. ' containing data on regression coefficients and regression models.
  28.  
  29. ' Dimensions
  30. DIM A(80),D(8,5),G%(12),S(1),S%(1),X(240),Y(8,400):' exptl data pairs, regression coeffs, flags, min/max, from/to,
  31. ' x-values are stored in a(2*i-1), y-values are stored in a(2*i), where i is the x,y-pair number
  32. ' regresion coefficients are in D(i,j), where
  33. '    1<=i<=8 is the model number, 1<=j<=5 is the coefficient's number, i.e. j=1 for coeff. A, =2 for coeff. B, =3 for
  34. '                                 coeff. C, =4 for coeff. of determination, =5 for the standard error of the estimate
  35. ' g%(i) are the flags
  36. ' s(i) are the x-min and x-max for the regression functions
  37. ' s%(i) are the first and the last data pair to be read
  38. ' x(i) are temporay variable array for data entry and calculation of coefficients
  39. ' y(i,j) are points for regression or spline functions, where i is the model number and 1<=j<=400 are function values
  40.  
  41. ' clear temporary data file (set all entries to 0, file is random access, data stored as REAL, i. e. 4 byte storage)
  42. OPEN"R",#1,"PLOT.$$$",4:FIELD#1,4 AS R$:J=0:Q$=MKS$(J):FOR I%=1 TO 640:LSET R$=Q$:PUT#1,I%:NEXT:CLOSE
  43. ' explanation of the records: A(2*i-1) is X value, A(2*i) is complimentary Y value for 1<=i<=40, where i is counter
  44. ' of X,Y data pairs.  The file contents is as follows
  45. ' record no. | content             | record no. | content            | record no. | content 
  46. '   1 to  80 - set no.1 x,y-pairs     81 to 160 - set no.2 x,y-pairs   161 to 240 - set no.3 x,y-pairs
  47. ' 241 to 320 - set no.4 x,y-pairs    321 to 400 - set no.5 x,y-pairs   401 to 480 - set no.6 x,y-pairs
  48. ' 481 to 560 - set no.5 x,y-pairs    561 to 640 - set no6. x,y-pairs
  49. ' Note that the starting position of a new set is always FIXED (multiples of 80) even when you use less than 40
  50. ' of x,y-pairs for each data set.
  51.  
  52. GOTO 2350:' to where the Main Section of this program is located
  53.  
  54. '                             Subroutine section
  55. '                             ==================
  56.  
  57. '    Subroutines to handle screen editing functions of the Televideo terminal (change to suit your terminal)
  58. ' 1<= M% <=24  is the row number and 1<= N% <=80 is the column number
  59. 1000 PRINT CHR$(27)CHR$(42):RETURN:' clear the whole screen
  60. 1100 PRINT CHR$(27)CHR$(61)CHR$(M%+31)CHR$(N%+31);:RETURN:' cursor addressing (set cursor position)
  61. 1110 GOSUB 1100:PRINT CHR$(27)CHR$(84):GOSUB 1100:RETURN:' clear to the end of line from the current cursor position
  62. 1120 GOSUB 1100:PRINT CHR$(27)CHR$(89):GOSUB 1100:RETURN:' clear to the end of screen from the current cursor position
  63.  
  64. '    Subroutine to open file for random access, FL%=rec length, F$=filename
  65. 1200 OPEN"R",#1,F$,FL%:FIELD#1,FL% AS R$:RETURN
  66.  
  67. '    Subroutine to record flag values (file PLOT.FLG, random access, data stored as INTEGERs, i. e. 2 bytes storage)
  68. 1230 FL%=2:F$="PLOT.FLG":GOSUB 1200:G%(0)=G%(3)+4+(FLG%=1):FOR J%=0 TO G%(3)+4:LSET R$=MKI$(G%(J%)):PUT#1,J%+1:NEXT
  69.      CLOSE:RETURN:' get filename and rec length, set no of recs to read, write all flag that have been used in routines
  70. ' explanation of records: 
  71. ' g(0) - no. of records to read from this file,
  72. ' g(1) - which option to plot:
  73. '        1 = points connected (if required) by a line (8 sets of experimental data allowed)
  74. '        2 = spline functions plot (only one set of experimental data allowed)
  75. '        3 = regression analysis (only one set of experimental data allowed)
  76. '        4 = user's function plot (8 sets of data experimental data and any nnumber of user defined functions points
  77. '            (calculated by the user and stored in a file)
  78. ' g(2) - where to plot:
  79. '        1 = graphics screen (Microangelo or compatible)
  80. '        2 = graphics plotter (Hewlett-Packard or any plotter using HPGL)
  81. '        3 = just experimental data entry for user defined function plot
  82. ' g(3) - number of experimental data sets (max 8)
  83. ' g(4) - number of x,y-pairs for set no.1 (max 40)
  84. ' g(5) .... g(11) - number of x,y-pairs for set no.2 ... 8
  85.  
  86. '    Subroutine to record the experimental data into a temporary file PLOT.$$$ (random access, rec length=4, 640 recs)
  87. 1250 FL%=4:F$="PLOT.$$$":GOSUB 1200:FOR I%=1 TO 2*G%(3+G%(3)):LSET R$=MKS$(X(I%)):PUT#1,I%+(G%(3)-1)*80:NEXT
  88.      CLOSE:RETURN:' get filename, set rec length, open file, write after each set of data entry
  89.  
  90. '    Subroutine to display the Main Menu (a bit of cosmetics but nothing spectracular - not enough memory)
  91. 1300 GOSUB 1000:M%=1:N%=28:GOSUB 1100:PRINT"<<<  Plotting Program  >>>":M%=3:N%=36:GOSUB 1100:PRINT"MAIN MENU":M%=5
  92.      N%=28:GOSUB 1100:PRINT"Select from the following:":M%=6:GOSUB 1100:PRINT STRING$((26),45):M%=8:N%=5:GOSUB 1100
  93.      PRINT"1. Plot experimental points":M%=10:GOSUB 1100:PRINT"2. Plot experimental points and draw a smooth":M%=11
  94.      N%=8:GOSUB 1100:PRINT"curve through the points using spline functions":M%=13:N%=5:GOSUB 1100
  95.      PRINT"3. Plot experimental points, do linear regression":M%=14:N%=8:GOSUB 1100
  96.      PRINT"analysis and draw a regression line or curve":M%=16:N%=5:GOSUB 1100
  97.      PRINT"4. Plot user defined function through":M%=17:N%=8:GOSUB 1100:PRINT"experimental points":M%=19:N%=5
  98.      GOSUB 1100:PRINT"5. Exit this program"
  99. 1310 M%=21:N%=1:GOSUB 1110:INPUT"Type the number of your selection: ",A$:A=VAL(A$):IF A<1 OR A>5 THEN 1310
  100.      IF A=5 THEN END ELSE G%(1)=INT(A):M%=5:N%=1:GOSUB 1120:M%=7:N%=5:GOSUB 1100:PRINT"1. Plot on the Graphics Screen"
  101.      M%=9:GOSUB 1100:PRINT"2. Plot on the plotter"
  102. 1320 M%=21:N%=1:GOSUB 1110:INPUT"Type the number of your selection: ",A$:A=VAL(A$):IF A<1 OR A>2 THEN 1320
  103.      G%(2)=INT(A):RETURN
  104.  
  105. '    Subroutine to display header for the data entry (a bit of cosmetics, too)
  106. '1500 RETURN:' if linker complaines about memory shortage then comment-out the lot below
  107. 1500 M%=1:N%=1:GOSUB 1110:PRINT"<";F1$;">":' clear first line, display header, F1$=filename, F2$=process name
  108.      IF G%(1)=1 THEN A$="Plot of experimental points":N%=INT(40-(LEN(A$))/2):GOSUB 1100:PRINT A$
  109.      IF G%(1)=2 THEN A$="Plot spline functions":N%=INT(40-(LEN(A$))/2):GOSUB 1100:PRINT A$
  110.      IF G%(1)=3 THEN A$="Plot regression function":N%=INT(40-(LEN(A$))/2):GOSUB 1100:PRINT A$
  111.      IF G%(1)=4 THEN A$="Plot user defined function":N%=INT(40-(LEN(A$))/2):GOSUB 1100:PRINT A$
  112.      N%=INT(78-LEN(F2$)):GOSUB 1100:PRINT"<";F2$;">":RETURN
  113.  
  114. '    Subroutine to get input (source) of experimental data
  115. 1510 F2$="data entry":GOSUB 1500:M%=2:N%=1:GOSUB 1120:' display header, FLG%=0 for data entry and =1 for no data entry
  116.      FLG%=0:M%=5:N%=28:GOSUB 1100:PRINT"Select from the following:":M%=6:N%=28:GOSUB 1100:PRINT STRING$((26),45)
  117.      M%=8:N%=8:GOSUB 1100:PRINT"1. Keyboard data entry":M%=9:GOSUB 1100:PRINT"2. Data entry from a data file"
  118.      IF G%(1)=4 THEN M%=11:GOSUB 1100:PRINT"Press ENTER for no data entry":' option for user defined function plot only
  119. 1520 M%=21:N%=1:GOSUB 1110:INPUT"The number of your selection: ",A$:A=VAL(A$):' get the user's response
  120.      IF A=0 AND G%(1)=4 THEN FLG%=1:G%(2)=3:RETURN:' for no data entry (user function plot)
  121.      IF A<1 OR A>2 THEN 1520:' get proper answer, DE%=1 for keyboard, =2 for data file
  122.      DE%=A:IF DE%=1 THEN F1$="PLOT.$$$":GOSUB 1500:' temp data file (to store entered exptl data)
  123.      G%(3)=G%(3)+1:ON DE% GOSUB 1600,1730:RETURN:' to keyboard or data-file routine
  124.  
  125. '    Subroutine to enter data from the keyboard
  126. 1600 M%=2:N%=1:GOSUB 1120:' cls, ask for no. of data pairs, G%(3)=no. of sets, G%(3+i)=no. of data pairs
  127. 1610 M%=21:N%=1:GOSUB 1110:PRINT"How many data pairs for plot no.";G%(3);:INPUT": ",A:IF A<1 OR A>40 THEN 1610
  128.      G%(G%(3)+3)=INT(A):IF G%(3)>1 THEN 1620 ELSE 1640:' if multiplot the x's could be the same
  129. 1620 G%=0:M%=21:N%=1:GOSUB 1110:INPUT"Are the x's from the previous set (Y/N)? ",A$:IF A$="Y" THEN G%=1:GOTO 1640
  130.      IF A$<>"Y" THEN IF A$<>"N" THEN 1620:' g%=+1 for the same x's, -1 for the same y's, 0 for new entry
  131. 1630 GOSUB 1110:INPUT"Are the y's from the previous set (Y/N)? ",A$:IF A$="Y" THEN G%=-1:GOTO 1640
  132.      IF A$<>"Y" THEN IF A$<>"N" THEN 1630
  133. 1640 FOR J%=1 TO G%(G%(3)+3):X%=2*J%-1:Y%=2*J%:IF G%=0 THEN 1650 ELSE 1660:' set the x and y entries counters
  134. 1650 PRINT"X(";J%;:INPUT") =",X(X%):PRINT"Y(";J%;:INPUT") =",X(Y%):GOTO 1700:' when both x and y entered
  135. 1660 IF G%=1 THEN 1670 ELSE 1680:' when x's would be the same, skip the entry but show the previous set x entry
  136. 1670 PRINT"X(";J%;") =";X(X%):PRINT"Y(";J%;:INPUT") =",X(Y%):GOTO 1700:' enter only y's
  137. 1680 IF G%=-1 THEN 1690 ELSE 1700:' when y's would be the same, skip the entry but show the previous set y entry
  138. 1690 PRINT"X(";J%;:INPUT") =",X(X%):PRINT"Y(";J%;") =";X(Y%):' enter only x's
  139. 1700 NEXT:FOR J%=1 TO 4:PRINT:NEXT:' make some room on the screen
  140. 1710 F1$="PLOT.$$$":F2$="keyboard entry":GOSUB 1500:' print header again
  141. 1720 M%=21:N%=1:GOSUB 1110:INPUT"Review and correct the data (Y/N)? ",A$:IF A$<>"Y" THEN IF A$<>"N" THEN 1720
  142.      IF A$="Y" THEN GOSUB 1780:GOTO 1710:' back to review and correct the data
  143.      GOSUB 1250:RETURN:' write data to the temporary file
  144.  
  145. '    Subroutine to enter data from user's data file
  146. 1730 M%=2:N%=1:GOSUB 1120:' clear screen, ask for the name of the data file
  147. 1740 M%=21:N%=1:GOSUB 1110:LINE INPUT"The name of the INPUT data file: ",A$:IF A$="" THEN 1740:' is there an entry?
  148. 1750 M%=21:N%=1:GOSUB 1120:INPUT"Start with the pair no.: ",S%(0):M%=22:N%=1:GOSUB 1110
  149.      INPUT"Finish with the pair no.: ",S%(1):S%=S%(1)-S%(0):G%(G%(3)+3)=S%+1:c%=1:IF S%<1 OR S%>40 THEN 1750
  150.      F$=A$:FL%=4:GOSUB 1200:FOR J%=S%(0) TO S%(1):GET#1,(2*j%-1):q$=r$:X(2*c%-1)=CVS(q$):' read rec for x,
  151.      GET#1,(2*j%):q$=r$:X(2*c%)=CVS(q$):c%=C%+1:NEXT:CLOSE:' increase counter, read rec. for y, go for next pair or close
  152. 1760 F1$=F$:F2$="data from file":GOSUB 1500:' display header
  153. 1770 M%=21:N%=1:GOSUB 1120:INPUT"Review and correct the data (Y/N)? ",A$:IF A$<>"Y" THEN IF A$<>"N"THEN 1770
  154.      IF A$="Y" THEN GOSUB 1780:FOR J%=1 TO 4:PRINT:NEXT:GOTO 1760:' view the data before recording into DF
  155.      GOSUB 1250:RETURN:' record the data into the program DF
  156.  
  157. '    Subroutine to correct the data entered
  158. 1780 F1$="PLOT.$$$":F2$="data review":GOSUB 1500:M%=2:N%=1:GOSUB 1120:PRINT"Type the correct data or press RETURN!"
  159.      FOR J%=1 TO G%(G%(3)+3):X%=2*J%-1:Y%=2*J%:PRINT"X(";G%(3);",";J%;") =";X(X%);:INPUT A$
  160.      IF A$="" THEN 1790 ELSE X(X%)=VAL(A$):' change x if wrong
  161. 1790 PRINT"Y(";G%(3);",";J%;") =";X(Y%);:INPUT A$
  162.      IF A$="" THEN 1800 ELSE X(Y%)=VAL(A$):' change y if wrong
  163. 1800 NEXT:FOR J%=1 TO 4:PRINT:NEXT:RETURN
  164.  
  165. '    Subroutine to calculate spline function coefficients
  166. ' Deatiled explanation of this sub is in following references:
  167. ' (usually found in undergraduate collections of Univ. libraries)
  168. ' A. Ralston and H. S. Wilf: Mathematical Methods for Digital Computers, John Wiley & Sons eds., N. Y. (1967)
  169. ' T. N. E. Greville; Theory and Applications of Spline Functions, Acad. Press publ., N. Y. (1969)
  170. ' NOTE should be taken here that the first and the second derivatives refer to the spline function; there are other
  171. '      methods how to draw a smooth line through a set of data which may differ from the spline curve.  Thus, for the
  172. '      data points fitting the equation y=x^2, when using spline function the first derivative is not necessarily
  173. '      y'=2x and the second one is not necessarily y''=2.  The user can compare the curves drawn by the regression
  174. '      and spline functions.
  175. 2100 X(0)=1E-3:A%=G%(4):X(4*A%+1)=0:X(5*A%)=0:B%=0:C=4*(2-SQR(3)):' set tolerance, no. of data pairs, counters
  176. 2110 B%=B%+1:Z%=3*A%+B%:X(Z%-A%)=X(2*B%+1)-X(2*B%-1):X(Z%)=(X(2*B%+2)-X(2*B%))/X(Z%-A%):IF B%>1 THEN 2120 ELSE 2110
  177. 2120 X(Z%+A%)=2*(X(Z%)-X(Z%-1))/(X(Z%-A%-1)+X(Z%-A%)):X(2*A%+Z%)=3*(X(A%+Z%))/2:IF (A%-1)>B% THEN 2110
  178. 2130 X=0:B%=2:' reset the counters (of points)
  179. 2140 Z%=4*A%+B%:Y=X(Z%+1):Y=C*(((Y-X(Z%-1))/(1+X(Z%-2*A%)/X(Z%-2*A%-1))-Y)/2-X(Z%)+X(Z%+A%)):IF ABS(Y)>X THEN X=ABS(Y)
  180. 2150 X(Z%)=X(Z%)+Y:B%=B%+1:IF B%<>A% THEN 2140 ELSE IF X(0)<=X THEN 2130:' count the area for the integral
  181. 2160 C=0:B%=1:' reset the counters for the calcn. of coeffs.
  182. 2170 Z%=4*A%+B%:Y=X(Z%-2*A%):C=Y*(X(2*B%)+X(2*B%+2))/2+(Y^3)*(X(Z%)+X(Z%+1))/24+C:B%=B%+1:IF A%=B% THEN 2180 ELSE 2170
  183. 2180 B%=1:' reset the counter again for other coeffs.
  184. 2190 Z%=4*A%+B%:X(Z%-2*A%)=(X(Z%+1)-X(Z%))/X(Z%-2*A%):Z%=Z%-2*A%:B%=B%+1:IF B%=A% THEN 2200 ELSE 2190
  185. 2200 B%=1:X=X(1):X(0)=(X(2*A%-1)-X(1))/400:H=C:FOR J%=1 TO 400:GOSUB 2240:B%=1:X=X+X(0):NEXT:RETURN:'for points to plot
  186. '    Subroutine to calculate spline functions points for plotting the spline curve
  187. ' this sub is called by sub 2100 for each one of 400 points - see line 2200
  188. 2240 IF B%>A% THEN 2290 ELSE IF A%<>B% THEN 2260 ELSE 2250
  189. 2250 Z=X-X(2*A%-1):IF A%=B% THEN 2290 ELSE IF Z<>0 THEN 2290 ELSE IF Z=0 THEN 2280
  190. 2260 Y=X-X(2*B%-1):Z=X-X(2*B%+1):IF Y<0 OR Z>=0 THEN 2270 ELSE 2280
  191. 2270 B%=B%+1:GOTO 2240
  192. 2280 C=X(4*A%+B%)+Y*X(2*A%+B%):Y(3,J%)=C:C=(X(4*A%+B%)+X(4*A%+B%+1)+C)/6:Y=X(2*B%)+Y*X(3*A%+B%)+C*Y*Z
  193.      Q=X(3*A%+B%)+(X-X(2*B%-1)+Z)*C+Z*(X-X(2*B%-1))*X(2*A%+B%)/6:Y(2,J%)=Q:Y(1,J%)=Y:Y(0,J%)=X
  194. 2290 RETURN
  195.  
  196. '    Subroutine to record all the calculated data; x, f(x), f'(x), f''(x), PLOT02.DAT, direct access, REAL numbers
  197. 2300 F$="PLOT02.DAT":FL%=4:GOSUB 1200:FOR I%=1 TO 400:LSET R$=MKS$(Y(0,I%)):PUT#1,I%:NEXT:L%=401:' x's are common
  198.      FOR I%=1 TO 3:FOR K%=1 TO 400:LSET R$=MKS$(Y(I%,K%)):PUT#1,L%:L%=L%+1:NEXT:NEXT:close:RETURN:' rec. set by set
  199.  
  200. '    Subroutine to read data from the experimental data file (temporary file)
  201. 2310 F$="PLOT.$$$":FL%=4:GOSUB 1200:C%=0:FOR I%=1 TO 2*G%(4):GET#1,I%:X(I%)=CVS(R$):IF X(I%)<=1E-37 THEN C%=1
  202.      NEXT:CLOSE:RETURN:' read the data, check for negative numbers or zeros
  203.  
  204. '    Subroutine to calculate all sums for all regression functions (takes a while for 40 x,y-pairs)
  205. ' This sub is based upon the following reference:
  206. ' N. R. Draper and H. Smith: Applied Regression Analysis, John Wiley & Sons publ., N.Y., 1st ed. (1966), 2nd ed. (1981)
  207. ' NOTE - any other book on the regression analysis would suit as well
  208. 2320 A=G%(4):FOR I%=1 TO G%(4):X%=2*I%-1:Y%=2*I%:A(14)=A(14)+X(X%)*X(Y%):A(15)=A(15)+X(Y%)/X(X%)
  209.      A(16)=A(16)+1/(X(X%)*X(Y%)):A(17)=A(17)+X(Y%)*SQR(X(X%)):A(18)=A(18)+X(X%)*LOG(X(Y%))
  210.      A(19)=A(19)+LOG(X(X%))*LOG(X(Y%)):A(20)=A(20)+X(Y%)*LOG(X(X%)):A(21)=A(21)+X(X%):A(22)=A(22)+1/X(X%)
  211.      A(23)=A(23)+SQR(X(X%)):A(24)=A(24)+LOG(X(X%)):A(25)=A(25)+X(X%)^2:A(26)=A(26)+1/(X(X%)^2)
  212.      A(27)=A(27)+(LOG(X(X%)))^2:A(28)=A(28)+X(Y%):A(29)=A(29)+1/X(Y%):A(30)=A(30)+LOG(X(Y%)):A(31)=A(31)+1/(X(Y%)^2)
  213.      A(32)=A(32)+X(Y%)^2:A(33)=A(33)+X(X%)^3:A(34)=A(34)+(LOG(X(Y%)))^2:A(35)=A(35)+X(X%)^4:A(36)=A(36)+X(Y%)*X(X%)^2
  214.      NEXT:FOR I%=1 TO 7:A(10)=A(13+I%):B=21-(I%=2)-(I%=3)-2*(I%=4)-3*(I%=6)-3*(I%=7):A(11)=A(B)
  215.      B=25-(I%=2)-(I%=3)-2*(I%=6)-2*(I%=7)+4*(I%=4):A(12)=A(B):B=28-(I%=3)-2*(I%=5)-2*(I%=6):A(13)=A(B)
  216.      B=32+(I%=3)-2*(I%=5)-2*(I%=6):A(7)=A(B):B=A*A(10)-A(11)*A(13):D(I%,2)=B/(A*A(12)-A(11)^2)
  217.      D(I%,1)=(A(13)-D(I%,2)*A(11))/A:D(I%,5)=SQR(ABS((A(12)-D(I%,1)*A(7)-D(I%,2)*A(10))/A))
  218.      D(I%,4)=(D(I%,1)*A(13)+D(I%,2)*A(10)-(A(13)^2)/A)/(A(7)-(A(13)^2)/A):IF I%=5 OR I%=6 THEN D(I%,1)=EXP(D(I%,1))
  219.      NEXT:A(37)=A*A(25)-A(21)^2:A(38)=A*A(36)-A(25)*A(28):A(39)=A*A(33)-A(21)*A(25):A(40)=A*A(14)-A(21)*A(28)
  220.      A(41)=A*A(35)-A(25)^2:A(42)=A(37)*A(38)-A(39)*A(40):D(8,3)=A(42)/(A(37)*A(41)-A(39)^2)
  221.      D(8,2)=(A(40)-D(8,3)*A(39))/A(37):D(8,1)=(A(28)-D(8,3)*A(25)-D(8,2)*A(21))/A
  222.      D(8,4)=(D(8,1)*A(28)+D(8,2)*A(14)+D(8,3)*A(36)-A(28)^2/A)/(A(32)-A(28)^2/A):' d(i,1)=A, d(i,2)=B, d(i,3)=C
  223.      D(8,5)=SQR(ABS((A(32)-D(8,1)*A(28)-D(8,2)*A(14)-D(8,3)*A(36))/(A-3))):' d(i,4)=dtmn coeff, d(i,5)=std err
  224.      RETURN
  225.  
  226. '    Subroutine to calculate regression function points to be plotted (using coeffs already calculated)
  227. 2330 X(0)=(S(1)-S(0))/400:X=S(0):FOR I%=1 TO 400:X=X+X(0):Y(0,I%)=X:Y(1,I%)=D(1,1)+D(1,2)*X:Y(5,I%)=D(5,1)*EXP(D(5,2)*X)
  228.      IF X>1E-37 THEN Y(2,I%)=D(2,1)+D(2,2)/X:' x must be > 10^-38
  229.      IF X>1E-37 THEN Y(3,I%)=1/(D(3,1)+D(3,2)/X):' x (and y) must be > 10^-38
  230.      IF X>1E-37 THEN Y(4,I%)=D(4,1)+D(4,2)*SQR(X):' x can be >=0
  231.      IF X>1E-37 THEN Y(7,I%)=D(7,1)+D(7,2)*LOG(X):' x must be > 10^-38
  232.      Y(6,I%)=D(6,1)*X^D(6,2):Y(8,I%)=D(8,1)+D(8,2)*X+D(8,3)*X*X
  233.      NEXT:RETURN
  234.  
  235. '    Subroutine to record all functions points and all coefficients
  236. '    Two data files are used here:
  237. '    PLOT03.DAT - direct access, binary, REAL numbers, 1600 records, 4 bytes each
  238. '    PLOT03.DOC - sequential, ascii, variable length, delimited by CR/LF sequence
  239. 2340 F$="PLOT03.DAT":FL%=4:GOSUB 1200:L%=1:FOR I%=0 TO 8:FOR J%=1 TO 400:LSET R$=MKS$(Y(I%,J%)):PUT#1,L%:L%=L%+1
  240.      NEXT:NEXT:FOR I%=1 TO 8:FOR J%=1 TO 5:LSET R$=MKS$(D(I%,J%)):PUT#1,L%:L%=L%+1:NEXT:NEXT:CLOSE:' write the coeffs
  241.      OPEN"O",#1,"PLOT03.DOC":' sequential file to record regression coeffs, coeff of detdn, std err of the estimate
  242.      A$="Y ="+STR$(D(1,1))+" + "+STR$(D(1,2))+" * X; r^2 ="+STR$(D(1,4))+", S ="+STR$(D(1,5)):WRITE#1,A$
  243.      A$="Y = "+STR$(D(2,1))+" + "+STR$(D(2,2))+" * (1/X); r^2 ="+STR$(D(2,4))+", S ="+STR$(D(2,5)):WRITE#1,A$
  244.      A$="(1/Y) = "+STR$(D(3,1))+" + "+STR$(D(3,2))+" * (1/X); r^2 ="+STR$(D(3,4))+", S ="+STR$(D(3,5)):WRITE#1,A$
  245.      A$="Y ="+STR$(D(4,1))+" + "+STR$(D(4,2))+" * X^(1/2); r^2 ="+STR$(D(4,4))+", S ="+STR$(D(4,5)):WRITE#1,A$
  246.      A$="Y ="+STR$(D(5,1))+" + EXP("+STR$(D(5,2))+" * X); r^2 ="+STR$(D(5,4))+", S ="+STR$(D(5,5)):WRITE#1,A$
  247.      A$="Y ="+STR$(D(6,1))+" * X^("+STR$(D(6,2))+"); r^2 ="+STR$(D(6,4))+", S ="+STR$(D(6,5)):WRITE#1,A$
  248.      A$="Y ="+STR$(D(7,1))+" + "+STR$(D(7,2))+" * LN(X); r^2 ="+STR$(D(7,4))+", S ="+STR$(D(7,5)):WRITE#1,A$
  249.      A$="Y ="+STR$(D(8,1))+" + "+STR$(D(8,2))+" * X +"+STR$(D(8,3))+" * X^2; r^2="+STR$(D(8,4))+", S ="+STR$(D(8,5))
  250.      WRITE#1,A$:CLOSE:RETURN
  251.  
  252.  
  253. '                      The Main Section of the program
  254. '                      ===============================
  255.  
  256. 2350 GOSUB 1300:' display the Main menu and get the options
  257. 2360 GOSUB 1510:IF G%(1)=1 THEN 2370:' max 8 sets of data for experimental points plot only
  258.      IF G%(1)=2 OR G%(1)=3 THEN 2380:' spline functions or regression analysis only one set of experimental data
  259.      IF G%(1)=4 THEN IF G%(2)<3 THEN 2370 ELSE 2380:' if user's function (up to 8 sets of experimental data allowed)
  260. 2370 IF G%(3)>7 THEN 2380 ELSE M%=21:N%=1:GOSUB 1120:PRINT"Another set (no.";G%(3)+1;:INPUT") to enter (Y/N)? ",A$
  261.      IF A$="Y" THEN 2360:' back for more data
  262.      IF A$="N" THEN 2380:' to record the data and then exit
  263.      GOTO 2370:' note that answers should be in CAPs
  264. 2380 GOSUB 1250:GOSUB 1230:' data entry finishes here, write flags and exptl data into data files
  265.      IF G%(1)=1 THEN CHAIN"PLOT01.COM":' plot points
  266.      IF G%(1)=2 THEN GOSUB 2500:CHAIN"PLOT01.COM":' do spline functions routine
  267.      IF G%(1)=3 THEN GOSUB 2600:CHAIN"PLOT01.COM":' do family regression routine
  268.      IF G%(1)=4 THEN CHAIN"PLOT01.COM":' plot user function
  269.  
  270. '    spline routine calls
  271. 2500 F2$="calculating":GOSUB 1500:M%=2:N%=1:GOSUB 1120:GOSUB 2310:GOSUB 2100:GOSUB 2300:RETURN
  272.  
  273. '    regression routine calls
  274. 2600 GOSUB 1000:FL%=4:F$="PLOT.$$$":F2$="regression models":GOSUB 2310:IF C%=1 THEN 2610 ELSE 2620:' cls, rd exptl data
  275. '    data should be larger then 10^-38 or else!
  276. 2610 M%=7:N%=1:GOSUB 1100:PRINT"Cannot proceed with calculations because your experimental data contain zero"
  277.      M%=8:GOSUB 1100:PRINT"or negative number. The program will return to the main menu, where you can"
  278.      M%=9:GOSUB 1100:PRINT"do corrections.":M%=21:GOSUB 1120:INPUT"Press ENTER to continue!",A$
  279. 2620 M%=1:N%=1:GOSUB 1500:M%=21:N%=1:GOSUB 1120:INPUT"X(MIN) = ",A$:S(0)=VAL(A$):IF S(0)<=1E-37 THEN 2620
  280.      M%=22:N%=1:GOSUB 1110:INPUT"X(MAX) = ",A$:S(1)=VAL(A$):IF S(1)-S(0)<=0 THEN 2620
  281.      M%=5:N%=1:GOSUB 1120:PRINT"Calculating all the sums ...":GOSUB 2320:M%=5:N%=1:GOSUB 1120
  282.      PRINT"Calculating regression coefficients ...":GOSUB 2330:M%=5:N%=1:GOSUB 1120:PRINT"Recording calculated data ..."
  283.      GOSUB 2340:RETURN:' write data (points and coeffs)
  284. END
  285.