Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4DataInterpolation Class Reference

#include <G4DataInterpolation.hh>

Public Member Functions

 G4DataInterpolation (G4double pX[], G4double pY[], G4int number)
 
 G4DataInterpolation (G4double pX[], G4double pY[], G4int number, G4double pFirstDerStart, G4double pFirstDerFinish)
 
 ~G4DataInterpolation ()
 
G4double PolynomInterpolation (G4double pX, G4double &deltaY) const
 
void PolIntCoefficient (G4double cof[]) const
 
G4double RationalPolInterpolation (G4double pX, G4double &deltaY) const
 
G4double CubicSplineInterpolation (G4double pX) const
 
G4double FastCubicSpline (G4double pX, G4int index) const
 
G4int LocateArgument (G4double pX) const
 
void CorrelatedSearch (G4double pX, G4int &index) const
 

Detailed Description

Definition at line 116 of file G4DataInterpolation.hh.

Constructor & Destructor Documentation

G4DataInterpolation::G4DataInterpolation ( G4double  pX[],
G4double  pY[],
G4int  number 
)

Definition at line 36 of file G4DataInterpolation.cc.

39  : fArgument(new G4double[number]),
40  fFunction(new G4double[number]),
41  fSecondDerivative(0),
42  fNumber(number)
43 {
44  for(G4int i=0;i<fNumber;i++)
45  {
46  fArgument[i] = pX[i] ;
47  fFunction[i] = pY[i] ;
48  }
49 }
int G4int
Definition: G4Types.hh:78
double G4double
Definition: G4Types.hh:76
G4DataInterpolation::G4DataInterpolation ( G4double  pX[],
G4double  pY[],
G4int  number,
G4double  pFirstDerStart,
G4double  pFirstDerFinish 
)

Definition at line 58 of file G4DataInterpolation.cc.

63  : fArgument(new G4double[number]),
64  fFunction(new G4double[number]),
65  fSecondDerivative(new G4double[number]),
66  fNumber(number)
67 {
68  G4int i=0 ;
69  G4double p=0.0, qn=0.0, sig=0.0, un=0.0 ;
70  const G4double maxDerivative = 0.99e30 ;
71  G4double* u = new G4double[fNumber - 1] ;
72 
73  for(i=0;i<fNumber;i++)
74  {
75  fArgument[i] = pX[i] ;
76  fFunction[i] = pY[i] ;
77  }
78  if(pFirstDerStart > maxDerivative)
79  {
80  fSecondDerivative[0] = 0.0 ;
81  u[0] = 0.0 ;
82  }
83  else
84  {
85  fSecondDerivative[0] = -0.5 ;
86  u[0] = (3.0/(fArgument[1]-fArgument[0]))
87  * ((fFunction[1]-fFunction[0])/(fArgument[1]-fArgument[0])
88  - pFirstDerStart) ;
89  }
90 
91  // Decomposition loop for tridiagonal algorithm. fSecondDerivative[i]
92  // and u[i] are used for temporary storage of the decomposed factors.
93 
94  for(i=1;i<fNumber-1;i++)
95  {
96  sig = (fArgument[i]-fArgument[i-1])/(fArgument[i+1]-fArgument[i-1]) ;
97  p = sig*fSecondDerivative[i-1] + 2.0 ;
98  fSecondDerivative[i] = (sig - 1.0)/p ;
99  u[i] = (fFunction[i+1]-fFunction[i])/(fArgument[i+1]-fArgument[i]) -
100  (fFunction[i]-fFunction[i-1])/(fArgument[i]-fArgument[i-1]) ;
101  u[i] =(6.0*u[i]/(fArgument[i+1]-fArgument[i-1]) - sig*u[i-1])/p ;
102  }
103  if(pFirstDerFinish > maxDerivative)
104  {
105  qn = 0.0 ;
106  un = 0.0 ;
107  }
108  else
109  {
110  qn = 0.5 ;
111  un = (3.0/(fArgument[fNumber-1]-fArgument[fNumber-2]))
112  * (pFirstDerFinish - (fFunction[fNumber-1]-fFunction[fNumber-2])
113  / (fArgument[fNumber-1]-fArgument[fNumber-2])) ;
114  }
115  fSecondDerivative[fNumber-1] = (un - qn*u[fNumber-2])/
116  (qn*fSecondDerivative[fNumber-2] + 1.0) ;
117 
118  // The backsubstitution loop for the triagonal algorithm of solving
119  // a linear system of equations.
120 
121  for(G4int k=fNumber-2;k>=0;k--)
122  {
123  fSecondDerivative[k] = fSecondDerivative[k]*fSecondDerivative[k+1] + u[k];
124  }
125  delete[] u ;
126 }
const char * p
Definition: xmltok.h:285
int G4int
Definition: G4Types.hh:78
double G4double
Definition: G4Types.hh:76
G4DataInterpolation::~G4DataInterpolation ( )

