TRadSpline3 Class Reference

#include <TRadSpline.h>

Inheritance diagram for TRadSpline3:

List of all members.

Public Member Functions
	TRadSpline3 ()
	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 ()
void	GetCoeff (int i, double &x, double &y, double &b, double &c, double &d)
void	GetKnot (int i, double &x, double &y) const
virtual void	SaveAs (const char *filename) const
virtual int	GetNpx () const
void	SetNpx (int n)
Static Public Member Functions
static void	Test ()
Protected Attributes
double	fDelta
double	fXmin
double	fXmax
int	fNp
bool	fKstep
int	fNpx
Private Member Functions
void	BuildCoeff ()
void	SetCond (const char *opt)
Private Attributes
TRadSplinePoly3 *	fPoly
double	fValBeg
double	fValEnd
int	fBegCond
int	fEndCond

---

Detailed Description

Definition at line 113 of file TRadSpline.h.

---

Constructor & Destructor Documentation

TRadSpline3::TRadSpline3

(

)

[inline]

Definition at line 125 of file TRadSpline.h.

Referenced by Test().

                : fPoly(0), fValBeg(0), fValEnd(0),
    fBegCond(-1), fEndCond(-1) {}

TRadSpline3::TRadSpline3	(	const char *	title,
		double	x[],
		double	y[],
		int	n,
		const char *	opt = 0,
		double	valbeg = 0,
		double	valend = 0
	)

Definition at line 24 of file TRadSpline.C.

References BuildCoeff(), fPoly, genRecEmupikp::i, SetCond(), and TRadSplinePoly::Y().

                                                 :
  TRadSpline(title,-1,x[0],x[n-1],n,false),
  fValBeg(valbeg), fValEnd(valend), fBegCond(0), fEndCond(0)
{
  //
  // Third spline creator given an array of
  // arbitrary knots in increasing abscissa order and
  // possibly end point conditions
  //

  //
  // Set endpoint conditions
  if(opt) SetCond(opt);
  //
  // Create the plynomial terms and fill
  // them with node information
  fPoly = new TRadSplinePoly3[n];
  for (int i=0; i<n; ++i) {
    fPoly[i].X() = x[i];
    fPoly[i].Y() = y[i];
  }
  //
  // Build the spline coefficients
  BuildCoeff();
}

TRadSpline3::TRadSpline3	(	const char *	title,
		double	xmin,
		double	xmax,
		double	y[],
		int	n,
		const char *	opt = 0,
		double	valbeg = 0,
		double	valend = 0
	)

Definition at line 53 of file TRadSpline.C.

References BuildCoeff(), TRadSpline::fDelta, fPoly, TRadSpline::fXmin, genRecEmupikp::i, SetCond(), and TRadSplinePoly::Y().

                                                 :
  TRadSpline(title,(xmax-xmin)/(n-1), xmin, xmax, n, true),
  fValBeg(valbeg), fValEnd(valend),
  fBegCond(0), fEndCond(0)
{
  //
  // Third spline creator given an array of
  // arbitrary function values on equidistant n abscissa
  // values from xmin to xmax and possibly end point conditions
  //

  //
  // Set endpoint conditions
  if(opt) SetCond(opt);
  //
  // Create the plynomial terms and fill
  // them with node information
  fPoly = new TRadSplinePoly3[n];
  for (int i=0; i<n; ++i) {
    fPoly[i].X() = fXmin+i*fDelta;
    fPoly[i].Y() = y[i];
  }
  //
  // Build the spline coefficients
  BuildCoeff();
}

TRadSpline3::TRadSpline3	(	const char *	title,
		double	xmin,
		double	xmax,
		double(*)(const double &)	func,
		int	n,
		const char *	opt = 0,
		double	valbeg = 0,
		double	valend = 0
	)

Definition at line 83 of file TRadSpline.C.

