1 | /* ------------------------------------------------------------------------- |
---|
2 | * spline.c |
---|
3 | * |
---|
4 | * --- cubic spline interpolation for non-periodic functions |
---|
5 | * "The latest algorithm collections in C" (1991) (in Japanese) |
---|
6 | * Haruhiko Okumura, Professor (Computer Science), Ph.D. |
---|
7 | * Matsusaka University |
---|
8 | * |
---|
9 | * modified by N.Sakaki |
---|
10 | * |
---|
11 | * $Id: spline.c,v 1.1 2003/03/27 07:38:45 sakaki Exp sakaki $ |
---|
12 | * ------------------------------------------------------------------------- |
---|
13 | */ |
---|
14 | /* |
---|
15 | */ |
---|
16 | #include <stdlib.h> |
---|
17 | #include "Mspline.hh" |
---|
18 | |
---|
19 | Int_t Mcsp_maketable(Double_t x[], Double_t y[], Double_t z[], size_t size) |
---|
20 | { |
---|
21 | Int_t i; |
---|
22 | Double_t t; |
---|
23 | Double_t *h, *d; |
---|
24 | |
---|
25 | h=(Double_t *)malloc(sizeof(Double_t)*size); |
---|
26 | d=(Double_t *)malloc(sizeof(Double_t)*size); |
---|
27 | if(h==NULL || d==NULL) |
---|
28 | return -1; |
---|
29 | |
---|
30 | z[0] = 0; z[size-1] = 0; /* y"(x) / 6 at both end points */ |
---|
31 | for(i=0; i<(Int_t)(size-1); i++) |
---|
32 | { |
---|
33 | h[i ] = x[i+1] - x[i]; |
---|
34 | d[i+1] = (y[i+1] - y[i])/h[i]; |
---|
35 | } |
---|
36 | |
---|
37 | z[1] = d[2] - d[1] - h[0]*z[0]; |
---|
38 | d[1] = 2.0*(x[2]-x[0]); |
---|
39 | |
---|
40 | for (i=1; i<(Int_t)(size-2); i++) |
---|
41 | { |
---|
42 | t = h[i]/d[i]; |
---|
43 | z[i+1] = d[i+2] - d[i+1] - z[i]*t; |
---|
44 | d[i+1] = 2.0*(x[i+2]-x[i]) - h[i]*t; |
---|
45 | } |
---|
46 | |
---|
47 | z[size-2] -= h[size-2]*z[size-1]; |
---|
48 | for (i=size-2; i>0; i--) |
---|
49 | z[i] = (z[i] - h[i]*z[i+1]) /d[i]; |
---|
50 | |
---|
51 | free(h); |
---|
52 | free(d); |
---|
53 | |
---|
54 | return 0; |
---|
55 | } |
---|
56 | |
---|
57 | Double_t Mcspline(Double_t t, Double_t x[], Double_t y[], Double_t z[], size_t size) |
---|
58 | { |
---|
59 | Int_t i, j, k; |
---|
60 | Double_t d, h; |
---|
61 | |
---|
62 | i = 0; j = size - 1; |
---|
63 | while (i<j) |
---|
64 | { |
---|
65 | k = (i+j)/2; |
---|
66 | if (x[k]<t) |
---|
67 | i = k+1; |
---|
68 | else |
---|
69 | j = k; |
---|
70 | } |
---|
71 | |
---|
72 | if (i>0) |
---|
73 | i--; |
---|
74 | h = x[i+1] - x[i]; |
---|
75 | d = t - x[i]; |
---|
76 | |
---|
77 | return (( (z[i+1] - z[i])*d/h + z[i]*3.0) * d |
---|
78 | + ( (y[i+1] - y[i])/h - (z[i]*2. + z[i+1])* h) )*d + y[i]; |
---|
79 | } |
---|
80 | |
---|
81 | #if 0 |
---|
82 | #include <stdio.h> |
---|
83 | |
---|
84 | #define N 5 |
---|
85 | Double_t x[N] = { 0, 10, 20, 30, 40 }, |
---|
86 | y[N] = { 610.66, 1227.4, 2338.1, 4244.9, 7381.2 }, |
---|
87 | z[N]; |
---|
88 | Int_t main() |
---|
89 | { |
---|
90 | Int_t i; |
---|
91 | |
---|
92 | csp_maketable(x, y, z, N); |
---|
93 | for (i = 10; i <= 30; i++) |
---|
94 | printf("%3d %6.1f\n", i, cspline(i, x, y, z, N)); |
---|
95 | return EXIT_SUCCESS; |
---|
96 | } |
---|
97 | #endif |
---|