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 31 of file EvtPolInt.hh.

Constructor & Destructor Documentation

EvtPolInt()

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

inline

Definition at line 33 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 172 of file EvtPolInt.cc.

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

getvalue()

double EvtPolInt::getvalue

(

)

Definition at line 165 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 59 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 102 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 141 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