References BuildCoeff(), TRadSpline::fDelta, fPoly, TRadSpline::fXmin, genRecEmupikp::i, SetCond(), TRadSplinePoly::X(), x, and TRadSplinePoly::Y().

                                                       :
  TRadSpline(title,(xmax-xmin)/(n-1), xmin, xmax, n, true),
  fValBeg(valbeg), fValEnd(valend),
  fBegCond(0), fEndCond(0)
{
  //
  // Third spline creator given a function to be
  // evaluated on n equidistand abscissa points between xmin
  // and xmax and possibly end point conditions
  //

  //
  // Set endpoint conditions
  if(opt) SetCond(opt);
  // Create the plynomial terms and fill
  // them with node information
  fPoly = new TRadSplinePoly3[n];
  for (int i=0; i<n; ++i) {
    double x=fXmin+i*fDelta;
    fPoly[i].X() = x;
    fPoly[i].Y() = func(x);
  }
  // Build the spline coefficients
  BuildCoeff();
}

virtual TRadSpline3::~TRadSpline3

(

)

[inline, virtual]

Definition at line 141 of file TRadSpline.h.

References fPoly.

{if (fPoly) delete [] fPoly;}

---

Member Function Documentation

void TRadSpline3::BuildCoeff

(

)

[private, virtual]

Implements TRadSpline.

Definition at line 492 of file TRadSpline.C.

References TRadSplinePoly3::B(), TRadSplinePoly3::C(), EvtCyclic3::C, TRadSplinePoly3::D(), fBegCond, fEndCond, TRadSpline::fNp, fPoly, fValBeg, fValEnd, genRecEmupikp::i, ganga-rec::j, and TRadSplinePoly::Y().

Referenced by TRadSpline3().

