BSD 4_3 release

[unix-history] / usr / src / usr.bin / spline.c

diff --git a/usr/src/usr.bin/spline.c b/usr/src/usr.bin/spline.c
index 09c12f0..df28f6e 100644 (file)
--- a/usr/src/usr.bin/spline.c
+++ b/usr/src/usr.bin/spline.c
@@ -1,4 +1,7 @@
-static char *sccsid = "@(#)spline.c    4.2 (Berkeley) 11/27/82";
+#ifndef lint
+static char *sccsid = "@(#)spline.c    4.3 (Berkeley) 9/21/85";
+#endif
+

 #include <stdio.h>
 #include <math.h>
 
 #include <stdio.h>
 #include <math.h>
 


@@ -17,49 +20,49 @@ float zero = 0.;
 
 /* Spline fit technique
 let x,y be vectors of abscissas and ordinates
 
 /* Spline fit technique
 let x,y be vectors of abscissas and ordinates

-    h   be vector of differences h\e9i\e8=x\e9i\e8-x\e9i-1\e\e9\e8\e8
+    h   be vector of differences hi=xi-xi-1

     y"  be vector of 2nd derivs of approx function
 If the points are numbered 0,1,2,...,n+1 then y" satisfies
 (R W Hamming, Numerical Methods for Engineers and Scientists,
 2nd Ed, p349ff)
     y"  be vector of 2nd derivs of approx function
 If the points are numbered 0,1,2,...,n+1 then y" satisfies
 (R W Hamming, Numerical Methods for Engineers and Scientists,
 2nd Ed, p349ff)

-       h\e9i\e8y"\b\e9i-1\e9\e8\e8+2(h\e9i\e8+h\e9i+1\e8)y"\b\e9i\e8+h\e9i+1\e8y"\b\e9i+1\e8
+       hiy"i-1+2(hi+hi+1)y"i+hi+1y"i+1

        
        

-       = 6[(y\e9i+1\e8-y\e9i\e8)/h\e9i+1\e8-(y\e9i\e8-y\e9i-1\e8)/h\e9i\e8]   i=1,2,...,n
+       = 6[(yi+1-yi)/hi+1-(yi-yi-1)/hi]   i=1,2,...,n

 
 

-where y"\b\e90\e8 = y"\b\e9n+1\e8 = 0
+where y"0 = y"n+1 = 0

 This is a symmetric tridiagonal system of the form
 
 This is a symmetric tridiagonal system of the form
 

-       | a\e91\e8 h\e92\e8               |  |y"\b\e91\e8|      |b\e91\e8|
-       | h\e92\e8 a\e92\e8 h\e93\e8            |  |y"\b\e92\e8|      |b\e92\e8|
-       |    h\e93\e8 a\e93\e8 h\e94\e8         |  |y"\b\e93\e8|  =   |b\e93\e8|
+       | a1 h2               |  |y"1|      |b1|
+       | h2 a2 h3            |  |y"2|      |b2|
+       |    h3 a3 h4         |  |y"3|  =   |b3|

        |         .           |  | .|      | .|
        |            .        |  | .|      | .|
 It can be triangularized into
        |         .           |  | .|      | .|
        |            .        |  | .|      | .|
 It can be triangularized into

-       | d\e91\e8 h\e92\e8               |  |y"\b\e91\e8|      |r\e91\e8|
-       |    d\e92\e8 h\e93\e8            |  |y"\b\e92\e8|      |r\e92\e8|
-       |       d\e93\e8 h\e94\e8         |  |y"\b\e93\e8|  =   |r\e93\e8|
+       | d1 h2               |  |y"1|      |r1|
+       |    d2 h3            |  |y"2|      |r2|
+       |       d3 h4         |  |y"3|  =   |r3|

        |          .          |  | .|      | .|
        |             .       |  | .|      | .|
 where
        |          .          |  | .|      | .|
        |             .       |  | .|      | .|
 where

-       d\e91\e8 = a\e91\e8
+       d1 = a1

 
 

-       r\e90\e8 = 0
+       r0 = 0

 
 

-       d\e9i\e8 = a\e9i\e8 - h\e9i\e8\b\e82\e9/d\e9i-1\e8 1<i<\b_n
+       di = ai - hi2/di-1      1<i<_n

 
 

-       r\e9i\e8 = b\e9i\e8 - h\e9i\e8r\e9i-1\e8/d\e9i-1\ei\e8     1<\b_i<\b_n
+       ri = bi - hiri-1/di-1\ei 1<_i<_n

 
 the back solution is
 
 the back solution is

-       y"\b\e9n\e8 = r\e9n\e8/d\e9n\e8
+       y"n = rn/dn

 
 

-       y"\b\e9i\e8 = (r\e9i\e8-h\e9i+1\e8y"\b\e9i+1\e8)/d\e9i\e8   1<\b_i<n
+       y"i = (ri-hi+1y"i+1)/di 1<_i<n

 
 

-superficially, d\e9i\e8 and r\e9i\e8 don't have to be stored for they can be
+superficially, di and ri don't have to be stored for they can be

 recalculated backward by the formulas
 
 recalculated backward by the formulas
 

-       d\e9i-1\e8 = h\e9i\e8\b\e82\e9/(a\e9i\e8-d\e9i\e8) 1<i<\b_n
+       di-1 = hi2/(ai-di)      1<i<_n

 
 

-       r\e9i-1\e8 = (b\e9i\e8-r\e9i\e8)d\e9i-1\e8/h\e9i\e8       1<i<\b_n
+       ri-1 = (bi-ri)di-1/hi   1<i<_n

 
 unhappily it turns out that the recursion forward for d
 is quite strongly geometrically convergent--and is wildly
 
 unhappily it turns out that the recursion forward for d
 is quite strongly geometrically convergent--and is wildly


@@ -69,74 +72,74 @@ results must be kept.
 
 Note that n-1 in the program below plays the role of n+1 in the theory
 
 
 Note that n-1 in the program below plays the role of n+1 in the theory
 

-Other boundary conditions\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b_________________________
+Other boundary conditions_________________________

 
 The boundary conditions are easily generalized to handle
 
 
 The boundary conditions are easily generalized to handle
 

-       y\e90\e8\b" = ky\e91\e8\b", y\e9n+1\e8\b\b\b"   = ky\e9n\e8\b"
+       y0" = ky1", yn+1"   = kyn"

 
 for some constant k.  The above analysis was for k = 0;
 k = 1 fits parabolas perfectly as well as stright lines;
 k = 1/2 has been recommended as somehow pleasant.
 
 
 for some constant k.  The above analysis was for k = 0;
 k = 1 fits parabolas perfectly as well as stright lines;
 k = 1/2 has been recommended as somehow pleasant.
 

-All that is necessary is to add h\e91\e8 to a\e91\e8 and h\e9n+1\e8 to a\e9n\e8.
+All that is necessary is to add h1 to a1 and hn+1 to an.

 
 
 
 

-Periodic case\b\b\b\b\b\b\b\b\b\b\b\b\b_____________
+Periodic case_____________

 
 To do this, add 1 more row and column thus
 
 
 To do this, add 1 more row and column thus
 

-       | a\e91\e8 h\e92\e8            h\e91\e8 |  |y\e91\e8\b"|     |b\e91\e8|
-       | h\e92\e8 a\e92\e8 h\e93\e8            |  |y\e92\e8\b"|     |b\e92\e8|
-       |    h\e93\e8 a\e94\e8 h\e94\e8         |  |y\e93\e8\b"|     |b\e93\e8|
+       | a1 h2            h1 |  |y1"|     |b1|
+       | h2 a2 h3            |  |y2"|     |b2|
+       |    h3 a4 h4         |  |y3"|     |b3|

        |                     |  | .|  =  | .|
        |             .       |  | .|     | .|
        |                     |  | .|  =  | .|
        |             .       |  | .|     | .|

-       | h\e91\e8            h\e90\e8 a\e90\e8 |  | .|     | .|
+       | h1            h0 a0 |  | .|     | .|

 
 

-where h\e90\e8=\b_ h\e9n+1\e8
+where h0=_ hn+1

 
 The same diagonalization procedure works, except for
 
 The same diagonalization procedure works, except for

-the effect of the 2 corner elements.  Let s\e9i\e8 be the part
-of the last element in the i\e8th\e9 "diagonalized" row that
+the effect of the 2 corner elements.  Let si be the part
+of the last element in the ith "diagonalized" row that

 arises from the extra top corner element.
 
 arises from the extra top corner element.
 

-               s\e91\e8 = h\e91\e8
+               s1 = h1

 
 

-               s\e9i\e8 = -s\e9i-1\e8h\e9i\e8/d\e9i-1\e8       2<\b_i<\b_n+1
+               si = -si-1hi/di-1       2<_i<_n+1

 
 After "diagonalizing", the lower corner element remains.
 
 After "diagonalizing", the lower corner element remains.

-Call t\e9i\e8 the bottom element that appears in the i\e8th\e9 colomn
+Call ti the bottom element that appears in the ith colomn

 as the bottom element to its left is eliminated
 
 as the bottom element to its left is eliminated
 

-               t\e91\e8 = h\e91\e8
+               t1 = h1

 
 

-               t\e9i\e8 = -t\e9i-1\e8h\e9i\e8/d\e9i-1\e8
+               ti = -ti-1hi/di-1

 
 

-Evidently t\e9i\e8 = s\e9i\e8.
+Evidently ti = si.

 Elimination along the bottom row
 introduces further corrections to the bottom right element
 and to the last element of the right hand side.
 Call these corrections u and v.
 
 Elimination along the bottom row
 introduces further corrections to the bottom right element
 and to the last element of the right hand side.
 Call these corrections u and v.
 

-       u\e91\e8 = v\e91\e8 = 0
+       u1 = v1 = 0

 
 

-       u\e9i\e8 = u\e9i-1\e8-s\e9i-1\e8*t\e9i-1\e8/d\e9i-1\e8
+       ui = ui-1-si-1*ti-1/di-1

 
 

-       v\e9i\e8 = v\e9i-1\e8-r\e9i-1\e8*t\e9i-1\e8/d\e9i-1\e8    2<\b_i<\b_n+1
+       vi = vi-1-ri-1*ti-1/di-1        2<_i<_n+1

 
 The back solution is now obtained as follows
 
 
 The back solution is now obtained as follows
 

-       y"\b\e9n+1\e8 = (r\e9n+1\e8+v\e9n+1\e8)/(d\e9n+1\e8+s\e9n+1\e8+t\e9n+1\e8+u\e9n+1\e8)
+       y"n+1 = (rn+1+vn+1)/(dn+1+sn+1+tn+1+un+1)

 
 

-       y"\b\e9i\e8 = (r\e9i\e8-h\e9i+1\e8*y\e9i+1\e8-s\e9i\e8*y\e9n+1\e8)/d\e9i\e8    1<\b_i<\b_n
+       y"i = (ri-hi+1*yi+1-si*yn+1)/di 1<_i<_n

 
 

-Interpolation in the interval x\e9i\e8<\b_x<\b_x\e9i+1\e8 is by the formula
+Interpolation in the interval xi<_x<_xi+1 is by the formula

 
 

-       y = y\e9i\e8x\e9+\e8 + y\e9i+1\e8x\e9-\e8 -(h\e82\e9\b\e9i+1\e8/6)[y"\b\e9i\e8(x\e9+\e8-x\e9+\e8\e8\b3\e9)+y"\b\e9i+1\e8(x\e9-\e8-x\e9-\e8\b\e83\e9)]
+       y = yix+ + yi+1x- -(h2i+1/6)[y"i(x+-x+3)+y"i+1(x--x-3)]

 where
 where

-       x\e9+\e8 = x\e9i+1\e8-x
+       x+ = xi+1-x

 
 

-       x\e9-\e8 = x-x\e9i\e8
+       x- = x-xi

 */
 
 float
 */
 
 float