CubicSpline.Interpolate()

        public override void Interpolate()
        {
            base.Interpolate();
            if (Points.Count == 0)
                return;
            if (SegmentCount < 1)
                a = b = c = d = null;
            else {
                if (a == null || a.Count() != SegmentCount) {
                    my = new float[SegmentCount - 1];
                    mc = new float[SegmentCount - 1];
                    mr = new float[Math.Max(0, SegmentCount - 2)];
                    a = new float[SegmentCount];
                    b = new float[SegmentCount];
                    c = new float[SegmentCount];
                    d = new float[SegmentCount];
                }

                float leftDeriv;
                bool leftNatural = !GetProperty("LeftDerivative", out leftDeriv);
                float rightDeriv;
                bool rightNatural = !GetProperty("RightDerivative", out rightDeriv);

                for (int i = 0; i < SegmentCount; ++i)
                    d[i] = Points[i].Y;

                for (int i = 1; i < SegmentCount; ++i) {
                    mc[i - 1] = 2f * (Points[i + 1].X - Points[i - 1].X);

                    my[i - 1] = 3f * ((Points[i + 1].Y - Points[i].Y) / (Points[i + 1].X - Points[i].X) -
                                  (Points[i].Y - Points[i - 1].Y) / (Points[i].X - Points[i - 1].X));

                    if (i < SegmentCount - 1)
                        mr[i - 1] = Points[i + 1].X - Points[i].X;
                }

                for (int i = 1; i < SegmentCount - 1; ++i) {
                    float k = mr[i - 1] / mc[i - 1];
                    mc[i] = mc[i] - mr[i - 1] * k;
                    my[i] = my[i] - my[i - 1] * k;
                }

                if (SegmentCount > 1) {
                    b[SegmentCount - 1] = my[SegmentCount - 2] / mc[SegmentCount - 2];

                    for (int i = SegmentCount - 2; i >= 1; --i)
                        b[i] = (my[i - 1] - b[i + 1] * mr[i - 1]) / mc[i - 1];
                }

                b[0] = 0f;

                for (int i = 0; i < SegmentCount - 1; ++i)
                    a[i] = 1f / (3f * (Points[i + 1].X - Points[i].X)) * (b[i + 1] - b[i]);
                a[SegmentCount - 1] = -1f / (3f * (Points[SegmentCount].X - Points[SegmentCount - 1].X)) * b[SegmentCount - 1];

                for (int i = 0; i < SegmentCount - 1; ++i)
                {
                    float h = Points[i + 1].X - Points[i].X;
                    c[i] = (Points[i + 1].Y - Points[i].Y) / h - h / 3f * (b[i + 1] + 2f * b[i]);
                }
                {
                    float h = Points[SegmentCount].X - Points[SegmentCount - 1].X;
                    c[SegmentCount - 1] = (Points[SegmentCount].Y - Points[SegmentCount - 1].Y) / h - h / 3f * (2f * b[SegmentCount - 1]);
                }
            }
        }

        public override float GetValue(float x)
        {
            if (SegmentCount < 1)
                return Points[0].Y;

            int s = FindSegment(x);
            if (x <= Points[0].X) {
                float xr = x - Points[0].X;
                return c[0] * xr + d[0];
            }
            else if (x >= Points[PointsCount - 1].X) {
                s = SegmentCount - 1;
                float xr = x - Points[SegmentCount].X;
                float h = Points[SegmentCount].X - Points[s].X;
                return ((3f * a[s] * h + 2f * b[s]) * h + c[s]) * xr + Points[SegmentCount].Y;
            }
            else {
                float xr = x - Points[s].X;
                return ((a[s] * xr + b[s]) * xr + c[s]) * xr + d[s];
            }
        }