{
//      subroutine cubspl ( tau, c, n, ibcbeg, ibcend )
//  from  * a practical guide to splines *  by c. de boor
//     ************************  input  ***************************
//     n = number of data points. assumed to be .ge. 2.
//     (tau(i), c(1,i), i=1,...,n) = abscissae and ordinates of the
//        data points. tau is assumed to be strictly increasing.
//     ibcbeg, ibcend = boundary condition indicators, and
//     c(2,1), c(2,n) = boundary condition information. specifically,
//        ibcbeg = 0  means no boundary condition at tau(1) is given.
//           in this case, the not-a-knot condition is used, i.e. the
//           jump in the third derivative across tau(2) is forced to
//           zero, thus the first and the second cubic polynomial pieces
//           are made to coincide.)
//        ibcbeg = 1  means that the slope at tau(1) is made to equal
//           c(2,1), supplied by input.
//        ibcbeg = 2  means that the second derivative at tau(1) is
//           made to equal c(2,1), supplied by input.
//        ibcend = 0, 1, or 2 has analogous meaning concerning the
//           boundary condition at tau(n), with the additional infor-
//           mation taken from c(2,n).
//     ***********************  output  **************************
//     c(j,i), j=1,...,4; i=1,...,l (= n-1) = the polynomial coefficients
//        of the cubic interpolating spline with interior knots (or
//        joints) tau(2), ..., tau(n-1). precisely, in the interval
//        (tau(i), tau(i+1)), the spline f is given by
//           f(x) = c(1,i)+h*(c(2,i)+h*(c(3,i)+h*c(4,i)/3.)/2.)
//        where h = x - tau(i). the function program *ppvalu* may be
//        used to evaluate f or its derivatives from tau,c, l = n-1,
//        and k=4.

  int i, j, l, m;
  double   divdf1,divdf3,dtau,g=0;
//***** a tridiagonal linear system for the unknown slopes s(i) of
//  f  at tau(i), i=1,...,n, is generated and then solved by gauss elim-
//  ination, with s(i) ending up in c(2,i), all i.
//     c(3,.) and c(4,.) are used initially for temporary storage.
  l = fNp-1;
// compute first differences of x sequence and store in C also,
// compute first divided difference of data and store in D.
  for (m=1; m<fNp ; ++m) {
    fPoly[m].C() = fPoly[m].X() - fPoly[m-1].X();
    fPoly[m].D() = (fPoly[m].Y() - fPoly[m-1].Y())/fPoly[m].C();
  }
// construct first equation from the boundary condition, of the form
//             D[0]*s[0] + C[0]*s[1] = B[0]
  if(fBegCond==0) {
    if(fNp == 2) {
//     no condition at left end and n = 2.
      fPoly[0].D() = 1.;
      fPoly[0].C() = 1.;
      fPoly[0].B() = 2.*fPoly[1].D();
    } else {
//     not-a-knot condition at left end and n .gt. 2.
      fPoly[0].D() = fPoly[2].C();
      fPoly[0].C() = fPoly[1].C() + fPoly[2].C();
      fPoly[0].B() =((fPoly[1].C()+2.*fPoly[0].C())*fPoly[1].D()*fPoly[2].C()+fPoly[1].C()*fPoly[1].C()*fPoly[2].D())/fPoly[0].C();
    }
  } else if (fBegCond==1) {
//     slope prescribed at left end.
    fPoly[0].B() = fValBeg;
    fPoly[0].D() = 1.;
    fPoly[0].C() = 0.;
  } else if (fBegCond==2) {
//     second derivative prescribed at left end.
    fPoly[0].D() = 2.;
    fPoly[0].C() = 1.;
    fPoly[0].B() = 3.*fPoly[1].D() - fPoly[1].C()/2.*fValBeg;
  }
  if(fNp > 2) {
//  if there are interior knots, generate the corresp. equations and car-
//  ry out the forward pass of gauss elimination, after which the m-th
//  equation reads    D[m]*s[m] + C[m]*s[m+1] = B[m].
    for (m=1; m<l; ++m) {
      g = -fPoly[m+1].C()/fPoly[m-1].D();
      fPoly[m].B() = g*fPoly[m-1].B() + 3.*(fPoly[m].C()*fPoly[m+1].D()+fPoly[m+1].C()*fPoly[m].D());
      fPoly[m].D() = g*fPoly[m-1].C() + 2.*(fPoly[m].C() + fPoly[m+1].C());
    }
// construct last equation from the second boundary condition, of the form
//           (-g*D[n-2])*s[n-2] + D[n-1]*s[n-1] = B[n-1]
//     if slope is prescribed at right end, one can go directly to back-
//     substitution, since c array happens to be set up just right for it
//     at this point.
    if(fEndCond == 0) {
      if (fNp > 3 || fBegCond != 0) {
//     not-a-knot and n .ge. 3, and either n.gt.3 or  also not-a-knot at
//     left end point.
        g = fPoly[fNp-2].C() + fPoly[fNp-1].C();
        fPoly[fNp-1].B() = ((fPoly[fNp-1].C()+2.*g)*fPoly[fNp-1].D()*fPoly[fNp-2].C()
                     + fPoly[fNp-1].C()*fPoly[fNp-1].C()*(fPoly[fNp-2].Y()-fPoly[fNp-3].Y())/fPoly[fNp-2].C())/g;
        g = -g/fPoly[fNp-2].D();
        fPoly[fNp-1].D() = fPoly[fNp-2].C();
      } else {
//     either (n=3 and not-a-knot also at left) or (n=2 and not not-a-
//     knot at left end point).
        fPoly[fNp-1].B() = 2.*fPoly[fNp-1].D();
        fPoly[fNp-1].D() = 1.;
        g = -1./fPoly[fNp-2].D();
      }
    } else if (fEndCond == 1) {
      fPoly[fNp-1].B() = fValEnd;
      goto L30;
    } else if (fEndCond == 2) {
//     second derivative prescribed at right endpoint.
      fPoly[fNp-1].B() = 3.*fPoly[fNp-1].D() + fPoly[fNp-1].C()/2.*fValEnd;
      fPoly[fNp-1].D() = 2.;
      g = -1./fPoly[fNp-2].D();
    }
   } else {
     if(fEndCond == 0) {
       if (fBegCond > 0) {
//     either (n=3 and not-a-knot also at left) or (n=2 and not not-a-
//     knot at left end point).
         fPoly[fNp-1].B() = 2.*fPoly[fNp-1].D();
         fPoly[fNp-1].D() = 1.;
         g = -1./fPoly[fNp-2].D();
       } else {
//     not-a-knot at right endpoint and at left endpoint and n = 2.
         fPoly[fNp-1].B() = fPoly[fNp-1].D();
         goto L30;
       }
     } else if(fEndCond == 1) {
       fPoly[fNp-1].B() = fValEnd;
       goto L30;
     } else if(fEndCond == 2) {
//     second derivative prescribed at right endpoint.
       fPoly[fNp-1].B() = 3.*fPoly[fNp-1].D() + fPoly[fNp-1].C()/2.*fValEnd;
       fPoly[fNp-1].D() = 2.;
       g = -1./fPoly[fNp-2].D();
     }
   }
// complete forward pass of gauss elimination.
  fPoly[fNp-1].D() = g*fPoly[fNp-2].C() + fPoly[fNp-1].D();
  fPoly[fNp-1].B() = (g*fPoly[fNp-2].B() + fPoly[fNp-1].B())/fPoly[fNp-1].D();
// carry out back substitution
 L30: j = l-1;
  do {
    fPoly[j].B() = (fPoly[j].B() - fPoly[j].C()*fPoly[j+1].B())/fPoly[j].D();
    --j;
  }  while (j>=0);
//****** generate cubic coefficients in each interval, i.e., the deriv.s
//  at its left endpoint, from value and slope at its endpoints.
  for (i=1; i<fNp; ++i) {
    dtau = fPoly[i].C();
    divdf1 = (fPoly[i].Y() - fPoly[i-1].Y())/dtau;
    divdf3 = fPoly[i-1].B() + fPoly[i].B() - 2.*divdf1;
    fPoly[i-1].C() = (divdf1 - fPoly[i-1].B() - divdf3)/dtau;
    fPoly[i-1].D() = (divdf3/dtau)/dtau;
  }
}

