/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// Spline interpolation                                                    // 
//                                                                         //
// Burkhard Militzer                                  Livermore 4-27-01    //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////

#ifndef _SPLINE_
#define _SPLINE_

#include <algorithm>
#include "Array.h"
#include "Vector.h"

class Spline {
 public:
  int n;
  Array1 <double> x;
  Array1 <double> y;
  Array1 <double> y2;
  //  static Array1 <double> u; // no longer static for openMP

 public:
  Spline():n(0){};
  Spline(const int m):n(0){ // just pre-allocate memory, do not fill in any elements yet
    Allocate(m); // just pre-allocate memory, do not fill in any elements yet
  };
  Spline(const Array1 <double> & xx, const Array1 <double> & yy):n(0) {
    PushBack(xx,yy);
    Init();
  }

  void Allocate(const int m) { // just pre-allocate memory, do not fill in any elements yet
    x.Allocate(m);
    y.Allocate(m);
    y2.Allocate(m);
  }

  inline static void CalculateY2(Array1 <double> & x,
			  Array1 <double> & y, 
			  const int n,
			  Array1 <double> & y2,
			  const bool natural=true, 
			  const double ypLow=0.0, 
			  const double ypHigh=0.0);
  inline void CalculateY2(const bool natural=true, 
		   const double ypLow=0.0, 
		   const double ypHigh=0.0) {
    CalculateY2(x,y,n,y2,natural,ypLow,ypHigh);
  }
  
  inline void Init(const bool natural=true, 
	    const double ypLow=0.0, 
	    const double ypHigh=0.0) {
    CalculateY2(x,y,n,y2,natural,ypLow,ypHigh);
  }
  
  // use finite difference to specify dy/dx at the boundaries
  void InitUsingFiniteDifferences() {
    if (n<=2) error("Spline:InitUsingFiniteDifferences is called with two few data points",n);
    Sort(x,y,n);
    MonotonicX(true);
    double ypLow =(y[1]-y[0])    /(x[1]-x[0]);
    double ypHigh=(y[n-1]-y[n-2])/(x[n-1]-x[n-2]);
    Init(false,ypLow,ypHigh);
  }

  void Sort() {
     Sort(x,y,n);
  }   

  static void Sort(Array1 <double> & x,
		   Array1 <double> & y, 
		   const int n);

  void Reset() {
    n=0;
    x.Resize(0);
    y.Resize(0);
  }
  void Resize(const int nn) {
    n = nn;
    x.Resize(n);
    y.Resize(n);
    y2.Resize(n);
  }
  void PushBack(const double xx, const double yy) {
    n++;
    x.PushBack(xx);
    y.PushBack(yy);
  }

  void PushBack(const Array1 <double> & xx, const Array1 <double> & yy) {
    if (xx.Size()!=yy.Size()) 
      error("Array size error ",xx.Size(),yy.Size());
    for(int i=0; i<xx.Size(); i++) {
      PushBack(xx[i],yy[i]);
    }
  }

  static double Interpolate(const Array1 <double> & x,
			    const Array1 <double> & y, 
			    const int n, 
			    const Array1 <double> & y2,
			    const double xx);
  static double Derivative(const Array1 <double> & x,
			   const Array1 <double> & y, 
			   const int n, 
			   const Array1 <double> & y2,
			   const double xx);
  static double SecondDerivative(const Array1 <double> & x,
				 const Array1 <double> & y, 
				 const int n, 
				 const Array1 <double> & y2,
				 const double xx);
  static double InterpolateLinearly(const Array1 <double> & x,
				    const Array1 <double> & y, 
				    const int n, 
				    const double xx);
  static double InterpolateLinearly(const Array1 <double> & x,
				    const Array1 <double> & y, 
				    const int n, 
				    const double xx,
				    int & k1, double & q1);
  
  double Interpolate(const double xx) const {
    return Interpolate(x,y,n,y2,xx);
  }
  double Derivative(const double xx) const {
    return Derivative(x,y,n,y2,xx);
  }
  double SecondDerivative(const double xx) const {
    return SecondDerivative(x,y,n,y2,xx);
  }
  double InterpolateLinearly(const double xx) const {
    return InterpolateLinearly(x,y,n,xx);
  }
  double InterpolateLinearly(const double xx, int & index, double & fraction1) const {
    return InterpolateLinearly(x,y,n,xx,index,fraction1);
  }

