Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
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() [1/2]

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}
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66

◆ G4DataInterpolation() [2/2]

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}

◆ ~G4DataInterpolation()

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

◆ CorrelatedSearch()

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}
G4bool
bool G4bool
Definition: G4Types.hh:67

◆ CubicSplineInterpolation()

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}
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41

◆ FastCubicSpline()

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}

◆ LocateArgument()

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}

◆ PolIntCoefficient()

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}

◆ PolynomInterpolation()

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}

◆ RationalPolInterpolation()

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}
The documentation for this class was generated from the following files: