home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l042 / 1.ddi / CHAP3.ARC / INTERP.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1987-12-30  |  38.3 KB  |  889 lines
  1. unit Interp;
  2.  
  3. {----------------------------------------------------------------------------}
  4. {-                                                                          -}
  5. {-     Turbo Pascal Numerical Methods Toolbox                               -}
  6. {-     Copyright (c) 1986, 87 by Borland International, Inc.                -}
  7. {-                                                                          -}
  8. {-  This unit provides procedures for performing interpolation.             -}
  9. {-                                                                          -}
  10. {----------------------------------------------------------------------------}
  11.  
  12. {$I Float.inc} { Determines the setting of the $N compiler directive }
  13.  
  14. interface
  15.  
  16. {$IFOPT N+}
  17. type
  18.   Float = Double; { 8 byte real, requires 8087 math chip }
  19.  
  20. const
  21.   TNNearlyZero = 1E-015;
  22. {$ELSE}
  23. type
  24.   Float = real;   { 6 byte real, no math chip required }
  25.  
  26. const
  27.   TNNearlyZero = 1E-07;
  28. {$ENDIF}
  29.  
  30.   TNArraySize = 50;            { Size of vectors }
  31.  
  32. type
  33.   TNvector = array[0..TNArraySize] of Float;
  34.   TNmatrix = array[0..TNArraySize] of TNvector;
  35.  
  36. procedure Lagrange(NumPoints : integer;
  37.                var XData     : TNvector;
  38.                var YData     : TNvector;
  39.                    NumInter  : integer;
  40.                var XInter    : TNvector;
  41.                var YInter    : TNvector;
  42.                var Poly      : TNvector;
  43.                var Error     : byte);
  44.  
  45. {--------------------------------------------------------------------------}
  46. {-                                                                        -}
  47. {-    Input:  NumPoints, XData, YData, NumInter, XInter                   -}
  48. {-    Output: YInter, Poly, Error                                         -}
  49. {-                                                                        -}
  50. {-           Purpose: Given a set of points (X,Y), use Lagrange           -}
  51. {-                    polynomials to construct a polynomial               -}
  52. {-                    to fit the points.                                  -}
  53. {-                    The degree of the polynomial is one less            -}
  54. {-                    than the number of points.                          -}
  55. {-                    Use this polynomial to interpolate between          -}
  56. {-                    the points.                                         -}
  57. {-                                                                        -}
  58. {- User-defined Types: TNvector = array[0..TNArraySize] of real;          -}
  59. {-                                                                        -}
  60. {-                                                                        -}
  61. {-  Global Variables: NumPoints : integer;  Number of points              -}
  62. {-                    XData     : TNvector; X-coordinate data points      -}
  63. {-                    YData     : TNvector; Y-coordinate data points      -}
  64. {-                    NumInter  : integer;  Number of interpolated points -}
  65. {-                    XInter    : TNvector; X-coordinate of points        -}
  66. {-                                          at which to interpolate       -}
  67. {-                    YInter    : TNvector; Calculated Y-coordinate of    -}
  68. {-                                          interpolated points           -}
  69. {-                    Poly      : TNvector; Coefficients of               -}
  70. {-                                          interpolating polynomial      -}
  71. {-                    Error     : byte;     Flags an error                -}
  72. {-                                                                        -}
  73. {-            Errors: 0: No errors                                        -}
  74. {-                    1: Points not unique                                -}
  75. {-                    2: NumPoints < 1                                    -}
  76. {-                                                                        -}
  77. {--------------------------------------------------------------------------}
  78.  
  79. procedure Divided_Difference(NumPoints : integer;
  80.                          var XData     : TNvector;
  81.                          var YData     : TNvector;
  82.                              NumInter  : integer;
  83.                          var XInter    : TNvector;
  84.                          var YInter    : TNvector;
  85.                          var Error     : byte);
  86.  
  87.  
  88. {----------------------------------------------------------------------------}
  89. {-                                                                          -}
  90. {-      Input:  NumPoints, XData, YData, NumInter, XInter                   -}
  91. {-      Output: YInter, Error                                               -}
  92. {-                                                                          -}
  93. {-            Purpose: This unit uses Newton's interpolary divided          -}
  94. {-                     difference equation as an interpolation algorithm.   -}
  95. {-                     This is a general divided difference equation and    -}
  96. {-                     so is not limited to uniform spacing; variable       -}
  97. {-                     spacing is allowed.  Furthermore, the data points    -}
  98. {-                     do not have to be in sequential order. The user      -}
  99. {-                     must supply the data points (X,Y) and the points     -}
  100. {-                     at which an interpolated value should be calculated. -}
  101. {-                                                                          -}
  102. {-  User-defined Types: TNvector = array[0..TNArraySize] OF real;           -}
  103. {-                      TNmatrix = array[0..TNArraySize] OF TNvector;       -}
  104. {-                                                                          -}
  105. {-   Global Variables: NumPoints : integer;   Number of data points         -}
  106. {-                     XData     : TNvector;  X-coordinate data points      -}
  107. {-                     YData     : TNvector;  Y-coordinate data points      -}
  108. {-                     NumInter  : integer;   Number of points at which     -}
  109. {-                                            to interpolate                -}
  110. {-                     XInter    : TNvector;  X-coordinate of               -}
  111. {-                                            interpolation points          -}
  112. {-                     YInter    : TNvector;  Calculated Y-coordinate       -}
  113. {-                                            of interpolation points       -}
  114. {-                     Error     : byte;      Flags if something went wrong -}
  115. {-                                                                          -}
  116. {-             Errors: 0: No Error                                          -}
  117. {-                     1: XData points not unique                           -}
  118. {-                     2: NumPoints < 1                                     -}
  119. {-                                                                          -}
  120. {----------------------------------------------------------------------------}
  121.  
  122. procedure CubicSplineFree(NumPoints : integer;
  123.                       var XData     : TNvector;
  124.                       var YData     : TNvector;
  125.                           NumInter  : integer;
  126.                       var XInter    : TNvector;
  127.                       var Coef0     : TNvector;
  128.                       var Coef1     : TNvector;
  129.                       var Coef2     : TNvector;
  130.                       var Coef3     : TNvector;
  131.                       var YInter    : TNvector;
  132.                       var Error     : byte);
  133.  
  134. {----------------------------------------------------------------------------}
  135. {-                                                                          -}
  136. {-    Input:  NumPoints, XData, YData, NumInter, XInter                     -}
  137. {-    Output: Coef0, Coef1, Coef2, Coef3, YInter, Error                     -}
  138. {-                                                                          -}
  139. {-            Purpose: To construct a cubic spline interpolant              -}
  140. {-                     to a set of points. The second derivative            -}
  141. {-                     of the interpolant is assumed to be zero             -}
  142. {-                     at the endpoints (free cubic spline)                 -}
  143. {-                     The spline is of the form                            -}
  144. {-                                                                          -}
  145. {-                        3               2                                 -}
  146. {-            D[I](X-X[I])  + C[I](X-X[I])  + B[I](X-X[I]) + A[I]           -}
  147. {-                                                                          -}
  148. {-                     where 1 < I < NumPoints.                             -}
  149. {-                                                                          -}
  150. {-  User-defined Types: TNvector = ARRAY[0..TNArraySize] OF real;           -}
  151. {-                                                                          -}
  152. {- Global Variables: NumPoints : integer;  Number of points                 -}
  153. {-                   XData     : TNvector; X-coordinate data points         -}
  154. {-                                         (must be in increasing order)    -}
  155. {-                   YData     : TNvector; Y-coordinate data points         -}
  156. {-                   Coef0     : TNvector; 0th coefficient of spline        -}
  157. {-                   Coef1     : TNvector; 1st coefficient of spline        -}
  158. {-                   Coef2     : TNvector; 2nd coefficient of spline        -}
  159. {-                   Coef3     : TNvector; 3rd coefficient of spline        -}
  160. {-                   NumInter  : integer;  Number of interpolated points    -}
  161. {-                   XInter    : TNvector; Points at which to interpolate   -}
  162. {-                   YInter    : TNvector; Interpolated values at XInter    -}
  163. {-                   Error     : byte;     0: flags if something goes wrong -}
  164. {-                                                                          -}
  165. {-              Errors: 0: No error                                         -}
  166. {-                      1: X-coordinate points not                          -}
  167. {-                         unique                                           -}
  168. {-                      2: X-coordinate points not                          -}
  169. {-                         in increasing order                              -}
  170. {-                      3: NumPoints < 2                                    -}
  171. {-                                                                          -}
  172. {----------------------------------------------------------------------------}
  173.  
  174. procedure CubicSplineClamped(NumPoints : integer;
  175.                          var XData     : TNvector;
  176.                          var YData     : TNvector;
  177.                              DerivLE   : Float;
  178.                              DerivRE   : Float;
  179.                              NumInter  : integer;
  180.                          var XInter    : TNvector;
  181.                          var Coef0     : TNvector;
  182.                          var Coef1     : TNvector;
  183.                          var Coef2     : TNvector;
  184.                          var Coef3     : TNvector;
  185.                          var YInter    : TNvector;
  186.                          var Error     : byte);
  187.  
  188. {----------------------------------------------------------------------------}
  189. {-                                                                          -}
  190. {-    Input:  NumPoints, XData, YData, DerivLE, DerivRE, NumInter, XInter   -}
  191. {-    Output: Coef0, Coef1, Coef2, Coef3, YInter, Error                     -}
  192. {-                                                                          -}
  193. {-            Purpose: To construct a cubic spline interpolant              -}
  194. {-                     to a set of points. The first derivative             -}
  195. {-                     of the interpolant is defined by the user            -}
  196. {-                     at the endpoints (clamped cubic spline)              -}
  197. {-                     The spline is of the form                            -}
  198. {-                                                                          -}
  199. {-                        3               2                                 -}
  200. {-            D[I](X-X[I])  + C[I](X-X[I])  + B[I](X-X[I]) + A[I]           -}
  201. {-                                                                          -}
  202. {-                     where 1 < I < NumPoints.                             -}
  203. {-                                                                          -}
  204. {-  User-defined Types: TNvector = ARRAY[0..TNArraySize] of real;           -}
  205. {-                                                                          -}
  206. {- Global Variables: NumPoints : integer;   Number of points                -}
  207. {-                   XData     : TNvector;  X-coordinate data points        -}
  208. {-                                          (must be in increasing order)   -}
  209. {-                   YData     : TNvector;  y-coordinate data points        -}
  210. {-                   DerivLE   : real;      derivative on the left end      -}
  211. {-                   DerivRE   : real;      derivative on the right end     -}
  212. {-                   Coef0     : TNvector;  0th coefficient of spline       -}
  213. {-                   Coef1     : TNvector;  1st coefficient of spline       -}
  214. {-                   Coef2     : TNvector;  2nd coefficient of spline       -}
  215. {-                   Coef3     : TNvector;  3rd coefficient of spline       -}
  216. {-                   NumInter  : integer;   Number of interpolated points   -}
  217. {-                   XInter    : TNvector;  Points at which to interpolate  -}
  218. {-                   YInter    : TNvector;  interpolated values at XInter   -}
  219. {-                   Error     : byte;      flags if something goes wrong   -}
  220. {-                                                                          -}
  221. {-             Errors: 0: No error                                          -}
  222. {-                     1: X-coordinate points not                           -}
  223. {-                        unique                                            -}
  224. {-                     2: X-coordinate points not                           -}
  225. {-                        in increasing order                               -}
  226. {-                     3: NumPoints < 2                                     -}
  227. {-                                                                          -}
  228. {----------------------------------------------------------------------------}
  229.  
  230. implementation
  231.  
  232. procedure Lagrange{(NumPoints : integer;
  233.                var XData      : TNvector;
  234.                var YData      : TNvector;
  235.                    NumInter   : integer;
  236.                var XInter     : TNvector;
  237.                var YInter     : TNvector;
  238.                var Poly       : TNvector;
  239.                var Error      : byte)};
  240.  
  241. var
  242.   Degree : integer;        { Degree of resulting polynomial }
  243.                            { one less than NumPoints        }
  244.  
  245. procedure SynDiv(Degree : integer;
  246.              var Poly   : TNvector;
  247.                  X      : Float;
  248.              var Y      : Float);
  249.  
  250. {----------------------------------------------------------------------}
  251. {- Input:  Degree, Poly, X                                            -}
  252. {- Output: Y                                                          -}
  253. {-                                                                    -}
  254. {- This procedure applies the technique of synthetic division         -}
  255. {- to a polynomial, Poly, at the value X.  The result is the value    -}
  256. {- of the polynomial at X.                                            -}
  257. {----------------------------------------------------------------------}
  258.  
  259. var
  260.   Index : integer;
  261.  
  262. begin
  263.   Y := Poly[Degree];
  264.   for Index := Degree-1 downto 0 do
  265.     Y := Y * X + Poly[Index];
  266. end; { procedure SynDiv }
  267.  
  268. procedure MakePolynomial(Degree : integer;
  269.                      var XData  : TNvector;
  270.                      var YData  : TNvector;
  271.                      var Poly   : TNvector;
  272.                      var Error  : byte);
  273.  
  274. {----------------------------------------------------------------}
  275. {- Input: Degree, XData, YData                                  -}
  276. {- Output: Poly, Error                                          -}
  277. {-                                                              -}
  278. {- This procedure constructs the polynomial which fits the data -}
  279. {- points (XData, YData).  The coefficients of the polynomial   -}
  280. {- are returned in Poly.  If the XData points are not unique,   -}
  281. {- then Error = 1 is returned.                                  -}
  282. {----------------------------------------------------------------}
  283.  
  284. var
  285.   T, D : TNvector;       { Iterative variables for computing }
  286.                          { the Lagrange Polynomial           }
  287.   Num, Denom, Coef : Float;
  288.   Term, Term2 : integer;
  289.  
  290. begin
  291.   for Term := 1 to Degree + 1 do
  292.   begin
  293.     T[Term] := 0;
  294.     D[Term] := 0
  295.   end;
  296.   T[0] := YData[1];
  297.   D[0] := -XData[1];
  298.   D[1] := 1;
  299.   for Term := 1 to Degree do
  300.   begin
  301.     { Evaluate D at XData[Term] }
  302.     SynDiv(Term, D, XData[Term + 1], Denom);
  303.     if ABS(Denom) < TNNearlyZero then
  304.       Error := 1    { Points not unique }
  305.     else
  306.       begin
  307.         { Evaluate T at XData[Term] }
  308.         SynDiv(Term - 1, T, XData[Term + 1], Num);
  309.         Coef := (YData[Term + 1] - Num) / Denom;
  310.         for Term2 := Term downto 0 do
  311.         begin
  312.           T[Term2] := Coef * D[Term2] + T[Term2];
  313.           D[Term2 + 1] := D[Term2] - XData[Term + 1] * D[Term2 + 1];
  314.         end;
  315.         D[0] := -XData[Term + 1] * D[0];
  316.       end;
  317.   end;
  318.   Poly := T;
  319. end; { procedure MakePoly }
  320.  
  321. procedure SolvePolynomial(Degree   : integer;
  322.                           NumInter : integer;
  323.                       var XInter   : TNvector;
  324.                       var Poly     : TNvector;
  325.                       var YInter   : TNvector);
  326.  
  327. {--------------------------------------------------------------}
  328. {- Input:  Degree, NumInter, XInter, Poly                     -}
  329. {- Output: YInter                                             -}
  330. {-                                                            -}
  331. {- This procedure uses the technique of synthetic division to -}
  332. {- solve the polynomial, Poly, at the X values contained in   -}
  333. {- XInter.  These interpolated values are returned in YInter. -}
  334. {--------------------------------------------------------------}
  335.  
  336. var
  337.   Term : integer;
  338.  
  339. begin
  340.   for Term := 1 to NumInter do
  341.     SynDiv(Degree, Poly, XInter[Term], YInter[Term]);
  342. end; { procedure SolvePolynomial }
  343.  
  344. begin { procedure Lagrange }
  345.   Degree := NumPoints - 1;
  346.   Error := 0;
  347.   if NumPoints < 1 then
  348.     Error := 2
  349.   else
  350.     begin
  351.       MakePolynomial(Degree, XData, YData, Poly, Error);
  352.       SolvePolynomial(Degree, NumInter, XInter, Poly, YInter);
  353.     end;
  354. end; { procedure Lagrange }
  355.  
  356. procedure Divided_Difference{(NumPoints : integer;
  357.                          var XData      : TNvector;
  358.                          var YData      : TNvector;
  359.                              NumInter   : integer;
  360.                          var XInter     : TNvector;
  361.                          var YInter     : TNvector;
  362.                          var Error      : byte)};
  363.  
  364.  
  365. var
  366.   DividedDif : TNmatrix;              { Divided difference table }
  367.  
  368. procedure MakeDivDifTable(NumPoints : integer;
  369.                       var XData     : TNvector;
  370.                       var YData     : TNvector;
  371.                       var DivDif    : TNmatrix;
  372.                       var Error     : byte);
  373.  
  374. {--------------------------------------------------------------}
  375. {- Input: NumPoints, XData, YData                             -}
  376. {- Output: DivDif, Error                                      -}
  377. {-                                                            -}
  378. {- This procedure constructs a NumPoints X NumPoints divided  -}
  379. {- difference table (DivDif). If the X values are not unique, -}
  380. {- then Error = 1 is returned.                                -}
  381. {--------------------------------------------------------------}
  382.  
  383. var
  384.   Row, Column : integer;
  385.   Den : Float;
  386.  
  387. begin
  388.   Error := 0;
  389.   for Row := 1 to NumPoints do
  390.     DivDif[Row, 1] := YData[Row];
  391.  
  392.   for Column := 2 to NumPoints do
  393.     for Row := Column to NumPoints do
  394.     begin
  395.       Den := XData[Row] - XData[Row - Column + 1];
  396.       if ABS(Den) < TNNearlyZero then
  397.         Error := 1
  398.       else
  399.         DivDif[Row, Column] :=
  400.                    (DivDif[Row, Column-1] - DivDif[Row-1, Column-1]) / Den;
  401.     end;
  402. end; { procedure MakeDivDifTable }
  403.  
  404. procedure EvaluateNewtonFormula(NumPoints : integer;
  405.                             var DivDif    : TNmatrix;
  406.                             var XData     : TNvector;
  407.                                 NumInter  : integer;
  408.                             var XInter    : TNvector;
  409.                             var YInter    : TNvector);
  410.  
  411. {---------------------------------------------------------------------}
  412. {- Input:  NumPoints, DivDif, XData, NumInter, XInter                -}
  413. {- Output: YInter                                                    -}
  414. {-                                                                   -}
  415. {- Newton's interpolary divided difference formula interpolates      -}
  416. {- between data points if a divided difference table of those points -}
  417. {- exists.  This procedure interpolates at the points contained in   -}
  418. {- XInter using Newton's interpolary  divided difference formula.    -}
  419. {- The interpolated values are returned in YInter.                   -}
  420. {---------------------------------------------------------------------}
  421.  
  422. var
  423.   Index, Coef : integer;
  424.   Factor, Unknown : Float;
  425.  
  426. begin { procedure EvaluateNewtonFormula }
  427.   for Index := 1 to NumInter do
  428.   begin
  429.     YInter[Index] := 0;
  430.     Unknown := XInter[Index];
  431.     for Coef := 1 to NumPoints do
  432.     begin
  433.       if Coef = 1 then
  434.         Factor := 1
  435.       else
  436.         Factor := Factor * (Unknown - XData[Coef-1]);
  437.       YInter[Index] := YInter[Index] + Factor * DivDif[Coef, Coef];
  438.     end;
  439.   end;
  440. end; { procedure EvaluateNewtonFormula }
  441.  
  442. begin { procedure Divided_Difference }
  443.   Error := 0;
  444.   if NumPoints < 1 then
  445.     Error := 2
  446.   else
  447.     MakeDivDifTable(NumPoints, XData, YData, DividedDif, Error);
  448.   if Error = 0 then
  449.     EvaluateNewtonFormula(NumPoints, DividedDif, XData,
  450.                           NumInter, XInter, YInter);
  451. end; { procedure Divided_Difference }
  452.  
  453. procedure CubicSplineFree{(NumPoints : integer;
  454.                       var XData      : TNvector;
  455.                       var YData      : TNvector;
  456.                           NumInter   : integer;
  457.                       var XInter     : TNvector;
  458.                       var Coef0      : TNvector;
  459.                       var Coef1      : TNvector;
  460.                       var Coef2      : TNvector;
  461.                       var Coef3      : TNvector;
  462.                       var YInter     : TNvector;
  463.                       var Error      : byte)};
  464.  
  465. type
  466.   TNcoefficients = record
  467.                      A, B, C, D : TNvector;
  468.                    end;
  469.  
  470. var
  471.   Interval : TNvector;        { Intervals between adjacent points }
  472.   Spline : TNcoefficients;    { All the cubics }
  473.  
  474. procedure CalculateIntervals(NumPoints : integer;
  475.                          var XData     : TNvector;
  476.                          var Interval  : TNvector;
  477.                          var Error     : byte);
  478.  
  479. {----------------------------------------------------------}
  480. {- Input: NumPoints, XData                                -}
  481. {- Output: Interval, Error                                -}
  482. {-                                                        -}
  483. {- This procedure calculates the length of the interval   -}
  484. {- between two adjacent X values, contained in XData. If  -}
  485. {- the X values are not sequential, Error = 2 is returned -}
  486. {- and if the X values are not unique, then Error = 1 is  -}
  487. {- returned.                                              -}
  488. {----------------------------------------------------------}
  489.  
  490. var
  491.   Index : integer;
  492.  
  493. begin
  494.   Error := 0;
  495.   for Index := 1 to NumPoints - 1 do
  496.   begin
  497.     Interval[Index] := XData[Index+1] - XData[Index];
  498.     if ABS(Interval[Index]) < TNNearlyZero then
  499.       Error := 1;     { Data not unique }
  500.     if Interval[Index] < 0 then
  501.       Error := 2;     { Data not in increasing order }
  502.   end;
  503. end; { procedure CalculateIntervals }
  504.  
  505. procedure CalculateCoefficients(NumPoints : integer;
  506.                             var XData     : TNvector;
  507.                             var YData     : TNvector;
  508.                             var Interval  : TNvector;
  509.                             var Spline    : TNcoefficients);
  510.  
  511. {---------------------------------------------------------------}
  512. {- Input: NumPoints, XData, YData, Interval                    -}
  513. {- Output: Spline                                              -}
  514. {-                                                             -}
  515. {- This procedure calculates the coefficients of each cubic    -}
  516. {- in the interpolating spline. A separate cubic is calculated -}
  517. {- for every interval between data points.  The coefficients   -}
  518. {- are returned in the variable Spline.                        -}
  519. {---------------------------------------------------------------}
  520.  
  521. procedure CalculateAs(NumPoints : integer;
  522.                   var YData     : TNvector;
  523.                   var Spline    : TNcoefficients);
  524.  
  525. {------------------------------------------}
  526. {- Input: NumPoints, YData                -}
  527. {- Ouput: Spline                          -}
  528. {-                                        -}
  529. {- This procedure calculates the constant -}
  530. {- Term in each cubic.                    -}
  531. {------------------------------------------}
  532.  
  533. var
  534.   Index : integer;
  535. begin
  536.   for Index := 1 to NumPoints do
  537.     Spline.A[Index] := YData[Index];
  538. end; { procedure CalculateAs }
  539.  
  540. procedure CalculateCs(NumPoints : integer;
  541.                   var XData     : TNvector;
  542.                   var Interval  : TNvector;
  543.                   var Spline    : TNcoefficients);
  544.  
  545. {---------------------------------------------}
  546. {- Input: NumPoints, XData, Interval         -}
  547. {- Ouput: Spline                             -}
  548. {-                                           -}
  549. {- This procedure calculates the coefficient -}
  550. {- of the squared Term in each cubic.        -}
  551. {---------------------------------------------}
  552.  
  553. var
  554.   Alpha, L, Mu, Z : TNvector;
  555.   Index : integer;
  556. begin
  557.   with Spline do
  558.   begin
  559.     { The next few lines solve a tridiagonal matrix }
  560.     for Index := 2 to NumPoints - 1 do
  561.       Alpha[Index] := 3 * ((A[Index+1] * Interval[Index-1])
  562.                          - (A[Index] * (XData[Index+1] - XData[Index-1]))
  563.                          + (A[Index-1] * Interval[Index]))
  564.                          / (Interval[Index-1] * Interval[Index]);
  565.     L[1] := 0;
  566.     Mu[1] := 0;
  567.     Z[1] := 0;
  568.     for Index := 2 to NumPoints - 1 do
  569.     begin
  570.       L[Index] := 2 * (XData[Index+1] - XData[Index-1])
  571.                     - Interval[Index-1] * Mu[Index-1];
  572.       Mu[Index] := Interval[Index]/L[Index];
  573.       Z[Index] := (Alpha[Index] - Interval[Index-1] * Z[Index-1]) / L[Index];
  574.     end;
  575.     { Now calculate the C's }
  576.     C[NumPoints] := 0;
  577.     for Index := NumPoints - 1 downto 1 do
  578.       C[Index] := Z[Index] - Mu[Index] * C[Index+1];
  579.   end; { with }
  580. end; { procedure CalculateCs }
  581.  
  582. procedure CalculateBandDs(NumPoints : integer;
  583.                       var Interval  : TNvector;
  584.                       var Spline    : TNcoefficients);
  585.  
  586. {------------------------------------------------}
  587. {- Input: NumPoints, Interval                   -}
  588. {- Ouput: Spline                                -}
  589. {-                                              -}
  590. {- This procedure calculates the coefficient of -}
  591. {- the linear and cubic terms in each cubic.    -}
  592. {------------------------------------------------}
  593.  
  594. var
  595.   Index : integer;
  596. begin
  597.   with Spline do
  598.     for Index := NumPoints - 1 downto 1 do
  599.     begin
  600.       B[Index] := (A[Index+1] - A[Index]) / Interval[Index]
  601.                    - Interval[Index] * (C[Index+1] + 2 * C[Index]) / 3;
  602.       D[Index] := (C[Index+1] - C[Index]) / (3 * Interval[Index]);
  603.     end;
  604. end; { procedure CalculateDs }
  605.  
  606. begin { procedure CalculateCoefficients }
  607.   CalculateAs(NumPoints, YData, Spline);
  608.   CalculateCs(NumPoints, XData, Interval, Spline);
  609.   CalculateBandDs(NumPoints, Interval, Spline);
  610. end; { procedure CalculateCoefficients }
  611.  
  612. procedure Interpolate(NumPoints : integer;
  613.                   var XData     : TNvector;
  614.                   var Spline    : TNcoefficients;
  615.                       NumInter  : integer;
  616.                   var XInter    : TNvector;
  617.                   var YInter    : TNvector);
  618.  
  619. {---------------------------------------------------------------}
  620. {- Input: NumPoints, XData, Spline, NumInter, XInter           -}
  621. {- Output: YInter                                              -}
  622. {-                                                             -}
  623. {- This procedure uses the interpolating cubic spline (Spline) -}
  624. {- to interpolate values for the X positions contained         -}
  625. {- in XInter.  The interpolated values are returned in YInter. -}
  626. {---------------------------------------------------------------}
  627.  
  628. var
  629.   Index, Location, Term : integer;
  630.   X : Float;
  631.  
  632. begin
  633.   for Index := 1 to NumInter do
  634.   begin
  635.     Location := 1;
  636.     for Term := 1 to NumPoints - 1 do
  637.       if XInter[Index] > XData[Term] then
  638.         Location := Term;
  639.     X := XInter[Index] - XData[Location];
  640.     with Spline do
  641.       YInter[Index] := ((D[Location] * X + C[Location]) * X +
  642.                          B[Location]) * X + A[Location];
  643.   end;
  644. end; { procedure Interpolate }
  645.  
  646. begin { procedure CubicSplineFree }
  647.   Error := 0;
  648.   if NumPoints < 2 then
  649.     Error := 3
  650.   else
  651.     CalculateIntervals(NumPoints, XData, Interval, Error);
  652.   if Error = 0 then
  653.   begin
  654.     CalculateCoefficients(NumPoints, XData, YData, Interval, Spline);
  655.     Interpolate(NumPoints, XData, Spline, NumInter, XInter, YInter);
  656.   end;
  657.   Coef0 := Spline.A;
  658.   Coef1 := Spline.B;
  659.   Coef2 := Spline.C;
  660.   Coef3 := Spline.D;
  661. end; { procedure CubicSplineFree }
  662.  
  663. procedure CubicSplineClamped{(NumPoints : integer;
  664.                          var XData      : TNvector;
  665.                          var YData      : TNvector;
  666.                              DerivLE    : Float;
  667.                              DerivRE    : Float;
  668.                              NumInter   : integer;
  669.                          var XInter     : TNvector;
  670.                          var Coef0      : TNvector;
  671.                          var Coef1      : TNvector;
  672.                          var Coef2      : TNvector;
  673.                          var Coef3      : TNvector;
  674.                          var YInter     : TNvector;
  675.                          var Error      : byte)};
  676.  
  677. type
  678.   TNcoefficient = record
  679.                     A, B, C, D : TNvector;
  680.                   end;
  681.  
  682. var
  683.   Interval : TNvector;
  684.   Spline : TNcoefficient;
  685.  
  686. procedure CalculateIntervals(NumPoints : integer;
  687.                          var XData     : TNvector;
  688.                          var Interval  : TNvector;
  689.                          var Error     : byte);
  690.  
  691. {----------------------------------------------------------}
  692. {- Input: NumPoints, XData                                -}
  693. {- Output: Interval, Error                                -}
  694. {-                                                        -}
  695. {- This procedure calculates the length of the interval   -}
  696. {- between two adjacent X values, contained in XData. If  -}
  697. {- the X values are not sequential, Error = 2 is returned -}
  698. {- and if the X values are not unique, then Error = 1 is  -}
  699. {- returned.                                              -}
  700. {----------------------------------------------------------}
  701.  
  702. var
  703.   Index : integer;
  704.  
  705. begin
  706.   Error := 0;
  707.   for Index := 1 to NumPoints - 1 do
  708.   begin
  709.     Interval[Index] := XData[Index+1] - XData[Index];
  710.     if ABS(Interval[Index]) < TNNearlyZero then
  711.       Error := 1;    { Points not unique }
  712.  
  713.     if Interval[Index] < 0 then
  714.       Error := 2;    { Points not in increasing order }
  715.   end;
  716. end; { procedure CalculateIntervals }
  717.  
  718. procedure CalculateCoefficients(NumPoints : integer;
  719.                             var XData     : TNvector;
  720.                             var YData     : TNvector;
  721.                             var Interval  : TNvector;
  722.                                 DerivLE   : Float;
  723.                                 DerivRE   : Float;
  724.                             var Spline    : TNcoefficient);
  725.  
  726. {---------------------------------------------------------------}
  727. {- Input: NumPoints, XData, YData, Interval, DerivLE, DerivRE  -}
  728. {- Output: Spline                                              -}
  729. {-                                                             -}
  730. {- This procedure calculates the coefficients of each cubic    -}
  731. {- in the interpolating spline. A separate cubic is calculated -}
  732. {- for every interval between data points.  The coefficients   -}
  733. {- are returned in the variable Spline.                        -}
  734. {---------------------------------------------------------------}
  735.  
  736. procedure CalculateAs(NumPoints : integer;
  737.                   var YData     : TNvector;
  738.                   var Spline    : TNcoefficient);
  739.  
  740. {------------------------------------------}
  741. {- Input: NumPoints, YData                -}
  742. {- Ouput: Spline                          -}
  743. {-                                        -}
  744. {- This procedure calculates the constant -}
  745. {- term in each cubic.                    -}
  746. {------------------------------------------}
  747.  
  748. var
  749.   Index : integer;
  750. begin
  751.   for Index := 1 to NumPoints do
  752.     Spline.A[Index] := YData[Index];
  753. end; { procedure CalculateAs }
  754.  
  755. procedure CalculateCs(NumPoints : integer;
  756.                   var XData     : TNvector;
  757.                   var Interval  : TNvector;
  758.                       DerivLE   : Float;
  759.                       DerivRE   : Float;
  760.                   var Spline    : TNcoefficient);
  761.  
  762. {---------------------------------------------}
  763. {- Input: NumPoints, XData, Interval,        -}
  764. {-        DerivLE, DerivRE                   -}
  765. {- Ouput: Spline                             -}
  766. {-                                           -}
  767. {- This procedure calculates the coefficient -}
  768. {- of the squared term in each cubic.        -}
  769. {---------------------------------------------}
  770.  
  771. var
  772.   Alpha, L, Mu, Z : TNvector;
  773.   Index : integer;
  774. begin
  775.   with Spline do
  776.   begin
  777.     { The next few lines solve a tridiagonal matrix }
  778.     Alpha[1] := 3*(A[2] - A[1]) / Interval[1] - 3 * DerivLE;
  779.     Alpha[NumPoints] := 3 * DerivRE - 3 * (A[NumPoints] - A[NumPoints - 1])
  780.                         / Interval[NumPoints - 1];
  781.     for Index := 2 to NumPoints - 1 do
  782.        Alpha[Index] := 3 * ((A[Index+1] * Interval[Index-1])
  783.                        - (A[Index] * (XData[Index+1] - XData[Index-1]))
  784.                        + (A[Index-1] * Interval[Index]))
  785.                        / (Interval[Index-1] * Interval[Index]);
  786.     L[1] := 2 * Interval[1];
  787.     Mu[1] := 0.5;
  788.     Z[1] := Alpha[1] / L[1];
  789.     for Index := 2 to NumPoints - 1 do
  790.     begin
  791.       L[Index] := 2 * (XData[Index+1] - XData[Index-1]) -
  792.                   Interval[Index-1] * Mu[Index-1];
  793.       Mu[Index] := Interval[Index] / L[Index];
  794.       Z[Index] := (Alpha[Index] - Interval[Index-1] * Z[Index-1]) / L[Index];
  795.     end;
  796.  
  797.     { Now calculate the C's }
  798.     C[NumPoints] := (Alpha[NumPoints] -
  799.                     Interval[NumPoints - 1] * Z[NumPoints - 1])
  800.                     / (Interval[NumPoints - 1] * (2-Mu[NumPoints - 1]));
  801.     for Index := NumPoints - 1 downto 1 do
  802.       C[Index] := Z[Index] - Mu[Index] * C[Index+1];
  803.   end; { with }
  804. end; { procedure CalculateCs }
  805.  
  806. procedure CalculateBandDs(NumPoints : integer;
  807.                       var Interval  : TNvector;
  808.                       var Spline    : TNcoefficient);
  809.  
  810. {------------------------------------------------}
  811. {- Input: NumPoints, Interval                   -}
  812. {- Ouput: Spline                                -}
  813. {-                                              -}
  814. {- This procedure calculates the coefficient of -}
  815. {- the linear and cubic terms in each cubic.    -}
  816. {------------------------------------------------}
  817.  
  818. var
  819.   Index : integer;
  820. begin
  821.   with Spline do
  822.     for Index := NumPoints - 1 downto 1 do
  823.     begin
  824.       B[Index] := (A[Index+1] - A[Index]) / Interval[Index]
  825.                   - Interval[Index] * (C[Index+1] + 2 * C[Index]) / 3;
  826.       D[Index] := (C[Index+1] - C[Index]) / (3 * Interval[Index]);
  827.     end;
  828. end; { procedure CalculateBandDs }
  829.  
  830. begin { procedure CalculateCoefficients }
  831.   CalculateAs(NumPoints, YData, Spline);
  832.   CalculateCs(NumPoints, XData, Interval, DerivLE, DerivRE, Spline);
  833.   CalculateBandDs(NumPoints, Interval, Spline);
  834. end; { procedure CalculateCoefficients }
  835.  
  836. procedure Interpolate(NumPoints : integer;
  837.                   var XData     : TNvector;
  838.                   var Spline    : TNcoefficient;
  839.                       NumInter  : integer;
  840.                   var XInter    : TNvector;
  841.                   var YInter    : TNvector);
  842.  
  843. {---------------------------------------------------------------}
  844. {- Input:  NumPoints, XData, Spline, NumInter, XInter          -}
  845. {- Output: YInter                                              -}
  846. {-                                                             -}
  847. {- This procedure uses the interpolating cubic spline (Spline) -}
  848. {- to interpolate values for the X positions contained         -}
  849. {- in XInter.  The interpolated values are returned in YInter. -}
  850. {---------------------------------------------------------------}
  851.  
  852. var
  853.   Index, Location, Term : integer;
  854.   X : Float;
  855.  
  856. begin
  857.   for Index := 1 to NumInter do
  858.   begin
  859.     Location := 1;
  860.     for Term := 1 to NumPoints - 1 do
  861.       if XInter[Index] > XData[Term] then
  862.         Location := Term;
  863.     X := XInter[Index] - XData[Location];
  864.     with Spline do
  865.       YInter[Index] := ((D[Location] * X + C[Location]) * X +
  866.                        B[Location]) * X + A[Location];
  867.   end;
  868. end; { procedure Interpolate }
  869.  
  870. begin { procedure CubicSplineClamped }
  871.   Error := 0;
  872.   if NumPoints < 2 then
  873.     Error := 3
  874.   else
  875.     CalculateIntervals(NumPoints, XData, Interval, Error);
  876.   if Error = 0 then
  877.   begin
  878.     CalculateCoefficients(NumPoints, XData, YData, Interval,
  879.                           DerivLE, DerivRE, Spline);
  880.     InterPolate(NumPoints, XData, Spline, NumInter, XInter, YInter);
  881.   end;
  882.   Coef0 := Spline.A;
  883.   Coef1 := Spline.B;
  884.   Coef2 := Spline.C;
  885.   Coef3 := Spline.D;
  886. end; { procedure CubicSplineClamped }
  887.  
  888. end. { Interp }
  889.