1 | /*************************************************************************
|
---|
2 | Copyright (c) 2007, Sergey Bochkanov (ALGLIB project).
|
---|
3 |
|
---|
4 | >>> SOURCE LICENSE >>>
|
---|
5 | This program is free software; you can redistribute it and/or modify
|
---|
6 | it under the terms of the GNU General Public License as published by
|
---|
7 | the Free Software Foundation (www.fsf.org); either version 2 of the
|
---|
8 | License, or (at your option) any later version.
|
---|
9 |
|
---|
10 | This program is distributed in the hope that it will be useful,
|
---|
11 | but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
12 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
13 | GNU General Public License for more details.
|
---|
14 |
|
---|
15 | A copy of the GNU General Public License is available at
|
---|
16 | http://www.fsf.org/licensing/licenses
|
---|
17 |
|
---|
18 | >>> END OF LICENSE >>>
|
---|
19 | *************************************************************************/
|
---|
20 |
|
---|
21 | using System;
|
---|
22 |
|
---|
23 | namespace alglib
|
---|
24 | {
|
---|
25 | public class spline3
|
---|
26 | {
|
---|
27 | public static void buildlinearspline(double[] x,
|
---|
28 | double[] y,
|
---|
29 | int n,
|
---|
30 | ref double[] c)
|
---|
31 | {
|
---|
32 | int i = 0;
|
---|
33 | int tblsize = 0;
|
---|
34 |
|
---|
35 | x = (double[])x.Clone();
|
---|
36 | y = (double[])y.Clone();
|
---|
37 |
|
---|
38 | System.Diagnostics.Debug.Assert(n>=2, "BuildLinearSpline: N<2!");
|
---|
39 |
|
---|
40 | //
|
---|
41 | // Sort points
|
---|
42 | //
|
---|
43 | heapsortpoints(ref x, ref y, n);
|
---|
44 |
|
---|
45 | //
|
---|
46 | // Fill C:
|
---|
47 | // C[0] - length(C)
|
---|
48 | // C[1] - type(C):
|
---|
49 | // 3 - general cubic spline
|
---|
50 | // C[2] - N
|
---|
51 | // C[3]...C[3+N-1] - x[i], i = 0...N-1
|
---|
52 | // C[3+N]...C[3+N+(N-1)*4-1] - coefficients table
|
---|
53 | //
|
---|
54 | tblsize = 3+n+(n-1)*4;
|
---|
55 | c = new double[tblsize-1+1];
|
---|
56 | c[0] = tblsize;
|
---|
57 | c[1] = 3;
|
---|
58 | c[2] = n;
|
---|
59 | for(i=0; i<=n-1; i++)
|
---|
60 | {
|
---|
61 | c[3+i] = x[i];
|
---|
62 | }
|
---|
63 | for(i=0; i<=n-2; i++)
|
---|
64 | {
|
---|
65 | c[3+n+4*i+0] = y[i];
|
---|
66 | c[3+n+4*i+1] = (y[i+1]-y[i])/(x[i+1]-x[i]);
|
---|
67 | c[3+n+4*i+2] = 0;
|
---|
68 | c[3+n+4*i+3] = 0;
|
---|
69 | }
|
---|
70 | }
|
---|
71 |
|
---|
72 |
|
---|
73 | public static void buildcubicspline(double[] x,
|
---|
74 | double[] y,
|
---|
75 | int n,
|
---|
76 | int boundltype,
|
---|
77 | double boundl,
|
---|
78 | int boundrtype,
|
---|
79 | double boundr,
|
---|
80 | ref double[] c)
|
---|
81 | {
|
---|
82 | double[] a1 = new double[0];
|
---|
83 | double[] a2 = new double[0];
|
---|
84 | double[] a3 = new double[0];
|
---|
85 | double[] b = new double[0];
|
---|
86 | double[] d = new double[0];
|
---|
87 | int i = 0;
|
---|
88 | int tblsize = 0;
|
---|
89 | double delta = 0;
|
---|
90 | double delta2 = 0;
|
---|
91 | double delta3 = 0;
|
---|
92 |
|
---|
93 | x = (double[])x.Clone();
|
---|
94 | y = (double[])y.Clone();
|
---|
95 |
|
---|
96 | System.Diagnostics.Debug.Assert(n>=2, "BuildCubicSpline: N<2!");
|
---|
97 | System.Diagnostics.Debug.Assert(boundltype==0 | boundltype==1 | boundltype==2, "BuildCubicSpline: incorrect BoundLType!");
|
---|
98 | System.Diagnostics.Debug.Assert(boundrtype==0 | boundrtype==1 | boundrtype==2, "BuildCubicSpline: incorrect BoundRType!");
|
---|
99 | a1 = new double[n-1+1];
|
---|
100 | a2 = new double[n-1+1];
|
---|
101 | a3 = new double[n-1+1];
|
---|
102 | b = new double[n-1+1];
|
---|
103 |
|
---|
104 | //
|
---|
105 | // Special case:
|
---|
106 | // * N=2
|
---|
107 | // * parabolic terminated boundary condition on both ends
|
---|
108 | //
|
---|
109 | if( n==2 & boundltype==0 & boundrtype==0 )
|
---|
110 | {
|
---|
111 |
|
---|
112 | //
|
---|
113 | // Change task type
|
---|
114 | //
|
---|
115 | boundltype = 2;
|
---|
116 | boundl = 0;
|
---|
117 | boundrtype = 2;
|
---|
118 | boundr = 0;
|
---|
119 | }
|
---|
120 |
|
---|
121 | //
|
---|
122 | //
|
---|
123 | // Sort points
|
---|
124 | //
|
---|
125 | heapsortpoints(ref x, ref y, n);
|
---|
126 |
|
---|
127 | //
|
---|
128 | // Left boundary conditions
|
---|
129 | //
|
---|
130 | if( boundltype==0 )
|
---|
131 | {
|
---|
132 | a1[0] = 0;
|
---|
133 | a2[0] = 1;
|
---|
134 | a3[0] = 1;
|
---|
135 | b[0] = 2*(y[1]-y[0])/(x[1]-x[0]);
|
---|
136 | }
|
---|
137 | if( boundltype==1 )
|
---|
138 | {
|
---|
139 | a1[0] = 0;
|
---|
140 | a2[0] = 1;
|
---|
141 | a3[0] = 0;
|
---|
142 | b[0] = boundl;
|
---|
143 | }
|
---|
144 | if( boundltype==2 )
|
---|
145 | {
|
---|
146 | a1[0] = 0;
|
---|
147 | a2[0] = 2;
|
---|
148 | a3[0] = 1;
|
---|
149 | b[0] = 3*(y[1]-y[0])/(x[1]-x[0])-0.5*boundl*(x[1]-x[0]);
|
---|
150 | }
|
---|
151 |
|
---|
152 | //
|
---|
153 | // Central conditions
|
---|
154 | //
|
---|
155 | for(i=1; i<=n-2; i++)
|
---|
156 | {
|
---|
157 | a1[i] = x[i+1]-x[i];
|
---|
158 | a2[i] = 2*(x[i+1]-x[i-1]);
|
---|
159 | a3[i] = x[i]-x[i-1];
|
---|
160 | 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]);
|
---|
161 | }
|
---|
162 |
|
---|
163 | //
|
---|
164 | // Right boundary conditions
|
---|
165 | //
|
---|
166 | if( boundrtype==0 )
|
---|
167 | {
|
---|
168 | a1[n-1] = 1;
|
---|
169 | a2[n-1] = 1;
|
---|
170 | a3[n-1] = 0;
|
---|
171 | b[n-1] = 2*(y[n-1]-y[n-2])/(x[n-1]-x[n-2]);
|
---|
172 | }
|
---|
173 | if( boundrtype==1 )
|
---|
174 | {
|
---|
175 | a1[n-1] = 0;
|
---|
176 | a2[n-1] = 1;
|
---|
177 | a3[n-1] = 0;
|
---|
178 | b[n-1] = boundr;
|
---|
179 | }
|
---|
180 | if( boundrtype==2 )
|
---|
181 | {
|
---|
182 | a1[n-1] = 1;
|
---|
183 | a2[n-1] = 2;
|
---|
184 | a3[n-1] = 0;
|
---|
185 | 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]);
|
---|
186 | }
|
---|
187 |
|
---|
188 | //
|
---|
189 | // Solve
|
---|
190 | //
|
---|
191 | solvetridiagonal(a1, a2, a3, b, n, ref d);
|
---|
192 |
|
---|
193 | //
|
---|
194 | // Now problem is reduced to the cubic Hermite spline
|
---|
195 | //
|
---|
196 | buildhermitespline(x, y, d, n, ref c);
|
---|
197 | }
|
---|
198 |
|
---|
199 |
|
---|
200 | public static void buildhermitespline(double[] x,
|
---|
201 | double[] y,
|
---|
202 | double[] d,
|
---|
203 | int n,
|
---|
204 | ref double[] c)
|
---|
205 | {
|
---|
206 | int i = 0;
|
---|
207 | int tblsize = 0;
|
---|
208 | double delta = 0;
|
---|
209 | double delta2 = 0;
|
---|
210 | double delta3 = 0;
|
---|
211 |
|
---|
212 | x = (double[])x.Clone();
|
---|
213 | y = (double[])y.Clone();
|
---|
214 | d = (double[])d.Clone();
|
---|
215 |
|
---|
216 | System.Diagnostics.Debug.Assert(n>=2, "BuildHermiteSpline: N<2!");
|
---|
217 |
|
---|
218 | //
|
---|
219 | // Sort points
|
---|
220 | //
|
---|
221 | heapsortdpoints(ref x, ref y, ref d, n);
|
---|
222 |
|
---|
223 | //
|
---|
224 | // Fill C:
|
---|
225 | // C[0] - length(C)
|
---|
226 | // C[1] - type(C):
|
---|
227 | // 3 - general cubic spline
|
---|
228 | // C[2] - N
|
---|
229 | // C[3]...C[3+N-1] - x[i], i = 0...N-1
|
---|
230 | // C[3+N]...C[3+N+(N-1)*4-1] - coefficients table
|
---|
231 | //
|
---|
232 | tblsize = 3+n+(n-1)*4;
|
---|
233 | c = new double[tblsize-1+1];
|
---|
234 | c[0] = tblsize;
|
---|
235 | c[1] = 3;
|
---|
236 | c[2] = n;
|
---|
237 | for(i=0; i<=n-1; i++)
|
---|
238 | {
|
---|
239 | c[3+i] = x[i];
|
---|
240 | }
|
---|
241 | for(i=0; i<=n-2; i++)
|
---|
242 | {
|
---|
243 | delta = x[i+1]-x[i];
|
---|
244 | delta2 = AP.Math.Sqr(delta);
|
---|
245 | delta3 = delta*delta2;
|
---|
246 | c[3+n+4*i+0] = y[i];
|
---|
247 | c[3+n+4*i+1] = d[i];
|
---|
248 | c[3+n+4*i+2] = (3*(y[i+1]-y[i])-2*d[i]*delta-d[i+1]*delta)/delta2;
|
---|
249 | c[3+n+4*i+3] = (2*(y[i]-y[i+1])+d[i]*delta+d[i+1]*delta)/delta3;
|
---|
250 | }
|
---|
251 | }
|
---|
252 |
|
---|
253 |
|
---|
254 | public static void buildakimaspline(double[] x,
|
---|
255 | double[] y,
|
---|
256 | int n,
|
---|
257 | ref double[] c)
|
---|
258 | {
|
---|
259 | int i = 0;
|
---|
260 | double[] d = new double[0];
|
---|
261 | double[] w = new double[0];
|
---|
262 | double[] diff = new double[0];
|
---|
263 |
|
---|
264 | x = (double[])x.Clone();
|
---|
265 | y = (double[])y.Clone();
|
---|
266 |
|
---|
267 | System.Diagnostics.Debug.Assert(n>=5, "BuildAkimaSpline: N<5!");
|
---|
268 |
|
---|
269 | //
|
---|
270 | // Sort points
|
---|
271 | //
|
---|
272 | heapsortpoints(ref x, ref y, n);
|
---|
273 |
|
---|
274 | //
|
---|
275 | // Prepare W (weights), Diff (divided differences)
|
---|
276 | //
|
---|
277 | w = new double[n-2+1];
|
---|
278 | diff = new double[n-2+1];
|
---|
279 | for(i=0; i<=n-2; i++)
|
---|
280 | {
|
---|
281 | diff[i] = (y[i+1]-y[i])/(x[i+1]-x[i]);
|
---|
282 | }
|
---|
283 | for(i=1; i<=n-2; i++)
|
---|
284 | {
|
---|
285 | w[i] = Math.Abs(diff[i]-diff[i-1]);
|
---|
286 | }
|
---|
287 |
|
---|
288 | //
|
---|
289 | // Prepare Hermite interpolation scheme
|
---|
290 | //
|
---|
291 | d = new double[n-1+1];
|
---|
292 | for(i=2; i<=n-3; i++)
|
---|
293 | {
|
---|
294 | if( (double)(Math.Abs(w[i-1])+Math.Abs(w[i+1]))!=(double)(0) )
|
---|
295 | {
|
---|
296 | d[i] = (w[i+1]*diff[i-1]+w[i-1]*diff[i])/(w[i+1]+w[i-1]);
|
---|
297 | }
|
---|
298 | else
|
---|
299 | {
|
---|
300 | d[i] = ((x[i+1]-x[i])*diff[i-1]+(x[i]-x[i-1])*diff[i])/(x[i+1]-x[i-1]);
|
---|
301 | }
|
---|
302 | }
|
---|
303 | d[0] = diffthreepoint(x[0], x[0], y[0], x[1], y[1], x[2], y[2]);
|
---|
304 | d[1] = diffthreepoint(x[1], x[0], y[0], x[1], y[1], x[2], y[2]);
|
---|
305 | 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]);
|
---|
306 | 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]);
|
---|
307 |
|
---|
308 | //
|
---|
309 | // Build Akima spline using Hermite interpolation scheme
|
---|
310 | //
|
---|
311 | buildhermitespline(x, y, d, n, ref c);
|
---|
312 | }
|
---|
313 |
|
---|
314 |
|
---|
315 | public static double splineinterpolation(ref double[] c,
|
---|
316 | double x)
|
---|
317 | {
|
---|
318 | double result = 0;
|
---|
319 | int n = 0;
|
---|
320 | int l = 0;
|
---|
321 | int r = 0;
|
---|
322 | int m = 0;
|
---|
323 |
|
---|
324 | System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineInterpolation: incorrect C!");
|
---|
325 | n = (int)Math.Round(c[2]);
|
---|
326 |
|
---|
327 | //
|
---|
328 | // Binary search in the [ x[0], ..., x[n-2] ] (x[n-1] is not included)
|
---|
329 | //
|
---|
330 | l = 3;
|
---|
331 | r = 3+n-2+1;
|
---|
332 | while( l!=r-1 )
|
---|
333 | {
|
---|
334 | m = (l+r)/2;
|
---|
335 | if( (double)(c[m])>=(double)(x) )
|
---|
336 | {
|
---|
337 | r = m;
|
---|
338 | }
|
---|
339 | else
|
---|
340 | {
|
---|
341 | l = m;
|
---|
342 | }
|
---|
343 | }
|
---|
344 |
|
---|
345 | //
|
---|
346 | // Interpolation
|
---|
347 | //
|
---|
348 | x = x-c[l];
|
---|
349 | m = 3+n+4*(l-3);
|
---|
350 | result = c[m]+x*(c[m+1]+x*(c[m+2]+x*c[m+3]));
|
---|
351 | return result;
|
---|
352 | }
|
---|
353 |
|
---|
354 |
|
---|
355 | public static void splinedifferentiation(ref double[] c,
|
---|
356 | double x,
|
---|
357 | ref double s,
|
---|
358 | ref double ds,
|
---|
359 | ref double d2s)
|
---|
360 | {
|
---|
361 | int n = 0;
|
---|
362 | int l = 0;
|
---|
363 | int r = 0;
|
---|
364 | int m = 0;
|
---|
365 |
|
---|
366 | System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineInterpolation: incorrect C!");
|
---|
367 | n = (int)Math.Round(c[2]);
|
---|
368 |
|
---|
369 | //
|
---|
370 | // Binary search
|
---|
371 | //
|
---|
372 | l = 3;
|
---|
373 | r = 3+n-2+1;
|
---|
374 | while( l!=r-1 )
|
---|
375 | {
|
---|
376 | m = (l+r)/2;
|
---|
377 | if( (double)(c[m])>=(double)(x) )
|
---|
378 | {
|
---|
379 | r = m;
|
---|
380 | }
|
---|
381 | else
|
---|
382 | {
|
---|
383 | l = m;
|
---|
384 | }
|
---|
385 | }
|
---|
386 |
|
---|
387 | //
|
---|
388 | // Differentiation
|
---|
389 | //
|
---|
390 | x = x-c[l];
|
---|
391 | m = 3+n+4*(l-3);
|
---|
392 | s = c[m]+x*(c[m+1]+x*(c[m+2]+x*c[m+3]));
|
---|
393 | ds = c[m+1]+2*x*c[m+2]+3*AP.Math.Sqr(x)*c[m+3];
|
---|
394 | d2s = 2*c[m+2]+6*x*c[m+3];
|
---|
395 | }
|
---|
396 |
|
---|
397 |
|
---|
398 | public static void splinecopy(ref double[] c,
|
---|
399 | ref double[] cc)
|
---|
400 | {
|
---|
401 | int s = 0;
|
---|
402 | int i_ = 0;
|
---|
403 |
|
---|
404 | s = (int)Math.Round(c[0]);
|
---|
405 | cc = new double[s-1+1];
|
---|
406 | for(i_=0; i_<=s-1;i_++)
|
---|
407 | {
|
---|
408 | cc[i_] = c[i_];
|
---|
409 | }
|
---|
410 | }
|
---|
411 |
|
---|
412 |
|
---|
413 | public static void splineunpack(ref double[] c,
|
---|
414 | ref int n,
|
---|
415 | ref double[,] tbl)
|
---|
416 | {
|
---|
417 | int i = 0;
|
---|
418 |
|
---|
419 | System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineUnpack: incorrect C!");
|
---|
420 | n = (int)Math.Round(c[2]);
|
---|
421 | tbl = new double[n-2+1, 5+1];
|
---|
422 |
|
---|
423 | //
|
---|
424 | // Fill
|
---|
425 | //
|
---|
426 | for(i=0; i<=n-2; i++)
|
---|
427 | {
|
---|
428 | tbl[i,0] = c[3+i];
|
---|
429 | tbl[i,1] = c[3+i+1];
|
---|
430 | tbl[i,2] = c[3+n+4*i];
|
---|
431 | tbl[i,3] = c[3+n+4*i+1];
|
---|
432 | tbl[i,4] = c[3+n+4*i+2];
|
---|
433 | tbl[i,5] = c[3+n+4*i+3];
|
---|
434 | }
|
---|
435 | }
|
---|
436 |
|
---|
437 |
|
---|
438 | public static void splinelintransx(ref double[] c,
|
---|
439 | double a,
|
---|
440 | double b)
|
---|
441 | {
|
---|
442 | int i = 0;
|
---|
443 | int n = 0;
|
---|
444 | double v = 0;
|
---|
445 | double dv = 0;
|
---|
446 | double d2v = 0;
|
---|
447 | double[] x = new double[0];
|
---|
448 | double[] y = new double[0];
|
---|
449 | double[] d = new double[0];
|
---|
450 |
|
---|
451 | System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineLinTransX: incorrect C!");
|
---|
452 | n = (int)Math.Round(c[2]);
|
---|
453 |
|
---|
454 | //
|
---|
455 | // Special case: A=0
|
---|
456 | //
|
---|
457 | if( (double)(a)==(double)(0) )
|
---|
458 | {
|
---|
459 | v = splineinterpolation(ref c, b);
|
---|
460 | for(i=0; i<=n-2; i++)
|
---|
461 | {
|
---|
462 | c[3+n+4*i] = v;
|
---|
463 | c[3+n+4*i+1] = 0;
|
---|
464 | c[3+n+4*i+2] = 0;
|
---|
465 | c[3+n+4*i+3] = 0;
|
---|
466 | }
|
---|
467 | return;
|
---|
468 | }
|
---|
469 |
|
---|
470 | //
|
---|
471 | // General case: A<>0.
|
---|
472 | // Unpack, X, Y, dY/dX.
|
---|
473 | // Scale and pack again.
|
---|
474 | //
|
---|
475 | x = new double[n-1+1];
|
---|
476 | y = new double[n-1+1];
|
---|
477 | d = new double[n-1+1];
|
---|
478 | for(i=0; i<=n-1; i++)
|
---|
479 | {
|
---|
480 | x[i] = c[3+i];
|
---|
481 | splinedifferentiation(ref c, x[i], ref v, ref dv, ref d2v);
|
---|
482 | x[i] = (x[i]-b)/a;
|
---|
483 | y[i] = v;
|
---|
484 | d[i] = a*dv;
|
---|
485 | }
|
---|
486 | buildhermitespline(x, y, d, n, ref c);
|
---|
487 | }
|
---|
488 |
|
---|
489 |
|
---|
490 | public static void splinelintransy(ref double[] c,
|
---|
491 | double a,
|
---|
492 | double b)
|
---|
493 | {
|
---|
494 | int i = 0;
|
---|
495 | int n = 0;
|
---|
496 | double v = 0;
|
---|
497 | double dv = 0;
|
---|
498 | double d2v = 0;
|
---|
499 | double[] x = new double[0];
|
---|
500 | double[] y = new double[0];
|
---|
501 | double[] d = new double[0];
|
---|
502 |
|
---|
503 | System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineLinTransX: incorrect C!");
|
---|
504 | n = (int)Math.Round(c[2]);
|
---|
505 |
|
---|
506 | //
|
---|
507 | // Special case: A=0
|
---|
508 | //
|
---|
509 | for(i=0; i<=n-2; i++)
|
---|
510 | {
|
---|
511 | c[3+n+4*i] = a*c[3+n+4*i]+b;
|
---|
512 | c[3+n+4*i+1] = a*c[3+n+4*i+1];
|
---|
513 | c[3+n+4*i+2] = a*c[3+n+4*i+2];
|
---|
514 | c[3+n+4*i+3] = a*c[3+n+4*i+3];
|
---|
515 | }
|
---|
516 | }
|
---|
517 |
|
---|
518 |
|
---|
519 | public static double splineintegration(ref double[] c,
|
---|
520 | double x)
|
---|
521 | {
|
---|
522 | double result = 0;
|
---|
523 | int n = 0;
|
---|
524 | int i = 0;
|
---|
525 | int l = 0;
|
---|
526 | int r = 0;
|
---|
527 | int m = 0;
|
---|
528 | double w = 0;
|
---|
529 |
|
---|
530 | System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineIntegration: incorrect C!");
|
---|
531 | n = (int)Math.Round(c[2]);
|
---|
532 |
|
---|
533 | //
|
---|
534 | // Binary search in the [ x[0], ..., x[n-2] ] (x[n-1] is not included)
|
---|
535 | //
|
---|
536 | l = 3;
|
---|
537 | r = 3+n-2+1;
|
---|
538 | while( l!=r-1 )
|
---|
539 | {
|
---|
540 | m = (l+r)/2;
|
---|
541 | if( (double)(c[m])>=(double)(x) )
|
---|
542 | {
|
---|
543 | r = m;
|
---|
544 | }
|
---|
545 | else
|
---|
546 | {
|
---|
547 | l = m;
|
---|
548 | }
|
---|
549 | }
|
---|
550 |
|
---|
551 | //
|
---|
552 | // Integration
|
---|
553 | //
|
---|
554 | result = 0;
|
---|
555 | for(i=3; i<=l-1; i++)
|
---|
556 | {
|
---|
557 | w = c[i+1]-c[i];
|
---|
558 | m = 3+n+4*(i-3);
|
---|
559 | result = result+c[m]*w;
|
---|
560 | result = result+c[m+1]*AP.Math.Sqr(w)/2;
|
---|
561 | result = result+c[m+2]*AP.Math.Sqr(w)*w/3;
|
---|
562 | result = result+c[m+3]*AP.Math.Sqr(AP.Math.Sqr(w))/4;
|
---|
563 | }
|
---|
564 | w = x-c[l];
|
---|
565 | m = 3+n+4*(l-3);
|
---|
566 | result = result+c[m]*w;
|
---|
567 | result = result+c[m+1]*AP.Math.Sqr(w)/2;
|
---|
568 | result = result+c[m+2]*AP.Math.Sqr(w)*w/3;
|
---|
569 | result = result+c[m+3]*AP.Math.Sqr(AP.Math.Sqr(w))/4;
|
---|
570 | return result;
|
---|
571 | }
|
---|
572 |
|
---|
573 |
|
---|
574 | public static void spline3buildtable(int n,
|
---|
575 | int diffn,
|
---|
576 | double[] x,
|
---|
577 | double[] y,
|
---|
578 | double boundl,
|
---|
579 | double boundr,
|
---|
580 | ref double[,] ctbl)
|
---|
581 | {
|
---|
582 | bool c = new bool();
|
---|
583 | int e = 0;
|
---|
584 | int g = 0;
|
---|
585 | double tmp = 0;
|
---|
586 | int nxm1 = 0;
|
---|
587 | int i = 0;
|
---|
588 | int j = 0;
|
---|
589 | double dx = 0;
|
---|
590 | double dxj = 0;
|
---|
591 | double dyj = 0;
|
---|
592 | double dxjp1 = 0;
|
---|
593 | double dyjp1 = 0;
|
---|
594 | double dxp = 0;
|
---|
595 | double dyp = 0;
|
---|
596 | double yppa = 0;
|
---|
597 | double yppb = 0;
|
---|
598 | double pj = 0;
|
---|
599 | double b1 = 0;
|
---|
600 | double b2 = 0;
|
---|
601 | double b3 = 0;
|
---|
602 | double b4 = 0;
|
---|
603 |
|
---|
604 | x = (double[])x.Clone();
|
---|
605 | y = (double[])y.Clone();
|
---|
606 |
|
---|
607 | n = n-1;
|
---|
608 | g = (n+1)/2;
|
---|
609 | do
|
---|
610 | {
|
---|
611 | i = g;
|
---|
612 | do
|
---|
613 | {
|
---|
614 | j = i-g;
|
---|
615 | c = true;
|
---|
616 | do
|
---|
617 | {
|
---|
618 | if( (double)(x[j])<=(double)(x[j+g]) )
|
---|
619 | {
|
---|
620 | c = false;
|
---|
621 | }
|
---|
622 | else
|
---|
623 | {
|
---|
624 | tmp = x[j];
|
---|
625 | x[j] = x[j+g];
|
---|
626 | x[j+g] = tmp;
|
---|
627 | tmp = y[j];
|
---|
628 | y[j] = y[j+g];
|
---|
629 | y[j+g] = tmp;
|
---|
630 | }
|
---|
631 | j = j-1;
|
---|
632 | }
|
---|
633 | while( j>=0 & c );
|
---|
634 | i = i+1;
|
---|
635 | }
|
---|
636 | while( i<=n );
|
---|
637 | g = g/2;
|
---|
638 | }
|
---|
639 | while( g>0 );
|
---|
640 | ctbl = new double[4+1, n+1];
|
---|
641 | n = n+1;
|
---|
642 | if( diffn==1 )
|
---|
643 | {
|
---|
644 | b1 = 1;
|
---|
645 | b2 = 6/(x[1]-x[0])*((y[1]-y[0])/(x[1]-x[0])-boundl);
|
---|
646 | b3 = 1;
|
---|
647 | b4 = 6/(x[n-1]-x[n-2])*(boundr-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
|
---|
648 | }
|
---|
649 | else
|
---|
650 | {
|
---|
651 | b1 = 0;
|
---|
652 | b2 = 2*boundl;
|
---|
653 | b3 = 0;
|
---|
654 | b4 = 2*boundr;
|
---|
655 | }
|
---|
656 | nxm1 = n-1;
|
---|
657 | if( n>=2 )
|
---|
658 | {
|
---|
659 | if( n>2 )
|
---|
660 | {
|
---|
661 | dxj = x[1]-x[0];
|
---|
662 | dyj = y[1]-y[0];
|
---|
663 | j = 2;
|
---|
664 | while( j<=nxm1 )
|
---|
665 | {
|
---|
666 | dxjp1 = x[j]-x[j-1];
|
---|
667 | dyjp1 = y[j]-y[j-1];
|
---|
668 | dxp = dxj+dxjp1;
|
---|
669 | ctbl[1,j-1] = dxjp1/dxp;
|
---|
670 | ctbl[2,j-1] = 1-ctbl[1,j-1];
|
---|
671 | ctbl[3,j-1] = 6*(dyjp1/dxjp1-dyj/dxj)/dxp;
|
---|
672 | dxj = dxjp1;
|
---|
673 | dyj = dyjp1;
|
---|
674 | j = j+1;
|
---|
675 | }
|
---|
676 | }
|
---|
677 | ctbl[1,0] = -(b1/2);
|
---|
678 | ctbl[2,0] = b2/2;
|
---|
679 | if( n!=2 )
|
---|
680 | {
|
---|
681 | j = 2;
|
---|
682 | while( j<=nxm1 )
|
---|
683 | {
|
---|
684 | pj = ctbl[2,j-1]*ctbl[1,j-2]+2;
|
---|
685 | ctbl[1,j-1] = -(ctbl[1,j-1]/pj);
|
---|
686 | ctbl[2,j-1] = (ctbl[3,j-1]-ctbl[2,j-1]*ctbl[2,j-2])/pj;
|
---|
687 | j = j+1;
|
---|
688 | }
|
---|
689 | }
|
---|
690 | yppb = (b4-b3*ctbl[2,nxm1-1])/(b3*ctbl[1,nxm1-1]+2);
|
---|
691 | i = 1;
|
---|
692 | while( i<=nxm1 )
|
---|
693 | {
|
---|
694 | j = n-i;
|
---|
695 | yppa = ctbl[1,j-1]*yppb+ctbl[2,j-1];
|
---|
696 | dx = x[j]-x[j-1];
|
---|
697 | ctbl[3,j-1] = (yppb-yppa)/dx/6;
|
---|
698 | ctbl[2,j-1] = yppa/2;
|
---|
699 | ctbl[1,j-1] = (y[j]-y[j-1])/dx-(ctbl[2,j-1]+ctbl[3,j-1]*dx)*dx;
|
---|
700 | yppb = yppa;
|
---|
701 | i = i+1;
|
---|
702 | }
|
---|
703 | for(i=1; i<=n; i++)
|
---|
704 | {
|
---|
705 | ctbl[0,i-1] = y[i-1];
|
---|
706 | ctbl[4,i-1] = x[i-1];
|
---|
707 | }
|
---|
708 | }
|
---|
709 | }
|
---|
710 |
|
---|
711 |
|
---|
712 | public static double spline3interpolate(int n,
|
---|
713 | ref double[,] c,
|
---|
714 | double x)
|
---|
715 | {
|
---|
716 | double result = 0;
|
---|
717 | int i = 0;
|
---|
718 | int l = 0;
|
---|
719 | int half = 0;
|
---|
720 | int first = 0;
|
---|
721 | int middle = 0;
|
---|
722 |
|
---|
723 | n = n-1;
|
---|
724 | l = n;
|
---|
725 | first = 0;
|
---|
726 | while( l>0 )
|
---|
727 | {
|
---|
728 | half = l/2;
|
---|
729 | middle = first+half;
|
---|
730 | if( (double)(c[4,middle])<(double)(x) )
|
---|
731 | {
|
---|
732 | first = middle+1;
|
---|
733 | l = l-half-1;
|
---|
734 | }
|
---|
735 | else
|
---|
736 | {
|
---|
737 | l = half;
|
---|
738 | }
|
---|
739 | }
|
---|
740 | i = first-1;
|
---|
741 | if( i<0 )
|
---|
742 | {
|
---|
743 | i = 0;
|
---|
744 | }
|
---|
745 | result = c[0,i]+(x-c[4,i])*(c[1,i]+(x-c[4,i])*(c[2,i]+c[3,i]*(x-c[4,i])));
|
---|
746 | return result;
|
---|
747 | }
|
---|
748 |
|
---|
749 |
|
---|
750 | private static void heapsortpoints(ref double[] x,
|
---|
751 | ref double[] y,
|
---|
752 | int n)
|
---|
753 | {
|
---|
754 | int i = 0;
|
---|
755 | int j = 0;
|
---|
756 | int k = 0;
|
---|
757 | int t = 0;
|
---|
758 | double tmp = 0;
|
---|
759 | bool isascending = new bool();
|
---|
760 | bool isdescending = new bool();
|
---|
761 |
|
---|
762 |
|
---|
763 | //
|
---|
764 | // Test for already sorted set
|
---|
765 | //
|
---|
766 | isascending = true;
|
---|
767 | isdescending = true;
|
---|
768 | for(i=1; i<=n-1; i++)
|
---|
769 | {
|
---|
770 | isascending = isascending & (double)(x[i])>(double)(x[i-1]);
|
---|
771 | isdescending = isdescending & (double)(x[i])<(double)(x[i-1]);
|
---|
772 | }
|
---|
773 | if( isascending )
|
---|
774 | {
|
---|
775 | return;
|
---|
776 | }
|
---|
777 | if( isdescending )
|
---|
778 | {
|
---|
779 | for(i=0; i<=n-1; i++)
|
---|
780 | {
|
---|
781 | j = n-1-i;
|
---|
782 | if( j<=i )
|
---|
783 | {
|
---|
784 | break;
|
---|
785 | }
|
---|
786 | tmp = x[i];
|
---|
787 | x[i] = x[j];
|
---|
788 | x[j] = tmp;
|
---|
789 | tmp = y[i];
|
---|
790 | y[i] = y[j];
|
---|
791 | y[j] = tmp;
|
---|
792 | }
|
---|
793 | return;
|
---|
794 | }
|
---|
795 |
|
---|
796 | //
|
---|
797 | // Special case: N=1
|
---|
798 | //
|
---|
799 | if( n==1 )
|
---|
800 | {
|
---|
801 | return;
|
---|
802 | }
|
---|
803 |
|
---|
804 | //
|
---|
805 | // General case
|
---|
806 | //
|
---|
807 | i = 2;
|
---|
808 | do
|
---|
809 | {
|
---|
810 | t = i;
|
---|
811 | while( t!=1 )
|
---|
812 | {
|
---|
813 | k = t/2;
|
---|
814 | if( (double)(x[k-1])>=(double)(x[t-1]) )
|
---|
815 | {
|
---|
816 | t = 1;
|
---|
817 | }
|
---|
818 | else
|
---|
819 | {
|
---|
820 | tmp = x[k-1];
|
---|
821 | x[k-1] = x[t-1];
|
---|
822 | x[t-1] = tmp;
|
---|
823 | tmp = y[k-1];
|
---|
824 | y[k-1] = y[t-1];
|
---|
825 | y[t-1] = tmp;
|
---|
826 | t = k;
|
---|
827 | }
|
---|
828 | }
|
---|
829 | i = i+1;
|
---|
830 | }
|
---|
831 | while( i<=n );
|
---|
832 | i = n-1;
|
---|
833 | do
|
---|
834 | {
|
---|
835 | tmp = x[i];
|
---|
836 | x[i] = x[0];
|
---|
837 | x[0] = tmp;
|
---|
838 | tmp = y[i];
|
---|
839 | y[i] = y[0];
|
---|
840 | y[0] = tmp;
|
---|
841 | t = 1;
|
---|
842 | while( t!=0 )
|
---|
843 | {
|
---|
844 | k = 2*t;
|
---|
845 | if( k>i )
|
---|
846 | {
|
---|
847 | t = 0;
|
---|
848 | }
|
---|
849 | else
|
---|
850 | {
|
---|
851 | if( k<i )
|
---|
852 | {
|
---|
853 | if( (double)(x[k])>(double)(x[k-1]) )
|
---|
854 | {
|
---|
855 | k = k+1;
|
---|
856 | }
|
---|
857 | }
|
---|
858 | if( (double)(x[t-1])>=(double)(x[k-1]) )
|
---|
859 | {
|
---|
860 | t = 0;
|
---|
861 | }
|
---|
862 | else
|
---|
863 | {
|
---|
864 | tmp = x[k-1];
|
---|
865 | x[k-1] = x[t-1];
|
---|
866 | x[t-1] = tmp;
|
---|
867 | tmp = y[k-1];
|
---|
868 | y[k-1] = y[t-1];
|
---|
869 | y[t-1] = tmp;
|
---|
870 | t = k;
|
---|
871 | }
|
---|
872 | }
|
---|
873 | }
|
---|
874 | i = i-1;
|
---|
875 | }
|
---|
876 | while( i>=1 );
|
---|
877 | }
|
---|
878 |
|
---|
879 |
|
---|
880 | private static void heapsortdpoints(ref double[] x,
|
---|
881 | ref double[] y,
|
---|
882 | ref double[] d,
|
---|
883 | int n)
|
---|
884 | {
|
---|
885 | int i = 0;
|
---|
886 | int j = 0;
|
---|
887 | int k = 0;
|
---|
888 | int t = 0;
|
---|
889 | double tmp = 0;
|
---|
890 | bool isascending = new bool();
|
---|
891 | bool isdescending = new bool();
|
---|
892 |
|
---|
893 |
|
---|
894 | //
|
---|
895 | // Test for already sorted set
|
---|
896 | //
|
---|
897 | isascending = true;
|
---|
898 | isdescending = true;
|
---|
899 | for(i=1; i<=n-1; i++)
|
---|
900 | {
|
---|
901 | isascending = isascending & (double)(x[i])>(double)(x[i-1]);
|
---|
902 | isdescending = isdescending & (double)(x[i])<(double)(x[i-1]);
|
---|
903 | }
|
---|
904 | if( isascending )
|
---|
905 | {
|
---|
906 | return;
|
---|
907 | }
|
---|
908 | if( isdescending )
|
---|
909 | {
|
---|
910 | for(i=0; i<=n-1; i++)
|
---|
911 | {
|
---|
912 | j = n-1-i;
|
---|
913 | if( j<=i )
|
---|
914 | {
|
---|
915 | break;
|
---|
916 | }
|
---|
917 | tmp = x[i];
|
---|
918 | x[i] = x[j];
|
---|
919 | x[j] = tmp;
|
---|
920 | tmp = y[i];
|
---|
921 | y[i] = y[j];
|
---|
922 | y[j] = tmp;
|
---|
923 | tmp = d[i];
|
---|
924 | d[i] = d[j];
|
---|
925 | d[j] = tmp;
|
---|
926 | }
|
---|
927 | return;
|
---|
928 | }
|
---|
929 |
|
---|
930 | //
|
---|
931 | // Special case: N=1
|
---|
932 | //
|
---|
933 | if( n==1 )
|
---|
934 | {
|
---|
935 | return;
|
---|
936 | }
|
---|
937 |
|
---|
938 | //
|
---|
939 | // General case
|
---|
940 | //
|
---|
941 | i = 2;
|
---|
942 | do
|
---|
943 | {
|
---|
944 | t = i;
|
---|
945 | while( t!=1 )
|
---|
946 | {
|
---|
947 | k = t/2;
|
---|
948 | if( (double)(x[k-1])>=(double)(x[t-1]) )
|
---|
949 | {
|
---|
950 | t = 1;
|
---|
951 | }
|
---|
952 | else
|
---|
953 | {
|
---|
954 | tmp = x[k-1];
|
---|
955 | x[k-1] = x[t-1];
|
---|
956 | x[t-1] = tmp;
|
---|
957 | tmp = y[k-1];
|
---|
958 | y[k-1] = y[t-1];
|
---|
959 | y[t-1] = tmp;
|
---|
960 | tmp = d[k-1];
|
---|
961 | d[k-1] = d[t-1];
|
---|
962 | d[t-1] = tmp;
|
---|
963 | t = k;
|
---|
964 | }
|
---|
965 | }
|
---|
966 | i = i+1;
|
---|
967 | }
|
---|
968 | while( i<=n );
|
---|
969 | i = n-1;
|
---|
970 | do
|
---|
971 | {
|
---|
972 | tmp = x[i];
|
---|
973 | x[i] = x[0];
|
---|
974 | x[0] = tmp;
|
---|
975 | tmp = y[i];
|
---|
976 | y[i] = y[0];
|
---|
977 | y[0] = tmp;
|
---|
978 | tmp = d[i];
|
---|
979 | d[i] = d[0];
|
---|
980 | d[0] = tmp;
|
---|
981 | t = 1;
|
---|
982 | while( t!=0 )
|
---|
983 | {
|
---|
984 | k = 2*t;
|
---|
985 | if( k>i )
|
---|
986 | {
|
---|
987 | t = 0;
|
---|
988 | }
|
---|
989 | else
|
---|
990 | {
|
---|
991 | if( k<i )
|
---|
992 | {
|
---|
993 | if( (double)(x[k])>(double)(x[k-1]) )
|
---|
994 | {
|
---|
995 | k = k+1;
|
---|
996 | }
|
---|
997 | }
|
---|
998 | if( (double)(x[t-1])>=(double)(x[k-1]) )
|
---|
999 | {
|
---|
1000 | t = 0;
|
---|
1001 | }
|
---|
1002 | else
|
---|
1003 | {
|
---|
1004 | tmp = x[k-1];
|
---|
1005 | x[k-1] = x[t-1];
|
---|
1006 | x[t-1] = tmp;
|
---|
1007 | tmp = y[k-1];
|
---|
1008 | y[k-1] = y[t-1];
|
---|
1009 | y[t-1] = tmp;
|
---|
1010 | tmp = d[k-1];
|
---|
1011 | d[k-1] = d[t-1];
|
---|
1012 | d[t-1] = tmp;
|
---|
1013 | t = k;
|
---|
1014 | }
|
---|
1015 | }
|
---|
1016 | }
|
---|
1017 | i = i-1;
|
---|
1018 | }
|
---|
1019 | while( i>=1 );
|
---|
1020 | }
|
---|
1021 |
|
---|
1022 |
|
---|
1023 | private static void solvetridiagonal(double[] a,
|
---|
1024 | double[] b,
|
---|
1025 | double[] c,
|
---|
1026 | double[] d,
|
---|
1027 | int n,
|
---|
1028 | ref double[] x)
|
---|
1029 | {
|
---|
1030 | int k = 0;
|
---|
1031 | double t = 0;
|
---|
1032 |
|
---|
1033 | a = (double[])a.Clone();
|
---|
1034 | b = (double[])b.Clone();
|
---|
1035 | c = (double[])c.Clone();
|
---|
1036 | d = (double[])d.Clone();
|
---|
1037 |
|
---|
1038 | x = new double[n-1+1];
|
---|
1039 | a[0] = 0;
|
---|
1040 | c[n-1] = 0;
|
---|
1041 | for(k=1; k<=n-1; k++)
|
---|
1042 | {
|
---|
1043 | t = a[k]/b[k-1];
|
---|
1044 | b[k] = b[k]-t*c[k-1];
|
---|
1045 | d[k] = d[k]-t*d[k-1];
|
---|
1046 | }
|
---|
1047 | x[n-1] = d[n-1]/b[n-1];
|
---|
1048 | for(k=n-2; k>=0; k--)
|
---|
1049 | {
|
---|
1050 | x[k] = (d[k]-c[k]*x[k+1])/b[k];
|
---|
1051 | }
|
---|
1052 | }
|
---|
1053 |
|
---|
1054 |
|
---|
1055 | private static double diffthreepoint(double t,
|
---|
1056 | double x0,
|
---|
1057 | double f0,
|
---|
1058 | double x1,
|
---|
1059 | double f1,
|
---|
1060 | double x2,
|
---|
1061 | double f2)
|
---|
1062 | {
|
---|
1063 | double result = 0;
|
---|
1064 | double a = 0;
|
---|
1065 | double b = 0;
|
---|
1066 |
|
---|
1067 | t = t-x0;
|
---|
1068 | x1 = x1-x0;
|
---|
1069 | x2 = x2-x0;
|
---|
1070 | a = (f2-f0-x2/x1*(f1-f0))/(AP.Math.Sqr(x2)-x1*x2);
|
---|
1071 | b = (f1-f0-a*AP.Math.Sqr(x1))/x1;
|
---|
1072 | result = 2*a*t+b;
|
---|
1073 | return result;
|
---|
1074 | }
|
---|
1075 | }
|
---|
1076 | }
|
---|