Spline.cpp

#include "Spline.h"
#include <cstring>
#include <algorithm>
using namespace std;

namespace SplineSpace{
	Spline::Spline(vector<double>& x0, vector<double>& y0, const int num,
				BoundaryCondition bc, const double leftBoundary, const double rightBoundary)
				:GivenX(x0), GivenY(y0), GivenNum(num), Bc(bc), LeftB(leftBoundary), RightB(rightBoundary)
	{
		if(x0.empty() || y0.empty() || (num<3)){
			throw SplineFailure("Points too few to create.");
		}
		PartialDerivative.resize(GivenNum, 0.0);	
		MaxX = *max_element(GivenX.begin(), GivenX.begin() + num);
		MinX = *min_element(GivenX.begin(), GivenX.begin() + num);
		if(Bc==GivenFirstOrder)	{
			PartialDerivative1();
		}
		else if(Bc == GivenSecondOrder){		
			PartialDerivative2();
		}
		else{
			PartialDerivative.clear();
			throw SplineFailure("Wrong parameters for boundary conditions.");
		}
	}

	//Type I Boundary Partial Derivative
	void Spline::PartialDerivative1(void)
	{
		// Using the tridiagonal matrix algorithm (Thomas algorithm) to solve for the second derivatives
		double *a=new double[GivenNum];                //  a: The bottom diagonal of the tridiagonal matrix
		double *b=new double[GivenNum];                //  b: The main diagonal of the tridiagonal matrix
		double *c=new double[GivenNum];                //  c: The top diagonal of the tridiagonal matrix
		double *d=new double[GivenNum];

		double* f = new double[GivenNum];

		double* bt = new double[GivenNum];
		double* gm = new double[GivenNum];

		double* h = new double[GivenNum];

		for(int i=0;i<GivenNum;i++)  b[i]=2;          //  middle nums are 2
		for(int i=0;i<GivenNum-1;i++)  h[i]=GivenX[i+1]-GivenX[i];                   // steps
		for(int i=1;i<GivenNum-1;i++)  a[i]=h[i-1]/(h[i-1]+h[i]);            
		a[GivenNum-1]=1;

		c[0] = 1;
		for (int i = 1; i < GivenNum - 1; i++)  c[i] = h[i] / (h[i - 1] + h[i]);

		for (int i = 0; i < GivenNum - 1; i++)
			f[i] = (GivenY[i + 1] - GivenY[i]) / (GivenX[i + 1] - GivenX[i]);

		d[0] = 6 * (f[0] - LeftB) / h[0];
		d[GivenNum - 1] = 6 * (RightB - f[GivenNum - 2]) / h[GivenNum - 2];

		for (int i = 1; i < GivenNum - 1; i++)  d[i] = 6 * (f[i] - f[i - 1]) / (h[i - 1] + h[i]);

		bt[0]=c[0]/b[0];                                            
		for(int i=1;i<GivenNum-1;i++)  bt[i]=c[i]/(b[i]-a[i]*bt[i-1]);


		gm[0] = d[0] / b[0];
		for (int i = 1; i <= GivenNum - 1; i++)  gm[i] = (d[i] - a[i] * gm[i - 1]) / (b[i] - a[i] * bt[i - 1]);

		PartialDerivative[GivenNum - 1] = gm[GivenNum - 1];
		for (int i = GivenNum - 2; i >= 0; i--)  PartialDerivative[i] = gm[i] - bt[i] * PartialDerivative[i + 1];

		delete[] a;
		delete[] b;
		delete[] c;
		delete[] d;
		delete[] gm;
		delete[] bt;
		delete[] f;
		delete[] h;
	}

