Geant4  10.00.p02
G4DataInterpolation.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4DataInterpolation.cc 67970 2013-03-13 10:10:06Z gcosmo $
28 //
29 #include "G4DataInterpolation.hh"
30 
32 //
33 // Constructor for initializing of fArgument, fFunction and fNumber
34 // data members
35 
37  G4double pY[],
38  G4int number )
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 }
50 
52 //
53 // Constructor for cubic spline interpolation. It creates the array
54 // fSecondDerivative[0,...fNumber-1] which is used in this interpolation by
55 // the function
56 
57 
59  G4double pY[],
60  G4int number,
61  G4double pFirstDerStart,
62  G4double pFirstDerFinish )
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  {
124  }
125  delete[] u ;
126 }
127 
129 //
130 // Destructor deletes dynamically created arrays for data members: fArgument,
131 // fFunction and fSecondDerivative, all have dimension of fNumber
132 
134 {
135  delete [] fArgument ;
136  delete [] fFunction ;
137  if(fSecondDerivative) { delete [] fSecondDerivative; }
138 }
139 
141 //
142 // This function returns the value P(pX), where P(x) is polynom of fNumber-1
143 // degree such that P(fArgument[i]) = fFunction[i], for i = 0, ..., fNumber-1.
144 // This is Lagrange's form of interpolation and it is based on Neville's
145 // algorithm
146 
147 G4double
149  G4double& deltaY ) const
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 }
192 
194 //
195 // Given arrays fArgument[0,..,fNumber-1] and fFunction[0,..,fNumber-1], this
196 // function calculates an array of coefficients. The coefficients don't provide
197 // usually (fNumber>10) better accuracy for polynom interpolation, as compared
198 // with PolynomInterpolation function. They could be used instead for derivate
199 // calculations and some other applications.
200 
201 void
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 }
240 
242 //
243 // The function returns diagonal rational function (Bulirsch and Stoer
244 // algorithm of Neville type) Pn(x)/Qm(x) where P and Q are polynoms.
245 // Tests showed the method is not stable and hasn't advantage if compared
246 // with polynomial interpolation ?!
247 
248 G4double
250  G4double& deltaY ) const
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 }
302 
304 //
305 // Cubic spline interpolation in point pX for function given by the table:
306 // fArgument, fFunction. The constructor, which creates fSecondDerivative,
307 // must be called before. The function works optimal, if sequential calls
308 // are in random values of pX.
309 
310 G4double
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 }
345 
347 //
348 // Return cubic spline interpolation in the point pX which is located between
349 // fArgument[index] and fArgument[index+1]. It is usually called in sequence
350 // of known from external analysis values of index.
351 
352 G4double
354  G4int index) const
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 }
371 
373 //
374 // Given argument pX, returns index k, so that pX bracketed by fArgument[k]
375 // and fArgument[k+1]
376 
377 G4int
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 }
404 
406 //
407 // Given a value pX, returns a value 'index' such that pX is between
408 // fArgument[index] and fArgument[index+1]. fArgument MUST BE MONOTONIC,
409 // either increasing or decreasing. If index = -1 or fNumber, this indicates
410 // that pX is out of range. The value index on input is taken as the initial
411 // approximation for index on output.
412 
413 void
415  G4int& index ) const
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 }
496 
497 //
498 //
G4double RationalPolInterpolation(G4double pX, G4double &deltaY) const
G4double FastCubicSpline(G4double pX, G4int index) const
void CorrelatedSearch(G4double pX, G4int &index) const
static const G4double cd
G4double a
Definition: TRTMaterials.hh:39
G4DataInterpolation(G4double pX[], G4double pY[], G4int number)
G4int LocateArgument(G4double pX) const
int G4int
Definition: G4Types.hh:78
G4double CubicSplineInterpolation(G4double pX) const
bool G4bool
Definition: G4Types.hh:79
void PolIntCoefficient(G4double cof[]) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double G4double
Definition: G4Types.hh:76
G4double PolynomInterpolation(G4double pX, G4double &deltaY) const