  double f(const double xx) const {
    return Interpolate(xx);
  }

  double operator()(const double xx) const {
    return Interpolate(xx);
  }

  double FirstX() const {
    return x[0];
  }

  double LastX() const {
    return x[x.Size()-1];
  }

  double FirstY() const {
    return y[0];
  }

  double LastY() const {
    return y[y.Size()-1];
  }

  int GetNumber() const {
    return x.Size();
  }

  int Size() const {
    return GetNumber();
  }

  void Print() const {
    for(int i=0; i<GetNumber(); i++) {
      //      Write4(i,x[i],y[i],y2[i]);
      Write3(i,x[i],y[i]);
    }
  }

  void PrintInterpolation(int n=-1) const {
    if (n<=0) n = 5*GetNumber();
    double x1 = FirstX();
    double x2 = LastX();
    for(int i=0; i<=n; i++) {
      double xx = x1+(x2-x1)*double(i)/double(n);
      double yy = Interpolate(xx);
      double yyp = Derivative(xx);
      double yypp = SecondDerivative(xx);
      cout << "spline x= " << xx << " y= " << yy << " y'= " << yyp << " y''= " << yypp << endl;
    }
  }

  // approach x[k1] from right x[k1]+eps but be careful: only good for k1<n-1
  double DerivativeAtPointFromRight(const int k1) const { 
    int k2=k1+1;
    double h=x[k2]-x[k1];
    //  double a=(x[k2]-xx)/h;  --> xx=x[k1] --> a=1
    //  double b=1.0-a;                      --> b=0
    double dfdx = (y[k2]-y[k1])/h - h/6.0 * ( 2.0*y2[k1] + y2[k2] );
    return dfdx;
  }
  // approach x[k2] from left x[k2]-eps but be careful: only good for k2>=1
  double DerivativeAtPointFromLeft(const int k2) const { 
    int k1=k2-1;
    double h=x[k2]-x[k1];
    //  double a=(x[k2]-xx)/h;  --> xx=x[k2] --> a=0
    //  double b=1.0-a;                      --> b=1
    double dfdx = (y[k2]-y[k1])/h + h/6.0 * ( y2[k1] + 2.0*y2[k2] );
    return dfdx;
  }
  double DerivativeAtPoint(const int k) const { 
    //    if (k>0 && k<GetNumber()-1) Write3(k,DerivativeAtPointFromLeft(k),DerivativeAtPointFromRight(k));
    return (k>0) ? DerivativeAtPointFromLeft(k) : DerivativeAtPointFromRight(k);
  }

  bool MonotonicX(const bool printAndStop) {
    for(int i=0; i<GetNumber()-1; i++) {
      if (x[i+1]<=x[i]) {
	if (printAndStop) {
	  error("Spline not monotonic, may have repeated entries",i,x[i],x[i+1]);
	} else {
	  return false;
	}
      }
    }
    return true;
  }

  bool Monotonic(const int nNotCheck=0, const bool print=false) { // must have been initialized
    bool monotonic=true;
    for(int i=nNotCheck; i<GetNumber()-1-nNotCheck; i++) {
      if (sign(DerivativeAtPoint(i))*sign(DerivativeAtPoint(i+1))<0.0) { // either one equal zero woud be ok
	if (!print) return false;
	cout << "Spline not monotonic: "; Write8(i,n,x[i],x[i+1],y[i],y[i+1],DerivativeAtPoint(i),DerivativeAtPoint(i+1));
	PrintInterpolation();
	Print();
	monotonic=false;
      }
    }
    return monotonic;
  }
  
  static double IntegralTerm1(const double x1, const double x2, const double xx);
  static double IntegralTerm2(const double x1, const double x2, const double xx);
  static double IntegralTerm3(const double x1, const double x2, const double xx);
  static double IntegralTerm4(const double x1, const double x2, const double xx);
  static double IntegralInterval1(const double x1, const double x2, const double xx, const double y1, const double y2, const double ys1, const double ys2);
  static double IntegralFullInterval(const double h, const double y1, const double y2, const double ys1, const double ys2);
  static double IntegralInterval2(const double x1, const double x2, const double xx, const double y1, const double y2, const double ys1, const double ys2);

