ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DataInterpolation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DataInterpolation.cc
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 //
28 #include "G4DataInterpolation.hh"
29 
31 //
32 // Constructor for initializing of fArgument, fFunction and fNumber
33 // data members
34 
36  G4double pY[],
37  G4int number )
38  : fArgument(new G4double[number]),
39  fFunction(new G4double[number]),
40  fSecondDerivative(0),
41  fNumber(number)
42 {
43  for(G4int i=0;i<fNumber;i++)
44  {
45  fArgument[i] = pX[i] ;
46  fFunction[i] = pY[i] ;
47  }
48 }
49 
51 //
52 // Constructor for cubic spline interpolation. It creates the array
53 // fSecondDerivative[0,...fNumber-1] which is used in this interpolation by
54 // the function
55 
56 
58  G4double pY[],
59  G4int number,
60  G4double pFirstDerStart,
61  G4double pFirstDerFinish )
62  : fArgument(new G4double[number]),
63  fFunction(new G4double[number]),
64  fSecondDerivative(new G4double[number]),
65  fNumber(number)
66 {
67  G4int i=0 ;
68  G4double p=0.0, qn=0.0, sig=0.0, un=0.0 ;
69  const G4double maxDerivative = 0.99e30 ;
70  G4double* u = new G4double[fNumber - 1] ;
71 
72  for(i=0;i<fNumber;i++)
73  {
74  fArgument[i] = pX[i] ;
75  fFunction[i] = pY[i] ;
76  }
77  if(pFirstDerStart > maxDerivative)
78  {
79  fSecondDerivative[0] = 0.0 ;
80  u[0] = 0.0 ;
81  }
82  else
83  {
84  fSecondDerivative[0] = -0.5 ;
85  u[0] = (3.0/(fArgument[1]-fArgument[0]))
86  * ((fFunction[1]-fFunction[0])/(fArgument[1]-fArgument[0])
87  - pFirstDerStart) ;
88  }
89 
90  // Decomposition loop for tridiagonal algorithm. fSecondDerivative[i]
91  // and u[i] are used for temporary storage of the decomposed factors.
92 
93  for(i=1;i<fNumber-1;i++)
94  {
95  sig = (fArgument[i]-fArgument[i-1])/(fArgument[i+1]-fArgument[i-1]) ;
96  p = sig*fSecondDerivative[i-1] + 2.0 ;
97  fSecondDerivative[i] = (sig - 1.0)/p ;
98  u[i] = (fFunction[i+1]-fFunction[i])/(fArgument[i+1]-fArgument[i]) -
99  (fFunction[i]-fFunction[i-1])/(fArgument[i]-fArgument[i-1]) ;
100  u[i] =(6.0*u[i]/(fArgument[i+1]-fArgument[i-1]) - sig*u[i-1])/p ;
101  }
102  if(pFirstDerFinish > maxDerivative)
103  {
104  qn = 0.0 ;
105  un = 0.0 ;
106  }
107  else
108  {
109  qn = 0.5 ;
110  un = (3.0/(fArgument[fNumber-1]-fArgument[fNumber-2]))
111  * (pFirstDerFinish - (fFunction[fNumber-1]-fFunction[fNumber-2])
112  / (fArgument[fNumber-1]-fArgument[fNumber-2])) ;
113  }
114  fSecondDerivative[fNumber-1] = (un - qn*u[fNumber-2])/
115  (qn*fSecondDerivative[fNumber-2] + 1.0) ;
116 
117  // The backsubstitution loop for the triagonal algorithm of solving
118  // a linear system of equations.
119 
120  for(G4int k=fNumber-2;k>=0;k--)
121  {
123  }
124  delete[] u ;
125 }
126 
128 //
129 // Destructor deletes dynamically created arrays for data members: fArgument,
130 // fFunction and fSecondDerivative, all have dimension of fNumber
131 
133 {
134  delete [] fArgument ;
135  delete [] fFunction ;
136  if(fSecondDerivative) { delete [] fSecondDerivative; }
137 }
138 
140 //
141 // This function returns the value P(pX), where P(x) is polynom of fNumber-1
142 // degree such that P(fArgument[i]) = fFunction[i], for i = 0, ..., fNumber-1.
143 // This is Lagrange's form of interpolation and it is based on Neville's
144 // algorithm
145 
146 G4double
148  G4double& deltaY ) const
149 {
150  G4int i=0, j=1, k=0 ;
151  G4double mult=0.0, difi=0.0, deltaLow=0.0, deltaUp=0.0, cd=0.0, y=0.0 ;
152  G4double* c = new G4double[fNumber] ;
153  G4double* d = new G4double[fNumber] ;
154  G4double diff = std::fabs(pX-fArgument[0]) ;
155  for(i=0;i<fNumber;i++)
156  {
157  difi = std::fabs(pX-fArgument[i]) ;
158  if(difi <diff)
159  {
160  k = i ;
161  diff = difi ;
162  }
163  c[i] = fFunction[i] ;
164  d[i] = fFunction[i] ;
165  }
166  y = fFunction[k--] ;
167  for(j=1;j<fNumber;j++)
168  {
169  for(i=0;i<fNumber-j;i++)
170  {
171  deltaLow = fArgument[i] - pX ;
172  deltaUp = fArgument[i+j] - pX ;
173  cd = c[i+1] - d[i] ;
174  mult = deltaLow - deltaUp ;
175  if (!(mult != 0.0))
176  {
177  G4Exception("G4DataInterpolation::PolynomInterpolation()",
178  "Error", FatalException, "Coincident nodes !") ;
179  }
180  mult = cd/mult ;
181  d[i] = deltaUp*mult ;
182  c[i] = deltaLow*mult ;
183  }
184  y += (deltaY = (2*k < (fNumber - j -1) ? c[k+1] : d[k--] )) ;
185  }
186  delete[] c ;
187  delete[] d ;
188 
189  return y ;
190 }
191 
193 //
194 // Given arrays fArgument[0,..,fNumber-1] and fFunction[0,..,fNumber-1], this
195 // function calculates an array of coefficients. The coefficients don't provide
196 // usually (fNumber>10) better accuracy for polynom interpolation, as compared
197 // with PolynomInterpolation function. They could be used instead for derivate
198 // calculations and some other applications.
199 
200 void
202 {
203  G4int i=0, j=0 ;
204  G4double factor;
205  G4double reducedY=0.0, mult=1.0 ;
206  G4double* tempArgument = new G4double[fNumber] ;
207 
208  for(i=0;i<fNumber;i++)
209  {
210  tempArgument[i] = cof[i] = 0.0 ;
211  }
212  tempArgument[fNumber-1] = -fArgument[0] ;
213 
214  for(i=1;i<fNumber;i++)
215  {
216  for(j=fNumber-1-i;j<fNumber-1;j++)
217  {
218  tempArgument[j] -= fArgument[i]*tempArgument[j+1] ;
219  }
220  tempArgument[fNumber-1] -= fArgument[i] ;
221  }
222  for(i=0;i<fNumber;i++)
223  {
224  factor = fNumber ;
225  for(j=fNumber-1;j>=1;j--)
226  {
227  factor = j*tempArgument[j] + factor*fArgument[i] ;
228  }
229  reducedY = fFunction[i]/factor ;
230  mult = 1.0 ;
231  for(j=fNumber-1;j>=0;j--)
232  {
233  cof[j] += mult*reducedY ;
234  mult = tempArgument[j] + mult*fArgument[i] ;
235  }
236  }
237  delete[] tempArgument ;
238 }
239 
241 //
242 // The function returns diagonal rational function (Bulirsch and Stoer
243 // algorithm of Neville type) Pn(x)/Qm(x) where P and Q are polynoms.
244 // Tests showed the method is not stable and hasn't advantage if compared
245 // with polynomial interpolation ?!
246 
247 G4double
249  G4double& deltaY ) const
250 {
251  G4int i=0, j=1, k=0 ;
252  const G4double tolerance = 1.6e-24 ;
253  G4double mult=0.0, difi=0.0, cd=0.0, y=0.0, cof=0.0 ;
254  G4double* c = new G4double[fNumber] ;
255  G4double* d = new G4double[fNumber] ;
256  G4double diff = std::fabs(pX-fArgument[0]) ;
257  for(i=0;i<fNumber;i++)
258  {
259  difi = std::fabs(pX-fArgument[i]) ;
260  if (!(difi != 0.0))
261  {
262  y = fFunction[i] ;
263  deltaY = 0.0 ;
264  delete[] c ;
265  delete[] d ;
266  return y ;
267  }
268  else if(difi < diff)
269  {
270  k = i ;
271  diff = difi ;
272  }
273  c[i] = fFunction[i] ;
274  d[i] = fFunction[i] + tolerance ; // to prevent rare zero/zero cases
275  }
276  y = fFunction[k--] ;
277  for(j=1;j<fNumber;j++)
278  {
279  for(i=0;i<fNumber-j;i++)
280  {
281  cd = c[i+1] - d[i] ;
282  difi = fArgument[i+j] - pX ;
283  cof = (fArgument[i] - pX)*d[i]/difi ;
284  mult = cof - c[i+1] ;
285  if (!(mult != 0.0)) // function to be interpolated has pole at pX
286  {
287  G4Exception("G4DataInterpolation::RationalPolInterpolation()",
288  "Error", FatalException, "Coincident nodes !") ;
289  }
290  mult = cd/mult ;
291  d[i] = c[i+1]*mult ;
292  c[i] = cof*mult ;
293  }
294  y += (deltaY = (2*k < (fNumber - j - 1) ? c[k+1] : d[k--] )) ;
295  }
296  delete[] c ;
297  delete[] d ;
298 
299  return y ;
300 }
301 
303 //
304 // Cubic spline interpolation in point pX for function given by the table:
305 // fArgument, fFunction. The constructor, which creates fSecondDerivative,
306 // must be called before. The function works optimal, if sequential calls
307 // are in random values of pX.
308 
309 G4double
311 {
312  G4int kLow=0, kHigh=fNumber-1, k=0 ;
313 
314  // Searching in the table by means of bisection method.
315  // fArgument must be monotonic, either increasing or decreasing
316 
317  while((kHigh - kLow) > 1)
318  {
319  k = (kHigh + kLow) >> 1 ; // compute midpoint 'bisection'
320  if(fArgument[k] > pX)
321  {
322  kHigh = k ;
323  }
324  else
325  {
326  kLow = k ;
327  }
328  } // kLow and kHigh now bracket the input value of pX
329  G4double deltaHL = fArgument[kHigh] - fArgument[kLow] ;
330  if (!(deltaHL != 0.0))
331  {
332  G4Exception("G4DataInterpolation::CubicSplineInterpolation()",
333  "Error", FatalException, "Bad fArgument input !") ;
334  }
335  G4double a = (fArgument[kHigh] - pX)/deltaHL ;
336  G4double b = (pX - fArgument[kLow])/deltaHL ;
337 
338  // Final evaluation of cubic spline polynomial for return
339 
340  return a*fFunction[kLow] + b*fFunction[kHigh] +
341  ((a*a*a - a)*fSecondDerivative[kLow] +
342  (b*b*b - b)*fSecondDerivative[kHigh])*deltaHL*deltaHL/6.0 ;
343 }
344 
346 //
347 // Return cubic spline interpolation in the point pX which is located between
348 // fArgument[index] and fArgument[index+1]. It is usually called in sequence
349 // of known from external analysis values of index.
350 
351 G4double
353  G4int index) const
354 {
355  G4double delta = fArgument[index+1] - fArgument[index] ;
356  if (!(delta != 0.0))
357  {
358  G4Exception("G4DataInterpolation::FastCubicSpline()",
359  "Error", FatalException, "Bad fArgument input !") ;
360  }
361  G4double a = (fArgument[index+1] - pX)/delta ;
362  G4double b = (pX - fArgument[index])/delta ;
363 
364  // Final evaluation of cubic spline polynomial for return
365 
366  return a*fFunction[index] + b*fFunction[index+1] +
367  ((a*a*a - a)*fSecondDerivative[index] +
368  (b*b*b - b)*fSecondDerivative[index+1])*delta*delta/6.0 ;
369 }
370 
372 //
373 // Given argument pX, returns index k, so that pX bracketed by fArgument[k]
374 // and fArgument[k+1]
375 
376 G4int
378 {
379  G4int kLow=-1, kHigh=fNumber, k=0 ;
380  G4bool ascend=(fArgument[fNumber-1] >= fArgument[0]) ;
381  while((kHigh - kLow) > 1)
382  {
383  k = (kHigh + kLow) >> 1 ; // compute midpoint 'bisection'
384  if( (pX >= fArgument[k]) == ascend)
385  {
386  kLow = k ;
387  }
388  else
389  {
390  kHigh = k ;
391  }
392  }
393  if (!(pX != fArgument[0]))
394  {
395  return 1 ;
396  }
397  else if (!(pX != fArgument[fNumber-1]))
398  {
399  return fNumber - 2 ;
400  }
401  else return kLow ;
402 }
403 
405 //
406 // Given a value pX, returns a value 'index' such that pX is between
407 // fArgument[index] and fArgument[index+1]. fArgument MUST BE MONOTONIC,
408 // either increasing or decreasing. If index = -1 or fNumber, this indicates
409 // that pX is out of range. The value index on input is taken as the initial
410 // approximation for index on output.
411 
412 void
414  G4int& index ) const
415 {
416  G4int kHigh=0, k=0, Increment=0 ;
417  // ascend = true for ascending order of table, false otherwise
418  G4bool ascend = (fArgument[fNumber-1] >= fArgument[0]) ;
419  if(index < 0 || index > fNumber-1)
420  {
421  index = -1 ;
422  kHigh = fNumber ;
423  }
424  else
425  {
426  Increment = 1 ; // What value would be the best ?
427  if((pX >= fArgument[index]) == ascend)
428  {
429  if(index == fNumber -1)
430  {
431  index = fNumber ;
432  return ;
433  }
434  kHigh = index + 1 ;
435  while((pX >= fArgument[kHigh]) == ascend)
436  {
437  index = kHigh ;
438  Increment += Increment ; // double the Increment
439  kHigh = index + Increment ;
440  if(kHigh > (fNumber - 1))
441  {
442  kHigh = fNumber ;
443  break ;
444  }
445  }
446  }
447  else
448  {
449  if(index == 0)
450  {
451  index = -1 ;
452  return ;
453  }
454  kHigh = index-- ;
455  while((pX < fArgument[index]) == ascend)
456  {
457  kHigh = index ;
458  Increment <<= 1 ; // double the Increment
459  if(Increment >= kHigh)
460  {
461  index = -1 ;
462  break ;
463  }
464  else
465  {
466  index = kHigh - Increment ;
467  }
468  }
469  } // Value bracketed
470  }
471  // final bisection searching
472 
473  while((kHigh - index) != 1)
474  {
475  k = (kHigh + index) >> 1 ;
476  if((pX >= fArgument[k]) == ascend)
477  {
478  index = k ;
479  }
480  else
481  {
482  kHigh = k ;
483  }
484  }
485  if (!(pX != fArgument[fNumber-1]))
486  {
487  index = fNumber - 2 ;
488  }
489  if (!(pX != fArgument[0]))
490  {
491  index = 0 ;
492  }
493  return ;
494 }
495 
496 //
497 //