TRadSpline.h

Go to the documentation of this file.

#ifndef ROOT_TRadSpline
#define ROOT_TRadSpline
 
class TRadSpline {
protected:
  double  fDelta;     // Distance between equidistant knots
  double  fXmin;      // Minimum value of abscissa
  double  fXmax;      // Maximum value of abscissa
  int     fNp;        // Number of knots
  bool    fKstep;     // True of equidistant knots
  int     fNpx;       // Number of points used for graphical representation
  
  virtual void BuildCoeff()=0;
  
public:
  TRadSpline() : fDelta(-1), fXmin(0), fXmax(0),
          fNp(0), fKstep(false), fNpx(100) {}
  TRadSpline(const char *title, double delta, double xmin,
      double xmax, int np, bool step) :
    fDelta(delta), fXmin(xmin),
    fXmax(xmax), fNp(np), fKstep(step),
    fNpx(100) {}
  virtual ~TRadSpline() {}
  virtual  void     GetKnot(int i, double &x, double &y) const =0;
  //  virtual  void     Draw(const char *option="");
  virtual  int    GetNpx() const {return fNpx;}
  //  virtual  void     Paint(const char *option="");
  virtual double  Eval(double x) const=0;
  virtual  void     SaveAs(const char * /*filename*/) const {;}
           void     SetNpx(int n) {fNpx=n;}
 
};
 
 
//________________________________________________________________________
class TRadSplinePoly {
protected:
  double fX;     // abscissa
  double fY;     // constant term
public:
  TRadSplinePoly() :
    fX(0), fY(0) {}
  TRadSplinePoly(double x, double y) :
    fX(x), fY(y) {}
  virtual ~TRadSplinePoly() {}
  double &X() {return fX;}
  double &Y() {return fY;}
  void GetKnot(double &x, double &y) const {x=fX; y=fY;}
 
  virtual double Eval(double) const {return fY;}
 
};
 
 
//________________________________________________________________________
class TRadSplinePoly3 : public TRadSplinePoly {
private:
  double fB; // first order expansion coefficient :  fB*1! is the first derivative at x
  double fC; // second order expansion coefficient : fC*2! is the second derivative at x
  double fD; // third order expansion coefficient :  fD*3! is the third derivative at x
public:
  TRadSplinePoly3() :
    fB(0), fC(0), fD(0) {}
  TRadSplinePoly3(double x, double y, double b, double c, double d) :
    TRadSplinePoly(x,y), fB(b), fC(c), fD(d) {}
  double &B() {return fB;}
  double &C() {return fC;}
  double &D() {return fD;}
  double Eval(double x) const {
    double dx=x-fX;
    return (fY+dx*(fB+dx*(fC+dx*fD)));
  }
  double Derivative(double x) const {
    double dx=x-fX;
    return (fB+2*fC*dx+3*fD*dx*dx);
  }
 
};
 
 
//________________________________________________________________________
class TRadSplinePoly5 : public TRadSplinePoly {
private:
  double fB; // first order expansion coefficient :  fB*1! is the first derivative at x
  double fC; // second order expansion coefficient : fC*2! is the second derivative at x
  double fD; // third order expansion coefficient :  fD*3! is the third derivative at x
  double fE; // fourth order expansion coefficient : fE*4! is the fourth derivative at x
  double fF; // fifth order expansion coefficient :  fF*5! is the fifth derivative at x
public:
  TRadSplinePoly5() :
    fB(0), fC(0), fD(0), fE(0), fF(0) {}
  TRadSplinePoly5(double x, double y, double b, double c,
     double d, double e, double f) :
    TRadSplinePoly(x,y), fB(b), fC(c), fD(d), fE(e), fF(f) {}
  double &B() {return fB;}
  double &C() {return fC;}
  double &D() {return fD;}
  double &E() {return fE;}
  double &F() {return fF;}
  double Eval(double x) const {
    double dx=x-fX;
    return (fY+dx*(fB+dx*(fC+dx*(fD+dx*(fE+dx*fF)))));
  }
  double Derivative(double x) const{
    double dx=x-fX;
    return (fB+2*fC*dx+3*fD*dx*dx+4*fE*dx*dx*dx+5*fF*dx*dx*dx*dx);
  }
 
};
 
 
//________________________________________________________________________
class TRadSpline3 : public TRadSpline {
private:
  TRadSplinePoly3  *fPoly;       //[fNp] Array of polynomial terms
  double       fValBeg;     // Initial value of first or second derivative
  double       fValEnd;     // End value of first or second derivative
  int          fBegCond;    // 0=no beg cond, 1=first derivative, 2=second derivative
  int          fEndCond;    // 0=no end cond, 1=first derivative, 2=second derivative
 