double TRadSpline3::Derivative

(

double

x

)

const

Definition at line 325 of file TRadSpline.C.

References TRadSplinePoly3::Derivative(), FindX(), and fPoly.

{
  int klow=FindX(x);
  return fPoly[klow].Derivative(x);
}

double TRadSpline3::Eval

(

double

x

)

const [virtual]

Implements TRadSpline.

Definition at line 318 of file TRadSpline.C.

References TRadSplinePoly3::Eval(), FindX(), and fPoly.

Referenced by TDFun::EvalSpline().

{
  int klow=FindX(x);
  return fPoly[klow].Eval(x);
}

int TRadSpline3::FindX

(

double

x

)

const

Definition at line 285 of file TRadSpline.C.

References TRadSpline::fDelta, TRadSpline::fKstep, TRadSpline::fNp, fPoly, TRadSpline::fXmax, and TRadSpline::fXmin.

Referenced by Derivative(), and Eval().

{
  int klow=0;
  //
  // If out of boundaries, extrapolate
  // It may be badly wrong
  if(x<=fXmin) klow=0;
  else if(x>=fXmax) klow=fNp-1;
  else {
    if(fKstep) {
      //
      // Equidistant knots, use histogramming
      klow = (int((x-fXmin)/fDelta)<fNp-1)?int((x-fXmin)/fDelta):fNp-1;
    } else {
      int khig=fNp-1, khalf;
      //
      // Non equidistant knots, binary search
      while(khig-klow>1)
        if(x>fPoly[khalf=(klow+khig)/2].X())
          klow=khalf;
        else
          khig=khalf;
    }
    //
    // This could be removed, sanity check
    if(!(fPoly[klow].X()<=x && x<=fPoly[klow+1].X()))
      printf("Binary search failed x(%d) = %f < %f < x(%d) = %f\n",
             klow, fPoly[klow].X(), x, klow+1, fPoly[klow+1].X());
  }
  return klow;
}

void TRadSpline3::GetCoeff	(	int	i,
		double &	x,
		double &	y,
		double &	b,
		double &	c,
		double &	d
	)			[inline]

Definition at line 142 of file TRadSpline.h.

References TRadSplinePoly3::B(), TRadSplinePoly3::C(), TRadSplinePoly3::D(), fPoly, TRadSplinePoly::X(), and TRadSplinePoly::Y().

Referenced by Test().

                                      {x=fPoly[i].X();y=fPoly[i].Y();
                b=fPoly[i].B();c=fPoly[i].C();d=fPoly[i].D();}

void TRadSpline3::GetKnot	(	int	i,
		double &	x,
		double &	y
	)			const [inline, virtual]

Implements TRadSpline.

Definition at line 145 of file TRadSpline.h.

References fPoly, TRadSplinePoly::X(), and TRadSplinePoly::Y().

    {x=fPoly[i].X(); y=fPoly[i].Y();}

virtual int TRadSpline::GetNpx

(

)

const [inline, virtual, inherited]

Definition at line 26 of file TRadSpline.h.

References TRadSpline::fNpx.

