CLHEP 2.4.7.1
C++ Class Library for High Energy Physics
CubicSplinePolynomial.icc
Go to the documentation of this file.
1// -*- C++ -*-
2// $Id:
5#include <cassert>
6#include <cmath>
7#include <cfloat>
8namespace Genfun {
9 FUNCTION_OBJECT_IMP(CubicSplinePolynomial)
10
12 public:
13
14 bool stale;
18 std::vector<std::pair<double,double> > xPoints;
19 inline void computeSlopes() const {
20
21 unsigned int N=xPoints.size()-1;
22 A=CLHEP::HepMatrix(N+1,N+1,0);
23 V=CLHEP::HepVector(N+1,0);
24
25 // First take care of the "normal elements, i=1,N-1;
26 for (unsigned int i=1;i<N;i++) {
27
28 double dxPlus = xPoints[i+1].first -xPoints[i].first;
29 double dyPlus = xPoints[i+1].second-xPoints[i].second;
30 double mPlus = dyPlus/dxPlus;
31
32 double dx = xPoints[i].first -xPoints[i-1].first;
33 double dy = xPoints[i].second-xPoints[i-1].second;
34 double m = dy/dx;
35
36 A[i][i-1] = 1/dx;
37 A[i][i] = 2*(1/dxPlus + 1/dx);
38 A[i][i+1] = 1/dxPlus;
39
40 V[i] = 3*(m/dx+mPlus/dxPlus);
41
42 }
43 // Special treatment for i=0;
44 {
45 double dx = xPoints[1].first -xPoints[0].first;
46 double dy = xPoints[1].second-xPoints[0].second;
47 double m = dy/dx;
48 A[0][0] = 2.0;
49 A[0][1] = 1;
50 V[0] = 3*m;
51 }
52 // Special treatment for i=N-1;
53 {
54 double dx = xPoints[N].first -xPoints[N-1].first;
55 double dy = xPoints[N].second-xPoints[N-1].second;
56 double m = dy/dx;
57 A[N][N-1] = 1.0;
58 A[N][N] = 2.0;
59 V[N] = 3*m;
60 }
61 int err;
62 Q=A.inverse(err)*V;
63 }
64 };
65
67 :Genfun::AbsFunction(),c(new Clockwork())
68 {
69 c->stale=true;
70 }
71
73 :Genfun::AbsFunction(),c(new Clockwork)
74 {
75 (*c) = (*right.c);
76 }
77
81
82 inline double CubicSplinePolynomial::operator() (double x) const {
83 unsigned int N=c->xPoints.size()-1;
84
85 if (c->xPoints.size()==0) return 0;
86 if (c->xPoints.size()==1) return c->xPoints[0].second;
87 if (c->stale) {
88 c->computeSlopes();
89 c->stale=false;
90 }
91 std::pair<double,double> fk(x,0);
92 std::vector<std::pair<double,double> >::const_iterator
93 n=std::lower_bound(c->xPoints.begin(),c->xPoints.end(),fk);
94 unsigned int i = n-c->xPoints.begin();
95 if (i==0) {
96 double x0=c->xPoints[0].first;
97 double y0=c->xPoints[0].second;
98 double m = c->Q[0];
99 return y0 + m*(x-x0);
100 }
101 else if (i==c->xPoints.size()) {
102 double x0=c->xPoints[N].first;
103 double y0=c->xPoints[N].second;
104 double m = c->Q[N];
105 return y0 + m*(x-x0);
106 }
107 else {
108 double x0=c->xPoints[i-1].first;
109 double x1=c->xPoints[i].first;
110 double y0=c->xPoints[i-1].second;
111 double y1=c->xPoints[i].second;
112 double dx = x1-x0;
113 double dy = y1-y0;
114 double m = dy/dx;
115 double Q0 = c->Q[i-1];
116 double Q1 = c->Q[i];
117 double a = (Q0-m)*dx;
118 double b = (m-Q1)*dx;
119 double t = (x-x0)/dx;
120 double u = (1-t);
121 return u*y0+t*y1 + t*u*(u*a + t*b);
122 }
123 }
124
125 inline void CubicSplinePolynomial::addPoint( double x, double y) {
126 c->xPoints.push_back(std::make_pair(x,y));
127 std::sort(c->xPoints.begin(),c->xPoints.end());
128 c->stale=true;
129 }
130
131 inline void CubicSplinePolynomial::getRange( double & min, double & max) const {
132 min=DBL_MAX, max=-DBL_MAX;
133 for (unsigned int i=0;i<c->xPoints.size();i++) {
134 min = std::min(min,c->xPoints[i].first);
135 max = std::max(max,c->xPoints[i].first);
136 }
137 }
138} // namespace Genfun
#define FUNCTION_OBJECT_IMP(classname)
HepMatrix inverse(int &ierr) const
Definition Matrix.icc:78
std::vector< std::pair< double, double > > xPoints
virtual double operator()(double argument) const override
void getRange(double &min, double &max) const
Definition Abs.hh:14