Definition at line 133 of file G4DataInterpolation.cc.

134 {
135  delete [] fArgument ;
136  delete [] fFunction ;
137  if(fSecondDerivative) { delete [] fSecondDerivative; }
138 }

Member Function Documentation

void G4DataInterpolation::CorrelatedSearch ( G4double  pX,
G4int index 
) const

Definition at line 414 of file G4DataInterpolation.cc.

416 {
417  G4int kHigh=0, k=0, Increment=0 ;
418  // ascend = true for ascending order of table, false otherwise
419  G4bool ascend = (fArgument[fNumber-1] >= fArgument[0]) ;
420  if(index < 0 || index > fNumber-1)
421  {
422  index = -1 ;
423  kHigh = fNumber ;
424  }
425  else
426  {
427  Increment = 1 ; // What value would be the best ?
428  if((pX >= fArgument[index]) == ascend)
429  {
430  if(index == fNumber -1)
431  {
432  index = fNumber ;
433  return ;
434  }
435  kHigh = index + 1 ;
436  while((pX >= fArgument[kHigh]) == ascend)
437  {
438  index = kHigh ;
439  Increment += Increment ; // double the Increment
440  kHigh = index + Increment ;
441  if(kHigh > (fNumber - 1))
442  {
443  kHigh = fNumber ;
444  break ;
445  }
446  }
447  }
448  else
449  {
450  if(index == 0)
451  {
452  index = -1 ;
453  return ;
454  }
455  kHigh = index-- ;
456  while((pX < fArgument[index]) == ascend)
457  {
458  kHigh = index ;
459  Increment <<= 1 ; // double the Increment
460  if(Increment >= kHigh)
461  {
462  index = -1 ;
463  break ;
464  }
465  else
466  {
467  index = kHigh - Increment ;
468  }
469  }
470  } // Value bracketed
471  }
472  // final bisection searching
473 
474  while((kHigh - index) != 1)
475  {
476  k = (kHigh + index) >> 1 ;
477  if((pX >= fArgument[k]) == ascend)
478  {
479  index = k ;
480  }
481  else
482  {
483  kHigh = k ;
484  }
485  }
486  if (!(pX != fArgument[fNumber-1]))
487  {
488  index = fNumber - 2 ;
489  }
490  if (!(pX != fArgument[0]))
491  {
492  index = 0 ;
493  }
494  return ;
495 }
int G4int
Definition: G4Types.hh:78
bool G4bool
Definition: G4Types.hh:79
G4double G4DataInterpolation::CubicSplineInterpolation ( G4double  pX) const

Definition at line 311 of file G4DataInterpolation.cc.