{return fNpx;}

void TRadSpline3::SaveAs

(

const char *

filename

)

const [virtual]

Reimplemented from TRadSpline.

Definition at line 332 of file TRadSpline.C.

References EvtCyclic3::B, EvtCyclic3::C, TRadSpline::fDelta, TRadSpline::fKstep, TRadSpline::fNp, fPoly, TRadSpline::fXmax, TRadSpline::fXmin, and genRecEmupikp::i.

{
  // write this spline as a C++ function that can be executed without ROOT
  // the name of the function is the name of the file up to the "." if any
    
   //open the file
   ofstream *f = new ofstream(filename,ios::out);
   if (f == 0) {
      printf("Cannot open file:%s\n",filename);
      return;
   }
   
   //write the function name and the spline constants
   char buffer[512];
   int nch = strlen(filename);
   sprintf(buffer,"double %s",filename);
   char *dot = strstr(buffer,".");
   if (dot) *dot = 0;
   strcat(buffer,"(double x) {\n");
   nch = strlen(buffer); f->write(buffer,nch);
   sprintf(buffer,"   const int fNp = %d, fKstep = %d;\n",fNp,fKstep);
   nch = strlen(buffer); f->write(buffer,nch);
   sprintf(buffer,"   const double fDelta = %g, fXmin = %g, fXmax = %g;\n",fDelta,fXmin,fXmax);
   nch = strlen(buffer); f->write(buffer,nch);

   //write the spline coefficients
   //array fX
   sprintf(buffer,"   const double fX[%d] = {",fNp);   
   nch = strlen(buffer); f->write(buffer,nch);
   buffer[0] = 0;
   int i;
   char numb[20];
   for (i=0;i<fNp;i++) {
      sprintf(numb," %g,",fPoly[i].X());
      nch = strlen(numb);
      if (i == fNp-1) numb[nch-1]=0;
      strcat(buffer,numb);
      if (i%5 == 4 || i == fNp-1) {
         nch = strlen(buffer); f->write(buffer,nch);
         if (i != fNp-1) sprintf(buffer,"\n                       ");
      }
   }
   sprintf(buffer," };\n");
   nch = strlen(buffer); f->write(buffer,nch);
   //array fY
   sprintf(buffer,"   const double fY[%d] = {",fNp);   
   nch = strlen(buffer); f->write(buffer,nch);
   buffer[0] = 0;
   for (i=0;i<fNp;i++) {
      sprintf(numb," %g,",fPoly[i].Y());
      nch = strlen(numb);
      if (i == fNp-1) numb[nch-1]=0;
      strcat(buffer,numb);
      if (i%5 == 4 || i == fNp-1) {
         nch = strlen(buffer); f->write(buffer,nch);
         if (i != fNp-1) sprintf(buffer,"\n                       ");
      }
   }
   sprintf(buffer," };\n");
   nch = strlen(buffer); f->write(buffer,nch);
   //array fB
   sprintf(buffer,"   const double fB[%d] = {",fNp);   
   nch = strlen(buffer); f->write(buffer,nch);
   buffer[0] = 0;
   for (i=0;i<fNp;i++) {
      sprintf(numb," %g,",fPoly[i].B());
      nch = strlen(numb);
      if (i == fNp-1) numb[nch-1]=0;
      strcat(buffer,numb);
      if (i%5 == 4 || i == fNp-1) {
         nch = strlen(buffer); f->write(buffer,nch);
         if (i != fNp-1) sprintf(buffer,"\n                       ");
      }
   }
   sprintf(buffer," };\n");
   nch = strlen(buffer); f->write(buffer,nch);
   //array fC
   sprintf(buffer,"   const double fC[%d] = {",fNp);   
   nch = strlen(buffer); f->write(buffer,nch);
   buffer[0] = 0;
   for (i=0;i<fNp;i++) {
      sprintf(numb," %g,",fPoly[i].C());
      nch = strlen(numb);
      if (i == fNp-1) numb[nch-1]=0;
      strcat(buffer,numb);
      if (i%5 == 4 || i == fNp-1) {
         nch = strlen(buffer); f->write(buffer,nch);
         if (i != fNp-1) sprintf(buffer,"\n                       ");
      }
   }
   sprintf(buffer," };\n");
   nch = strlen(buffer); f->write(buffer,nch);
    //array fD
   sprintf(buffer,"   const double fD[%d] = {",fNp);   
   nch = strlen(buffer); f->write(buffer,nch);
   buffer[0] = 0;
   for (i=0;i<fNp;i++) {
      sprintf(numb," %g,",fPoly[i].D());
      nch = strlen(numb);
      if (i == fNp-1) numb[nch-1]=0;
      strcat(buffer,numb);
      if (i%5 == 4 || i == fNp-1) {
         nch = strlen(buffer); f->write(buffer,nch);
         if (i != fNp-1) sprintf(buffer,"\n                       ");
      }
   }
   sprintf(buffer," };\n");
   nch = strlen(buffer); f->write(buffer,nch);
    
   //generate code for the spline evaluation
   sprintf(buffer,"   int klow=0;\n");
   nch = strlen(buffer); f->write(buffer,nch);

   sprintf(buffer,"   // If out of boundaries, extrapolate. It may be badly wrong\n");
   sprintf(buffer,"   if(x<=fXmin) klow=0;\n");
   nch = strlen(buffer); f->write(buffer,nch);
   sprintf(buffer,"   else if(x>=fXmax) klow=fNp-1;\n");
   nch = strlen(buffer); f->write(buffer,nch);
   sprintf(buffer,"   else {\n");
   nch = strlen(buffer); f->write(buffer,nch);
   sprintf(buffer,"     if(fKstep) {\n");
   nch = strlen(buffer); f->write(buffer,nch);

   sprintf(buffer,"       // Equidistant knots, use histogramming\n");
   nch = strlen(buffer); f->write(buffer,nch);
   sprintf(buffer,"       klow = int((x-fXmin)/fDelta);\n");
   nch = strlen(buffer); f->write(buffer,nch);
   sprintf(buffer,"       if (klow < fNp-1) klow = fNp-1;\n");
   nch = strlen(buffer); f->write(buffer,nch);
   sprintf(buffer,"     } else {\n");
   nch = strlen(buffer); f->write(buffer,nch);
   sprintf(buffer,"       int khig=fNp-1, khalf;\n");
   nch = strlen(buffer); f->write(buffer,nch);

   sprintf(buffer,"       // Non equidistant knots, binary search\n");
   nch = strlen(buffer); f->write(buffer,nch);
   sprintf(buffer,"       while(khig-klow>1)\n");
   nch = strlen(buffer); f->write(buffer,nch);
   sprintf(buffer,"      if(x>fX[khalf=(klow+khig)/2]) klow=khalf;\n");
   nch = strlen(buffer); f->write(buffer,nch);
   sprintf(buffer,"      else khig=khalf;\n");
   nch = strlen(buffer); f->write(buffer,nch);
   sprintf(buffer,"     }\n");
   nch = strlen(buffer); f->write(buffer,nch);
   sprintf(buffer,"   }\n");
   nch = strlen(buffer); f->write(buffer,nch);
   sprintf(buffer,"   // Evaluate now\n");
   nch = strlen(buffer); f->write(buffer,nch);
   sprintf(buffer,"   double dx=x-fX[klow];\n");
   nch = strlen(buffer); f->write(buffer,nch);
   sprintf(buffer,"   return (fY[klow]+dx*(fB[klow]+dx*(fC[klow]+dx*fD[klow])));\n");
   nch = strlen(buffer); f->write(buffer,nch);
    
   //close file
   f->write("}\n",2);
      
   if (f) { f->close(); delete f;}
}

