#pragma once

#include "./numerics/interpolation1d/interpolation1d_base.h"

//#include "./numerics/abs.h"
#include "./utils/vector.h"


namespace numerics{
	
	template <typename T>
	struct interp_cubic_spline : Base_interp<T> {
		using Base   = Base_interp<T>;
	  	// bring base data members into scope (or use this->xx / this->yy below)
	  	using Base::xx;
	  	using Base::yy;
	  	//using Base::mm;

	  	utils::Vector<T> y2;


		interp_cubic_spline(utils::Vector<T> &xv, utils::Vector<T> &yv, T yp1=T{1.e99}, T ypn=T{1.e99})
		: Base_interp<T>(xv, &yv[0], 2), y2(xv.size(),T{0}) {
			sety2(&xv[0], &yv[0], yp1, ypn);
		}

		interp_cubic_spline(utils::Vector<T> &xv, const T *yv, T yp1=T{1.e99}, T ypn=T{1.e99})
		: Base_interp<T>(xv, yv, 2), y2(xv.size(),T{0}) {
			sety2(&xv[0], yv, yp1, ypn);
		}


		void sety2(const T *xv, const T *yv, T yp1, T ypn){

			T p, qn, sig, un;
			uint64_t n = y2.size();
			utils::Vector<T> u(n-1, T{0});

			if (yp1 > static_cast<T>(0.99e99)){						// The lower boundary condition is set either to be “natural”
				y2[0] = u[0] = T{0};
			}else{													// or else to have a specified first derivative.
				y2[0] = T{-0.5};
				u[0] = (3.0/(xv[1]-xv[0]))*(((yv[1]-yv[0])/(xv[1]-xv[0]))-yp1);
			}
			for (uint64_t i = 1; i < n-1; ++i){						// This is the decomposition loop of the tridiagonal algorithm
				sig = (xv[i]-xv[i-1])/(xv[i+1]-xv[i-1]);
				p = sig*y2[i-1]+T{2};
				y2[i] = (sig - T{1})/p;								// y2 and u are used for temporary storage of the decomposed factors.
				u[i]=((yv[i+1]-yv[i])/(xv[i+1]-xv[i])) - ((yv[i]-yv[i-1])/(xv[i]-xv[i-1]));
				u[i]=((T{6}*u[i]/(xv[i+1]-xv[i-1])) - sig*u[i-1])/p;
			}
			if (ypn > static_cast<T>(0.99e99)){						// The upper boundary condition is set either to be “natural”
				qn = un = T{0};
			}else{													// or else to have a specified first derivative.
				qn = T{0.5};
				un = (T{3}/(xv[n-1]-xv[n-2]))*(ypn-((yv[n-1]-yv[n-2])/(xv[n-1]-xv[n-2])));
			}
			y2[n-1] = (un-(qn*u[n-2]))/((qn*y2[n-2])+T{1});
			for (int64_t k = n-2; k >= 0; --k){
				y2[k] = y2[k] * y2[k+1]+u[k];
			}
		}


		T rawinterp(int64_t jl, T x) override{

			int64_t klo=jl, khi=jl+1;
			T y, h, b, a;

			h = xx[khi] - xx[klo];
			if (h == T{0}){											// The xa’s must be distinct.
				throw std::invalid_argument("interp_cubic_spline: Bad input to routine splint");
			}

			a = (xx[khi] - x)/h;									// Cubic spline polynomial is now evaluated.
			b = (x - xx[klo])/h;
			y = a*yy[klo] + b*yy[khi] + ( ((a*a*a) - a)*y2[klo] + ((b*b*b) - b)*y2[khi] ) * (h*h) / T{6};


			return y;
		}

	};

} // namespace numerics