//SEE: http://www.tol-project.org:8081/tolapi/spline3_8h.html ************************************************************************* Copyright (c) 2007, Sergey Bochkanov (ALGLIB project). Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: - Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. - Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer listed in this license in the documentation and/or other materials provided with the distribution. - Neither the name of the copyright holders nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *************************************************************************/ #ifndef _spline3_h #define _spline3_h #include "ap.h" /************************************************************************* This subroutine builds linear spline coefficients table. Input parameters: X - spline nodes, array[0..N-1] Y - function values, array[0..N-1] N - points count, N>=2 Output parameters: C - coefficients table. Used by SplineInterpolation and other subroutines from this file. -- ALGLIB PROJECT -- Copyright 24.06.2007 by Bochkanov Sergey *************************************************************************/ void buildlinearspline(ap::real_1d_array x, ap::real_1d_array y, int n, ap::real_1d_array& c); /************************************************************************* This subroutine builds cubic spline coefficients table. Input parameters: X - spline nodes, array[0..N-1] Y - function values, array[0..N-1] N - points count, N>=2 BoundLType - boundary condition type for the left boundary BoundL - left boundary condition (first or second derivative, depending on the BoundLType) BoundRType - boundary condition type for the right boundary BoundR - right boundary condition (first or second derivative, depending on the BoundRType) Output parameters: C - coefficients table. Used by SplineInterpolation and other subroutines from this file. The BoundLType/BoundRType parameters can have the following values: * 0, which corresponds to the parabolically terminated spline (BoundL/BoundR are ignored). * 1, which corresponds to the first derivative boundary condition * 2, which corresponds to the second derivative boundary condition -- ALGLIB PROJECT -- Copyright 23.06.2007 by Bochkanov Sergey *************************************************************************/ void buildcubicspline(ap::real_1d_array x, ap::real_1d_array y, int n, int boundltype, double boundl, int boundrtype, double boundr, ap::real_1d_array& c); /************************************************************************* This subroutine builds cubic Hermite spline coefficients table. Input parameters: X - spline nodes, array[0..N-1] Y - function values, array[0..N-1] D - derivatives, array[0..N-1] N - points count, N>=2 Output parameters: C - coefficients table. Used by SplineInterpolation and other subroutines from this file. -- ALGLIB PROJECT -- Copyright 23.06.2007 by Bochkanov Sergey *************************************************************************/ void buildhermitespline(ap::real_1d_array x, ap::real_1d_array y, ap::real_1d_array d, int n, ap::real_1d_array& c); /************************************************************************* This subroutine builds Akima spline coefficients table. Input parameters: X - spline nodes, array[0..N-1] Y - function values, array[0..N-1] N - points count, N>=5 Output parameters: C - coefficients table. Used by SplineInterpolation and other subroutines from this file. -- ALGLIB PROJECT -- Copyright 24.06.2007 by Bochkanov Sergey *************************************************************************/ void buildakimaspline(ap::real_1d_array x, ap::real_1d_array y, int n, ap::real_1d_array& c); /************************************************************************* This subroutine calculates the value of the spline at the given point X. Input parameters: C - coefficients table. Built by BuildLinearSpline, BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline. X - point Result: S(x) -- ALGLIB PROJECT -- Copyright 23.06.2007 by Bochkanov Sergey *************************************************************************/ double splineinterpolation(const ap::real_1d_array& c, double x); /************************************************************************* This subroutine differentiates the spline. Input parameters: C - coefficients table. Built by BuildLinearSpline, BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline. X - point Result: S - S(x) DS - S'(x) D2S - S''(x) -- ALGLIB PROJECT -- Copyright 24.06.2007 by Bochkanov Sergey *************************************************************************/ void splinedifferentiation(const ap::real_1d_array& c, double x, double& s, double& ds, double& d2s); /************************************************************************* This subroutine makes the copy of the spline. Input parameters: C - coefficients table. Built by BuildLinearSpline, BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline. Result: CC - spline copy -- ALGLIB PROJECT -- Copyright 29.06.2007 by Bochkanov Sergey *************************************************************************/ void splinecopy(const ap::real_1d_array& c, ap::real_1d_array& cc); /************************************************************************* This subroutine unpacks the spline into the coefficients table. Input parameters: C - coefficients table. Built by BuildLinearSpline, BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline. X - point Result: Tbl - coefficients table, unpacked format, array[0..N-2, 0..5]. For I = 0...N-2: Tbl[I,0] = X[i] Tbl[I,1] = X[i+1] Tbl[I,2] = C0 Tbl[I,3] = C1 Tbl[I,4] = C2 Tbl[I,5] = C3 On [x[i], x[i+1]] spline is equals to: S(x) = C0 + C1*t + C2*t^2 + C3*t^3 t = x-x[i] -- ALGLIB PROJECT -- Copyright 29.06.2007 by Bochkanov Sergey *************************************************************************/ void splineunpack(const ap::real_1d_array& c, int& n, ap::real_2d_array& tbl); /************************************************************************* This subroutine performs linear transformation of the spline argument. Input parameters: C - coefficients table. Built by BuildLinearSpline, BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline. A, B- transformation coefficients: x = A*t + B Result: C - transformed spline -- ALGLIB PROJECT -- Copyright 30.06.2007 by Bochkanov Sergey *************************************************************************/ void splinelintransx(ap::real_1d_array& c, double a, double b); /************************************************************************* This subroutine performs linear transformation of the spline. Input parameters: C - coefficients table. Built by BuildLinearSpline, BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline. A, B- transformation coefficients: S2(x) = A*S(x) + B Result: C - transformed spline -- ALGLIB PROJECT -- Copyright 30.06.2007 by Bochkanov Sergey *************************************************************************/ void splinelintransy(ap::real_1d_array& c, double a, double b); /************************************************************************* This subroutine integrates the spline. Input parameters: C - coefficients table. Built by BuildLinearSpline, BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline. X - right bound of the integration interval [a, x] Result: integral(S(t)dt,a,x) -- ALGLIB PROJECT -- Copyright 23.06.2007 by Bochkanov Sergey *************************************************************************/ double splineintegration(const ap::real_1d_array& c, double x); /************************************************************************* Obsolete subroutine, left for backward compatibility. *************************************************************************/ void spline3buildtable(int n, const int& diffn, ap::real_1d_array x, ap::real_1d_array y, const double& boundl, const double& boundr, ap::real_2d_array& ctbl); /************************************************************************* Obsolete subroutine, left for backward compatibility. *************************************************************************/ double spline3interpolate(int n, const ap::real_2d_array& c, const double& x); #endif MORE /************************************************************************* 00002 Copyright (c) 2007, Sergey Bochkanov (ALGLIB project). 00003 00004 Redistribution and use in source and binary forms, with or without 00005 modification, are permitted provided that the following conditions are 00006 met: 00007 00008 - Redistributions of source code must retain the above copyright 00009 notice, this list of conditions and the following disclaimer. 00010 00011 - Redistributions in binary form must reproduce the above copyright 00012 notice, this list of conditions and the following disclaimer listed 00013 in this license in the documentation and/or other materials 00014 provided with the distribution. 00015 00016 - Neither the name of the copyright holders nor the names of its 00017 contributors may be used to endorse or promote products derived from 00018 this software without specific prior written permission. 00019 00020 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 00021 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 00022 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 00023 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 00024 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 00025 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 00026 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 00027 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 00028 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 00029 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 00030 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00031 *************************************************************************/ 00032 00033 #include 00034 #include "spline3.h" 00035 00036 void heapsortpoints(ap::real_1d_array& x, ap::real_1d_array& y, int n); 00037 void heapsortdpoints(ap::real_1d_array& x, 00038 ap::real_1d_array& y, 00039 ap::real_1d_array& d, 00040 int n); 00041 void solvetridiagonal(ap::real_1d_array a, 00042 ap::real_1d_array b, 00043 ap::real_1d_array c, 00044 ap::real_1d_array d, 00045 int n, 00046 ap::real_1d_array& x); 00047 double diffthreepoint(double t, 00048 double x0, 00049 double f0, 00050 double x1, 00051 double f1, 00052 double x2, 00053 double f2); 00054 00055 /************************************************************************* 00056 This subroutine builds linear spline coefficients table. 00057 00058 Input parameters: 00059 X - spline nodes, array[0..N-1] 00060 Y - function values, array[0..N-1] 00061 N - points count, N>=2 00062 00063 Output parameters: 00064 C - coefficients table. Used by SplineInterpolation and other 00065 subroutines from this file. 00066 00067 -- ALGLIB PROJECT -- 00068 Copyright 24.06.2007 by Bochkanov Sergey 00069 *************************************************************************/ 00070 void buildlinearspline(ap::real_1d_array x, 00071 ap::real_1d_array y, 00072 int n, 00073 ap::real_1d_array& c) 00074 { 00075 int i; 00076 int tblsize; 00077 00078 ap::ap_error::make_assertion(n>=2); 00079 00080 // 00081 // Sort points 00082 // 00083 heapsortpoints(x, y, n); 00084 00085 // 00086 // Fill C: 00087 // C[0] - length(C) 00088 // C[1] - type(C): 00089 // 3 - general cubic spline 00090 // C[2] - N 00091 // C[3]...C[3+N-1] - x[i], i = 0...N-1 00092 // C[3+N]...C[3+N+(N-1)*4-1] - coefficients table 00093 // 00094 tblsize = 3+n+(n-1)*4; 00095 c.setbounds(0, tblsize-1); 00096 c(0) = tblsize; 00097 c(1) = 3; 00098 c(2) = n; 00099 for(i = 0; i <= n-1; i++) 00100 { 00101 c(3+i) = x(i); 00102 } 00103 for(i = 0; i <= n-2; i++) 00104 { 00105 c(3+n+4*i+0) = y(i); 00106 c(3+n+4*i+1) = (y(i+1)-y(i))/(x(i+1)-x(i)); 00107 c(3+n+4*i+2) = 0; 00108 c(3+n+4*i+3) = 0; 00109 } 00110 } 00111 00112 00113 /************************************************************************* 00114 This subroutine builds cubic spline coefficients table. 00115 00116 Input parameters: 00117 X - spline nodes, array[0..N-1] 00118 Y - function values, array[0..N-1] 00119 N - points count, N>=2 00120 BoundLType - boundary condition type for the left boundary 00121 BoundL - left boundary condition (first or second derivative, 00122 depending on the BoundLType) 00123 BoundRType - boundary condition type for the right boundary 00124 BoundR - right boundary condition (first or second derivative, 00125 depending on the BoundRType) 00126 00127 Output parameters: 00128 C - coefficients table. Used by SplineInterpolation and 00129 other subroutines from this file. 00130 00131 The BoundLType/BoundRType parameters can have the following values: 00132 * 0, which corresponds to the parabolically terminated spline 00133 (BoundL/BoundR are ignored). 00134 * 1, which corresponds to the first derivative boundary condition 00135 * 2, which corresponds to the second derivative boundary condition 00136 00137 -- ALGLIB PROJECT -- 00138 Copyright 23.06.2007 by Bochkanov Sergey 00139 *************************************************************************/ 00140 void buildcubicspline(ap::real_1d_array x, 00141 ap::real_1d_array y, 00142 int n, 00143 int boundltype, 00144 double boundl, 00145 int boundrtype, 00146 double boundr, 00147 ap::real_1d_array& c) 00148 { 00149 ap::real_1d_array a1; 00150 ap::real_1d_array a2; 00151 ap::real_1d_array a3; 00152 ap::real_1d_array b; 00153 ap::real_1d_array d; 00154 int i; 00155 int tblsize; 00156 double delta; 00157 double delta2; 00158 double delta3; 00159 00160 ap::ap_error::make_assertion(n>=2); 00161 ap::ap_error::make_assertion(boundltype==0||boundltype==1||boundltype==2); 00162 ap::ap_error::make_assertion(boundrtype==0||boundrtype==1||boundrtype==2); 00163 a1.setbounds(0, n-1); 00164 a2.setbounds(0, n-1); 00165 a3.setbounds(0, n-1); 00166 b.setbounds(0, n-1); 00167 00168 // 00169 // Special case: 00170 // * N=2 00171 // * parabolic terminated boundary condition on both ends 00172 // 00173 if( n==2&&boundltype==0&&boundrtype==0 ) 00174 { 00175 00176 // 00177 // Change task type 00178 // 00179 boundltype = 2; 00180 boundl = 0; 00181 boundrtype = 2; 00182 boundr = 0; 00183 } 00184 00185 // 00186 // 00187 // Sort points 00188 // 00189 heapsortpoints(x, y, n); 00190 00191 // 00192 // Left boundary conditions 00193 // 00194 if( boundltype==0 ) 00195 { 00196 a1(0) = 0; 00197 a2(0) = 1; 00198 a3(0) = 1; 00199 b(0) = 2*(y(1)-y(0))/(x(1)-x(0)); 00200 } 00201 if( boundltype==1 ) 00202 { 00203 a1(0) = 0; 00204 a2(0) = 1; 00205 a3(0) = 0; 00206 b(0) = boundl; 00207 } 00208 if( boundltype==2 ) 00209 { 00210 a1(0) = 0; 00211 a2(0) = 2; 00212 a3(0) = 1; 00213 b(0) = 3*(y(1)-y(0))/(x(1)-x(0))-0.5*boundl*(x(1)-x(0)); 00214 } 00215 00216 // 00217 // Central conditions 00218 // 00219 for(i = 1; i <= n-2; i++) 00220 { 00221 a1(i) = x(i+1)-x(i); 00222 a2(i) = 2*(x(i+1)-x(i-1)); 00223 a3(i) = x(i)-x(i-1); 00224 b(i) = 3*(y(i)-y(i-1))/(x(i)-x(i-1))*(x(i+1)-x(i))+3*(y(i+1)-y(i))/(x(i+1)-x(i))*(x(i)-x(i-1)); 00225 } 00226 00227 // 00228 // Right boundary conditions 00229 // 00230 if( boundrtype==0 ) 00231 { 00232 a1(n-1) = 1; 00233 a2(n-1) = 1; 00234 a3(n-1) = 0; 00235 b(n-1) = 2*(y(n-1)-y(n-2))/(x(n-1)-x(n-2)); 00236 } 00237 if( boundrtype==1 ) 00238 { 00239 a1(n-1) = 0; 00240 a2(n-1) = 1; 00241 a3(n-1) = 0; 00242 b(n-1) = boundr; 00243 } 00244 if( boundrtype==2 ) 00245 { 00246 a1(n-1) = 1; 00247 a2(n-1) = 2; 00248 a3(n-1) = 0; 00249 b(n-1) = 3*(y(n-1)-y(n-2))/(x(n-1)-x(n-2))+0.5*boundr*(x(n-1)-x(n-2)); 00250 } 00251 00252 // 00253 // Solve 00254 // 00255 solvetridiagonal(a1, a2, a3, b, n, d); 00256 00257 // 00258 // Now problem is reduced to the cubic Hermite spline 00259 // 00260 buildhermitespline(x, y, d, n, c); 00261 } 00262 00263 00264 /************************************************************************* 00265 This subroutine builds cubic Hermite spline coefficients table. 00266 00267 Input parameters: 00268 X - spline nodes, array[0..N-1] 00269 Y - function values, array[0..N-1] 00270 D - derivatives, array[0..N-1] 00271 N - points count, N>=2 00272 00273 Output parameters: 00274 C - coefficients table. Used by SplineInterpolation and 00275 other subroutines from this file. 00276 00277 -- ALGLIB PROJECT -- 00278 Copyright 23.06.2007 by Bochkanov Sergey 00279 *************************************************************************/ 00280 void buildhermitespline(ap::real_1d_array x, 00281 ap::real_1d_array y, 00282 ap::real_1d_array d, 00283 int n, 00284 ap::real_1d_array& c) 00285 { 00286 int i; 00287 int tblsize; 00288 double delta; 00289 double delta2; 00290 double delta3; 00291 00292 ap::ap_error::make_assertion(n>=2); 00293 00294 // 00295 // Sort points 00296 // 00297 heapsortdpoints(x, y, d, n); 00298 00299 // 00300 // Fill C: 00301 // C[0] - length(C) 00302 // C[1] - type(C): 00303 // 3 - general cubic spline 00304 // C[2] - N 00305 // C[3]...C[3+N-1] - x[i], i = 0...N-1 00306 // C[3+N]...C[3+N+(N-1)*4-1] - coefficients table 00307 // 00308 tblsize = 3+n+(n-1)*4; 00309 c.setbounds(0, tblsize-1); 00310 c(0) = tblsize; 00311 c(1) = 3; 00312 c(2) = n; 00313 for(i = 0; i <= n-1; i++) 00314 { 00315 c(3+i) = x(i); 00316 } 00317 for(i = 0; i <= n-2; i++) 00318 { 00319 delta = x(i+1)-x(i); 00320 delta2 = ap::sqr(delta); 00321 delta3 = delta*delta2; 00322 c(3+n+4*i+0) = y(i); 00323 c(3+n+4*i+1) = d(i); 00324 c(3+n+4*i+2) = (3*(y(i+1)-y(i))-2*d(i)*delta-d(i+1)*delta)/delta2; 00325 c(3+n+4*i+3) = (2*(y(i)-y(i+1))+d(i)*delta+d(i+1)*delta)/delta3; 00326 } 00327 } 00328 00329 00330 /************************************************************************* 00331 This subroutine builds Akima spline coefficients table. 00332 00333 Input parameters: 00334 X - spline nodes, array[0..N-1] 00335 Y - function values, array[0..N-1] 00336 N - points count, N>=5 00337 00338 Output parameters: 00339 C - coefficients table. Used by SplineInterpolation and 00340 other subroutines from this file. 00341 00342 -- ALGLIB PROJECT -- 00343 Copyright 24.06.2007 by Bochkanov Sergey 00344 *************************************************************************/ 00345 void buildakimaspline(ap::real_1d_array x, 00346 ap::real_1d_array y, 00347 int n, 00348 ap::real_1d_array& c) 00349 { 00350 int i; 00351 ap::real_1d_array d; 00352 ap::real_1d_array w; 00353 ap::real_1d_array diff; 00354 00355 ap::ap_error::make_assertion(n>=5); 00356 00357 // 00358 // Sort points 00359 // 00360 heapsortpoints(x, y, n); 00361 00362 // 00363 // Prepare W (weights), Diff (divided differences) 00364 // 00365 w.setbounds(1, n-2); 00366 diff.setbounds(0, n-2); 00367 for(i = 0; i <= n-2; i++) 00368 { 00369 diff(i) = (y(i+1)-y(i))/(x(i+1)-x(i)); 00370 } 00371 for(i = 1; i <= n-2; i++) 00372 { 00373 w(i) = fabs(diff(i)-diff(i-1)); 00374 } 00375 00376 // 00377 // Prepare Hermite interpolation scheme 00378 // 00379 d.setbounds(0, n-1); 00380 for(i = 2; i <= n-3; i++) 00381 { 00382 if( fabs(w(i-1))+fabs(w(i+1))!=0 ) 00383 { 00384 d(i) = (w(i+1)*diff(i-1)+w(i-1)*diff(i))/(w(i+1)+w(i-1)); 00385 } 00386 else 00387 { 00388 d(i) = ((x(i+1)-x(i))*diff(i-1)+(x(i)-x(i-1))*diff(i))/(x(i+1)-x(i-1)); 00389 } 00390 } 00391 d(0) = diffthreepoint(x(0), x(0), y(0), x(1), y(1), x(2), y(2)); 00392 d(1) = diffthreepoint(x(1), x(0), y(0), x(1), y(1), x(2), y(2)); 00393 d(n-2) = diffthreepoint(x(n-2), x(n-3), y(n-3), x(n-2), y(n-2), x(n-1), y(n-1)); 00394 d(n-1) = diffthreepoint(x(n-1), x(n-3), y(n-3), x(n-2), y(n-2), x(n-1), y(n-1)); 00395 00396 // 00397 // Build Akima spline using Hermite interpolation scheme 00398 // 00399 buildhermitespline(x, y, d, n, c); 00400 } 00401 00402 00403 /************************************************************************* 00404 This subroutine calculates the value of the spline at the given point X. 00405 00406 Input parameters: 00407 C - coefficients table. Built by BuildLinearSpline, 00408 BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline. 00409 X - point 00410 00411 Result: 00412 S(x) 00413 00414 -- ALGLIB PROJECT -- 00415 Copyright 23.06.2007 by Bochkanov Sergey 00416 *************************************************************************/ 00417 double splineinterpolation(const ap::real_1d_array& c, double x) 00418 { 00419 double result; 00420 int n; 00421 int l; 00422 int r; 00423 int m; 00424 00425 ap::ap_error::make_assertion(ap::round(c(1))==3); 00426 n = ap::round(c(2)); 00427 00428 // 00429 // Binary search in the [ x[0], ..., x[n-2] ] (x[n-1] is not included) 00430 // 00431 l = 3; 00432 r = 3+n-2+1; 00433 while(l!=r-1) 00434 { 00435 m = (l+r)/2; 00436 if( c(m)>=x ) 00437 { 00438 r = m; 00439 } 00440 else 00441 { 00442 l = m; 00443 } 00444 } 00445 00446 // 00447 // Interpolation 00448 // 00449 x = x-c(l); 00450 m = 3+n+4*(l-3); 00451 result = c(m)+x*(c(m+1)+x*(c(m+2)+x*c(m+3))); 00452 return result; 00453 } 00454 00455 00456 /************************************************************************* 00457 This subroutine differentiates the spline. 00458 00459 Input parameters: 00460 C - coefficients table. Built by BuildLinearSpline, 00461 BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline. 00462 X - point 00463 00464 Result: 00465 S - S(x) 00466 DS - S'(x) 00467 D2S - S''(x) 00468 00469 -- ALGLIB PROJECT -- 00470 Copyright 24.06.2007 by Bochkanov Sergey 00471 *************************************************************************/ 00472 void splinedifferentiation(const ap::real_1d_array& c, 00473 double x, 00474 double& s, 00475 double& ds, 00476 double& d2s) 00477 { 00478 int n; 00479 int l; 00480 int r; 00481 int m; 00482 00483 ap::ap_error::make_assertion(ap::round(c(1))==3); 00484 n = ap::round(c(2)); 00485 00486 // 00487 // Binary search 00488 // 00489 l = 3; 00490 r = 3+n-2+1; 00491 while(l!=r-1) 00492 { 00493 m = (l+r)/2; 00494 if( c(m)>=x ) 00495 { 00496 r = m; 00497 } 00498 else 00499 { 00500 l = m; 00501 } 00502 } 00503 00504 // 00505 // Differentiation 00506 // 00507 x = x-c(l); 00508 m = 3+n+4*(l-3); 00509 s = c(m)+x*(c(m+1)+x*(c(m+2)+x*c(m+3))); 00510 ds = c(m+1)+2*x*c(m+2)+3*ap::sqr(x)*c(m+3); 00511 d2s = 2*c(m+2)+6*x*c(m+3); 00512 } 00513 00514 00515 /************************************************************************* 00516 This subroutine makes the copy of the spline. 00517 00518 Input parameters: 00519 C - coefficients table. Built by BuildLinearSpline, 00520 BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline. 00521 00522 Result: 00523 CC - spline copy 00524 00525 -- ALGLIB PROJECT -- 00526 Copyright 29.06.2007 by Bochkanov Sergey 00527 *************************************************************************/ 00528 void splinecopy(const ap::real_1d_array& c, ap::real_1d_array& cc) 00529 { 00530 int s; 00531 00532 s = ap::round(c(0)); 00533 cc.setbounds(0, s-1); 00534 ap::vmove(cc.getvector(0, s-1), c.getvector(0, s-1)); 00535 } 00536 00537 00538 /************************************************************************* 00539 This subroutine unpacks the spline into the coefficients table. 00540 00541 Input parameters: 00542 C - coefficients table. Built by BuildLinearSpline, 00543 BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline. 00544 X - point 00545 00546 Result: 00547 Tbl - coefficients table, unpacked format, array[0..N-2, 0..5]. 00548 For I = 0...N-2: 00549 Tbl[I,0] = X[i] 00550 Tbl[I,1] = X[i+1] 00551 Tbl[I,2] = C0 00552 Tbl[I,3] = C1 00553 Tbl[I,4] = C2 00554 Tbl[I,5] = C3 00555 On [x[i], x[i+1]] spline is equals to: 00556 S(x) = C0 + C1*t + C2*t^2 + C3*t^3 00557 t = x-x[i] 00558 00559 -- ALGLIB PROJECT -- 00560 Copyright 29.06.2007 by Bochkanov Sergey 00561 *************************************************************************/ 00562 void splineunpack(const ap::real_1d_array& c, int& n, ap::real_2d_array& tbl) 00563 { 00564 int i; 00565 00566 ap::ap_error::make_assertion(ap::round(c(1))==3); 00567 n = ap::round(c(2)); 00568 tbl.setbounds(0, n-2, 0, 5); 00569 00570 // 00571 // Fill 00572 // 00573 for(i = 0; i <= n-2; i++) 00574 { 00575 tbl(i,0) = c(3+i); 00576 tbl(i,1) = c(3+i+1); 00577 tbl(i,2) = c(3+n+4*i); 00578 tbl(i,3) = c(3+n+4*i+1); 00579 tbl(i,4) = c(3+n+4*i+2); 00580 tbl(i,5) = c(3+n+4*i+3); 00581 } 00582 } 00583 00584 00585 /************************************************************************* 00586 This subroutine performs linear transformation of the spline argument. 00587 00588 Input parameters: 00589 C - coefficients table. Built by BuildLinearSpline, 00590 BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline. 00591 A, B- transformation coefficients: x = A*t + B 00592 Result: 00593 C - transformed spline 00594 00595 -- ALGLIB PROJECT -- 00596 Copyright 30.06.2007 by Bochkanov Sergey 00597 *************************************************************************/ 00598 void splinelintransx(ap::real_1d_array& c, double a, double b) 00599 { 00600 int i; 00601 int n; 00602 double v; 00603 double dv; 00604 double d2v; 00605 ap::real_1d_array x; 00606 ap::real_1d_array y; 00607 ap::real_1d_array d; 00608 00609 ap::ap_error::make_assertion(ap::round(c(1))==3); 00610 n = ap::round(c(2)); 00611 00612 // 00613 // Special case: A=0 00614 // 00615 if( a==0 ) 00616 { 00617 v = splineinterpolation(c, b); 00618 for(i = 0; i <= n-2; i++) 00619 { 00620 c(3+n+4*i) = v; 00621 c(3+n+4*i+1) = 0; 00622 c(3+n+4*i+2) = 0; 00623 c(3+n+4*i+3) = 0; 00624 } 00625 return; 00626 } 00627 00628 // 00629 // General case: A<>0. 00630 // Unpack, X, Y, dY/dX. 00631 // Scale and pack again. 00632 // 00633 x.setbounds(0, n-1); 00634 y.setbounds(0, n-1); 00635 d.setbounds(0, n-1); 00636 for(i = 0; i <= n-1; i++) 00637 { 00638 x(i) = c(3+i); 00639 splinedifferentiation(c, x(i), v, dv, d2v); 00640 x(i) = (x(i)-b)/a; 00641 y(i) = v; 00642 d(i) = a*dv; 00643 } 00644 buildhermitespline(x, y, d, n, c); 00645 } 00646 00647 00648 /************************************************************************* 00649 This subroutine performs linear transformation of the spline. 00650 00651 Input parameters: 00652 C - coefficients table. Built by BuildLinearSpline, 00653 BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline. 00654 A, B- transformation coefficients: S2(x) = A*S(x) + B 00655 Result: 00656 C - transformed spline 00657 00658 -- ALGLIB PROJECT -- 00659 Copyright 30.06.2007 by Bochkanov Sergey 00660 *************************************************************************/ 00661 void splinelintransy(ap::real_1d_array& c, double a, double b) 00662 { 00663 int i; 00664 int n; 00665 double v; 00666 double dv; 00667 double d2v; 00668 ap::real_1d_array x; 00669 ap::real_1d_array y; 00670 ap::real_1d_array d; 00671 00672 ap::ap_error::make_assertion(ap::round(c(1))==3); 00673 n = ap::round(c(2)); 00674 00675 // 00676 // Special case: A=0 00677 // 00678 for(i = 0; i <= n-2; i++) 00679 { 00680 c(3+n+4*i) = a*c(3+n+4*i)+b; 00681 c(3+n+4*i+1) = a*c(3+n+4*i+1); 00682 c(3+n+4*i+2) = a*c(3+n+4*i+2); 00683 c(3+n+4*i+3) = a*c(3+n+4*i+3); 00684 } 00685 } 00686 00687 00688 /************************************************************************* 00689 This subroutine integrates the spline. 00690 00691 Input parameters: 00692 C - coefficients table. Built by BuildLinearSpline, 00693 BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline. 00694 X - right bound of the integration interval [a, x] 00695 Result: 00696 integral(S(t)dt,a,x) 00697 00698 -- ALGLIB PROJECT -- 00699 Copyright 23.06.2007 by Bochkanov Sergey 00700 *************************************************************************/ 00701 double splineintegration(const ap::real_1d_array& c, double x) 00702 { 00703 double result; 00704 int n; 00705 int i; 00706 int l; 00707 int r; 00708 int m; 00709 double w; 00710 00711 ap::ap_error::make_assertion(ap::round(c(1))==3); 00712 n = ap::round(c(2)); 00713 00714 // 00715 // Binary search in the [ x[0], ..., x[n-2] ] (x[n-1] is not included) 00716 // 00717 l = 3; 00718 r = 3+n-2+1; 00719 while(l!=r-1) 00720 { 00721 m = (l+r)/2; 00722 if( c(m)>=x ) 00723 { 00724 r = m; 00725 } 00726 else 00727 { 00728 l = m; 00729 } 00730 } 00731 00732 // 00733 // Integration 00734 // 00735 result = 0; 00736 for(i = 3; i <= l-1; i++) 00737 { 00738 w = c(i+1)-c(i); 00739 m = 3+n+4*(i-3); 00740 result = result+c(m)*w; 00741 result = result+c(m+1)*ap::sqr(w)/2; 00742 result = result+c(m+2)*ap::sqr(w)*w/3; 00743 result = result+c(m+3)*ap::sqr(ap::sqr(w))/4; 00744 } 00745 w = x-c(l); 00746 m = 3+n+4*(l-3); 00747 result = result+c(m)*w; 00748 result = result+c(m+1)*ap::sqr(w)/2; 00749 result = result+c(m+2)*ap::sqr(w)*w/3; 00750 result = result+c(m+3)*ap::sqr(ap::sqr(w))/4; 00751 return result; 00752 } 00753 00754 00755 /************************************************************************* 00756 Obsolete subroutine, left for backward compatibility. 00757 *************************************************************************/ 00758 void spline3buildtable(int n, 00759 const int& diffn, 00760 ap::real_1d_array x, 00761 ap::real_1d_array y, 00762 const double& boundl, 00763 const double& boundr, 00764 ap::real_2d_array& ctbl) 00765 { 00766 bool c; 00767 int e; 00768 int g; 00769 double tmp; 00770 int nxm1; 00771 int i; 00772 int j; 00773 double dx; 00774 double dxj; 00775 double dyj; 00776 double dxjp1; 00777 double dyjp1; 00778 double dxp; 00779 double dyp; 00780 double yppa; 00781 double yppb; 00782 double pj; 00783 double b1; 00784 double b2; 00785 double b3; 00786 double b4; 00787 00788 n = n-1; 00789 g = (n+1)/2; 00790 do 00791 { 00792 i = g; 00793 do 00794 { 00795 j = i-g; 00796 c = true; 00797 do 00798 { 00799 if( x(j)<=x(j+g) ) 00800 { 00801 c = false; 00802 } 00803 else 00804 { 00805 tmp = x(j); 00806 x(j) = x(j+g); 00807 x(j+g) = tmp; 00808 tmp = y(j); 00809 y(j) = y(j+g); 00810 y(j+g) = tmp; 00811 } 00812 j = j-1; 00813 } 00814 while(j>=0&&c); 00815 i = i+1; 00816 } 00817 while(i<=n); 00818 g = g/2; 00819 } 00820 while(g>0); 00821 ctbl.setbounds(0, 4, 0, n); 00822 n = n+1; 00823 if( diffn==1 ) 00824 { 00825 b1 = 1; 00826 b2 = 6/(x(1)-x(0))*((y(1)-y(0))/(x(1)-x(0))-boundl); 00827 b3 = 1; 00828 b4 = 6/(x(n-1)-x(n-2))*(boundr-(y(n-1)-y(n-2))/(x(n-1)-x(n-2))); 00829 } 00830 else 00831 { 00832 b1 = 0; 00833 b2 = 2*boundl; 00834 b3 = 0; 00835 b4 = 2*boundr; 00836 } 00837 nxm1 = n-1; 00838 if( n>=2 ) 00839 { 00840 if( n>2 ) 00841 { 00842 dxj = x(1)-x(0); 00843 dyj = y(1)-y(0); 00844 j = 2; 00845 while(j<=nxm1) 00846 { 00847 dxjp1 = x(j)-x(j-1); 00848 dyjp1 = y(j)-y(j-1); 00849 dxp = dxj+dxjp1; 00850 ctbl(1,j-1) = dxjp1/dxp; 00851 ctbl(2,j-1) = 1-ctbl(1,j-1); 00852 ctbl(3,j-1) = 6*(dyjp1/dxjp1-dyj/dxj)/dxp; 00853 dxj = dxjp1; 00854 dyj = dyjp1; 00855 j = j+1; 00856 } 00857 } 00858 ctbl(1,0) = -b1/2; 00859 ctbl(2,0) = b2/2; 00860 if( n!=2 ) 00861 { 00862 j = 2; 00863 while(j<=nxm1) 00864 { 00865 pj = ctbl(2,j-1)*ctbl(1,j-2)+2; 00866 ctbl(1,j-1) = -ctbl(1,j-1)/pj; 00867 ctbl(2,j-1) = (ctbl(3,j-1)-ctbl(2,j-1)*ctbl(2,j-2))/pj; 00868 j = j+1; 00869 } 00870 } 00871 yppb = (b4-b3*ctbl(2,nxm1-1))/(b3*ctbl(1,nxm1-1)+2); 00872 i = 1; 00873 while(i<=nxm1) 00874 { 00875 j = n-i; 00876 yppa = ctbl(1,j-1)*yppb+ctbl(2,j-1); 00877 dx = x(j)-x(j-1); 00878 ctbl(3,j-1) = (yppb-yppa)/dx/6; 00879 ctbl(2,j-1) = yppa/2; 00880 ctbl(1,j-1) = (y(j)-y(j-1))/dx-(ctbl(2,j-1)+ctbl(3,j-1)*dx)*dx; 00881 yppb = yppa; 00882 i = i+1; 00883 } 00884 for(i = 1; i <= n; i++) 00885 { 00886 ctbl(0,i-1) = y(i-1); 00887 ctbl(4,i-1) = x(i-1); 00888 } 00889 } 00890 } 00891 00892 00893 /************************************************************************* 00894 Obsolete subroutine, left for backward compatibility. 00895 *************************************************************************/ 00896 double spline3interpolate(int n, const ap::real_2d_array& c, const double& x) 00897 { 00898 double result; 00899 int i; 00900 int l; 00901 int half; 00902 int first; 00903 int middle; 00904 00905 n = n-1; 00906 l = n; 00907 first = 0; 00908 while(l>0) 00909 { 00910 half = l/2; 00911 middle = first+half; 00912 if( c(4,middle)x(i-1); 00954 isdescending = isdescending&&x(i)=x(t-1) ) 00998 { 00999 t = 1; 01000 } 01001 else 01002 { 01003 tmp = x(k-1); 01004 x(k-1) = x(t-1); 01005 x(t-1) = tmp; 01006 tmp = y(k-1); 01007 y(k-1) = y(t-1); 01008 y(t-1) = tmp; 01009 t = k; 01010 } 01011 } 01012 i = i+1; 01013 } 01014 while(i<=n); 01015 i = n-1; 01016 do 01017 { 01018 tmp = x(i); 01019 x(i) = x(0); 01020 x(0) = tmp; 01021 tmp = y(i); 01022 y(i) = y(0); 01023 y(0) = tmp; 01024 t = 1; 01025 while(t!=0) 01026 { 01027 k = 2*t; 01028 if( k>i ) 01029 { 01030 t = 0; 01031 } 01032 else 01033 { 01034 if( kx(k-1) ) 01037 { 01038 k = k+1; 01039 } 01040 } 01041 if( x(t-1)>=x(k-1) ) 01042 { 01043 t = 0; 01044 } 01045 else 01046 { 01047 tmp = x(k-1); 01048 x(k-1) = x(t-1); 01049 x(t-1) = tmp; 01050 tmp = y(k-1); 01051 y(k-1) = y(t-1); 01052 y(t-1) = tmp; 01053 t = k; 01054 } 01055 } 01056 } 01057 i = i-1; 01058 } 01059 while(i>=1); 01060 } 01061 01062 01063 /************************************************************************* 01064 Internal subroutine. Heap sort. 01065 *************************************************************************/ 01066 void heapsortdpoints(ap::real_1d_array& x, 01067 ap::real_1d_array& y, 01068 ap::real_1d_array& d, 01069 int n) 01070 { 01071 int i; 01072 int j; 01073 int k; 01074 int t; 01075 double tmp; 01076 bool isascending; 01077 bool isdescending; 01078 01079 01080 // 01081 // Test for already sorted set 01082 // 01083 isascending = true; 01084 isdescending = true; 01085 for(i = 1; i <= n-1; i++) 01086 { 01087 isascending = isascending&&x(i)>x(i-1); 01088 isdescending = isdescending&&x(i)=x(t-1) ) 01135 { 01136 t = 1; 01137 } 01138 else 01139 { 01140 tmp = x(k-1); 01141 x(k-1) = x(t-1); 01142 x(t-1) = tmp; 01143 tmp = y(k-1); 01144 y(k-1) = y(t-1); 01145 y(t-1) = tmp; 01146 tmp = d(k-1); 01147 d(k-1) = d(t-1); 01148 d(t-1) = tmp; 01149 t = k; 01150 } 01151 } 01152 i = i+1; 01153 } 01154 while(i<=n); 01155 i = n-1; 01156 do 01157 { 01158 tmp = x(i); 01159 x(i) = x(0); 01160 x(0) = tmp; 01161 tmp = y(i); 01162 y(i) = y(0); 01163 y(0) = tmp; 01164 tmp = d(i); 01165 d(i) = d(0); 01166 d(0) = tmp; 01167 t = 1; 01168 while(t!=0) 01169 { 01170 k = 2*t; 01171 if( k>i ) 01172 { 01173 t = 0; 01174 } 01175 else 01176 { 01177 if( kx(k-1) ) 01180 { 01181 k = k+1; 01182 } 01183 } 01184 if( x(t-1)>=x(k-1) ) 01185 { 01186 t = 0; 01187 } 01188 else 01189 { 01190 tmp = x(k-1); 01191 x(k-1) = x(t-1); 01192 x(t-1) = tmp; 01193 tmp = y(k-1); 01194 y(k-1) = y(t-1); 01195 y(t-1) = tmp; 01196 tmp = d(k-1); 01197 d(k-1) = d(t-1); 01198 d(t-1) = tmp; 01199 t = k; 01200 } 01201 } 01202 } 01203 i = i-1; 01204 } 01205 while(i>=1); 01206 } 01207 01208 01209 /************************************************************************* 01210 Internal subroutine. Tridiagonal solver. 01211 *************************************************************************/ 01212 void solvetridiagonal(ap::real_1d_array a, 01213 ap::real_1d_array b, 01214 ap::real_1d_array c, 01215 ap::real_1d_array d, 01216 int n, 01217 ap::real_1d_array& x) 01218 { 01219 int k; 01220 double t; 01221 01222 x.setbounds(0, n-1); 01223 a(0) = 0; 01224 c(n-1) = 0; 01225 for(k = 1; k <= n-1; k++) 01226 { 01227 t = a(k)/b(k-1); 01228 b(k) = b(k)-t*c(k-1); 01229 d(k) = d(k)-t*d(k-1); 01230 } 01231 x(n-1) = d(n-1)/b(n-1); 01232 for(k = n-2; k >= 0; k--) 01233 { 01234 x(k) = (d(k)-c(k)*x(k+1))/b(k); 01235 } 01236 } 01237 01238 01239 /************************************************************************* 01240 Internal subroutine. Three-point differentiation 01241 *************************************************************************/ 01242 double diffthreepoint(double t, 01243 double x0, 01244 double f0, 01245 double x1, 01246 double f1, 01247 double x2, 01248 double f2) 01249 { 01250 double result; 01251 double a; 01252 double b; 01253 01254 t = t-x0; 01255 x1 = x1-x0; 01256 x2 = x2-x0; 01257 a = (f2-f0-x2/x1*(f1-f0))/(ap::sqr(x2)-x1*x2); 01258 b = (f1-f0-a*ap::sqr(x1))/x1; 01259 result = 2*a*t+b; 01260 return result; 01261 } 01262 01263 01264