312 {
313  G4int kLow=0, kHigh=fNumber-1, k=0 ;
314 
315  // Searching in the table by means of bisection method.
316  // fArgument must be monotonic, either increasing or decreasing
317 
318  while((kHigh - kLow) > 1)
319  {
320  k = (kHigh + kLow) >> 1 ; // compute midpoint 'bisection'
321  if(fArgument[k] > pX)
322  {
323  kHigh = k ;
324  }
325  else
326  {
327  kLow = k ;
328  }
329  } // kLow and kHigh now bracket the input value of pX
330  G4double deltaHL = fArgument[kHigh] - fArgument[kLow] ;
331  if (!(deltaHL != 0.0))
332  {
333  G4Exception("G4DataInterpolation::CubicSplineInterpolation()",
334  "Error", FatalException, "Bad fArgument input !") ;
335  }
336  G4double a = (fArgument[kHigh] - pX)/deltaHL ;
337  G4double b = (pX - fArgument[kLow])/deltaHL ;
338 
339  // Final evaluation of cubic spline polynomial for return
340 
341  return a*fFunction[kLow] + b*fFunction[kHigh] +
342  ((a*a*a - a)*fSecondDerivative[kLow] +
343  (b*b*b - b)*fSecondDerivative[kHigh])*deltaHL*deltaHL/6.0 ;
344 }
int G4int
Definition: G4Types.hh:78
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4DataInterpolation::FastCubicSpline ( G4double  pX,
G4int  index 
) const

Definition at line 353 of file G4DataInterpolation.cc.

355 {
356  G4double delta = fArgument[index+1] - fArgument[index] ;
357  if (!(delta != 0.0))
358  {
359  G4Exception("G4DataInterpolation::FastCubicSpline()",
360  "Error", FatalException, "Bad fArgument input !") ;
361  }
362  G4double a = (fArgument[index+1] - pX)/delta ;
363  G4double b = (pX - fArgument[index])/delta ;
364 
365  // Final evaluation of cubic spline polynomial for return
366 
367  return a*fFunction[index] + b*fFunction[index+1] +
368  ((a*a*a - a)*fSecondDerivative[index] +
369  (b*b*b - b)*fSecondDerivative[index+1])*delta*delta/6.0 ;
370 }
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4int G4DataInterpolation::LocateArgument ( G4double  pX) const

Definition at line 378 of file G4DataInterpolation.cc.

379 {
380  G4int kLow=-1, kHigh=fNumber, k=0 ;
381  G4bool ascend=(fArgument[fNumber-1] >= fArgument[0]) ;
382  while((kHigh - kLow) > 1)
383  {
384  k = (kHigh + kLow) >> 1 ; // compute midpoint 'bisection'
385  if( (pX >= fArgument[k]) == ascend)
386  {
387  kLow = k ;
388  }
389  else
390  {
391  kHigh = k ;
392  }
393  }
394  if (!(pX != fArgument[0]))
395  {
396  return 1 ;
397  }
398  else if (!(pX != fArgument[fNumber-1]))
399  {
400  return fNumber - 2 ;
401  }
402  else return kLow ;
403 }
int G4int
Definition: G4Types.hh:78
bool G4bool
Definition: G4Types.hh:79
void G4DataInterpolation::PolIntCoefficient ( G4double  cof[]) const

Definition at line 202 of file G4DataInterpolation.cc.

203 {
204  G4int i=0, j=0 ;
205  G4double factor;
206  G4double reducedY=0.0, mult=1.0 ;
207  G4double* tempArgument = new G4double[fNumber] ;
208 
209  for(i=0;i<fNumber;i++)
210  {
211  tempArgument[i] = cof[i] = 0.0 ;
212  }
213  tempArgument[fNumber-1] = -fArgument[0] ;
214 
215  for(i=1;i<fNumber;i++)
216  {
217  for(j=fNumber-1-i;j<fNumber-1;j++)
218  {
219  tempArgument[j] -= fArgument[i]*tempArgument[j+1] ;
220  }
221  tempArgument[fNumber-1] -= fArgument[i] ;
222  }
223  for(i=0;i<fNumber;i++)
224  {
225  factor = fNumber ;
226  for(j=fNumber-1;j>=1;j--)
227  {
228  factor = j*tempArgument[j] + factor*fArgument[i] ;
229  }
230  reducedY = fFunction[i]/factor ;
231  mult = 1.0 ;
232  for(j=fNumber-1;j>=0;j--)
233  {
234  cof[j] += mult*reducedY ;
235  mult = tempArgument[j] + mult*fArgument[i] ;
236  }
237  }
238  delete[] tempArgument ;
239 }
int G4int
Definition: G4Types.hh:78
double G4double
Definition: G4Types.hh:76
G4double G4DataInterpolation::PolynomInterpolation ( G4double  pX,
G4double deltaY 
) const

