EvtPolInt Class Reference

#include <EvtPolInt.hh>

Public Member Functions
	EvtPolInt (vector< double > xi, vector< double > yi, double x)
virtual	~EvtPolInt ()
void	polynomial ()
void	ratint ()
vector< double >	spline ()
void	splint ()
double	getvalue ()
double	geterror ()

Detailed Description

Definition at line 30 of file EvtPolInt.hh.

Constructor & Destructor Documentation

EvtPolInt()

EvtPolInt::EvtPolInt ( vector< double > xi, vector< double > yi, double x )

inline

Definition at line 32 of file EvtPolInt.hh.

                                                              {
    xx   = x;
    size = xi.size();
    vx.clear();
    vy.clear();
    for ( int i = 0; i < size; i++ )
    {
      vx.push_back( xi[i] );
      vy.push_back( yi[i] );
    }
  }

~EvtPolInt()

virtual EvtPolInt::~EvtPolInt ( )

inlinevirtual

Definition at line 43 of file EvtPolInt.hh.

{}

Member Function Documentation

geterror()

double EvtPolInt::geterror

(

)

Definition at line 186 of file EvtPolInt.cc.

                           {
  polynomial();
  // ratint();
  return error;
}

getvalue()

double EvtPolInt::getvalue

(

)

Definition at line 179 of file EvtPolInt.cc.

                           {
  polynomial();
  //  ratint();
  //  splint();
  return value;
}

polynomial()

void EvtPolInt::polynomial

(

)

Definition at line 23 of file EvtPolInt.cc.

                           {
  double y, dy;
  int    i, m, ns = 0;
  double den, dif, dift, ho, hp, w;
 
  int            n = vx.size();
  vector<double> c( n ), d( n );
  dif = fabs( xx - vx[0] );
  for ( i = 0; i < n; i++ )
  {
    if ( ( dift = fabs( xx - vx[i] ) ) < dif )
    {
      ns  = i;
      dif = dift;
    }
    c[i] = vy[i];
    d[i] = vy[i];
  }
  y = vy[ns--];
  for ( m = 1; m < n; m++ )
  {
    for ( i = 0; i < n - m; i++ )
    {
      ho = vx[i] - xx;
      hp = vx[i + m] - xx;
      w  = c[i + 1] - d[i];
      if ( ( den = ho - hp ) == 0.0 ) std::cout << "Error in routine polint" << std::endl;
      den  = w / den;
      d[i] = hp * den;
      c[i] = ho * den;
    }
    y += ( dy = ( 2 * ( ns + 1 ) < ( n - m ) ? c[ns + 1] : d[ns--] ) );
  }
  value = y;
  error = dy;
  // std::cout<<"value= "<<value<<std::endl;
  if ( value < 0 ) value = 0;
  return;
}

Referenced by geterror(), and getvalue().

ratint()

void EvtPolInt::ratint

(

)

Definition at line 63 of file EvtPolInt.cc.

                       {
  double       y, dy;
  const double TINY = 1.0e-25;
  int          m, i, ns = 0;
  double       w, t, hh, h, dd;
 
  int            n = vx.size();
  vector<double> c( n ), d( n );
  hh = fabs( xx - vx[0] );
  for ( i = 0; i < n; i++ )
  {
    h = fabs( xx - vx[i] );
    if ( h == 0.0 )
    {
      y  = vy[i];
      dy = 0.0;
      return;
    }
    else if ( h < hh )
    {
      ns = i;
      hh = h;
    }
    c[i] = vy[i];
    d[i] = vy[i] + TINY;
  }
  y = vy[ns--];
  for ( m = 1; m < n; m++ )
  {
    for ( i = 0; i < n - m; i++ )
    {
      w  = c[i + 1] - d[i];
      h  = vx[i + m] - xx;
      t  = ( vx[i] - xx ) * d[i] / h;
      dd = t - c[i + 1];
      if ( dd == 0.0 ) std::cout << "Error in routine ratint" << std::endl;
      dd   = w / dd;
      d[i] = c[i + 1] * dd;
      c[i] = t * dd;
    }
    y += ( dy = ( 2 * ( ns + 1 ) < ( n - m ) ? c[ns + 1] : d[ns--] ) );
  }
  value = y;
  error = dy;
  if ( value < 0 ) value = 0;
  return;
}

spline()

vector< double > EvtPolInt::spline

(

)

Definition at line 111 of file EvtPolInt.cc.

                                 {
  vector<double> y2;
  y2.clear();
  for ( int i = 0; i < size; i++ ) { y2.push_back( 0.0 ); }
  double yp1 = 0;
  double ypn = 0;
  int    i, k;
  double p, qn, sig, un;
 
  int            n = size;
  vector<double> u; //(n-1);
  if ( yp1 > 0.99e30 )
  {
    u[0] = 0.0;
    y2.push_back( 0.0 );
  }
  else
  {
    y2.push_back( -0.5 );
    u[0] = ( 3.0 / ( vx[1] - vx[0] ) ) * ( ( vy[1] - vy[0] ) / ( vx[1] - vx[0] ) - yp1 );
  }
  for ( i = 1; i < n - 1; i++ )
  {
    sig        = ( vx[i] - vx[i - 1] ) / ( vx[i + 1] - vx[i - 1] );
    p          = sig * y2[i - 1] + 2.0;
    double yy2 = ( sig - 1.0 ) / p;
    y2.push_back( yy2 );
    u[i] = ( vy[i + 1] - vy[i] ) / ( vx[i + 1] - vx[i] ) -
           ( vy[i] - vy[i - 1] ) / ( vx[i] - vx[i - 1] );
    u[i] = ( 6.0 * u[i] / ( vx[i + 1] - vx[i - 1] ) - sig * u[i - 1] ) / p;
  }
  if ( ypn > 0.99e30 ) qn = un = 0.0;
  else
  {
    qn = 0.5;
    un = ( 3.0 / ( vx[n - 1] - vx[n - 2] ) ) *
         ( ypn - ( vy[n - 1] - vy[n - 2] ) / ( vx[n - 1] - vx[n - 2] ) );
  }
  y2[n - 1] = ( un - qn * u[n - 2] ) / ( qn * y2[n - 2] + 1.0 );
  for ( k = n - 2; k >= 0; k-- ) y2[k] = y2[k] * y2[k + 1] + u[k];
  return y2;
}

Referenced by splint().

splint()

void EvtPolInt::splint

(

)

Definition at line 154 of file EvtPolInt.cc.

                       {
  vector<double> y2a = spline();
  double         y;
  int            k;
  double         h, b, a;
 
  int n   = vx.size();
  int klo = 0;
  int khi = n - 1;
  while ( khi - klo > 1 )
  {
    k = ( khi + klo ) >> 1;
    if ( vx[k] > xx ) khi = k;
    else klo = k;
  }
  h = vx[khi] - vx[klo];
  if ( h == 0.0 ) std::cout << "Bad xa input to routine splint" << std::endl;
  a = ( vx[khi] - xx ) / h;
  b = ( xx - vx[klo] ) / h;
  y = a * vy[klo] + b * vy[khi] +
      ( ( a * a * a - a ) * y2a[klo] + ( b * b * b - b ) * y2a[khi] ) * ( h * h ) / 6.0;
  value = y;
  return;
}

---

The documentation for this class was generated from the following files:

• EvtPolInt.hh
• EvtPolInt.cc