  void BuildCoeff();
  void SetCond(const char *opt);
 
public:
  TRadSpline3() : fPoly(0), fValBeg(0), fValEnd(0),
    fBegCond(-1), fEndCond(-1) {}
  TRadSpline3(const char *title,
       double x[], double y[], int n, const char *opt=0,
       double valbeg=0, double valend=0);
  TRadSpline3(const char *title,
       double xmin, double xmax,
       double y[], int n, const char *opt=0,
       double valbeg=0, double valend=0);
  TRadSpline3(const char *title,
          double xmin, double xmax,
          double (*func)(const double&), int n, const char *opt=0,
          double valbeg=0, double valend=0);
  int    FindX(double x) const;
  double Eval(double x) const;
  double Derivative(double x) const;
  virtual ~TRadSpline3() {if (fPoly) delete [] fPoly;}
  void GetCoeff(int i, double &x, double &y, double &b,
        double &c, double &d) {x=fPoly[i].X();y=fPoly[i].Y();
        b=fPoly[i].B();c=fPoly[i].C();d=fPoly[i].D();}
  void GetKnot(int i, double &x, double &y) const
    {x=fPoly[i].X(); y=fPoly[i].Y();}
  virtual  void     SaveAs(const char *filename) const;
  static void Test();
 
};
 
 
//________________________________________________________________________
class TRadSpline5 : public TRadSpline {
private:
  TRadSplinePoly5  *fPoly;     //[fNp] Array of polynomial terms
 
  void BuildCoeff();
  void BoundaryConditions(const char *opt, int &beg, int &end,
              const char *&cb1, const char *&ce1, const char *&cb2,
                          const char *&ce2);
  void SetBoundaries(double b1, double e1, double b2, double e2,
             const char *cb1, const char *ce1, const char *cb2,
                     const char *ce2);
 
public:
  TRadSpline5() : fPoly(0) {}
  TRadSpline5(const char *title,
       double x[], double y[], int n,
       const char *opt=0, double b1=0, double e1=0,
       double b2=0, double e2=0);
  TRadSpline5(const char *title,
       double xmin, double xmax,
       double y[], int n,
       const char *opt=0, double b1=0, double e1=0,
       double b2=0, double e2=0);
  TRadSpline5(const char *title,
          double xmin, double xmax,
          double(*func)(const double&), 
          int n, const char *opt=0,
          double b1=0, double e1=0,
          double b2=0, double e2=0);
  int    FindX(double x) const;
  double Eval(double x) const;
  double Derivative(double x) const;
  ~TRadSpline5() {if (fPoly) delete [] fPoly;}
  void GetCoeff(int i, double &x, double &y, double &b,
        double &c, double &d, double &e, double &f)
    {x=fPoly[i].X();y=fPoly[i].Y();b=fPoly[i].B();
    c=fPoly[i].C();d=fPoly[i].D();
    e=fPoly[i].E();f=fPoly[i].F();}
  void GetKnot(int i, double &x, double &y) const
    {x=fPoly[i].X(); y=fPoly[i].Y();}
  virtual  void     SaveAs(const char *filename) const;
  static void Test();
 
};
 
#endif