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