#ifndef CSPLINE_H_SEEN
#define CSPLINE_H_SEEN

#include "machdefs.h"
#include "pexceptions.h"

///////////////////////////////////////////////////////////////////
// Class CSpline et CSpline2 selon Numerical Receipes (cmv 17/02/97)
///////////////////////////////////////////////////////////////////
namespace SOPHYA {

//! Spline fit to a set of points Y=f(X) (as in Numerical Receipes).
class CSpline {

public:

  friend class CSpline2;

  enum { NaturalAll = 3, NoNatural = 0, Natural1 = 1, NaturalN = 2 };

  CSpline(int n,double* x,double* y,double yp1=0.,double ypn=0.
         ,int natural=NaturalAll,bool order=true);

  CSpline(double yp1=0.,double ypn=0.,int natural=NaturalAll);

  virtual ~CSpline();

  void SetNewTab(int n,double* x,double* y
                ,bool order=true,bool force=false);

  void SetBound1er(double yp1 = 0.,double yp2 = 0.);

  inline void SetNaturalCSpline(int type = NaturalAll)
                 { Natural = type;}

  inline void Free_Tmp()
                { if(tmp != NULL) delete [] tmp; tmp=NULL;}

  void ComputeCSpline();

  double CSplineInt(double x);

protected:

  void DelTab();

  // nombre d elements dans les tableaux X et Y
  int Nel;
  // true si les tableaux ont ete changes
  // et qu il faut recalculer ComputeCSpline()
  bool corrupt_Y2;
  // true si les tableaux X,Y ont ete alloues
  bool XY_Created;
  // type de contraintes sur la derivee 2sd
  int Natural;
  // valeurs imposees de la derivee 1ere aux limites
  double YP1, YPn;

  // tableaux rellement alloues si "order=true" ou seulement
  // connectes aux tableaux externes si "order=false"
  double* X;
  double* Y;

  // tableau des coeff permettant l interpolation,
  // remplis par ComputeCSpline()
  double* Y2;

  // tableau tampon utilise dans ComputeCSpline()
  double* tmp;
  int_4*  ind;

};


///////////////////////////////////////////////////////////////////
//! 2D Spline fit to a set of points Y=f(X1,X2) (as in Numerical Receipes).
class CSpline2  {

public:
  CSpline2(int n1,double* x1,int n2,double* x2,double* y
          ,int natural=CSpline::NaturalAll,bool order=true);

  CSpline2(int natural=CSpline::NaturalAll);

  virtual ~CSpline2();

  void SetNewTab(int n1,double* x1,int n2,double* x2,double* y
                ,bool order=true,bool force=false);

  inline void SetNaturalCSpline(int type = CSpline::NaturalAll)
                 { Natural = type;}

  void ComputeCSpline();

  double CSplineInt(double x1,double x2);

protected:
  void DelTab();

  int Nel1, Nel2;
  bool corrupt_Y2;
  bool XY_Created;
  int Natural;

  double* X1;
  double* X2;
  double* Y;
  double* Y2;

  // nombre d elements alloues dans S
  int Nel_S;
  // tableau de CSpline pour interpolation selon X1
  CSpline** S;       // S[0->n2]
  CSpline*  Sint;

  // tableau tampon utilise dans CSplineInt()
  double* tmp;  // tmp[max(n1,n2)]
  int_4*  ind;

};

} // Fin du namespace

#endif /* CSPLINE_H_SEEN */