home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l042 / 1.ddi / CHAP3.ARC / CUBE_CLA.PAS next >
Encoding:
Pascal/Delphi Source File  |  1987-12-30  |  11.5 KB  |  323 lines
  1. program CubicSplineClamped_Prog;
  2.  
  3. {------------------------------------------------------------------------}
  4. {-                                                                      -}
  5. {-     Turbo Pascal Numerical Methods Toolbox                           -}
  6. {-     Copyright (c) 1986, 87 by Borland International, Inc.            -}
  7. {-                                                                      -}
  8. {-       Purpose: This program demonstrates interpolation with          -}
  9. {-                a clamped cubic spline.                               -}
  10. {-                                                                      -}
  11. {-       Unit   : Interp    procedure CubicSplineClamped                -}
  12. {-                                                                      -}
  13. {------------------------------------------------------------------------}
  14.  
  15. {$I-}           { Disable I/O error trapping }
  16. {$R+}           { Enable range checking }
  17.  
  18. uses
  19.   Interp, Dos, Crt, Common;
  20.  
  21. var
  22.   XData, YData : TNvector;                  { Data points (X,Y) }
  23.   NumPoints : integer;                      { Number of points }
  24.   DerivLE, DerivRE : Float;                 { Derivative at the endpoints }
  25.   Coef0, Coef1, Coef2, Coef3 : TNvector;    { Coefficients of the spline }
  26.   NumInter : integer;                       { Number interpolating points }
  27.   XInter, YInter : TNvector;                { Interpolating points }
  28.   Error : byte;                             { Flags an error }
  29.  
  30. procedure Initialize(var Coef0     : TNvector;
  31.                      var Coef1     : TNvector;
  32.                      var Coef2     : TNvector;
  33.                      var Coef3     : TNvector;
  34.                      var XData     : TNvector;
  35.                      var YData     : TNvector;
  36.                      var DerivLE   : Float;
  37.                      var DerivRE   : Float;
  38.                      var NumPoints : integer;
  39.                      var NumInter  : integer;
  40.                      var XInter    : TNvector;
  41.                      var YInter    : TNvector;
  42.                      var Error     : byte);
  43.  
  44. {----------------------------------------------------------------}
  45. {- Output: Coef0, Coef1, Coef2, Coef3, XData, YData, DerivLE,   -}
  46. {-         DerivRe, NumPoints, NumInter, XInter, YInter, Error  -}
  47. {-                                                              -}
  48. {- This procedure initializes the above variables to zero.      -}
  49. {----------------------------------------------------------------}
  50.  
  51. begin
  52.   FillChar(Coef0, SizeOf(Coef0), 0);
  53.   FillChar(Coef1, SizeOf(Coef1), 0);
  54.   FillChar(Coef2, SizeOf(Coef2), 0);
  55.   FillChar(Coef3, SizeOf(Coef3), 0);
  56.   FillChar(XData, SizeOf(XData), 0);
  57.   FillChar(YData, SizeOf(YData), 0);
  58.   FillChar(XInter, SizeOf(XInter), 0);
  59.   FillChar(YInter, SizeOf(YInter), 0);
  60.   NumPoints := 0;
  61.   NumInter := 0;
  62.   DerivLE := 0;
  63.   DerivRE := 0;
  64.   Error := 0;
  65.   Writeln;
  66. end; { procedure Initialize }
  67.  
  68. procedure GetData(var NumPoints : integer;
  69.                   var XData     : TNvector;
  70.                   var YData     : TNvector;
  71.                   var DerivLE   : Float;
  72.                   var DerivRE   : Float;
  73.                   var NumInter  : integer;
  74.                   var XInter    : TNvector);
  75.  
  76. {----------------------------------------------------------------}
  77. {- Output: NumPoints, XData, YData, DerivLE, DerivRE,           -}
  78. {-         NumInter, XInter                                     -}
  79. {-                                                              -}
  80. {- This procedure reads in data from either the keyboard        -}
  81. {- or a data file.  The number of data points (NumPoints),      -}
  82. {- the data points (XData, YData), the second derivatives at    -}
  83. {- the endpoints (DerivLE, DerivRE), the number of interpolated -}
  84. {- points (NumInter) and the X values at which to interpolate   -}
  85. {- (XInter) are all read in here.                               -}
  86. {----------------------------------------------------------------}
  87.  
  88. var
  89.   Ch : char;
  90.  
  91. procedure GetTwoVectorsFromFile(var NumPoints : integer;
  92.                                 var XData     : TNvector;
  93.                                 var YData     : TNvector);
  94.  
  95. {-------------------------------------------------------------}
  96. {- Output: NumPoints, XData, YData                           -}
  97. {-                                                           -}
  98. {- This procedure reads in the data points from a data file. -}
  99. {-------------------------------------------------------------}
  100.  
  101. var
  102.   Filename : string[255];
  103.   InFile : text;
  104. begin
  105.   Writeln;
  106.   repeat
  107.     Write('File name? ');
  108.     Readln(Filename);
  109.     Assign(InFile, Filename);
  110.     Reset(InFile);
  111.     IOCheck;
  112.   until not IOerr;
  113.   NumPoints := 0;
  114.   while not EOF(InFile) do
  115.   begin
  116.     NumPoints := Succ(NumPoints);
  117.     Readln(InFile, XData[NumPoints], YData[NumPoints]);
  118.     IOCheck;
  119.   end;
  120.   Close(InFile);
  121. end; { procedure GetTwoVectorsFromFile }
  122.  
  123. procedure GetTwoVectorsFromKeyboard(var NumPoints : integer;
  124.                                     var XData     : TNvector;
  125.                                     var YData     : TNvector);
  126.  
  127. {--------------------------------------------------------------}
  128. {- Output: NumPoints, XData, YData                            -}
  129. {-                                                            -}
  130. {- This procedure reads in the data points from the keyboard. -}
  131. {--------------------------------------------------------------}
  132.  
  133. var
  134.   Term : integer;
  135.  
  136. begin
  137.   NumPoints := 0;
  138.   Writeln;
  139.   repeat
  140.     Write('Number of points (0-', TNArraySize, ')? ');
  141.     Readln(NumPoints);
  142.     IOCheck;
  143.   until ((NumPoints >= 0) and (NumPoints <= TNArraySize) and not IOerr);
  144.   Writeln;
  145.   Write('Type in the X ');
  146.   Writeln('and Y values, separated by a space (not a comma):');
  147.   for Term := 1 to NumPoints do
  148.   repeat
  149.     Write('X[',Term,'], Y[',Term,']:');
  150.     Read(XData[Term], YData[Term]);
  151.     Writeln;
  152.     IOCheck;
  153.   until not IOerr;
  154. end; { procedure GetTwoVectorsFromKeyboard }
  155.  
  156. procedure GetOneVectorFromFile(var NumInter : integer;
  157.                                var XInter   : TNvector);
  158.  
  159. {------------------------------------------}
  160. {- Output: NumInter, XInter               -}
  161. {-                                        -}
  162. {- This procedure reads in the points at  -}
  163. {- which to interpolate from a data file. -}
  164. {------------------------------------------}
  165.  
  166. var
  167.   Filename : string[255];
  168.   InFile : text;
  169.  
  170. begin
  171.   Writeln;
  172.   repeat
  173.     Write('File name? ');
  174.     Readln(Filename);
  175.     Assign(InFile, Filename);
  176.     Reset(InFile);
  177.     IOCheck;
  178.   until not IOerr;
  179.   NumInter := 0;
  180.   while not EOF(InFile) do
  181.   begin
  182.     NumInter := Succ(NumInter);
  183.     Readln(InFile, XInter[NumInter]);
  184.     IOCheck;
  185.   end;
  186.   Close(InFile);
  187. end; { procedure GetOneVectorFromFile }
  188.  
  189. procedure GetOneVectorFromKeyboard(var NumInter : integer;
  190.                                    var XInter   : TNvector);
  191.  
  192. {-------------------------------------------}
  193. {- Output: NumInter, XInter                -}
  194. {-                                         -}
  195. {- This procedure reads in the points at   -}
  196. {- which to interpolate from the keyboard. -}
  197. {-------------------------------------------}
  198.  
  199. var
  200.   Term : integer;
  201.  
  202. begin
  203.   NumInter := 0;
  204.   Writeln;
  205.   repeat
  206.     Write('Number of points (0-',TNArraySize,')?');
  207.     Readln(NumInter);
  208.     IOCheck;
  209.   until((NumInter >= 0) and (NumInter <= TNArraySize) and not IOerr);
  210.   Writeln;
  211.   for Term := 1 to NumInter do
  212.   repeat
  213.     Write('Point ',Term,':');
  214.     Readln(XInter[Term]);
  215.     IOCheck;
  216.   until not IOerr;
  217. end; { procedure GetOneVectorFromKeyboard }
  218.  
  219. begin { procedure GetData }
  220.   case InputChannel('Input Data Points From') of
  221.     'K' : begin
  222.             GetTwoVectorsFromKeyboard(NumPoints, XData, YData);
  223.             repeat
  224.               Writeln;
  225.               Write('Derivative of the function at the left endpoint: ');
  226.               Readln(DerivLE);
  227.               IOCheck;
  228.             until not IOerr;
  229.             repeat
  230.               Writeln;
  231.               Write('Derivative of the function at the right endpoint: ');
  232.               Readln(DerivRE);
  233.               IOCheck;
  234.             until not IOerr;
  235.           end;
  236.     'F' : begin
  237.             GetTwoVectorsFromFile(NumPoints, XData, YData);
  238.             { The endpoint derivatives are the last data point }
  239.             NumPoints := Pred(NumPoints);
  240.             DerivLE := XData[NumPoints+1];
  241.             DerivRE := YData[NumPoints+1];
  242.           end;
  243.   end;
  244.   Writeln;
  245.   case InputChannel('Input Interpolated Points From') of
  246.     'K' : GetOneVectorFromKeyboard(NumInter, XInter);
  247.     'F' : GetOneVectorFromFile(NumInter, XInter);
  248.   end;
  249.   GetOutputFile(OutFile);
  250. end; { procedure GetData }
  251.  
  252. procedure Results(NumPoints : integer;
  253.               var XData     : TNvector;
  254.               var YData     : TNvector;
  255.                   DerivLE   : Float;
  256.                   DerivRE   : Float;
  257.               var Coef0     : TNvector;
  258.               var Coef1     : TNvector;
  259.               var Coef2     : TNvector;
  260.               var Coef3     : TNvector;
  261.                   NumInter  : integer;
  262.               var XInter    : TNvector;
  263.               var YInter    : TNvector;
  264.                   Error     : byte);
  265.  
  266. {------------------------------------------------------------}
  267. {- This procedure outputs the results to the device OutFile -}
  268. {------------------------------------------------------------}
  269.  
  270. var
  271.   Index : integer;
  272.  
  273. begin
  274.   Writeln(OutFile);
  275.   Writeln(OutFile);
  276.   Writeln(OutFile, 'Data :           X                   Y');
  277.   for Index := 1 to NumPoints do
  278.     Writeln(OutFile, Index:3, ':    ', XData[Index]:15:10, '      ',
  279.             YData[Index] : 15 : 10);
  280.   Writeln(OutFile);
  281.   Writeln(OutFile, 'Derivative at X=', XData[1], ' :  ', DerivLE);
  282.   Writeln(OutFile, 'Derivative at X=', XData[NumPoints], ' :  ', DerivRE);
  283.   Writeln(OutFile);
  284.   if Error >= 1 then
  285.     DisplayError;
  286.   case Error of
  287.     0 : begin
  288.           Writeln(OutFile, 'Splines:',  ' ':7, 'Coef0', ' ':13,
  289.                            'Coef1', ' ':14, 'Coef2', ' ':13, 'Coef3');
  290.           for Index := 1 to NumPoints-1 do
  291.             Writeln(OutFile, '  ', Index : 3, ':', '  ', Coef0[Index]:15:10,
  292.                              ' ':3, Coef1[Index]:15:10, ' ':3,
  293.                              Coef2[Index]:15:10, ' ':3, Coef3[Index]:15:10);
  294.           Writeln(OutFile);
  295.           Writeln(OutFile, 'Interpolated Points:    X                   Y');
  296.           for Index := 1 to NumInter do
  297.             Writeln(OutFile, Index:10, ':    ',
  298.                              XInter[Index] : 15 : 10, '      ',
  299.                              YInter[Index] : 15 : 10);
  300.         end;
  301.  
  302.     1 : Writeln(OutFile, 'The X points must be unique.');
  303.  
  304.     2 : Writeln(OutFile,
  305.                        'The X points must be in increasing sequential order.');
  306.  
  307.     3 : Writeln(OutFile, 'There must be at least two data points.');
  308.  
  309.   end; { case }
  310. end; { procedure Results }
  311.  
  312. begin { program CubicSplineClamped }
  313.   ClrScr;
  314.   Initialize(Coef0, Coef1, Coef2, Coef3, XData, YData, DerivLE, DerivRE,
  315.              NumPoints, NumInter, XInter, YInter, Error);
  316.   GetData(NumPoints, XData, YData, DerivLE, DerivRE, NumInter, XInter);
  317.   CubicSplineClamped(NumPoints, XData, YData, DerivLE, DerivRE, NumInter,
  318.                      XInter, Coef0, Coef1, Coef2, Coef3, YInter, Error);
  319.   Results(NumPoints, XData, YData, DerivLE, DerivRE,
  320.           Coef0, Coef1, Coef2, Coef3, NumInter, XInter, YInter, Error);
  321.   Close(OutFile);
  322. end. { program CubicSplineClamped }
  323.