void TRadSpline3::SetCond

(

const char *

opt

)

[private]

Definition at line 112 of file TRadSpline.C.

References fBegCond, and fEndCond.

Referenced by TRadSpline3().

{
  //
  // Check the boundary conditions
  //
  const char *b1 = strstr(opt,"b1");
  const char *e1 = strstr(opt,"e1");
  const char *b2 = strstr(opt,"b2");
  const char *e2 = strstr(opt,"e2");
  if (b1 && b2)
    printf("%s %s\n","SetCond","Cannot specify first and second derivative at first point");
  if (e1 && e2)
    printf("%s %s\n","SetCond","Cannot specify first and second derivative at last point");
  if (b1) fBegCond=1;
  else if (b2) fBegCond=2;
  if (e1) fEndCond=1;
  else if (e2) fEndCond=2;
}

void TRadSpline::SetNpx

(

int

n

)

[inline, inherited]

Definition at line 30 of file TRadSpline.h.

References TRadSpline::fNpx.

{fNpx=n;}

void TRadSpline3::Test

(

)

[static]

Definition at line 132 of file TRadSpline.C.

References GetCoeff(), genRecEmupikp::i, ganga-rec::j, TRadSpline3(), and x.

{
  //
  // Test method for TRadSpline5
  //
  //   n          number of data points.
  //   m          2*m-1 is order of spline.
  //                 m = 2 always for third spline.
  //   nn,nm1,mm,
  //   mm1,i,k,
  //   j,jj       temporary integer variables.
  //   z,p        temporary double precision variables.
  //   x[n]       the sequence of knots.
  //   y[n]       the prescribed function values at the knots.
  //   a[200][4]  two dimensional array whose columns are
  //                 the computed spline coefficients
  //   diff[3]    maximum values of differences of values and
  //                 derivatives to right and left of knots.
  //   com[3]     maximum values of coefficients.
  //
  //
  //   test of TRadSpline3 with nonequidistant knots and
  //      equidistant knots follows.
  //
  //
  double hx;
  double diff[3];
  double a[800], c[4];
  int i, j, k, m, n;
  double x[200], y[200], z;
  int jj, mm;
  int mm1, nm1;
  double com[3];
  printf("1         TEST OF TRadSpline3 WITH NONEQUIDISTANT KNOTS\n");
  n = 5;
  x[0] = -3;
  x[1] = -1;
  x[2] = 0;
  x[3] = 3;
  x[4] = 4;
  y[0] = 7;
  y[1] = 11;
  y[2] = 26;
  y[3] = 56;
  y[4] = 29;
  m = 2;
  mm = m << 1;
  mm1 = mm-1;
  printf("\n-N = %3d    M =%2d\n",n,m);
  TRadSpline3 *spline = new TRadSpline3("Test",x,y,n);
  for (i = 0; i < n; ++i)
    spline->GetCoeff(i,hx, a[i],a[i+200],a[i+400],a[i+600]);
  delete spline;
  for (i = 0; i < mm1; ++i) diff[i] = com[i] = 0;
  for (k = 0; k < n; ++k) {
    for (i = 0; i < mm; ++i) c[i] = a[k+i*200];
    printf(" ---------------------------------------%3d --------------------------------------------\n",k+1);
    printf("%12.8f\n",x[k]);
    if (k == n-1) {
      printf("%16.8f\n",c[0]);
    } else {
      for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
      printf("\n");
      for (i = 0; i < mm1; ++i)
        if ((z=fabs(a[k+i*200])) > com[i]) com[i] = z;
      z = x[k+1]-x[k];
      for (i = 1; i < mm; ++i)
        for (jj = i; jj < mm; ++jj) {
          j = mm+i-jj;
          c[j-2] = c[j-1]*z+c[j-2];
        }
      for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
      printf("\n");
      for (i = 0; i < mm1; ++i)
        if (!(k >= n-2 && i != 0))
          if((z = fabs(c[i]-a[k+1+i*200]))
             > diff[i]) diff[i] = z;
    }
  }
  printf("  MAXIMUM ABSOLUTE VALUES OF DIFFERENCES \n");
  for (i = 0; i < mm1; ++i) printf("%18.9E",diff[i]);
  printf("\n");
  printf("  MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS \n");
  if (fabs(c[0]) > com[0])
    com[0] = fabs(c[0]);
  for (i = 0; i < mm1; ++i) printf("%16.8f",com[i]);
  printf("\n");
  m = 2;
  for (n = 10; n <= 100; n += 10) {
    mm = m << 1;
    mm1 = mm-1;
    nm1 = n-1;
    for (i = 0; i < nm1; i += 2) {
      x[i] = i+1;
      x[i+1] = i+2;
      y[i] = 1;
      y[i+1] = 0;
    }
    if (n % 2 != 0) {
      x[n-1] = n;
      y[n-1] = 1;
    }
    printf("\n-N = %3d    M =%2d\n",n,m);
    spline = new TRadSpline3("Test",x,y,n);
    for (i = 0; i < n; ++i)
      spline->GetCoeff(i,hx,a[i],a[i+200],a[i+400],a[i+600]);
    delete spline;
    for (i = 0; i < mm1; ++i)
      diff[i] = com[i] = 0;
    for (k = 0; k < n; ++k) {
      for (i = 0; i < mm; ++i)
        c[i] = a[k+i*200];
      if (n < 11) {
        printf(" ---------------------------------------%3d --------------------------------------------\n",k+1);
        printf("%12.8f\n",x[k]);
        if (k == n-1) printf("%16.8f\n",c[0]);
      }
      if (k == n-1) break;
      if (n <= 10) {
        for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
        printf("\n");
      }
      for (i = 0; i < mm1; ++i)
        if ((z=fabs(a[k+i*200])) > com[i])
          com[i] = z;
      z = x[k+1]-x[k];
      for (i = 1; i < mm; ++i)
        for (jj = i; jj < mm; ++jj) {
          j = mm+i-jj;
          c[j-2] = c[j-1]*z+c[j-2];
        }
      if (n <= 10) {
        for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
        printf("\n");
      }
      for (i = 0; i < mm1; ++i)
        if (!(k >= n-2 && i != 0))
          if ((z = fabs(c[i]-a[k+1+i*200]))
              > diff[i]) diff[i] = z;
    }
    printf("  MAXIMUM ABSOLUTE VALUES OF DIFFERENCES \n");
    for (i = 0; i < mm1; ++i) printf("%18.9E",diff[i]);
    printf("\n");
    printf("  MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS \n");
    if (fabs(c[0]) > com[0])
      com[0] = fabs(c[0]);
    for (i = 0; i < mm1; ++i) printf("%16.8E",com[i]);
    printf("\n");
  }
}

