Free cookie consent management tool by TermsFeed Policy Generator

Ignore:
Timestamp:
10/15/09 13:22:07 (15 years ago)
Author:
gkronber
Message:

Moved ALGLIB code into a separate plugin. #783

Location:
trunk/sources/ALGLIB
Files:
1 added
1 copied

Legend:

Unmodified
Added
Removed

  • trunk/sources/ALGLIB/spline3.cs

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