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>
7namespace Genfun {
8 FUNCTION_OBJECT_IMP(InterpolatingPolynomial)
9
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
#define FUNCTION_OBJECT_IMP(classname)
void getRange(double &min, double &max) const
virtual double operator()(double argument) const