Free cookie consent management tool by TermsFeed Policy Generator

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

Last change on this file since 2643 was 2563, checked in by gkronber, 15 years ago

Updated ALGLIB to latest version. #751 (Plugin for for data-modeling with ANN (integrated into CEDMA))
File size: 31.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        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}
Note: See TracBrowser for help on using the repository browser.