home *** CD-ROM | disk | FTP | other *** search
/ CP/M / CPM_CDROM.iso / mbug / mbug113.arc / CURVEFIT.BAS < prev    next >
BASIC Source File  |  1979-12-31  |  6KB  |  149 lines
  1. 100   REM  GENERAL CURVE FITTING PROGRAM
  2. 110   REM  CURVEFIT.BAS
  3. 120   REM  Written in Microsoft BASIC (Ver. 5.21)
  4. 130   REM
  5. 140   REM  Andrew Waltho, 1986
  6. 500   REM  ****  Program Data
  7. 510   DATA "Linear","Y=A+B*X","X=(Y-A)/B"
  8. 520   DATA "Exponential","Y=A*EXP(B*X)","X=LOG(Y/A)/B"
  9. 530   DATA "Logarithmic","Y=A+B*LOG(X)","X=EXP((Y-A)/B)"
  10. 540   DATA "Power","Y=A*X^B","X=(Y/A)^(1/B)"
  11. 1000  REM  ****  Define variable types, dimension arrays
  12. 1010  OPTION BASE 1:REM  Arrays do not have zeroth elements
  13. 1020  DEFINT F,I-N:REM  All other variables single precision real or strings
  14. 1030  DIM POINTS(2,300), EQUATION$(4,3)
  15. 1040  FLAG1=1:REM  Set initial value of flag for main calculation loop
  16. 1500  REM  ****  Main Program
  17. 1510  GOSUB 10000:REM                  Print opening message
  18. 1520  WHILE FLAG1
  19. 1540    GOSUB 11000:REM                Enter standard data for calculation
  20. 1550    GOSUB 12000:REM                Clear variables used in curve fitting equation
  21. 1560    GOSUB 12500:REM                Read data for curve fitting equations
  22. 1570    GOSUB 13000:REM                Determine curve fitting relationships
  23. 1580    GOSUB 13500:REM                Determine most suitable curve type
  24. 1610    PRINT "Another data set?  Y/N"
  25. 1620    GOSUB 17500:REM                Yes/No selection
  26. 1630    IF KEY=78 THEN FLAG1=0
  27. 1640    PRINT CHR$(26)
  28. 1650  WEND
  29. 1660  SYSTEM
  30. 1670  END
  31. 10000 REM  ****  Print opening message
  32. 10010 PRINT CHR$(26);TAB(27);"***************************"
  33. 10020 PRINT TAB(27);"*                         *"
  34. 10030 PRINT TAB(27);"*  CURVE FITTING PROGRAM  *"
  35. 10040 PRINT TAB(27);"*                         *"
  36. 10050 PRINT TAB(27);"***************************":PRINT:PRINT
  37. 10060 RETURN
  38. 11000 REM  ****  Standard data entry subroutine
  39. 11010 PRINT "Enter standard data":PRINT
  40. 11020 FLAG2=1
  41. 11030 WHILE FLAG2
  42. 11040   INPUT "Number of points (300 max.): ";NPOINTS
  43. 11050   IF (NPOINTS<=1) OR (NPOINTS>300) THEN PRINT "Data for between 1 and 300 points required" ELSE FLAG2=0:REM  Exit loop if data for more than one standard entered
  44. 11060 WEND
  45. 11080 FOR K=1 TO NPOINTS
  46. 11090   FLAG2=1
  47. 11100   WHILE FLAG2
  48. 11110     PRINT CHR$(26);"Point ";K
  49. 11120     INPUT "X - value ";POINTS(1,K)
  50. 11130     INPUT "Y - value ";POINTS(2,K):PRINT
  51. 11140     PRINT CHR$(26);"X - value: ";POINTS(1,K);TAB(25);"Y - value: "POINTS(2,K);TAB(50);"Data correct? Y/N":PRINT
  52. 11150     KEY$=INKEY$:IF KEY$="" THEN 11150
  53. 11160     KEY=ASC(KEY$):KEY=(KEY AND 95):IF (KEY<>89) AND (KEY<>78) THEN 11150
  54. 11170     IF KEY=89 THEN FLAG2=0:REM  Exit loop if data correct
  55. 11180   WEND
  56. 11190   IF K=1 THEN XMIN=POINTS(1,K):YMIN=POINTS(2,K):GOTO 11220
  57. 11200   IF POINTS(1,K)<XMIN THEN XMIN=POINTS(1,K)
  58. 11210   IF POINTS(2,K)<YMIN THEN YMIN=POINTS(2,K)
  59. 11220 NEXT K
  60. 11225 PRINT CHR$(26):REM  Clear screen to indicate calculations in progress
  61. 11230 RETURN
  62. 12000 REM  ****  Clear variables used in curve fitting equations
  63. 12010 XTOTAL=0:YTOTAL=0
  64. 12020 XYTOTAL=0:XXTOTAL=0:YYTOTAL=0
  65. 12030 TXLOG=0:TYLOG=0
  66. 12040 SX=0:SY=0
  67. 12050 X2=0:Y2=0
  68. 12060 SLOG=0
  69. 12500 REM  ****  Read data for curve fitting routine
  70. 12510 RESTORE 500
  71. 12520 FOR K=1 TO 4
  72. 12530   FOR L=1 TO 3
  73. 12540     READ EQUATION$(K,L)
  74. 12550   NEXT L
  75. 12560 NEXT K
  76. 12570 RETURN
  77. 13000 REM  ****  Determine curve fitting relationships
  78. 13010 FOR K=1 TO NPOINTS
  79. 13020   X=POINTS(1,K):Y=POINTS(2,K)
  80. 13030   IF X>0 THEN XLOG=LOG(X) ELSE XLOG=0
  81. 13040   IF Y>0 THEN YLOG=LOG(Y) ELSE YLOG=0
  82. 13050   XTOTAL=XTOTAL+X
  83. 13060   YTOTAL=YTOTAL+Y
  84. 13070   XYTOTAL=XYTOTAL+X*Y
  85. 13080   XXTOTAL=XXTOTAL+X*X
  86. 13090   YYTOTAL=YYTOTAL+Y*Y
  87. 13100   TXLOG=TXLOG+XLOG
  88. 13110   TYLOG=TYLOG+YLOG
  89. 13120   SX=SX+X*YLOG
  90. 13130   SY=SY+Y*XLOG
  91. 13140   X2=X2+XLOG*XLOG
  92. 13150   Y2=Y2+YLOG*YLOG
  93. 13160   SLOG=SLOG+XLOG*YLOG
  94. 13170 NEXT K
  95. 13180 S1=XTOTAL*YTOTAL
  96. 13190 S2=XTOTAL*XTOTAL
  97. 13200 S3=YTOTAL*YTOTAL
  98. 13210 S4=TXLOG*TXLOG
  99. 13220 S5=TYLOG*TYLOG
  100. 13230 RETURN
  101. 13500 REM  ****  Determine most suitable curve type
  102. 13510 RBEST=0
  103. 13520 PRINT CHR$(26);"Regression Coefficients":PRINT
  104. 13530 FOR K=1 TO 4
  105. 13540   ON K GOSUB 14000,14500,15000,15500:REM  Determine curve equations
  106. 13550   R$=STR$(R2):IF R2=0 THEN R$="NOT APPLICABLE"
  107. 13560   PRINT EQUATION$(K,1);TAB(15);R$
  108. 13565   LPRINT EQUATION$(K,1);TAB(15);R$:REM  List results on printer
  109. 13570   IF R2>RBEST THEN RBEST=R2:ABEST=A:BBEST=B:KBEST=K         
  110. 13580 NEXT K
  111. 13590 R2=RBEST:A=ABEST:B=BBEST
  112. 13600 PRINT:PRINT EQUATION$(KBEST,1);" curve fit is best."
  113. 13605 LPRINT:LPRINT EQUATION$(KBEST,1);" curve fit is best.":LPRINT:LPRINT:LPRINT
  114. 13610 PRINT "The formulae are:"
  115. 13620 PRINT TAB(10);"1)  ";EQUATION$(KBEST,2)
  116. 13630 PRINT TAB(10);"2)  ";EQUATION$(KBEST,3)
  117. 13640 PRINT "The values of A and B are:"
  118. 13650 PRINT TAB(10);"A=";A
  119. 13660 PRINT TAB(10);"B=";B
  120. 13740 RETURN
  121. 14000 REM  ****  Determine linear equation
  122. 14010 B=(XYTOTAL-S1/NPOINTS)/(XXTOTAL-S2/NPOINTS)
  123. 14020 A=YTOTAL/NPOINTS-B*XTOTAL/NPOINTS
  124. 14030 R2=(XYTOTAL-S1/NPOINTS)^2/((XXTOTAL-S2/NPOINTS)*(YYTOTAL-S3/NPOINTS))
  125. 14040 RETURN
  126. 14500 REM  ****  Determine expotential equation
  127. 14510 IF YMIN<=0 THEN R2=0:RETURN:REM  No solution if YMIN is zero or negative
  128. 14520 B=(SX-XTOTAL*TYLOG/NPOINTS)/(XXTOTAL-S2/NPOINTS)
  129. 14530 A=EXP((TYLOG-B*XTOTAL)/NPOINTS)
  130. 14540 R2=(SX-XTOTAL*TYLOG/NPOINTS)^2/((XXTOTAL-S2/NPOINTS)*(Y2-S5/NPOINTS))
  131. 14550 RETURN
  132. 15000 REM  ****  Determine logarithmic equation
  133. 15010 IF XMIN<=0 THEN R2=0:RETURN:REM  No solution if XMIN is zero or negative
  134. 15020 B=(SY-TXLOG*YTOTAL/NPOINTS)/(X2-S4/NPOINTS)
  135. 15030 A=(YTOTAL-B*TXLOG)/NPOINTS
  136. 15040 R2=(SY-TXLOG*YTOTAL/NPOINTS)^2/((X2-S4/NPOINTS)*(YYTOTAL-S3/NPOINTS))
  137. 15050 RETURN
  138. 15500 REM  ****  Determine power equation
  139. 15510 IF (XMIN<=0) OR (YMIN<=0) THEN R2=0:RETURN:REM  No solution if either XMIN or YMIN is zero or negative
  140. 15520 B=(SLOG-TXLOG*TYLOG/NPOINTS)/(X2-S4/NPOINTS)
  141. 15530 A=EXP(TYLOG/NPOINTS-B*TXLOG/NPOINTS)
  142. 15540 R2=(SLOG-TXLOG*TYLOG/NPOINTS)^2/((X2-S4/NPOINTS)*(Y2-S5/NPOINTS))
  143. 15550 RETURN
  144. 17500 REM  ****  Y/N Key entry subroutine
  145. 17510 KEY$=INKEY$:IF KEY$<>"" THEN KEY=ASC(KEY$) ELSE 17510
  146. 17520 KEY=(KEY AND 95):REM  Convert any lowercase to uppercase
  147. 17530 IF (KEY<>78) AND (KEY<>89) THEN 17510
  148. 17540 RETURN
  149.