1
2 /*
3 * spline.c
4 *
5 * Routines from numerical recipes for cubic spline interpolation
6 *
7 */
8
9 #include "nr.h"
10
11 /*
12 * Given arrays x[1..n] and y[1..n] containing a tabulated function, i.e., yi = f(xi), with
13 * x1 <x2 < :: : < xN, and given values yp1 and ypn for the first derivative of the interpolating
14 * function at points 1 and n, respectively, this routine returns an array y2[1..n] that contains
15 * the second derivatives of the interpolating function at the tabulated points xi. If yp1 and/or
16 * ypn are equal to 1 x 10^30 or larger, the routine is signaled to set the corresponding boundary
17 * condition for a natural spline, with zero second derivative on that boundary.
18 */
19
20 void spline(float x[], float y[], int n, float yp1, float ypn, float y2[]) {
21 int i,k;
22 float p,qn,sig,un,*u;
23
24 u=vector(1,n-1);
25 if (yp1 > 0.99e30) // The lower boundary condition is set either to be "natural
26 y2[1]=u[1]=0.0;
27 else { // or else to have a specified first derivative.
28 y2[1] = -0.5;
29 u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);
30 }
31
32 // This is the decomposition loop of the tridiagonal algorithm.
33 // y2 and u are used for temporary storage of the decomposed factors.
34 for (i=2;i<=n-1;i++) {
35 sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
36 p=sig*y2[i-1]+2.0;
37 y2[i]=(sig-1.0)/p;
38 u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
39 u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
40 }
41
42 if (ypn > 0.99e30) // The upper boundary condition is set either to be "natural"
43 qn=un=0.0;
44 else { // or else to have a speci ed rst derivative.
45 qn=0.5;
46 un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));
47 }
48
49 y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
50 for (k=n-1;k>=1;k--) // This is the backsubstitution loop of the tridiagonal
51 y2[k]=y2[k]*y2[k+1]+u[k]; // algorithm.
52 free_vector(u,1,n-1);
53 }
54
55
56 /*
57 * Given the arrays xa[1..n] and ya[1..n], which tabulate a function (with the xai's in order),
58 * and given the array y2a[1..n], which is the output from spline above, and given a value of
59 * x, this routine returns a cubic-spline interpolated value y.
60 */
61 void splint(float xa[], float ya[], float y2a[], int n, float x, float *y)
62 {
63 void nrerror(char error_text[]);
64 int klo,khi,k;
65 float h,b,a;
66
67 /*
68 * We will find the right place in the table by means of
69 * bisection. This is optimal if sequential calls to this
70 * routine are at random values of x. If sequential calls
71 * are in order, and closely spaced, one would do better
72 * to store previous values of klo and khi and test if
73 * they remain appropriate on the next call.
74 */
75
76 klo=1;
77 khi=n;
78 while (khi-klo > 1) {
79 k=(khi+klo) >> 1;
80 if (xa[k] > x) khi=k;
81 else klo=k;
82 } // klo and khi now bracket the input value of x.
83
84 h=xa[khi]-xa[klo];
85 if (h == 0.0) nrerror("Bad xa input to routine splint"); //The xa's must be
86 a=(xa[khi]-x)/h; // distinct.
87 b=(x-xa[klo])/h; // Cubic spline polynomial is now evaluated.
88 *y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
89 }
90
91
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.