	//Type II Boundary Partial Derivative
	void Spline::PartialDerivative2(void)
	{

		double *a=new double[GivenNum];             
		double *b=new double[GivenNum];        
		double *c=new double[GivenNum];               
		double *d=new double[GivenNum];

		double* f = new double[GivenNum];

		double* bt = new double[GivenNum];
		double* gm = new double[GivenNum];

		double* h = new double[GivenNum];

		for (int i = 0; i < GivenNum; i++)  b[i] = 2;
		for (int i = 0; i < GivenNum - 1; i++)  h[i] = GivenX[i + 1] - GivenX[i];
		for (int i = 1; i < GivenNum - 1; i++)  a[i] = h[i - 1] / (h[i - 1] + h[i]);
		a[GivenNum - 1] = 1;

		c[0] = 1;
		for (int i = 1; i < GivenNum - 1; i++)  c[i] = h[i] / (h[i - 1] + h[i]);

		for (int i = 0; i < GivenNum - 1; i++)
			f[i] = (GivenY[i + 1] - GivenY[i]) / (GivenX[i + 1] - GivenX[i]);

		for (int i = 1; i < GivenNum - 1; i++)  d[i] = 6 * (f[i] - f[i - 1]) / (h[i - 1] + h[i]);

		d[1] = d[1] - a[1] * LeftB;
		d[GivenNum - 2] = d[GivenNum - 2] - c[GivenNum - 2] * RightB;

		bt[1] = c[1] / b[1];
		for (int i = 2; i < GivenNum - 2; i++)  bt[i] = c[i] / (b[i] - a[i] * bt[i - 1]);

		gm[1] = d[1] / b[1];
		for (int i = 2; i <= GivenNum - 2; i++)  gm[i] = (d[i] - a[i] * gm[i - 1]) / (b[i] - a[i] * bt[i - 1]);

		PartialDerivative[GivenNum - 2] = gm[GivenNum - 2];//
		for (int i = GivenNum - 3; i >= 1; i--)  PartialDerivative[i] = gm[i] - bt[i] * PartialDerivative[i + 1];

		PartialDerivative[0] = LeftB;
		PartialDerivative[GivenNum - 1] = RightB;

		delete[] a;
		delete[] b;
		delete[] c;
		delete[] d;
		delete[] gm;
		delete[] bt;
		delete[] f;
		delete[] h;
	}

	//Single Point interpretation 
	bool Spline::SinglePointInterp(const double x, double& y)throw(SplineFailure){
		if( (x<MinX) || (x>MaxX) )
			throw SplineFailure("No support for outer interoperation.");
		int klo,khi,k;
		klo=0; khi=GivenNum-1;
		double hh,bb,aa;

		// Using the binary search method 
		while(khi-klo>1)         
		{
			k = (khi + klo) >> 1;
			if (GivenX[k] > x)  khi = k;
			else klo = k;
		}
		hh = GivenX[khi] - GivenX[klo];

		aa = (GivenX[khi] - x) / hh;
		bb = (x - GivenX[klo]) / hh;

		y = aa * GivenY[klo] + bb * GivenY[khi] + ((aa * aa * aa - aa) * PartialDerivative[klo] + (bb * bb * bb - bb) * PartialDerivative[khi]) * hh * hh / 6.0;
		return true;
	}

	//Multiple Points Interpretation.
	bool Spline::MultiPointInterp(vector<double>& x, const int num, vector<double>& y)throw(SplineFailure){
		for(int i = 0;i < num;i++)
		{
			SinglePointInterp(x[i], y[i]);
		}
		return true;
	}

	// Auto Interpretation
	bool Spline::AutoInterp(const int num, vector<double>& x, vector<double>& y)throw(SplineFailure){
		if(num < 2)
			throw SplineFailure("At least 2 points.");

		double perStep = (MaxX-MinX)/(num-1);
		for(int i = 0;i < num;i++){
			x[i] = MinX+i*perStep;
			SinglePointInterp(x[i],y[i]);
		}
		return true;
	}

	Spline::~Spline(){
		PartialDerivative.clear();	
	}

	SplineFailure::SplineFailure(const char* msg):Message(msg){};
	const char* SplineFailure::GetMessage(){return Message;}

}