home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l042 / 1.ddi / CHAP4.ARC / DIFFER.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1987-12-30  |  59.9 KB  |  1,334 lines
  1. unit Differ;
  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 numerical differentiation. -}
  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. {$ELSE}
  20. type
  21.   Float = real;   { 6 byte real, no math chip required }
  22. {$ENDIF}
  23.  
  24. const
  25.   TNArraySize = 100;              { Size of the vectors }
  26.  
  27. type
  28.   TNvector = array[1..TNArraySize] of Float;
  29.  
  30. procedure First_Derivative(NumPoints : integer;
  31.                        var XData     : TNvector;
  32.                        var YData     : TNvector;
  33.                            Point     : byte;
  34.                            NumDeriv  : integer;
  35.                        var XDeriv    : TNvector;
  36.                        var YDeriv    : TNvector;
  37.                        var Error     : byte);
  38.  
  39. {----------------------------------------------------------------------------}
  40. {-                                                                          -}
  41. {-    Input: NumPoints, XData, YData, Point, NumDeriv, XDeriv               -}
  42. {-   Output: YDeriv, Error                                                  -}
  43. {-                                                                          -}
  44. {-  Purpose: Given a set of points (X,Y) to a function f, determine the     -}
  45. {-           first derivative of f at the points XDeriv.  All the XDeriv    -}
  46. {-           must be in the data set (X,Y).  Either 2, 3 or 5 point         -}
  47. {-           differentiation may be used.                                   -}
  48. {-                                                                          -}
  49. {-  User-defined Types:  TNvector = array[1..TNArraySize] of real;          -}
  50. {-                                                                          -}
  51. {-    Global Variables:  NumPoints : integer;  Number of points             -}
  52. {-                       XData     : TNvector; X-coordinate data points     -}
  53. {-                       YData     : TNvector; Y-coordinate data points     -}
  54. {-                       NumDeriv  : integer;  Number points at which to    -}
  55. {-                                             take the derivative          -}
  56. {-                       XDeriv    : TNvector; Point at which to take       -}
  57. {-                                             the derivative               -}
  58. {-                       YDeriv    : TNvector; Derivative at                -}
  59. {-                                             points XDeriv                -}
  60. {-                       Point     : byte;     2, 3, 5 point diff.          -}
  61. {-                       Error     : byte;     Flags if something goes      -}
  62. {-                                             wrong                        -}
  63. {-                                                                          -}
  64. {-              Errors:  0: No errors                                       -}
  65. {-                       1: Warning - not all the 2nd derivatives           -}
  66. {-                          were calculated; NOT A FATAL ERROR! See Note    -}
  67. {-                       2: Data not unique                                 -}
  68. {-                       3: Data not in increasing order                    -}
  69. {-                       4: Not enough data                                 -}
  70. {-                       5: NOT(Point in [2, 3, 5])                         -}
  71. {-                       6: Data points not evenly spaced for 5 point       -}
  72. {-                                                                          -}
  73. {-                Note:  If the point at which to take the derivative       -}
  74. {-                       is not among the data points, the value of         -}
  75. {-                       -9.9999999E35 is arbitrarily assigned as the       -}
  76. {-                       derivative at that point. Error = 1 is returned.   -}
  77. {-                       When trying to do 5 point differentiation          -}
  78. {-                       with only 5 points, there is not enough            -}
  79. {-                       information to calculate the 2nd derivative at     -}
  80. {-                       the 1st, 2nd, 4th or 5th point. Likewise, if there -}
  81. {-                       are only 6 points, there is not enough information -}
  82. {-                       to calculate the derivative at the 2nd or 5th      -}
  83. {-                       point. Should an attempt be made to evaluate the   -}
  84. {-                       derivative at any of these points, the             -}
  85. {-                       value of 9.9999999E35 is arbitrarily assigned      -}
  86. {-                       as the derivative at that point. Error = 1 is      -}
  87. {-                       returned.                                          -}
  88. {-                                                                          -}
  89. {----------------------------------------------------------------------------}
  90.  
  91. procedure Second_Derivative(NumPoints : integer;
  92.                         var XData     : TNvector;
  93.                         var YData     : TNvector;
  94.                             Point     : byte;
  95.                             NumDeriv  : integer;
  96.                         var XDeriv    : TNvector;
  97.                         var YDeriv    : TNvector;
  98.                         var Error     : byte);
  99.  
  100. {----------------------------------------------------------------------------}
  101. {-                                                                          -}
  102. {-    Input: NumPoints, XData, YData, Point, NumDeriv, XDeriv               -}
  103. {-   Output: YDeriv, Error                                                  -}
  104. {-                                                                          -}
  105. {-  Purpose: Given a set of points (X,Y) to a function f, determine the     -}
  106. {-           second derivative of f at the points XDeriv.  All the XDeriv   -}
  107. {-           must be in the data set (X,Y).  Either 3 or 5 point            -}
  108. {-           second differentiation may be used.                            -}
  109. {-                                                                          -}
  110. {-  User-defined Types:  TNvector = array[1..TNArraySize] of real;          -}
  111. {-                                                                          -}
  112. {-    Global Variables:  NumPoints : integer;  Number of points             -}
  113. {-                       XData     : TNvector; X-coordinate data points     -}
  114. {-                       YData     : TNvector; Y-coordinate data points     -}
  115. {-                       NumDeriv  : integer;  Number points at which to    -}
  116. {-                                             take the 2nd derivative      -}
  117. {-                       XDeriv    : TNvector; Point at which to take       -}
  118. {-                                             the 2nd derivative           -}
  119. {-                       YDeriv    : TNvector; 2nd derivative at            -}
  120. {-                                             points XDeriv                -}
  121. {-                       Point     : byte;     3 or 5 point diff.           -}
  122. {-                       Error     : byte;     Flags if something goes      -}
  123. {-                                             wrong                        -}
  124. {-                                                                          -}
  125. {-              Errors:  0: No errors                                       -}
  126. {-                       1: Warning, at least one point was                 -}
  127. {-                          not calculated. NOT A FATAL ERROR! (See Note)   -}
  128. {-                       2: Data not unique                                 -}
  129. {-                       3: Data not in increasing order                    -}
  130. {-                       4: Not enough data                                 -}
  131. {-                       5: not(Point IN [3, 5])                            -}
  132. {-                       6: Data points not evenly spaced                   -}
  133. {-                                                                          -}
  134. {-                Note:  If the point at which to take the 2nd derivative   -}
  135. {-                       is not among the data points, the value of         -}
  136. {-                       NotADataPoint is arbitrarily assigned as the       -}
  137. {-                       derivative at that point. Error = 1 is returned.   -}
  138. {-                       When trying to do 5 point 2nd differentiation      -}
  139. {-                       with only 5 points, there is not enough            -}
  140. {-                       information to calculate the 2nd derivative at     -}
  141. {-                       the 2nd or 4th point.  Should an attempt be made   -}
  142. {-                       to evaluate the 2nd derivative at these points,    -}
  143. {-                       the value of NotEnoughPoints is arbitrarily        -}
  144. {-                       assigned as the 2nd derivative at that point.      -}
  145. {-                       Error = 1 is returned.                             -}
  146. {-                                                                          -}
  147. {----------------------------------------------------------------------------}
  148.  
  149. procedure Interpolate_Derivative(NumPoints : integer;
  150.                              var XData     : TNvector;
  151.                              var YData     : TNvector;
  152.                                  NumDeriv  : integer;
  153.                              var XDeriv    : TNvector;
  154.                              var YInter    : TNvector;
  155.                              var YDeriv    : TNvector;
  156.                              var YDeriv2   : TNvector;
  157.                              var Error     : byte);
  158.  
  159.  
  160. {----------------------------------------------------------------------------}
  161. {-                                                                          -}
  162. {-    Input:  NumPoints, XData, YData, NumDeriv, XDeriv                     -}
  163. {-    Output: YInter, YDeriv, YDeriv2, Error                                -}
  164. {-                                                                          -}
  165. {-            Purpose: To construct a cubic spline interpolant              -}
  166. {-                     and its derivative given a set of points.            -}
  167. {-                     Use these interpolants to find the value and         -}
  168. {-                     first and second derivative at a set of points       -}
  169. {-                     XDeriv.  The second derivative of the interpolant    -}
  170. {-                     is assumed to be zero at the endpoints.              -}
  171. {-                     The spline is of the form                            -}
  172. {-                                                                          -}
  173. {-                                                2              3          -}
  174. {-            A[I]  + B[I](X-X[I])  + C[I](X-X[I]) + D[I](X-X[I])           -}
  175. {-                                                                          -}
  176. {-                     where 1 < I < NumPoints - 1.                         -}
  177. {-                                                                          -}
  178. {-  User-defined Types: TNvector = array[0..TNArraySize] of real;           -}
  179. {-                                                                          -}
  180. {- Global Variables: NumPoints : integer;  Number of points                 -}
  181. {-                   XData     : TNvector; X-coordinate data points         -}
  182. {-                                         (must be in increasing order)    -}
  183. {-                   YData     : TNvector; Y-coordinate data points         -}
  184. {-                   NumDeriv  : integer;  Number of interpolated points    -}
  185. {-                   XDeriv    : TNvector; Points at which to interpolate   -}
  186. {-                   YInter    : TNvector; Interpolated values at XDeriv    -}
  187. {-                   YDeriv    : TNvector; Interpolated first derivative    -}
  188. {-                                         at XDeriv                        -}
  189. {-                   YDeriv2   : TNvector; Interpolated second derivative   -}
  190. {-                                         at XDeriv                        -}
  191. {-                   Error     : byte;     Flags if something goes wrong    -}
  192. {-                                                                          -}
  193. {-             Errors: 0: No error                                          -}
  194. {-                     1: X-coordinate points not                           -}
  195. {-                        unique                                            -}
  196. {-                     2: X-coordinate points not                           -}
  197. {-                        in increasing order                               -}
  198. {-                     3: NumPoints < 2                                     -}
  199. {-                                                                          -}
  200. {----------------------------------------------------------------------------}
  201.  
  202. procedure FirstDerivative(NumDeriv : integer;
  203.                      var XDeriv    : TNvector;
  204.                      var YDeriv    : TNvector;
  205.                          Tolerance : Float;
  206.                      var Error     : byte;
  207.                          FuncPtr   : Pointer);
  208.  
  209. {----------------------------------------------------------------------------}
  210. {-                                                                          -}
  211. {-    Input: NumDeriv, XDeriv, Tolerance                                    -}
  212. {-   Output: YDeriv, Error                                                  -}
  213. {-                                                                          -}
  214. {-  Purpose: Given a function TNTargetF(X), approximate the first           -}
  215. {-           derivative of TNTargetF(X) at the points XDeriv. The algorithm -}
  216. {-           uses a three point formula and Richardson extrapolation        -}
  217. {-                                                                          -}
  218. {- User-defined Function:  FUNCTION TNTargetF(X : real) : real;             -}
  219. {-                                                                          -}
  220. {-  User-defined Types:  TNvector = array[1..TNArraySize] of real;          -}
  221. {-                                                                          -}
  222. {-    Global Variables:  NumDeriv  : integer;   Number points at which to   -}
  223. {-                                              take the derivative         -}
  224. {-                       XDeriv    : TNvector;  Point at which to take      -}
  225. {-                                              the derivative              -}
  226. {-                       YDeriv    : TNvector;  Derivative at               -}
  227. {-                                              points XDeriv               -}
  228. {-                       Tolerance : real;      Tolerance in answer         -}
  229. {-                       Error     : byte;      Flags if something goes     -}
  230. {-                                              wrong                       -}
  231. {-                                                                          -}
  232. {-              Errors:  0: No errors                                       -}
  233. {-                       1: Tolerance < TNNearlyZero                        -}
  234. {-                                                                          -}
  235. {----------------------------------------------------------------------------}
  236.  
  237. procedure SecondDerivative(NumDeriv  : integer;
  238.                        var XDeriv    : TNvector;
  239.                        var YDeriv    : TNvector;
  240.                            Tolerance : Float;
  241.                        var Error     : byte;
  242.                            FuncPtr   : Pointer);
  243.  
  244. {----------------------------------------------------------------------------}
  245. {-                                                                          -}
  246. {-    Input: NumDeriv, XDeriv, Tolerance                                    -}
  247. {-   Output: YDeriv, Error                                                  -}
  248. {-                                                                          -}
  249. {-  Purpose: Given a function TNTargetF(X), approximate the second          -}
  250. {-           derivative of TNTargetF(X) at the points XDeriv. The algorithm -}
  251. {-           uses a three point formula and Richardson extrapolation.       -}
  252. {-                                                                          -}
  253. {- User-defined function:  function TNTargetF(X : real) : real;             -}
  254. {-                                                                          -}
  255. {-  User-defined Types:  TNvector = array[1..TNArraySize] of real;          -}
  256. {-                                                                          -}
  257. {-    Global Variables:  NumDeriv  : integer;  Number points at which to    -}
  258. {-                                             take the derivative          -}
  259. {-                       XDeriv    : TNvector; Point at which to take       -}
  260. {-                                             the derivative               -}
  261. {-                       YDeriv    : TNvector; Derivative at                -}
  262. {-                                             points XDeriv                -}
  263. {-                       Tolerance : real;     Tolerance in answer          -}
  264. {-                       Error     : byte;     Flags if something goes      -}
  265. {-                                             wrong                        -}
  266. {-                                                                          -}
  267. {-              Errors:  0: No errors                                       -}
  268. {-                       1: Tolerance < TNNearlyZero                        -}
  269. {-                                                                          -}
  270. {----------------------------------------------------------------------------}
  271.  
  272. implementation
  273.  
  274. {$F+}
  275. {$L Differ.OBJ}
  276. function UserFunction(X : Float; ProcAddr : Pointer) : Float; external;
  277. {$F-}
  278.  
  279. procedure First_Derivative{(NumPoints : integer;
  280.                        var XData      : TNvector;
  281.                        var YData      : TNvector;
  282.                            Point      : byte;
  283.                            NumDeriv   : integer;
  284.                        var XDeriv     : TNvector;
  285.                        var YDeriv     : TNvector;
  286.                        var Error      : byte)};
  287.  
  288. {$IFOPT N+}
  289. const
  290.   TNNearlyZero = 1E-013;
  291. {$ELSE}
  292. const
  293.   TNNearlyZero = 1E-06;
  294. {$ENDIF}
  295.  
  296.   NotADataPoint = -9.9999999E35;    { See Note above }
  297.   NotEnoughPoints = 9.9999999E35;   { See Note above }
  298.  
  299.   procedure CheckData(NumPoints : integer;
  300.                   var XData     : TNvector;
  301.                       Point     : byte;
  302.                   var Error     : byte);
  303.  
  304. {---------------------------------------------------------}
  305. {- Input: NumPoints, XData, Point                        -}
  306. {- Output: Error                                         -}
  307. {-                                                       -}
  308. {- This procedure checks the data for errors:            -}
  309. {-                                                       -}
  310. {- Errors:  0: No errors                                 -}
  311. {-          2: Data not unique                           -}
  312. {-          3: Data not in increasing order              -}
  313. {-          4: Not enough data                           -}
  314. {-          5: Not(Point in [2, 3, 5])                   -}
  315. {-          6: Data points not evenly spaced for 5 point -}
  316. {---------------------------------------------------------}
  317.  
  318. var
  319.   Spacing, CurrentSpacing : Float;
  320.   Index : integer;
  321.  
  322. begin { procedure CheckData }
  323.   Error := 0;
  324.   if not(Point in [2, 3, 5]) then
  325.     Error := 5
  326.   else
  327.     if (NumPoints < Point) then
  328.       Error := 4;       { Not enough data }
  329.   Spacing := XData[2] - XData[1];
  330.   for Index := 2 to NumPoints do
  331.   begin
  332.     CurrentSpacing := XData[Index] - XData[Index - 1];
  333.     if CurrentSpacing < 0 then
  334.       Error := 3;  { Data not increasing }
  335.     if ABS(CurrentSpacing) < TNNearlyZero then
  336.       Error := 2;  { Data not unique }
  337.     if (Point = 5) and (ABS(Spacing - CurrentSpacing) > TNNearlyZero) then
  338.       Error := 6; { Data not evenly spaced }
  339.   end;
  340. end; { procedure CheckData }
  341.  
  342. procedure Differentiate(Point     : byte;
  343.                         NumPoints : integer;
  344.                     var XData     : TNvector;
  345.                     var YData     : TNvector;
  346.                         NumDeriv  : integer;
  347.                     var XDeriv    : TNvector;
  348.                     var YDeriv    : TNvector;
  349.                     var Error     : byte);
  350.  
  351. {------------------------------------------------------------}
  352. {- Input:  Point, NumPoints, XData, YData, NumDeriv, XDeriv -}
  353. {- Output: YDeriv, Error                                    -}
  354. {-                                                          -}
  355. {- This procedure approximates the derivative at the        -}
  356. {- XDeriv points using a 2, 3 or 5 point formula and stores -}
  357. {- them in YDeriv. If any of the XDeriv points are not in   -}
  358. {- XData then Error = 1 (a warning - not a fatal error) is  -}
  359. {- returned.                                                -}
  360. {------------------------------------------------------------}
  361.  
  362. var
  363.   Term, Index : integer;
  364.  
  365. procedure FindIndex(NumPoints : integer;
  366.                 var XData     : TNvector;
  367.                     XDeriv    : Float;
  368.                 var Index     : integer);
  369.  
  370. {------------------------------------------------------------}
  371. {- Input:  NumPoints, XData, XDeriv                         -}
  372. {- Output: Index                                            -}
  373. {-                                                          -}
  374. {- This procedure searches the vector XData for the value   -}
  375. {- XDeriv. The position of XDeriv withing XData is returned -}
  376. {- as Index.                                                -}
  377. {------------------------------------------------------------}
  378.  
  379. begin
  380.   Index := 1;
  381.   while (Index <= NumPoints) and (ABS(XData[Index] - XDeriv) > TNNearlyZero) do
  382.     Index := Succ(Index);
  383. end; { procedure FindIndex }
  384.  
  385. procedure TwoPoint(var XData  : TNvector;
  386.                    var YData  : TNvector;
  387.                        Index  : integer;
  388.                    var YDeriv : Float);
  389.  
  390. {----------------------------------------------------}
  391. {- Input:  XData, YData, Index                      -}
  392. {- Output: YDeriv                                   -}
  393. {-                                                  -}
  394. {- This procedure applies two point differentiation -}
  395. {- to the function represented by (XData, YData) at -}
  396. {- the point XData[Index].  The value of the        -}
  397. {- derivative is returned in YDeriv.                -}
  398. {- Due to the nature of the two point formula, a    -}
  399. {- different (though equivalent) formula is used    -}
  400. {- for the leftmost point than for the rest of the  -}
  401. {- points.                                          -}
  402. {----------------------------------------------------}
  403.  
  404. begin { procedure TwoPoint }
  405.   if Index = 1 then
  406.     YDeriv := (YData[Index + 1] - YData[Index]) /
  407.               (XData[Index + 1] - XData[Index])
  408.   else
  409.     YDeriv := (YData[Index - 1] - YData[Index]) /
  410.               (XData[Index - 1] - XData[Index]);
  411. end; { procedure TwoPoint }
  412.  
  413. procedure ThreePoint(NumPoints : integer;
  414.                  var XData     : TNvector;
  415.                  var YData     : TNvector;
  416.                      Index     : integer;
  417.                  var YDeriv    : Float);
  418.  
  419. {------------------------------------------------------}
  420. {- Input:  NumPoints, XData, YData, Index             -}
  421. {- Output: YDeriv                                     -}
  422. {-                                                    -}
  423. {- This procedure applies three point differentiation -}
  424. {- to the function represented by (XData, YData) at   -}
  425. {- the point XData[Index].  The value of the          -}
  426. {- derivative is returned in YDeriv.                  -}
  427. {- Due to the nature of the three point formula, a    -}
  428. {- different (less accurate) formula is used for the  -}
  429. {- endpoints than for the rest of the points.         -}
  430. {------------------------------------------------------}
  431.  
  432. var
  433.   Dif1, Dif2, Dif3 : Float;  { Differences }
  434.  
  435. begin { procedure ThreePoint }
  436.   if Index = 1 then   { Leftmost point }
  437.     begin
  438.       Dif1 := XData[1] - XData[2];
  439.       Dif2 := XData[1] - XData[3];
  440.       Dif3 := XData[2] - Xdata[3];
  441.       YDeriv := YData[Index] * (Dif1 + Dif2) / (Dif1 * Dif2)
  442.                 - YData[Index + 1] * Dif2 / (Dif1 * Dif3)
  443.                 + YData[Index + 2] * Dif1 / (Dif2 * Dif3);
  444.     end
  445.   else
  446.     if Index = NumPoints then { Rightmost point }
  447.       begin
  448.         Dif1 := XData[NumPoints - 2] - XData[NumPoints - 1];
  449.         Dif2 := XData[NumPoints - 2] - XData[NumPoints];
  450.         Dif3 := XData[NumPoints - 1] - Xdata[NumPoints];
  451.         YDeriv := - YData[Index - 2] * Dif3 / (Dif1 * Dif2)
  452.                   + YData[Index - 1] * Dif2 / (Dif1 * Dif3)
  453.                   - YData[Index] * (Dif2 + Dif3) / (Dif2 * Dif3);
  454.       end
  455.     else { Middle points }
  456.       begin
  457.         Dif1 := XData[Index - 1] - XData[Index];
  458.         Dif2 := XData[Index - 1] - XData[Index + 1];
  459.         Dif3 := XData[Index] - XData[Index + 1];
  460.         YDeriv := YData[Index - 1] * Dif3 / (Dif1 * Dif2)
  461.                   - YData[Index] * (Dif3 - Dif1) / (Dif1 * Dif3)
  462.                   - YData[Index + 1] * Dif1 / (Dif2 * Dif3);
  463.       end;
  464. end; { procedure ThreePoint }
  465.  
  466. procedure FivePoint(NumPoints : integer;
  467.                 var XData     : TNvector;
  468.                 var YData     : TNvector;
  469.                     Index     : integer;
  470.                 var YDeriv    : Float;
  471.                 var Error     : byte);
  472.  
  473. {------------------------------------------------------}
  474. {- Input:  NumPoints, XData, YData, Index             -}
  475. {- Output: YDeriv, Error                              -}
  476. {-                                                    -}
  477. {- This procedure applies five point differentiation  -}
  478. {- to the function represented by (XData, YData) at   -}
  479. {- the point XData[Index].  The value of the          -}
  480. {- derivative is returned in YDeriv.                  -}
  481. {- Due to the nature of the five point formula,       -}
  482. {- different (less accurate) formulas are used for    -}
  483. {- the two leftmost and two rightmost points than for -}
  484. {- the rest of the points.                            -}
  485. {- If there are only 5 data points, the only point    -}
  486. {- at which a 5 point formula may be applied is point -}
  487. {- 3.  If there are only 6 data points, the only      -}
  488. {- points at which a 5 point formula may be applied   -}
  489. {- are points 1, 3, 4 and 6.  If these conditions are -}
  490. {- not met, then Error = 1 (a warning - not a fatal   -}
  491. {- error is returned.  The value of the derivative at -}
  492. {- such an "error point" is arbitrarily assigned to   -}
  493. {- be NotEnoughPoints.                                -}
  494. {------------------------------------------------------}
  495.  
  496. var
  497.   Spacing : Float;
  498.  
  499. begin  { procedure FivePoint }
  500.   Spacing := XData[2] - XData[1];
  501.   if ((NumPoints = 6) and (Index in [2, 5])) or
  502.      ((NumPoints = 5) and not(Index = 3))    then
  503.     begin
  504.       YDeriv := NotEnoughPoints;
  505.       Error := 1;    { Warning about this point }
  506.     end
  507.   else
  508.     if Index <= 2 then { Leftmost points }
  509.       YDeriv := (-25 * YData[Index] + 48 * YData[Index + 1]
  510.                  - 36 * YData[Index + 2] + 16 * YData[Index + 3]
  511.                  - 3 * YData[Index + 4]) / (12 * Spacing)
  512.     else
  513.       if Index >= NumPoints-1  then { Rightmost points }
  514.         YDeriv := (-25*YData[Index] + 48 * YData[Index - 1]
  515.                    - 36 * YData[Index - 2] + 16 * YData[Index - 3]
  516.                    - 3 * YData[Index - 4]) / (-12 * Spacing)
  517.       else { Middle points }
  518.          YDeriv := (YData[Index - 2] - 8 * YData[Index - 1]
  519.                    + 8 * YData[Index + 1] - YData[Index + 2]) / (12 * Spacing);
  520. end; { procedure FivePoint }
  521.  
  522. begin { procedure Differentiate }
  523.   for Term := 1 to NumDeriv do
  524.   begin
  525.     FindIndex(NumPoints, XData, XDeriv[Term], Index);
  526.     if Index > NumPoints then
  527.       begin
  528.         YDeriv[Term] := NotADataPoint; { YDeriv[Term] isn't }
  529.                                        { in XData           }
  530.          Error := 1;                   { Warning about this point }
  531.        end
  532.      else
  533.        case Point of
  534.          2 : TwoPoint(XData, YData, Index, YDeriv[Term]);
  535.          3 : ThreePoint(NumPoints, XData, YData, Index, YDeriv[Term]);
  536.          5 : FivePoint(NumPoints, XData, YData, Index, YDeriv[Term], Error);
  537.        end;
  538.   end;
  539. end; { procedure Differentiate }
  540.  
  541. begin { procedure First_Derivative }
  542.   CheckData(NumPoints, XData, Point, Error);
  543.   if Error = 0 then
  544.     Differentiate(Point, NumPoints, XData, YData,
  545.                   NumDeriv, XDeriv, YDeriv, Error);
  546. end; { procedure First_Derivative }
  547.  
  548. procedure Second_Derivative{(NumPoints : integer;
  549.                         var XData      : TNvector;
  550.                         var YData      : TNvector;
  551.                             Point      : byte;
  552.                             NumDeriv   : integer;
  553.                         var XDeriv     : TNvector;
  554.                         var YDeriv     : TNvector;
  555.                         var Error      : byte)};
  556.  
  557. {$IFOPT N+}
  558. const
  559.   NearlyZero = 1E-013;
  560. {$ELSE}
  561. const
  562.   NearlyZero = 1E-06;
  563. {$ENDIF}
  564.  
  565.   NotADataPoint = -9.9999999E35;    { See Note above }
  566.   NotEnoughPoints = 9.9999999E35;   { See Note above }
  567.  
  568. var
  569.   SpacingSquared : Float;   { Square of the spacing between data points }
  570.  
  571. procedure CheckData(NumPoints      : integer;
  572.                 var XData          : TNvector;
  573.                     Point          : byte;
  574.                 var SpacingSquared : Float;
  575.                 var Error          : byte);
  576.  
  577. {---------------------------------------------------------}
  578. {- Input:  NumPoints, XData, Point                       -}
  579. {- Output: Error                                         -}
  580. {-                                                       -}
  581. {- This procedure checks the data for errors             -}
  582. {-                                                       -}
  583. {- Errors:  0: No errors                                 -}
  584. {-          2: Data not unique                           -}
  585. {-          3: Data not in increasing order              -}
  586. {-          4: Not enough data                           -}
  587. {-          5: not(Point IN [2, 3, 5])                   -}
  588. {-          6: Data points not evenly spaced for 5 point -}
  589. {---------------------------------------------------------}
  590.  
  591. var
  592.   Spacing, CurrentSpacing : Float;
  593.   Index : integer;
  594.  
  595. begin { procedure CheckData }
  596.   if not(Point in [3, 5]) then
  597.     Error := 5      { Point <> 3 or 5 }
  598.   else
  599.     if (NumPoints < Point)  then
  600.       Error := 4;   { Too few points }
  601.   Spacing := XData[2] - XData[1];
  602.   for Index := 2 to NumPoints do
  603.   begin
  604.     CurrentSpacing := XData[Index] - XData[Index - 1];
  605.     if ABS(CurrentSpacing) < NearlyZero then
  606.       Error := 2;  { Data not unique }
  607.     if CurrentSpacing < 0 then
  608.       Error := 3;  { Data not increasing }
  609.     if ABS(Spacing - CurrentSpacing) > NearlyZero then
  610.       Error := 6;  { Data not evenly spaced }
  611.   end;
  612.   SpacingSquared := Sqr(Spacing);
  613. end; { procedure CheckData }
  614.  
  615. procedure Differentiate(Point          : byte;
  616.                         NumPoints      : integer;
  617.                     var XData          : TNvector;
  618.                     var YData          : TNvector;
  619.                         SpacingSquared : Float;
  620.                         NumDeriv       : integer;
  621.                     var XDeriv         : TNvector;
  622.                     var YDeriv         : TNvector;
  623.                     var Error          : byte);
  624.  
  625. {------------------------------------------------------------}
  626. {- Input: Point, NumPoints, XData, YData, SpacingSquared,   -}
  627. {-        NumDeriv, XDeriv                                  -}
  628. {- Output: YDeriv, Error                                    -}
  629. {-                                                          -}
  630. {- This procedure approximates the second derivative of the -}
  631. {- function represented by (XData, YData) at the XDeriv     -}
  632. {- points using a 3 or 5 point formula and stores them in   -}
  633. {- YDeriv. If any of the XDeriv points are not in XData,    -}
  634. {- then Error = 1 (a warning - not a fatal error) is        -}
  635. {- returned. The value of the second derivative at one of   -}
  636. {- these "error points" is arbitrarily assigned to be       -}
  637. {- NotADataPoint.                                           -}
  638. {------------------------------------------------------------}
  639.  
  640. var
  641.   Term, Index : integer;
  642.  
  643. procedure FindIndex(NumPoints : integer;
  644.                 var XData     : TNvector;
  645.                 var XDeriv    : Float;
  646.                 var Index     : integer);
  647.  
  648. {------------------------------------------------------------}
  649. {- Input: NumPoints, XData, XDeriv                          -}
  650. {- Output: Index                                            -}
  651. {-                                                          -}
  652. {- This procedure searches the vector XData for the value   -}
  653. {- XDeriv. The position of XDeriv withing XData is returned -}
  654. {- as Index.                                                -}
  655. {------------------------------------------------------------}
  656.  
  657. begin
  658.   Index := 1;
  659.   while (Index <= NumPoints) and (XData[Index] <> XDeriv) do
  660.     Index := Succ(Index);
  661. end; { procedure FindIndex }
  662.  
  663. procedure ThreePoint(NumPoints      : integer;
  664.                  var YData          : TNvector;
  665.                      Index          : integer;
  666.                      SpacingSquared : Float;
  667.                  var YDeriv         : Float);
  668.  
  669. {------------------------------------------------------}
  670. {- Input: NumPoints, YData, Index, SpacingSquared     -}
  671. {- Output: YDeriv                                     -}
  672. {-                                                    -}
  673. {- This procedure applies three point second          -}
  674. {- differentiation to the function represented by     -}
  675. {- (XData, YData) at the point XData[Index].  The     -}
  676. {- value of the derivative is returned in YDeriv.     -}
  677. {- Due to the nature of the three point formula, a    -}
  678. {- different (less accurate) formula is used for the  -}
  679. {- endpoints than for the rest of the points.         -}
  680. {------------------------------------------------------}
  681.  
  682. begin { procedure ThreePoint }
  683.   if Index = 1 then  { Leftmost point }
  684.     YDeriv := (YData[1] - 2 * YData[2] + YData[3]) / SpacingSquared
  685.     { The error in this Term is O(Spacing) }
  686.   else
  687.     if Index = NumPoints then { Rightmost point }
  688.       YDeriv := (YData[NumPoints] - 2 * YData[NumPoints - 1]
  689.                 + YData[NumPoints - 2]) / SpacingSquared
  690.       { The error in this Term is O(Spacing) }
  691.     else  { Middle points }
  692.       YDeriv := (YData[Index - 1] - 2 * YData[Index]
  693.                 + YData[Index + 1]) / SpacingSquared;
  694.       { The error in this Term is O(SpacingSquared) }
  695. end; { procedure ThreePoint }
  696.  
  697. procedure FivePoint(NumPoints      : integer;
  698.                 var YData          : TNvector;
  699.                     Index          : integer;
  700.                     SpacingSquared : Float;
  701.                 var YDeriv         : Float;
  702.                 var Error          : byte);
  703.  
  704. {------------------------------------------------------}
  705. {- Input: NumPoints, YData, Index                     -}
  706. {- Output: YDeriv, Error                              -}
  707. {-                                                    -}
  708. {- This procedure applies five point differentiation  -}
  709. {- to the function represented by (XData, YData) at   -}
  710. {- the point XData[Index].  The value of the          -}
  711. {- derivative is returned in YDeriv.                  -}
  712. {- Due to the nature of the five point formula,       -}
  713. {- different (less accurate) formulas are used for    -}
  714. {- the two leftmost and two rightmost points than for -}
  715. {- the rest of the points.                            -}
  716. {- If there are only 5 data points, the only points   -}
  717. {- at which a 5 point formula may be applied are      -}
  718. {- points 1, 3 or 5. If this condition is not met,    -}
  719. {- then Error = 1 (a warning - not a fatal error is   -}
  720. {- returned.  The value of the derivative at such an  -}
  721. {- "error point" is arbitrarily assigned to be        -}
  722. {- NotEnoughPoints.                                   -}
  723. {------------------------------------------------------}
  724.  
  725. begin { procedure FivePoint }
  726.   if (NumPoints = 5) and (Index in [2, 4]) then
  727.     begin
  728.       YDeriv := NotEnoughPoints;    { Too few points to       }
  729.                                     { compute derivative here }
  730.       Error := 1;     { Warning about this point }
  731.     end
  732.   else
  733.     if Index <= 2 then { Leftmost points }
  734.       YDeriv := (2 * YData[Index] - 5 * YData[Index + 1]
  735.                 + 4 * YData[Index + 2] - YData[Index + 3]) / SpacingSquared
  736.                 { The error in this Term is O(SpacingSquared) }
  737.     else
  738.       if Index >= NumPoints - 1 then { Rightmost points }
  739.         YDeriv := (2 * YData[Index] - 5 * YData[Index - 1]
  740.                   + 4 * YData[Index - 2] - YData[Index - 3]) / SpacingSquared
  741.                   { The error in this Term is O(SpacingSquared) }
  742.       else { Middle points }
  743.         YDeriv := (-YData[Index - 2] + 16 * YData[Index - 1]
  744.         - 30 * YData[Index] + 16 * YData[Index + 1]
  745.         - YData[Index + 2]) / (12 * SpacingSquared);
  746.         { The error in this Term is O(SQR(SpacingSquared)) }
  747. end; { procedure FivePoint }
  748.  
  749. begin  { procedure Differentiate }
  750.   for Term := 1 to NumDeriv do
  751.   begin
  752.     FindIndex(NumPoints, XData, XDeriv[Term], Index);
  753.     if Index > NumPoints then
  754.       begin
  755.         YDeriv[Term] := NotADataPoint;    { YDeriv[Term] isn't in XData }
  756.         Error := 1;   { Warning about this point }
  757.       end
  758.     else
  759.       if Point = 3 then
  760.         ThreePoint(NumPoints, YData, Index, SpacingSquared, YDeriv[Term])
  761.       else
  762.         FivePoint(NumPoints, YData, Index, SpacingSquared, YDeriv[Term], Error);
  763.   end;
  764. end; { procedure Differentiate }
  765.  
  766. begin  { procedure Second_Derivative }
  767.   CheckData(NumPoints, XData, Point, SpacingSquared, Error);
  768.   if Error = 0 then
  769.     Differentiate(Point, NumPoints, XData, YData, SpacingSquared,
  770.                   NumDeriv, XDeriv, YDeriv, Error);
  771. end; { procedure Second_Derivative }
  772.  
  773. procedure Interpolate_Derivative{(NumPoints : integer;
  774.                              var XData      : TNvector;
  775.                              var YData      : TNvector;
  776.                                  NumDeriv   : integer;
  777.                              var XDeriv     : TNvector;
  778.                              var YInter     : TNvector;
  779.                              var YDeriv     : TNvector;
  780.                              var YDeriv2    : TNvector;
  781.                              var Error      : byte)};
  782.  
  783. {$IFOPT N+}
  784. const
  785.   TNNearlyZero = 1E-015;
  786. {$ELSE}
  787. const
  788.   TNNearlyZero = 1E-08;
  789. {$ENDIF}
  790.  
  791. type
  792.   TNcoefficients = record
  793.                      A, B, C, D : TNvector;
  794.                    end;
  795. var
  796.   Spline : TNcoefficients;     { Cubics defining the interpolated function }
  797.  
  798. {------------------------------------------------------------------------}
  799. procedure CubicSplineFree(NumPoints : integer;
  800.                       var XData     : TNvector;
  801.                       var YData     : TNvector;
  802.                       var CoefA     : TNvector;
  803.                       var CoefB     : TNvector;
  804.                       var CoefC     : TNvector;
  805.                       var CoefD     : TNvector;
  806.                           NumInter  : integer;
  807.                       var XInter    : TNvector;
  808.                       var YInter    : TNvector;
  809.                       var Error     : byte);
  810.  
  811. {----------------------------------------------------------------------------}
  812. {-                                                                          -}
  813. {-    Input: NumPoints, XData, YData, NumInter, XInter                      -}
  814. {-    Output: CoefA, CoefB, CoefC, CoefD, YInter, Error                     -}
  815. {-                                                                          -}
  816. {-            Purpose: To construct a cubicspline interpolant               -}
  817. {-                     to a set of points. The second derivative            -}
  818. {-                     of the interpolant is assumed to be zero             -}
  819. {-                     at the endpoints (free cubic spline)                 -}
  820. {-                     The spline is of the form                            -}
  821. {-                                                                          -}
  822. {-                        3               2                                 -}
  823. {-            D[I](X-X[I])  + C[I](X-X[I])  + B[I](X-X[I]) + A[I]           -}
  824. {-                                                                          -}
  825. {-                     Where 1 < I < NumPoints.                             -}
  826. {-                                                                          -}
  827. {-  User-defined Types: TNvector = array[0..TNArraySize] of real;           -}
  828. {-                                                                          -}
  829. {- Global Variables: NumPoints : integer;  Number of points                 -}
  830. {-                   XData     : TNvector; X-coordinate data points         -}
  831. {-                                         (must be in increasing order)    -}
  832. {-                   YData     : TNvector; Y-coordinate data points         -}
  833. {-                   CoefA     : TNvector; 0th coefficient of spline        -}
  834. {-                   CoefB     : TNvector; 1st coefficient of spline        -}
  835. {-                   CoefC     : TNvector; 2nd coefficient of spline        -}
  836. {-                   CoefD     : TNvector; 3rd coefficient of spline        -}
  837. {-                   NumInter  : integer;  Number of interpolated points    -}
  838. {-                   XInter    : TNvector; Points at which to interpolate   -}
  839. {-                   YInter    : TNvector; Interpolated values at XInter    -}
  840. {-                   Error     : byte;     Flags if something goes wrong    -}
  841. {-                                                                          -}
  842. {-              Errors: 0: No error                                         -}
  843. {-                      1: X-coordinate points not                          -}
  844. {-                         unique                                           -}
  845. {-                      2: X-coordinate points not                          -}
  846. {-                         in increasing order                              -}
  847. {-                      3: NumPoints < 2                                    -}
  848. {-                                                                          -}
  849. {----------------------------------------------------------------------------}
  850.  
  851. var
  852.   Interval : TNvector;        { Intervals between adjacent points }
  853.   Spline : TNcoefficients;    { All the cubics }
  854.  
  855. procedure CalculateIntervals(NumPoints : integer;
  856.                          var XData     : TNvector;
  857.                          var Interval  : TNvector;
  858.                          var Error     : byte);
  859.  
  860. {----------------------------------------------------------}
  861. {- Input:  NumPoints, XData                               -}
  862. {- Output: Interval, Error                                -}
  863. {-                                                        -}
  864. {- This procedure calculates the length of the interval   -}
  865. {- between two adjacent X values, contained in XData. If  -}
  866. {- the X values are not sequential, Error = 2 is returned -}
  867. {- and if the X values are not unique, then Error = 1 is  -}
  868. {- returned.                                              -}
  869. {----------------------------------------------------------}
  870.  
  871. var
  872.   Index : integer;
  873.  
  874. begin
  875.   Error := 0;
  876.   for Index := 1 to NumPoints - 1 do
  877.   begin
  878.     Interval[Index] := XData[Index+1] - XData[Index];
  879.     if ABS(Interval[Index]) < TNNearlyZero then
  880.       Error := 1;     { Data not unique }
  881.     if Interval[Index] < 0 then
  882.       Error := 2;     { Data not in increasing order }
  883.   end;
  884. end; { procedure CalculateIntervals }
  885.  
  886. procedure CalculateCoefficients(NumPoints : integer;
  887.                             var XData     : TNvector;
  888.                             var YData     : TNvector;
  889.                             var Interval  : TNvector;
  890.                             var Spline    : TNcoefficients);
  891.  
  892. {---------------------------------------------------------------}
  893. {- Input:  NumPoints, XData, YData, Interval                   -}
  894. {- Output: Spline                                              -}
  895. {-                                                             -}
  896. {- This procedure calculates the coefficients of each cubic    -}
  897. {- in the interpolating spline. A separate cubic is calculated -}
  898. {- for every interval between data points.  The coefficients   -}
  899. {- are returned in the variable Spline.                        -}
  900. {---------------------------------------------------------------}
  901.  
  902. procedure CalculateAs(NumPoints : integer;
  903.                   var YData     : TNvector;
  904.                   var Spline    : TNcoefficients);
  905.  
  906. {------------------------------------------}
  907. {- Input: NumPoints, YData                -}
  908. {- Ouput: Spline                          -}
  909. {-                                        -}
  910. {- This procedure calculates the constant -}
  911. {- term in each cubic.                    -}
  912. {------------------------------------------}
  913.  
  914. var
  915.   Index : integer;
  916.  
  917. begin
  918.   for Index := 1 to NumPoints do
  919.     Spline.A[Index] := YData[Index];
  920. end; { procedure CalculateAs }
  921.  
  922.   procedure CalculateCs(NumPoints : integer;
  923.                     var XData     : TNvector;
  924.                     var Interval  : TNvector;
  925.                     var Spline    : TNcoefficients);
  926.  
  927. {---------------------------------------------}
  928. {- Input: NumPoints, XData, Interval         -}
  929. {- Ouput: Spline                             -}
  930. {-                                           -}
  931. {- This procedure calculates the coefficient -}
  932. {- of the squared term in each cubic.        -}
  933. {---------------------------------------------}
  934.  
  935. var
  936.   Alpha, L, Mu, Z : TNvector;
  937.   Index : integer;
  938.  
  939. begin
  940.   with Spline do
  941.   begin
  942.     { The next few lines solve a tridiagonal matrix }
  943.     for Index := 2 to NumPoints - 1 do
  944.       Alpha[Index] := 3 * ((A[Index+1] * Interval[Index-1])
  945.                       - (A[Index] * (XData[Index+1] - XData[Index-1]))
  946.                       + (A[Index-1] * Interval[Index]))
  947.                       / (Interval[Index-1] * Interval[Index]);
  948.     L[1] := 0;
  949.     Mu[1] := 0;
  950.     Z[1] := 0;
  951.     for Index := 2 to NumPoints - 1 do
  952.     begin
  953.       L[Index] := 2 * (XData[Index+1] - XData[Index-1]) -
  954.                   Interval[Index-1] * Mu[Index-1];
  955.       Mu[Index] := Interval[Index] / L[Index];
  956.       Z[Index] := (Alpha[Index] - Interval[Index-1] * Z[Index-1]) / L[Index];
  957.     end;
  958.  
  959.     { Now calculate the C's }
  960.     C[NumPoints] := 0;
  961.     for Index := NumPoints - 1 downto 1 do
  962.       C[Index] := Z[Index] - Mu[Index] * C[Index+1];
  963.   end; { with }
  964. end; { procedure CalculateCs }
  965.  
  966. procedure CalculateBandDs(NumPoints : integer;
  967.                       var Interval  : TNvector;
  968.                       var Spline    : TNcoefficients);
  969.  
  970. {------------------------------------------------}
  971. {- Input: NumPoints, Interval                   -}
  972. {- Ouput: Spline                                -}
  973. {-                                              -}
  974. {- This procedure calculates the coefficient of -}
  975. {- the linear and cubic terms in each cubic.    -}
  976. {------------------------------------------------}
  977.  
  978. var
  979.   Index : integer;
  980.  
  981. begin
  982.   with Spline do
  983.     for Index := NumPoints - 1 downto 1 do
  984.     begin
  985.       B[Index] := (A[Index+1] - A[Index]) / Interval[Index]
  986.                   - Interval[Index] * (C[Index+1] + 2 * C[Index]) / 3;
  987.       D[Index] := (C[Index+1] - C[Index]) / (3 * Interval[Index]);
  988.     end;
  989. end; { procedure CalculateDs }
  990.  
  991. begin { procedure CalculateCoefficients }
  992.   CalculateAs(NumPoints, YData, Spline);
  993.   CalculateCs(NumPoints, XData, Interval, Spline);
  994.   CalculateBandDs(NumPoints, Interval, Spline);
  995. end; { procedure CalculateCoefficients }
  996.  
  997. procedure Interpolate(NumPoints : integer;
  998.                   var XData     : TNvector;
  999.                   var Spline    : TNcoefficients;
  1000.                       NumInter  : integer;
  1001.                   var XInter    : TNvector;
  1002.                   var YInter    : TNvector);
  1003.  
  1004. {---------------------------------------------------------------}
  1005. {- Input:  NumPoints, XData, Spline, NumInter, XInter          -}
  1006. {- Output: YInter                                              -}
  1007. {-                                                             -}
  1008. {- This procedure uses the interpolating cubic spline (Spline) -}
  1009. {- to interpolate values for the X positions contained         -}
  1010. {- in XInter.  The interpolated values are returned in YInter. -}
  1011. {---------------------------------------------------------------}
  1012.  
  1013. var
  1014.   Index, Location, Term : integer;
  1015.   X : Float;
  1016.  
  1017. begin
  1018.   for Index := 1 to NumInter do
  1019.   begin
  1020.     Location := 1;
  1021.     for Term := 1 to NumPoints - 1 do
  1022.       if XInter[Index] > XData[Term] then
  1023.         Location := Term;
  1024.     X := XInter[Index] - XData[Location];
  1025.     with Spline do
  1026.       YInter[Index] := ((D[Location] * X + C[Location]) * X +
  1027.                          B[Location]) * X + A[Location];
  1028.   end;
  1029. end; { procedure Interpolate }
  1030.  
  1031. begin { procedure CubicSplineFree }
  1032.   Error := 0;
  1033.   if NumPoints < 2 then
  1034.     Error := 3
  1035.   else
  1036.     CalculateIntervals(NumPoints, XData, Interval, Error);
  1037.   if Error = 0 then
  1038.   begin
  1039.     CalculateCoefficients(NumPoints, XData, YData, Interval, Spline);
  1040.     Interpolate(NumPoints, XData, Spline, NumInter, XInter, YInter);
  1041.   end;
  1042.   CoefA := Spline.A;
  1043.   CoefB := Spline.B;
  1044.   CoefC := Spline.C;
  1045.   CoefD := Spline.D;
  1046. end;  { procedure CubicSplineFree }
  1047.  
  1048. {------------------------------------------------------------------------}
  1049.  
  1050. procedure Differentiate(NumPoints : integer;
  1051.                     var XData     : TNvector;
  1052.                     var Spline    : TNcoefficients;
  1053.                         NumDeriv  : integer;
  1054.                     var XDeriv    : TNvector;
  1055.                     var YDeriv    : TNvector;
  1056.                     var YDeriv2   : TNvector);
  1057.  
  1058. {-----------------------------------------------------}
  1059. {- Input: NumPoints, XData, Spline, NumDeriv, XDeriv -}
  1060. {- Output: YDeriv, YDeriv2                           -}
  1061. {-                                                   -}
  1062. {-                                                   -}
  1063. {- This procedure uses the cubic spline interpolant, -}
  1064. {- Spline, to approximate the first and second       -}
  1065. {- derivatives.  The derivatives are returned in     -}
  1066. {- YDeriv, YDeriv2.                                  -}
  1067. {-----------------------------------------------------}
  1068.  
  1069. var
  1070.   Index, Location, Term : integer;
  1071.   X : Float;
  1072.  
  1073. begin
  1074.   for Index := 1 to NumDeriv do
  1075.   begin
  1076.     Location := 1;
  1077.     for Term := 1 to NumPoints-1 do
  1078.       if XDeriv[Index] > XData[Term] then
  1079.         Location := Term;
  1080.     X := XDeriv[Index] - XData[Location];
  1081.  
  1082.     with Spline do    { Interpolate first derivative }
  1083.       YDeriv[Index] := (3 * D[Location] * X + 2 * C[Location]) * X
  1084.                        + B[Location];
  1085.     with Spline do  { Interpolate second derivative }
  1086.       YDeriv2[Index] := 6 * D[Location] * X + 2 * C[Location];
  1087.   end;
  1088. end; { procedure Differentiate }
  1089.  
  1090. begin { procedure Interpolate_Derivative }
  1091.   Error := 0;
  1092.   if NumPoints < 2 then
  1093.     Error := 3
  1094.   else
  1095.     CubicSplineFree(NumPoints, XData, YData, Spline.A, Spline.B, Spline.C,
  1096.                     Spline.D, NumDeriv, XDeriv, YInter, Error);
  1097.   if Error = 0 then
  1098.     Differentiate(NumPoints, XData, Spline, NumDeriv, XDeriv, YDeriv, YDeriv2);
  1099. end; { procedure Interpolate_Derivative }
  1100.  
  1101. procedure FirstDerivative{(NumDeriv : integer;
  1102.                      var XDeriv     : TNvector;
  1103.                      var YDeriv     : TNvector;
  1104.                          Tolerance  : Float;
  1105.                      var Error      : byte;
  1106.                          FuncPtr    : Pointer)};
  1107.  
  1108. {$IFOPT N+}
  1109. const
  1110.   TNNearlyZero = 1E-010;
  1111. {$ELSE}
  1112. const
  1113.   TNNearlyZero = 1E-05;
  1114. {$ENDIF}
  1115.  
  1116. procedure Differentiate(Tolerance : Float;
  1117.                         NumDeriv  : integer;
  1118.                     var XDeriv    : TNvector;
  1119.                     var YDeriv    : TNvector);
  1120.  
  1121. {------------------------------------------------------------}
  1122. {- Input: Tolerance, NumDeriv, XDeriv                       -}
  1123. {- Output: YDeriv                                           -}
  1124. {-                                                          -}
  1125. {- This procedure uses a three point formula and Richardson -}
  1126. {- extrapolation to approximate the value of the derivative -}
  1127. {- of the function TNTargetF at the points XDeriv.  The     -}
  1128. {- derivatives are returned in YDeriv.                      -}
  1129. {------------------------------------------------------------}
  1130.  
  1131. var
  1132.   Term, Iter : integer;
  1133.   TwoToTheIterMinus2 : integer;
  1134.   X, DeltaX : Float;
  1135.   OldEstimate, NewEstimate : TNvector;
  1136.  
  1137. function EvaluateFirstDeriv(X      : Float;
  1138.                             DeltaX : Float) : Float;
  1139.  
  1140. {------------------------------------------------------}
  1141. {- Input: X, DeltaX                                   -}
  1142. {- Output: EvaluateFirstDeriv                         -}
  1143. {-                                                    -}
  1144. {- This procedure uses a three point formula to       -}
  1145. {- evaluate the derivative of the function TNTargetF  -}
  1146. {- at the point X.  The function returns the value of -}
  1147. {- the derivative.                                    -}
  1148. {------------------------------------------------------}
  1149.  
  1150. var
  1151.   LeftPoint, RightPoint : Float;   { Points on either side of X }
  1152.  
  1153. begin  { function EvaluateFirstDeriv }
  1154.   LeftPoint := UserFunction(X - DeltaX, FuncPtr);
  1155.   RightPoint := UserFunction(X + DeltaX, FuncPtr);
  1156.   EvaluateFirstDeriv := (RightPoint - LeftPoint)/(2 * DeltaX);
  1157. end; { function EvaluateFirstDeriv }
  1158.  
  1159. procedure Extrapolate(Iter        : integer;
  1160.                   var OldEstimate : TNvector;
  1161.                   var NewEstimate : TNvector);
  1162.  
  1163. {--------------------------------------------------------}
  1164. {- Input:  Iter, OldEstimate                            -}
  1165. {- Output: NewEstimate                                  -}
  1166. {-                                                      -}
  1167. {- This procedure uses Richardson extrapolation         -}
  1168. {- to improve the current approximation to the          -}
  1169. {- derivative (OldEstimate).  The result is returned    -}
  1170. {- as the Iter element in the array NewEstimate         -}
  1171. {--------------------------------------------------------}
  1172.  
  1173. var
  1174.   Extrap : integer;
  1175.   FourToTheExtrapMinus1 : Float;
  1176.  
  1177. begin
  1178.   FourToTheExtrapMinus1 := 1;
  1179.   for Extrap := 2 to Iter do
  1180.   begin
  1181.     FourToTheExtrapMinus1 := FourToTheExtrapMinus1 * 4;
  1182.     NewEstimate[Extrap] := (FourToTheExtrapMinus1 * NewEstimate[Extrap - 1] -
  1183.                         OldEstimate[Extrap - 1]) / (FourToTheExtrapMinus1 - 1);
  1184.   end;
  1185. end; { procedure Extrapolate }
  1186.  
  1187. begin { procedure Differentiate }
  1188.   for Term := 1 to NumDeriv do
  1189.   begin
  1190.     X := XDeriv[Term];
  1191.     if ABS(X) < TNNearlyZero then
  1192.       DeltaX := Sqrt(Tolerance)
  1193.     else
  1194.       DeltaX := ABS(X * Sqrt(Tolerance));
  1195.     OldEstimate[1] := EvaluateFirstDeriv(X, DeltaX);
  1196.     Iter := 1;
  1197.     TwoToTheIterMinus2 := 1;
  1198.     repeat
  1199.       Iter := Succ(Iter);
  1200.       DeltaX := DeltaX / 2;
  1201.       NewEstimate[1] := EvaluateFirstDeriv(X, DeltaX);
  1202.       TwoToTheIterMinus2 := TwoToTheIterMinus2 * 2;
  1203.       Extrapolate(Iter, OldEstimate, NewEstimate);
  1204.       OldEstimate := NewEstimate;
  1205.     until { The fractional difference between }
  1206.           { iterations is within Tolerance    }
  1207.           (ABS(NewEstimate[Iter - 1] - NewEstimate[Iter]) <=
  1208.            ABS(Tolerance * NewEstimate[Iter])) or
  1209.           (ABS(DeltaX) < TNNearlyZero);  { Prevents infinite loops }
  1210.     YDeriv[Term] := NewEstimate[Iter];
  1211.   end;
  1212. end; { procedure Differentiate }
  1213.  
  1214. begin  { procedure FirstDerivative }
  1215.   if Tolerance < TNNearlyZero then
  1216.     Error := 1
  1217.   else
  1218.     Differentiate(Tolerance, NumDeriv, XDeriv, YDeriv);
  1219. end; { procedure FirstDerivative }
  1220.  
  1221. procedure SecondDerivative{(NumDeriv : integer;
  1222.                        var XDeriv    : TNvector;
  1223.                        var YDeriv    : TNvector;
  1224.                            Tolerance : Float;
  1225.                        var Error     : byte;
  1226.                            FuncPtr   : Pointer)};
  1227.  
  1228. {$IFOPT N+}
  1229. const
  1230.   TNNearlyZero = 1E-004;
  1231. {$ELSE}
  1232. const
  1233.   TNNearlyZero = 1E-02;
  1234. {$ENDIF}
  1235.  
  1236. procedure Differentiate(Tolerance : Float;
  1237.                         NumDeriv  : integer;
  1238.                     var XDeriv    : TNvector;
  1239.                     var YDeriv    : TNvector);
  1240.  
  1241. var
  1242.   Term, Iter : integer;
  1243.   TwoToTheIterMinus2 : integer;
  1244.   X, DeltaX : Float;
  1245.   OldEstimate, NewEstimate : TNvector;
  1246.  
  1247.  
  1248. function EvaluateSecondDeriv(X      : Float;
  1249.                              DeltaX : Float) : Float;
  1250.  
  1251. {------------------------------------------------------}
  1252. {- Input: X, DeltaX                                   -}
  1253. {- Output: EvaluateSecondDeriv                        -}
  1254. {-                                                    -}
  1255. {- This procedure uses a three point formula to       -}
  1256. {- evaluate the second derivative of the function     -}
  1257. {- TNTargetF at the point X.  The function returns    -}
  1258. {- the value of the second derivative.                -}
  1259. {------------------------------------------------------}
  1260.  
  1261. var
  1262.   Point, LeftPoint, RightPoint : Float;   { Value of function at X  }
  1263.                                          { and on either side of X }
  1264.  
  1265. begin { function EvaluateSecondDeriv }
  1266.   LeftPoint := UserFunction(X - DeltaX, FuncPtr);
  1267.   RightPoint := UserFunction(X + DeltaX, FuncPtr);
  1268.   Point := UserFunction(X, FuncPtr);
  1269.   EvaluateSecondDeriv := (RightPoint - 2 * Point + LeftPoint) / Sqr(DeltaX);
  1270. end; { function EvaluateSecondDeriv }
  1271.  
  1272. procedure Extrapolate(Iter        : integer;
  1273.                   var OldEstimate : TNvector;
  1274.                   var NewEstimate : TNvector);
  1275.  
  1276. {--------------------------------------------------------}
  1277. {- Input: Iter, OldEstimate                             -}
  1278. {- Output: NewEstimate                                  -}
  1279. {-                                                      -}
  1280. {- This procedure uses Richardson extrapolation         -}
  1281. {- to improve the current approximation to the          -}
  1282. {- derivative (OldEstimate).  The result is returned    -}
  1283. {- as the Iter element in the array NewEstimate         -}
  1284. {--------------------------------------------------------}
  1285.  
  1286. var
  1287.   Extrap : integer;
  1288.   FourToTheExtrapMinus1 : Float;
  1289.  
  1290. begin
  1291.   FourToTheExtrapMinus1 := 1;
  1292.   for Extrap := 2 to Iter do
  1293.   begin
  1294.     FourToTheExtrapMinus1 := FourToTheExtrapMinus1 * 4;
  1295.     NewEstimate[Extrap] := (FourToTheExtrapMinus1 * NewEstimate[Extrap - 1] -
  1296.                         OldEstimate[Extrap - 1]) / (FourToTheExtrapMinus1 - 1);
  1297.   end;
  1298. end; { procedure Extrapolate }
  1299.  
  1300. begin { procedure Differentiate }
  1301.   for Term := 1 to NumDeriv do
  1302.   begin
  1303.     X := XDeriv[Term];
  1304.     if ABS(X) < TNNearlyZero then
  1305.       DeltaX := Tolerance
  1306.     else
  1307.       DeltaX := ABS(X * Tolerance);
  1308.     OldEstimate[1] := EvaluateSecondDeriv(X, DeltaX);
  1309.     Iter := 1;
  1310.     TwoToTheIterMinus2 := 1;
  1311.     repeat
  1312.       Iter := Succ(Iter);
  1313.       DeltaX := DeltaX / 2;
  1314.       NewEstimate[1] := EvaluateSecondDeriv(X, DeltaX);
  1315.       TwoToTheIterMinus2 := TwoToTheIterMinus2 * 2;
  1316.       Extrapolate(Iter, OldEstimate, NewEstimate);
  1317.       OldEstimate := NewEstimate;
  1318.     until (ABS(NewEstimate[Iter - 1] - NewEstimate[Iter]) <=
  1319.            ABS(Tolerance * NewEstimate[Iter])) or
  1320.            (ABS(DeltaX) < TNNearlyZero);
  1321.     { Prevents excessive roundoff }
  1322.     YDeriv[Term] := NewEstimate[Iter];
  1323.   end;
  1324. end; { procedure Differentiate }
  1325.  
  1326. begin { procedure SecondDerivative }
  1327.     if Tolerance < TNNearlyZero then
  1328.       Error := 1
  1329.     else
  1330.       Differentiate(Tolerance, NumDeriv, XDeriv, YDeriv);
  1331. end; { procedure SecondDerivative }
  1332.  
  1333. end. { Differ }
  1334.