Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/spline3.cs @ 2465

Last change on this file since 2465 was 2430, checked in by gkronber, 15 years ago

Moved ALGLIB code into a separate plugin. #783

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