//////////////////////////////////////////////////////////////////////////////
//                                                                          //
// Modified Spline interpolation with linear extrapolation                  //
//                                                                          //
// Burkhard Militzer                                 Washington 08-03-04    //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////

#ifndef _LSPLINE_
#define _LSPLINE_

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

class LSpline {
 public:
  
  int n;
  Array1 <double> x;
  Array1 <double> y;
  Array1 <double> y2;
  static Array1 <double> u;
  
  LSpline():n(0){};
  LSpline(const int m):n(0) { 
    Allocate(m); // just pre-allocate memory, do not fill in any elements yet
  };
  LSpline(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);
  void CalculateY2(const bool natural=true, 
		   const double ypLow=0.0, 
		   const double ypHigh=0.0) {
    CalculateY2(x,y,n,y2,natural,ypLow,ypHigh);
  }
  
  // this class is supposed to be initialized InitUsingFiniteDifferences() 
  // could make this private so that we do not call it by mistake
  void Init(const bool natural=true, 
	    const double ypLow=0.0, 
	    const double ypHigh=0.0) {
    CalculateY2(x,y,n,y2,natural,ypLow,ypHigh);
  }
  
  public:
  // 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);
    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);
  }

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

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

  void Reset() {
    n=0;
    x.Resize(0);
    y.Resize(0);
  }

  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 InterpolateLinearly(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);
  
  double Interpolate(const double xx) const {
    return Interpolate(x,y,n,y2,xx);
  }
  static double LinearInterpolation(const double x1, const double x2, 
				    const double y1, const double y2, 
				    const double xx) {
    double q2 = (xx-x1)/(x2-x1);
    double q1 = (x2-xx)/(x2-x1);
    return q1*y1+q2*y2;
  }

  double InterpolateLinearly(const double xx) const {
    return InterpolateLinearly(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 f(const double xx) const {
    return Interpolate(xx);
  }

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

  bool Inside(const double xx, const double eps=0.0) const {
    //    Write4(xx,FirstX(),LastX(),eps);
    if (xx*(1.0+sign(xx)*eps)<FirstX()) return false;
    if (xx*(1.0-sign(xx)*eps)>LastX()) return false;
    return true;
  }

  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 Size() const {
    return GetNumber();
  }

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

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

  void WriteToFile(const string & fileName, const int p=0) const {
    ofstream of(fileName.c_str());
    if (p>0) of.precision(p);
    if (!of) error("Could not open file",fileName);
    for(int i=0; i<GetNumber(); i++) {
      of << "i= " << i << " x= " << x[i] << " y= " << y[i] << endl;
      //      Write3(i,x[i],y[i]);
    }
    of.close();
  }
  
  void PrintInterpolation(int nn=-1) const { // argument 'nn' should not shadow variable 'n'
    double x1 = FirstX();
    double x2 = LastX();
    if (nn<=0) nn = 5*GetNumber();
    int iMin =   -nn/5;
    int iMax = nn+nn/5;
    //    for(int i=0; i<=nn; i++) {
    for(int i=iMin; i<=iMax; i++) {
      double xx   = x1+(x2-x1)*double(i)/double(nn);
      double yy   = Interpolate(xx);
      double yyp  = Derivative(xx);
      double yypp = SecondDerivative(xx);
      cout << "LSpline x= " << xx << " y= " << yy << " y'= " << yyp << " y''= " << yypp << endl;
    }
  }

  void WriteInterpolationToFile(const string & fileName, int nn=-1) {
    ofstream of(fileName.c_str());
    if (!of) error("Could not open file",fileName);
    double x1 = FirstX();
    double x2 = LastX();
    if (nn<=0) nn = 5*GetNumber();
    int iMin =   -nn/5;
    int iMax = nn+nn/5;
    //    for(int i=0; i<=nn; i++) {
    for(int i=iMin; i<=iMax; i++) {
      double xx   = x1+(x2-x1)*double(i)/double(nn);
      double yy   = Interpolate(xx);
      double yyp  = Derivative(xx);
      double yypp = SecondDerivative(xx);
      of << "i= " << i << " x= " << xx << " y= " << yy << " y'= " << yyp << " y''= " << yypp << endl;
      //      Write3(i,x[i],y[i]);
    }
    of.close();
  }
  // protected:

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

  friend ostream& operator<<(ostream & os, const LSpline & s ) {
    os << "lSpline.x= " << s.x << endl;
    os << "lSpline.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.");
    }
  }
  */

  // 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++) {
      double s1 = sign(DerivativeAtPoint(i)  );
      double s2 = sign(DerivativeAtPoint(i+1));
      if (s1*s2<0.0) { // either one equal zero would 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;
  }

  bool PointsMonotonic(const int nNotCheck=0, const bool print=false) { // must have been initialized
    bool monotonic=true;
    for(int i=nNotCheck; i<GetNumber()-2-nNotCheck; i++) {
      double s1 = sign(y[i+1]-y[i]  );
      double s2 = sign(y[i+2]-y[i+1]);
      if (s1*s2<0.0) { // either one equal zero would be ok
	if (!print) return false;
	cout << "Spline points 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 IntegralOfLinearFunction(const double x1, const double x2, 
					 const double y1, const double y2, 
					 const double xx) {
    double h = x2-x1;
    return (0.5*xx*xx*(y2-y1) + xx*(x2*y1-x1*y2))/h;
  }

  static double IntegralOfLinearFunction(const double x1, const double x2, 
					 const double y1, const double y2, 
					 const double xx1, const double xx2) {
    return IntegralOfLinearFunction(x1,x2,y1,y2,xx2) - IntegralOfLinearFunction(x1,x2,y1,y2,xx1);
  }
  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);
  }

  double SolveForXLinearly(const double x1, const double y1, const double x2, const double y2, const double yy) const {
    const double dy1 = y1-yy;
    const double dy2 = y2-yy;
    return (dy2*x1-dy1*x2)/(dy2-dy1);
  }
  double SolveForXAtTheBeginning(const double yy) const {
    return SolveForXLinearly(x[0],y[0],x[1],y[1],yy);
  }
  double SolveForXAtTheEnd(const double yy) const {
    return SolveForXLinearly(x[n-2],y[n-2],x[n-1],y[n-1],yy);
  }
};

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

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