Definition at line 148 of file G4DataInterpolation.cc.

150 {
151  G4int i=0, j=1, k=0 ;
152  G4double mult=0.0, difi=0.0, deltaLow=0.0, deltaUp=0.0, cd=0.0, y=0.0 ;
153  G4double* c = new G4double[fNumber] ;
154  G4double* d = new G4double[fNumber] ;
155  G4double diff = std::fabs(pX-fArgument[0]) ;
156  for(i=0;i<fNumber;i++)
157  {
158  difi = std::fabs(pX-fArgument[i]) ;
159  if(difi <diff)
160  {
161  k = i ;
162  diff = difi ;
163  }
164  c[i] = fFunction[i] ;
165  d[i] = fFunction[i] ;
166  }
167  y = fFunction[k--] ;
168  for(j=1;j<fNumber;j++)
169  {
170  for(i=0;i<fNumber-j;i++)
171  {
172  deltaLow = fArgument[i] - pX ;
173  deltaUp = fArgument[i+j] - pX ;
174  cd = c[i+1] - d[i] ;
175  mult = deltaLow - deltaUp ;
176  if (!(mult != 0.0))
177  {
178  G4Exception("G4DataInterpolation::PolynomInterpolation()",
179  "Error", FatalException, "Coincident nodes !") ;
180  }
181  mult = cd/mult ;
182  d[i] = deltaUp*mult ;
183  c[i] = deltaLow*mult ;
184  }
185  y += (deltaY = (2*k < (fNumber - j -1) ? c[k+1] : d[k--] )) ;
186  }
187  delete[] c ;
188  delete[] d ;
189 
190  return y ;
191 }
static const G4double cd
int G4int
Definition: G4Types.hh:78
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4DataInterpolation::RationalPolInterpolation ( G4double  pX,
G4double deltaY 
) const

Definition at line 249 of file G4DataInterpolation.cc.

251 {
252  G4int i=0, j=1, k=0 ;
253  const G4double tolerance = 1.6e-24 ;
254  G4double mult=0.0, difi=0.0, cd=0.0, y=0.0, cof=0.0 ;
255  G4double* c = new G4double[fNumber] ;
256  G4double* d = new G4double[fNumber] ;
257  G4double diff = std::fabs(pX-fArgument[0]) ;
258  for(i=0;i<fNumber;i++)
259  {
260  difi = std::fabs(pX-fArgument[i]) ;
261  if (!(difi != 0.0))
262  {
263  y = fFunction[i] ;
264  deltaY = 0.0 ;
265  delete[] c ;
266  delete[] d ;
267  return y ;
268  }
269  else if(difi < diff)
270  {
271  k = i ;
272  diff = difi ;
273  }
274  c[i] = fFunction[i] ;
275  d[i] = fFunction[i] + tolerance ; // to prevent rare zero/zero cases
276  }
277  y = fFunction[k--] ;
278  for(j=1;j<fNumber;j++)
279  {
280  for(i=0;i<fNumber-j;i++)
281  {
282  cd = c[i+1] - d[i] ;
283  difi = fArgument[i+j] - pX ;
284  cof = (fArgument[i] - pX)*d[i]/difi ;
285  mult = cof - c[i+1] ;
286  if (!(mult != 0.0)) // function to be interpolated has pole at pX
287  {
288  G4Exception("G4DataInterpolation::RationalPolInterpolation()",
289  "Error", FatalException, "Coincident nodes !") ;
290  }
291  mult = cd/mult ;
292  d[i] = c[i+1]*mult ;
293  c[i] = cof*mult ;
294  }
295  y += (deltaY = (2*k < (fNumber - j - 1) ? c[k+1] : d[k--] )) ;
296  }
297  delete[] c ;
298  delete[] d ;
299 
300  return y ;
301 }
static const G4double cd
int G4int
Definition: G4Types.hh:78
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:


The documentation for this class was generated from the following files: