1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2008 The Regents of the University of California
4 //
5 // This file is part of Qbox
6 //
7 // Qbox is distributed under the terms of the GNU General Public License
8 // as published by the Free Software Foundation, either version 2 of
9 // the License, or (at your option) any later version.
10 // See the file COPYING in the root directory of this distribution
11 // or <http://www.gnu.org/licenses/>.
12 //
13 ///////////////////////////////////////////////////////////////////////////////
14 //
15 // spline.cpp
16 //
17 ///////////////////////////////////////////////////////////////////////////////
18 #include "spline.h"
19 #include <cassert>
20
tridsolve(int n,double * d,double * e,double * f,double * x)21 void tridsolve(int n, double* d, double* e, double* f, double* x)
22 {
23 // solve the tridiagonal system Ax=b
24 // d[i] = a(i,i)
25 // e[i] = a(i,i+1) (superdiagonal of A, e[n-1] not defined)
26 // f[i] = a(i,i-1) (subdiagonal of A, f[0] not defined)
27 // x[i] = right-hand side b as input
28 // x[i] = solution as output
29
30 for ( int i = 1; i < n; i++ )
31 {
32 f[i] /= d[i-1];
33 d[i] -= f[i]*e[i-1];
34 }
35
36 for ( int i = 1; i < n; i++ )
37 x[i] -= f[i]*x[i-1];
38
39 x[n-1] /= d[n-1];
40
41 for ( int i = n-2; i >= 0; i-- )
42 x[i] = (x[i]-e[i]*x[i+1])/d[i];
43 }
44
spline(int n,double * x,double * y,double yp_left,double yp_right,int bcnat_left,int bcnat_right,double * y2)45 void spline(int n, double *x, double *y, double yp_left, double yp_right,
46 int bcnat_left, int bcnat_right, double *y2)
47 {
48 const double third = 1.0/3.0;
49 const double sixth = 1.0/6.0;
50 double *d = new double[n];
51 double *e = new double[n];
52 double *f = new double[n];
53 if ( bcnat_left == 0 )
54 {
55 // use derivative yp_left at x[0]
56 const double h = x[1]-x[0];
57 assert(h>0.0);
58 d[0] = third*h;
59 e[0] = sixth*h;
60 f[0] = 0.0;
61 y2[0] = (y[1]-y[0])/h - yp_left;
62 }
63 else
64 {
65 // use natural spline at x[0]
66 d[0] = 1.0;
67 e[0] = 0.0;
68 f[0] = 0.0;
69 y2[0] = 0.0;
70 }
71 if ( bcnat_right == 0 )
72 {
73 // use derivative yp_right at x[n-1]
74 const double h = x[n-1]-x[n-2];
75 assert(h>0.0);
76 d[n-1] = third*h;
77 e[n-1] = 0.0;
78 f[n-1] = sixth*h;
79 y2[n-1] = yp_right - (y[n-1]-y[n-2])/h;
80 }
81 else
82 {
83 // use natural spline at x[n-1]
84 d[n-1] = 1.0;
85 e[n-1] = 0.0;
86 f[n-1] = 0.0;
87 y2[n-1] = 0.0;
88 }
89
90 // tridiagonal matrix
91 for ( int i = 1; i < n-1; i++ )
92 {
93 const double hp = x[i+1]-x[i];
94 const double hm = x[i]-x[i-1];
95 assert(hp>0.0);
96 assert(hm>0.0);
97 d[i] = third * (hp+hm);
98 e[i] = sixth * hp;
99 f[i] = sixth * hm;
100 y2[i] = (y[i+1]-y[i])/hp - (y[i]-y[i-1])/hm;
101 }
102
103 tridsolve(n,d,e,f,y2);
104
105 delete [] d;
106 delete [] e;
107 delete [] f;
108 }
109
splint(int n,double * xa,double * ya,double * y2a,double x,double * y)110 void splint (int n, double *xa, double *ya, double *y2a, double x, double *y)
111 {
112 int k;
113 double a,b,h;
114
115 int kl = 0;
116 int kh = n-1;
117
118 while ( kh - kl > 1 )
119 {
120 k = ( kh + kl ) / 2;
121 if ( xa[k] > x )
122 kh = k;
123 else
124 kl = k;
125 }
126
127 h = xa[kh] - xa[kl];
128 assert ( h > 0.0 );
129
130 a = ( xa[kh] - x ) / h;
131 b = ( x - xa[kl] ) / h;
132
133 *y = a * ya[kl] + b * ya[kh] + h * h * (1.0/6.0) *
134 ( (a*a*a-a) * y2a[kl] + (b*b*b-b) * y2a[kh] );
135
136 }
137
splintd(int n,double * xa,double * ya,double * y2a,double x,double * y,double * dy)138 void splintd (int n, double *xa, double *ya, double *y2a,
139 double x, double *y, double *dy)
140 {
141 int k;
142 double a,b,h;
143
144 int kl = 0;
145 int kh = n-1;
146
147 while ( kh - kl > 1 )
148 {
149 k = ( kh + kl ) / 2;
150 if ( xa[k] > x )
151 kh = k;
152 else
153 kl = k;
154 }
155
156 h = xa[kh] - xa[kl];
157 assert ( h > 0.0 );
158
159 a = ( xa[kh] - x ) / h;
160 b = ( x - xa[kl] ) / h;
161
162 *y = a * ya[kl] + b * ya[kh] + h * h * (1.0/6.0) *
163 ( (a*a*a-a) * y2a[kl] + (b*b*b-b) * y2a[kh] );
164
165 *dy = ( ya[kh] - ya[kl] ) / h +
166 h * ( ( (1.0/6.0) - 0.5 * a * a ) * y2a[kl] +
167 ( 0.5 * b * b - (1.0/6.0) ) * y2a[kh] );
168 }
169