Cubic Spline CSplineA

Spline is a sufficiently smooth polynomial function that is piecewise-defined, and possesses a high degree of smoothness at the places where the polynomial pieces connect (wikipedia). There is some of method to get Spline graphics. The most commonly used splines are cubic spline, B-Spline, Spline interpolation, etc.

I get one of cubic spline method with name Cubic Spline CSplineA method. I get this method from excel VBA file. This is a free software, so you can get the complete Excel computation file for this Cubic Spline CSplineA method. I have success to convert this excel VBA file to C#. So, you can merge this code if you want to use this code in your program.

This is the complete code for Cubic Spline CSplineA code in C# :

using System;
using System.Collections.Generic;
using System.Text;

namespace consoleSplineCubic
{
    class CSplineA
    {
        public double[,] CSplineAFunc(double[] XVal, double[] YVal, long Numrows, 
            double[] Xint, long NumOut, long Out = 1)
        {
            double[,] ResA;
            double[,] FactA;
            double X, Chord, ArcAng, AvR;
            long i, j;

            AvR = 1.0;
            FactA = new double[Numrows,8+1];
            ResA = new double[NumOut, 6+1];

            for (i = 0; i<Numrows - 1; i++)
            {
                FactA[i, 1] = XVal[i + 1] - XVal[i];
            }

            for (i = 1; i<Numrows - 1; i++) 
            {
                FactA[i, 2] = (3 / FactA[i, 1]) * (YVal[i + 1] - YVal[i]) - (3 / FactA[i - 1, 1]) * (YVal[i] - YVal[i - 1]);
            }

            FactA[0, 3] = 1;
            FactA[0, 4] = 0;
            FactA[0, 5] = 0;

            for (i = 1; i<Numrows - 1; i++)
            {
                FactA[i, 3] = 2 * (XVal[i + 1] - XVal[i - 1]) - FactA[i - 1, 1] * FactA[i - 1, 4];
                FactA[i, 4] = FactA[i, 1] / FactA[i, 3];
                FactA[i, 5] = (FactA[i, 2] - FactA[i - 1, 1] * FactA[i - 1, 5]) / FactA[i, 3];
            }

            FactA[i, 3] = 1;
            FactA[i, 5] = 0;
            FactA[i, 7] = 0;

            for (i = Numrows - 2; i>= 0; i--) 
            {
                FactA[i, 7] = FactA[i, 5] - FactA[i, 4] * FactA[i + 1, 7];
                FactA[i, 6] = (YVal[i + 1] - YVal[i]) / FactA[i, 1] - FactA[i, 1] * (FactA[i + 1, 7] + 2 * FactA[i, 7]) / 3;
                FactA[i, 8] = (FactA[i + 1, 7] - FactA[i, 7]) / (3 * FactA[i, 1]);
            }

            i = Numrows-1;
            X = XVal[i] - XVal[i - 1];
            //' FactA(i, 6) = 3 * X ^ 2 * FactA(i - 1, 8) + 2 * X * FactA(i - 1, 7) + FactA(i - 1, 6)

            if (Out == 2)
            {
                return FactA;
            }

            j = 0;
            for (i = 0; i<NumOut; i++)
            {
                X = Xint[i];
                if (X >= XVal[0] && X <= XVal[Numrows-1]) 
                {
                    if (X > XVal[j])  j = 0;
                    while (j < Numrows-1) 
                    {
                        if (X < XVal[j]) break;
                        j = j + 1;
                    }

                    j = j - 1;
                    X = X - XVal[j];
                    ResA[i, 1] = Math.Pow(X, 3) * FactA[j, 8] + Math.Pow(X, 2) * FactA[j, 7] + X * FactA[j, 6] + YVal[j];
                    ResA[i, 2] = 3 * Math.Pow(X, 2) * FactA[j, 8] + 2 * X * FactA[j, 7] + FactA[j, 6];
                    ResA[i, 3] = 6 * X * FactA[j, 8] + 2 * FactA[j, 7];
                    if (Math.Abs(ResA[i, 3]) > 0.0000000001)
                        ResA[i, 4] = Math.Pow((1 + Math.Pow(ResA[i, 2], 2)), (3 / 2)) / ResA[i, 3];
                    else
                        ResA[i, 4] = 0.0; //"-"

                    if (i > 0) 
                    {
                        Chord = Math.Pow( Math.Pow(Xint[i] - Xint[i - 1], 2) + Math.Pow(ResA[i, 1] - ResA[i - 1, 1], 2) , 0.5);
                        ArcAng = Math.Atan(ResA[i, 2]) - Math.Atan(ResA[i - 1, 2]);
                        if (ArcAng != 0) AvR = Chord / 2 / Math.Sin(ArcAng / 2);
                        ResA[i, 5] = ArcAng * AvR;
                        ResA[i, 6] = Chord;
                    }
                    else 
                    {
                        ResA[i, 5] = 0.0;//"-"
                        ResA[i, 6] = 0.0; //"-"
                    }
                }
                else 
                {
                    if (X < XVal[0]) j = 0; 
                    else j = Numrows-1;

                    X = X - XVal[j];
                    ResA[i, 1] = X * FactA[j, 6] + YVal[j];
                    ResA[i, 2] = FactA[j, 6];
                    ResA[i, 3] = 0;
                }
            }

            return ResA;

        }
    }
}

If you want to get the complete excel file for this program, you can get the complete excel file at here :

http://interactiveds.com.au/software/

You can get the demo about output from this Cubic Spline CSplineA at here :

http://newtonexcelbach.wordpress.com/2009/07/02/cubic-splines/

If you have a question about this Cubic Spline CSplineA code above, please send question to me.

2 Comments

Leave a Reply to Przemek Cancel reply

Your email address will not be published. Required fields are marked *