  static double Integrate(const Array1 <double> & x,
			  const Array1 <double> & y, 
			  const int n, 
			  const Array1 <double> & y2,
			  const double xx1, 
			  const double xx2);
  double Integrate(const double xx1, const double xx2) const {
    return Integrate(x,y,n,y2,xx1,xx2);
  }
  static double Average(const Array1 <double> & x,
			const Array1 <double> & y, 
			const int n, 
			const Array1 <double> & y2,
			const double xx1, 
			const double xx2) {
    double dx = xx2 - xx1;
    double sInt = Integrate(x,y,n,y2,xx2,xx1);
    return sInt / dx;
  }
  double Average(const double xx1, const double xx2) const {
    return Average(x,y,n,y2,xx1,xx2);
  }
  // Add "BetweenPoints" to the following function name because
  // because the wrong functions above would be called and
  // int arguments would be converted to double, I think.
  static double IntegrateBetweenPoints(const Array1 <double> & x,
				       const Array1 <double> & y, 
				       const int n, 
				       const Array1 <double> & y2,
				       const int k1, 
				       const int k2);
  double IntegrateBetweenPoints(const int k1, const int k2) const {
    return IntegrateBetweenPoints(x,y,n,y2,k1,k2);
  }
  double IntegrateOverEntireSpline() const {
    //    Write2(k1,k2);
    return IntegrateBetweenPoints(0,n-1);
  }
  static double AverageBetweenPoints(const Array1 <double> & x,
				     const Array1 <double> & y, 
				     const int n, 
				     const Array1 <double> & y2,
				     const int k1, 
				     const int k2);
  double AverageBetweenPoints(const int k1, const int k2) const {
    //    Write2(k1,k2);
    return AverageBetweenPoints(x,y,n,y2,k1,k2);
  }
  double AverageOverEntireSpline() const {
    //    Write2(k1,k2);
    return AverageBetweenPoints(0,n-1);
  }

  // protected:

  static int GetIndex(const Array1 <double> & x, const int n, const double xx);
 
  bool operator==(const Spline & s) const {
    return ((n==s.n) && (x==s.x) && (y==s.y) && (y2==s.y2));
  }
  bool operator!=(const Spline & s) const {
    return !(*this==s);
  }

  friend ostream& operator<<(ostream & os, const Spline & s ) {
    os << "spline.x= " << s.x << endl;
    os << "spline.y= " << s.y; // << endl;
    return os;
  }

  void Check() {
    if (n!=x.Size()) error("Spline::Check() x");
    if (n!=y.Size()) error("Spline::Check() y");
    if (n!=y2.Size()) error("Spline::Check() y2");
    Array1 <double> yCopy(y);
    Array1 <double> y2Copy(y2);
    CalculateY2(x,yCopy,n,y2Copy); // only works if spline was initialized earlier with default arguments of natural,ypLow,ypHigh
    if (y2Copy!=y2) { 
      Write(y2);
      Write(y2Copy);
      error("Spline::Check() y2 different.");
    }
  }

  bool Inside(const double xx, const double eps=0.0) const {
    return (FirstX()-eps<=xx) && (xx<=LastX()+eps); // assume x points have been ordered
  }

};

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

class XY {
public:
  double x,y;
  bool operator< (const XY & xy) const {
    return x<xy.x;
  }
};

inline void Spline::Sort(Array1 <double> & x,
			 Array1 <double> & y, 
			 const int n) {
  bool sortFlag=false;
  for (int i=0; i<n-1; i++) {
    if (x[i]>x[i+1]) {
      sortFlag=true;
      break;
    }
  }
  if (!sortFlag) return;

  Vector <XY> v(n);
  for(int i=0; i<n; i++) {
    v[i].x = x[i];
    v[i].y = y[i];
  }
  sort(v.begin(),v.end());
  for(int i=0; i<n; i++) {
    x[i] = v[i].x;
    y[i] = v[i].y;
  }
  
}