---

Member Data Documentation

int TRadSpline3::fBegCond [private]

Definition at line 118 of file TRadSpline.h.

Referenced by BuildCoeff(), and SetCond().

double TRadSpline::fDelta [protected, inherited]

Definition at line 6 of file TRadSpline.h.

Referenced by TRadSpline5::FindX(), FindX(), TRadSpline5::SaveAs(), SaveAs(), TRadSpline3(), and TRadSpline5::TRadSpline5().

int TRadSpline3::fEndCond [private]

Definition at line 119 of file TRadSpline.h.

Referenced by BuildCoeff(), and SetCond().

bool TRadSpline::fKstep [protected, inherited]

Definition at line 10 of file TRadSpline.h.

Referenced by TRadSpline5::FindX(), FindX(), TRadSpline5::SaveAs(), and SaveAs().

int TRadSpline::fNp [protected, inherited]

Definition at line 9 of file TRadSpline.h.

Referenced by TRadSpline5::BoundaryConditions(), TRadSpline5::BuildCoeff(), BuildCoeff(), TRadSpline5::FindX(), FindX(), TRadSpline5::SaveAs(), SaveAs(), TRadSpline5::SetBoundaries(), and TRadSpline5::TRadSpline5().

int TRadSpline::fNpx [protected, inherited]

Definition at line 11 of file TRadSpline.h.

Referenced by TRadSpline::GetNpx(), and TRadSpline::SetNpx().

TRadSplinePoly3* TRadSpline3::fPoly [private]

Definition at line 115 of file TRadSpline.h.

Referenced by BuildCoeff(), Derivative(), Eval(), FindX(), GetCoeff(), GetKnot(), SaveAs(), TRadSpline3(), and ~TRadSpline3().

double TRadSpline3::fValBeg [private]

Definition at line 116 of file TRadSpline.h.

Referenced by BuildCoeff().

double TRadSpline3::fValEnd [private]

Definition at line 117 of file TRadSpline.h.

Referenced by BuildCoeff().

double TRadSpline::fXmax [protected, inherited]

Definition at line 8 of file TRadSpline.h.

Referenced by TRadSpline5::FindX(), FindX(), TRadSpline5::SaveAs(), and SaveAs().

double TRadSpline::fXmin [protected, inherited]

Definition at line 7 of file TRadSpline.h.

Referenced by TRadSpline5::FindX(), FindX(), TRadSpline5::SaveAs(), SaveAs(), TRadSpline3(), and TRadSpline5::TRadSpline5().

---