inline void LSpline::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 <XYL> 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 LSpline::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);
  
  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 monotonic",i,x[i-1],x[i],x[i+1]);

    const double sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
    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;
  }
  
  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];
}

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

// #define NOT_INSIDE_WARNING_LSPLINE
#ifdef NOT_INSIDE_WARNING_LSPLINE
  double dx = 1e-6 * (x[k2] - x[k1]); // small step but not too small
  //  if (xx<x[k1]-dx || xx>x[k2]+dx) warning("LSpline argument outside",xx,x[k1],x[k2]);
  if (xx<x[k1]-dx || xx>x[k2]+dx) error("LSpline 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 LSpline::Interpolate(const Array1 <double> & x,
			   const Array1 <double> & y, 
			   const int n, 
			   const Array1 <double> & y2,
			   const double xx) {
  if (xx<x[0])   return LinearInterpolation(x[0],  x[1],  y[0],  y[1],  xx);
  if (xx>x[n-1]) return LinearInterpolation(x[n-2],x[n-1],y[n-2],y[n-1],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 LSpline::InterpolateLinearly(const Array1 <double> & x,
				    const Array1 <double> & y, 
				    const int n, 
				    const Array1 <double> & y2,
				    const double xx) {
  if (xx<x[0])   return LinearInterpolation(x[0],  x[1],  y[0],  y[1],  xx);
  if (xx>x[n-1]) return LinearInterpolation(x[n-2],x[n-1],y[n-2],y[n-1],xx);
  int k1=GetIndex(x,n,xx);
  int k2=k1+1;
  return LinearInterpolation(x[k1],  x[k2],  y[k1],  y[k2],  xx);
}

inline double LSpline::Derivative(const Array1 <double> & x,
			  const Array1 <double> & y, 
			  const int n, 
			  const Array1 <double> & y2,
			  const double xx) {
  if (xx<x[0])   return (y[1]  -y[0]  )/(x[1]  -x[0]  );
  if (xx>x[n-1]) return (y[n-1]-y[n-2])/(x[n-1]-x[n-2]);
  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 LSpline::SecondDerivative(const Array1 <double> & x,
				const Array1 <double> & y, 
				const int n, 
				const Array1 <double> & y2,
				const double xx) {
  if (xx<x[0])   return 0.0;
  if (xx>x[n-1]) return 0.0;

  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 LSpline::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 LSpline::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 LSpline::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 LSpline::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 LSpline::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 LSpline::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 LSpline::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 LSpline::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);

  if (xx2>x[n-1]) { // upper limit outside spline interval, linear extrapolation applies
    if (xx1>x[n-1]) { // both above upper spline limit
      return IntegralOfLinearFunction(x[n-2],x[n-1],y[n-2],y[n-1],xx1,xx2);
    }
    return IntegralOfLinearFunction(x[n-2],x[n-1],y[n-2],y[n-1],x[n-1],xx2)+Integrate(x,y,n,y2,xx1,x[n-1]); // takes care of xx1<x[0] also
  }

  if (xx1<x[0]) { // lower limit outside spline interval, linear extrapolation applies
    if (xx2<x[0]) { // both below lower spline limit
       return IntegralOfLinearFunction(x[0],x[1],y[0],y[1],xx1,xx2);
    }
    return IntegralOfLinearFunction(x[0],x[1],y[0],y[1],xx1,x[0])+Integrate(x,y,n,y2,x[0],xx2); // case xx2>x[n-1] should have been handle earlier
  }
     
  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;
}

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

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

  for(int i=0; i<n; i++) {
    x[i] = i;
    y[i] = sin(i)+cos(10*i);
    Write2(x[i],y[i]);
  }


  Array1 <double> y2;
  Spline s;
  //  s.CalculateY2(x,y,n,y2);
  s.CalculateY2(x,y,n,y2,false,5,5);

  const int m=10;
  for(int i=0; i<n*m; i++) {
    double xx = double(i)/double(m);
    double yy = s.Interpolate(x,y,n,y2,xx);
    double yyd = s.Derivative(x,y,n,y2,xx);
    double yydd = s.SecondDerivative(x,y,n,y2,xx);
    Write4(xx,yy,yyd,yydd);
  }
}
*/


#endif // _LSPLINE_