CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

InterpolatingPolynomial.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // $Id:
4 #include <cassert>
5 #include <cmath>
6 #include <cfloat>
7 namespace Genfun {
8  FUNCTION_OBJECT_IMP(InterpolatingPolynomial)
9 
11  :Genfun::AbsFunction()
12  {}
13 
15  :Genfun::AbsFunction(),xPoints(right.xPoints)
16  {}
17 
19  }
20 
21  double InterpolatingPolynomial::operator() (double x) const {
22  double y=0.0;
23  double deltay=0; // also gives error;
24  double dif = fabs(x-xPoints[0].first),dift;
25  const unsigned int _K=xPoints.size(),_KP=_K+1;
26  std::vector<double>c(_KP),d(_KP);
27  int ns=0;
28  for (unsigned int i=0;i<_K;i++) {
29  dift=fabs(x-xPoints[i].first);
30  if (dift<dif) {
31  ns=i;
32  dif=dift;
33  }
34  c[i]=d[i]=xPoints[i].second;
35  }
36  y = xPoints[ns--].second;
37  for (unsigned int m=0;m<_K-1;m++) {
38  for (unsigned int i=0;i<_K-m-1;i++) {
39  double ho = xPoints[i].first-x;
40  double hp= xPoints[i+m+1].first-x;
41  double w=c[i+1]-d[i];
42  double den=ho-hp;
43  if (den==0)
44  std::cerr
45  << "Error in polynomial extrapolation"
46  << std::endl;
47  den=w/den;
48  d[i]=hp*den;
49  c[i]=ho*den;
50  }
51  deltay = 2*(ns+1) < (int)(_K-m-1) ? c[ns+1]: d[ns--];
52  y += deltay;
53  }
54  return y;
55  }
56 
57  void InterpolatingPolynomial::addPoint( double x, double y) {
58  xPoints.push_back(std::make_pair(x,y));
59  }
60 
61  void InterpolatingPolynomial::getRange( double & min, double & max) const {
62  min=DBL_MAX, max=-DBL_MAX;
63  for (unsigned int i=0;i<xPoints.size();i++) {
64  min = std::min(min,xPoints[i].first);
65  max = std::max(max,xPoints[i].first);
66  }
67  }
68 } // namespace Genfun
void getRange(double &min, double &max) const
#define FUNCTION_OBJECT_IMP(classname)
virtual double operator()(double argument) const