3 /* natural cubic spline interpolation based on Numerical Recipes in C */
5 /* calculate second derivatives; assumes strictly increasing x values */
7 find_spline_derivs(double *x
, double *y
, int n
,
13 y2
[0] = u
[0] = 0.0; /* lower boundary condition for natural spline */
15 /* decomposition loop for the tridiagonal algorithm */
16 for (i
= 1; i
< n
- 1; i
++) {
17 y2
[i
] = u
[i
] = 0.0; /* set in case a zero test is failed */
18 if (x
[i
- 1] < x
[i
] && x
[i
] < x
[i
+ 1]) {
19 sig
= (x
[i
] - x
[i
- 1]) / (x
[i
+ 1] - x
[i
- 1]);
20 p
= sig
* y2
[i
- 1] + 2.0;
22 y2
[i
] = (sig
- 1.0) / p
;
23 u
[i
] = (y
[i
+ 1] - y
[i
]) / (x
[i
+ 1] - x
[i
])
24 - (y
[i
] - y
[i
- 1]) / (x
[i
] - x
[i
- 1]);
25 u
[i
] = (6.0 * u
[i
] / (x
[i
+ 1] - x
[i
- 1]) - sig
* u
[i
- 1]) / p
;
30 /* upper boundary condition for natural spline */
33 /* backsubstitution loop of the tridiagonal algorithm */
34 for (k
= n
- 2; k
>= 0; k
--)
35 y2
[k
] = y2
[k
] * y2
[k
+ 1] + u
[k
];
38 /* interpolate or extrapolate value at x using results of find_spline_derivs */
40 spline_interp(double *xa
, double *ya
, double *y2a
,
41 int n
, double x
, double *y
)
48 b
= (h
> 0.0) ? (ya
[1] - ya
[0]) / h
- h
* y2a
[1] / 6.0 : 0.0;
49 *y
= ya
[0] + b
* (x
- xa
[0]);
51 else if (x
>= xa
[n
- 1]) {
52 h
= xa
[n
- 1] - xa
[n
- 2];
53 b
= (h
> 0.0) ? (ya
[n
- 1] - ya
[n
- 2]) / h
+ h
* y2a
[n
- 2] / 6.0 : 0.0;
54 *y
= ya
[n
- 1] + b
* (x
- xa
[n
- 1]);
57 /* try a linear interpolation for equally spaced x values */
58 k
= (n
- 1) * (x
- xa
[0]) / (xa
[n
- 1] - xa
[0]);
60 /* make sure the range is right */
62 if (k
> n
- 2) k
= n
- 2;
64 /* bisect if necessary until the bracketing interval is found */
65 klo
= (x
>= xa
[k
]) ? k
: 0;
66 khi
= (x
< xa
[k
+ 1]) ? k
+ 1 : n
- 1;
67 while (khi
- klo
> 1) {
69 if (xa
[k
] > x
) khi
= k
;
74 h
= xa
[khi
] - xa
[klo
];
76 a
= (xa
[khi
] - x
) / h
;
77 b
= (x
- xa
[klo
]) / h
;
78 *y
= a
* ya
[klo
] + b
* ya
[khi
]
79 + ((a
* a
* a
- a
) * y2a
[klo
] + (b
* b
* b
- b
) * y2a
[khi
]) * (h
* h
) / 6.0;
81 else *y
= (ya
[klo
] + ya
[khi
]) / 2.0; /* should not be needed */
86 fit_spline(int n
, double *x
, double *y
,
87 int ns
, double *xs
, double *ys
, double *work
)
92 y2
= work
; u
= work
+ n
;
94 if (n
< 2 || ns
< 1) return (1); /* signal an error */
95 for (i
= 1; i
< n
; i
++)
96 if (x
[i
- 1] >= x
[i
]) return(1); /* signal an error */
98 find_spline_derivs(x
, y
, n
, y2
, u
);
100 for (i
= 0; i
< ns
; i
++)
101 spline_interp(x
, y
, y2
, n
, xs
[i
], &ys
[i
]);