inline void Spline::CalculateY2(Array1 <double> & x,
				Array1 <double> & y, 
				const int n,
				Array1 <double> & y2,
				const bool natural, 
				const double ypLow, 
				const double ypHigh) {
  if (n<2) 
    error("Spline::CalculateY2 call for empty spline",n);
  if (n!=x.Size() || n!=y.Size()) 
    error("Array size problem",n,x.Size(),y.Size());
  
  Sort(x,y,n);
  
  Array1 <double> u(n-1);
  //  u.Resize(n-1);
  
  y2.Resize(n);

  if (natural) {
    u[0] =0.0;
    y2[0]=0.0;
  } else {
    y2[0] = -0.5;
    u[0]  = 3.0/(x[1]-x[0]) * ((y[1]-y[0])/(x[1]-x[0])-ypLow);
  }
  
  for (int i=1; i<n-1; i++) {
    if (x[i-1]>=x[i] || x[i-1]>=x[i+1]) 
      error("Spline:x values are not monotonically increasing",i,x[i-1],x[i],x[i+1]);
    
    const double dpm1 = 1.0/(x[i+1]-x[i-1]);
    //    const double sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
    const double sig=(x[i]-x[i-1])*dpm1;
    const double p  =sig*y2[i-1]+2.0;
    y2[i]=(sig-1.0)/p;
    u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
    //    u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
    u[i]=(6.0*u[i]*dpm1-sig*u[i-1])/p;
  }
  
  double qn,un;
  if (natural) {
    qn= 0.0;
    un= 0.0;
  } else {
    qn= 0.5;
    un= 3.0/(x[n-1]-x[n-2]) * (ypHigh-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
  }

  y2[n-1] = (un-qn*u[n-2]) / (qn*y2[n-2]+1.0);
  for (int k=n-2; k>=0; k--)
    y2[k]=y2[k]*y2[k+1]+u[k];

  //  y2=0.0;
  //  warning("0.0=y2");
}

inline int Spline::GetIndex(const Array1 <double> & x, const int n, const double xx) {
  int k1=0;
  int k2=n-1;

// #define NOT_INSIDE_WARNING_SPLINE
#ifdef NOT_INSIDE_WARNING_SPLINE
  double dx = 1e-6 * (x[k2] - x[k1]); // small step but not too small
  //  if (xx<x[k1]-dx || xx>x[k2]+dx) warning("Spline argument outside",xx,x[k1],x[k2]);
  if (xx<x[k1]-dx || xx>x[k2]+dx) error("Spline argument outside",xx,x[k1],x[k2]);
#endif

  while (k2-k1 > 1) {
    int k=(k2+k1) >> 1;
    if (x[k] > xx) k2=k;
    else k1=k;
  }

  if (x[k2] == x[k1]) error("SplineInterpolation: Bad x-vector input",k1,k2,n,x[k1],x[k2],xx);
  return k1;
}

inline double Spline::Interpolate(const Array1 <double> & x,
				  const Array1 <double> & y, 
				  const int n, 
				  const Array1 <double> & y2,
				  const double xx) {
  int k1=GetIndex(x,n,xx);
  int k2=k1+1;
  double h=x[k2]-x[k1];
  double a=(x[k2]-xx)/h;
  double b=1.0-a;
  return a*y[k1] + b*y[k2] + 
    ( (a*a*a-a)*y2[k1]+(b*b*b-b)*y2[k2] ) * (h*h)/6.0;
}

inline double Spline::InterpolateLinearly(const Array1 <double> & x,
					  const Array1 <double> & y, 
					  const int n, 
					  const double xx) {
  int k1=GetIndex(x,n,xx);
  int k2=k1+1;
  double h=x[k2]-x[k1];
  double q1 = (x[k2] - xx   )/h;
  //  double q2 = (xx    - x[k1])/h;
  //  return q1*y[k1] + q2*y[k2];
  return y[k2] + q1*(y[k1]-y[k2]);
}

inline double Spline::InterpolateLinearly(const Array1 <double> & x,
					  const Array1 <double> & y, 
					  const int n, 
					  const double xx, int & k1, double & q1) {
  k1=GetIndex(x,n,xx);
  int k2=k1+1;
  double h=x[k2]-x[k1];
  q1 = (x[k2] - xx   )/h;
  return y[k2] + q1*(y[k1]-y[k2]);
}

inline double Spline::Derivative(const Array1 <double> & x,
				 const Array1 <double> & y, 
				 const int n, 
				 const Array1 <double> & y2,
				 const double xx) {
  int k1=GetIndex(x,n,xx);
  int k2=k1+1;
  double h=x[k2]-x[k1];
  double a=(x[k2]-xx)/h;
  double b=1.0-a;
  
  return (y[k2]-y[k1])/h + h/6.0* 
    ( (1.0-3.0*a*a)*y2[k1] + (3.0*b*b-1.0)*y2[k2] ); // bug fixed 6/18/06: y2[k2] instead of y2[k1]
}

inline double Spline::SecondDerivative(const Array1 <double> & x,
				       const Array1 <double> & y, 
				       const int n, 
				       const Array1 <double> & y2,
				       const double xx) {
  int k1=GetIndex(x,n,xx);
  int k2=k1+1;
  double h=x[k2]-x[k1];
  double a=(x[k2]-xx)/h;
  double b=(xx-x[k1])/h;
  
  return a*y2[k1]+b*y2[k2];
}

inline double Spline::IntegralTerm1(const double x1, const double x2, const double xx) {
  //  Simplify[Integrate[a*y1,{x,x1,xx}]]
  //         -((x1 - xx) (x1 - 2 x2 + xx) y1)
  //Out[37]= --------------------------------
  //                   2 (x1 - x2)
  return 0.5*(xx-x1)/(x1-x2) * (x1 - 2.0*x2 + xx);
}

inline double Spline::IntegralTerm2(const double x1, const double x2, const double xx) {
  //Simplify[Integrate[b*y2,{x,x1,xx}]]
  //                    2
  //         -((x1 - xx)  y2)
  //Out[38]= ----------------
  //           2 (x1 - x2)
  return sqr(xx-x1)/(2.0*(x2-x1));
}

inline double Spline::IntegralTerm3(const double x1, const double x2, const double xx) {
  //In[32]:= Integrate[(a^3-a)*ys1*h*h/6,{x,x1,xx}]
  //                  2                 2
  //         (x1 - xx)  (x1 - 2 x2 + xx)  ys1
  //Out[32]= --------------------------------
  //                   24 (x1 - x2)
  return sqr(xx-x1)*sqr(x1 - 2.0*x2 + xx) / ( 24.0 * (x1-x2) );
}

inline double Spline::IntegralTerm4(const double x1, const double x2, const double xx) {
  //In[33]:= Integrate[(b^3-b)*ys2*h*h/6,{x,x1,xx}]
  //                  2    2                 2               2
  //         (x1 - xx)  (x1  - 4 x1 x2 + 2 x2  + 2 x1 xx - xx ) ys2
  //Out[33]= ------------------------------------------------------
  //                              24 (x1 - x2)
  return sqr(xx-x1)*(x1*x1 - 4.0*x1*x2 + 2.0*x2*x2 + 2.0*x1*xx - xx*xx) / ( 24.0 * (x1-x2) );
}

// integral dx_{x1...xx}
inline double Spline::IntegralInterval1(const double x1, const double x2, const double xx, const double y1, const double y2, const double ys1, const double ys2) {
  return 
    y1 *IntegralTerm1(x1,x2,xx)+y2 *IntegralTerm2(x1,x2,xx)+
    ys1*IntegralTerm3(x1,x2,xx)+ys2*IntegralTerm4(x1,x2,xx);
}

// integral dx_{x1...x2}
inline double Spline::IntegralFullInterval(const double h, const double y1, const double y2, const double ys1, const double ys2) {
  //  In[39]:= Integrate[f,{x,x1,x2}]
  //                                              2
  //         (x1 - x2) (-12 y1 - 12 y2 + (x1 - x2)  (ys1 + ys2))
  //Out[39]= ---------------------------------------------------
  //                                 24
  return h * ( 12.0*(y1+y2) - h*h*(ys1+ys2) ) / 24.0;
}

// integral dx_{xx...x2}
inline double Spline::IntegralInterval2(const double x1, const double x2, const double xx, const double y1, const double y2, const double ys1, const double ys2) {
  return 
    IntegralFullInterval(x2-x1,y1,y2,ys1,ys2) - IntegralInterval1(x1,x2,xx,y1,y2,ys1,ys2);
}

inline double Spline::Integrate(const Array1 <double> & x,
				const Array1 <double> & y, 
				const int n, 
				const Array1 <double> & y2,
				const double xx1, 
				const double xx2) {
  if (xx2<xx1) return -Integrate(x,y,n,y2,xx2,xx1);

  int k1=GetIndex(x,n,xx1);
  int k2=GetIndex(x,n,xx2); // definition differs from above!
  int kk1 = k1+1;
  int kk2 = k2+1;

  //  Write2(k1,k2);

  if (k1==k2) {
    /*
    int nn=100000;
    double ss=0.0;
    for(int i=0; i<nn; i++) {
      ss+=Interpolate(x,y,n,y2,xx1+(xx2-xx1)*(0.5+i)/nn)*(xx2-xx1)/nn;
    }
    Write6(xx1,xx2,k1,k2,n,ss);
    */

    //    double f1 = IntegralInterval1(x[k1],x[kk1],xx2,y[k1],y[kk1],y2[k1],y2[kk1])
    //      -IntegralInterval1(x[k1],x[kk1],xx1,y[k1],y[kk1],y2[k1],y2[kk1]);
    //    double f2 = IntegralFullInterval(x[kk1]-x[k1],y[k1],y[kk1],y2[k1],y2[kk1]);
    //    Write2(f1,f2);

    return 
       IntegralInterval1(x[k1],x[kk1],xx2,y[k1],y[kk1],y2[k1],y2[kk1])
      -IntegralInterval1(x[k1],x[kk1],xx1,y[k1],y[kk1],y2[k1],y2[kk1]);
  }
  double sInt = IntegralInterval2(x[k1],x[kk1],xx1,y[k1],y[kk1],y2[k1],y2[kk1]);
  sInt += IntegralInterval1(x[k2],x[kk2],xx2,y[k2],y[kk2],y2[k2],y2[kk2]);

  for(int k=k1+1; k<k2; k++) {
    int kk = k+1;
    sInt += IntegralFullInterval(x[kk]-x[k],y[k],y[kk],y2[k],y2[kk]);
  }

  /*
  int nn=100000;
  double ss=0.0;
  for(int i=0; i<nn; i++) {
    ss+=Interpolate(x,y,n,y2,xx1+(xx2-xx1)*(0.5+i)/nn)*(xx2-xx1)/nn;
  }
  Write7(xx1,xx2,k1,k2,n,sInt,ss);
  */

  return sInt;
}

inline double Spline::IntegrateBetweenPoints(const Array1 <double> & x,
				const Array1 <double> & y, 
				const int n, 
				const Array1 <double> & y2,
				const int k1, 
				const int k2) {
  double sInt = 0.0;
  for(int k=k1; k<k2; k++) {
    int kk = k+1;
    sInt += IntegralFullInterval(x[kk]-x[k],y[k],y[kk],y2[k],y2[kk]);
  }
  return sInt;
}

inline double Spline::AverageBetweenPoints(const Array1 <double> & x,
			      const Array1 <double> & y, 
			      const int n, 
			      const Array1 <double> & y2,
			      const int k1, 
			      const int k2) {
  double dx = x[k2] - x[k1];
  double sInt = IntegrateBetweenPoints(x,y,n,y2,k1,k2);
  //  Write6(k1,k2,x[k1],x[k2],dx,sInt);
  return sInt / dx;
}

/*
int main() {
  const int n=10;
  Array1 <double> x(n);
  Array1 <double> y(n);

  Spline s;
  for(int i=0; i<n; i++) {
    x[i] = i*1.5;
    //    x[i] = i*1.0;
    y[i] = sin(i)+cos(10*i);
    Write2(x[i],y[i]);
    s.PushBack(x[i],y[i]);
  }


  Array1 <double> y2;
  Spline ss;
  s.CalculateY2(x,y,n,y2);
  //  ss.CalculateY2(x,y,n,y2,false,5,5);
  s.Init();
  
  //  Write(s.Integrate(x[1],x[2]*1.5));
  //  Write(s.Integrate(x[1]*0.1,1.5));
  Write(s.Integrate(0.5,13.0+1e-15));
  error("Q");

  const int m=10;
  for(int i=0; i<n*m; i++) {
    double xx = double(i)/double(m);
    double yy   = ss.Interpolate(x,y,n,y2,xx);
    double yyd  = ss.Derivative(x,y,n,y2,xx);
    double yydd = ss.SecondDerivative(x,y,n,y2,xx);
    double y   = s.Interpolate(xx);
    double yd  = s.Derivative(xx);
    double ydd = s.SecondDerivative(xx);
    double yII = s.Integrate(x[1],x[n/2]);
    Write4(xx,yy,yyd,yydd);
    Write5(xx,y,yd,ydd,yII);
  }
}
*/


#endif // _SPLINE_