Geant4 10.7.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 ()
 
 G4DataInterpolation (const G4DataInterpolation &)=delete
 
G4DataInterpolation & operator= (const G4DataInterpolation &)=delete
 
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 42 of file G4DataInterpolation.hh.

Constructor & Destructor Documentation

◆ G4DataInterpolation() [1/3]

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

Definition at line 38 of file G4DataInterpolation.cc.

40 : fArgument(new G4double[number])
41 , fFunction(new G4double[number])
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:83
G4int
int G4int
Definition: G4Types.hh:85

◆ G4DataInterpolation() [2/3]

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

Definition at line 57 of file G4DataInterpolation.cc.

60 : fArgument(new G4double[number])
61 , fFunction(new G4double[number])
62 , fSecondDerivative(new G4double[number])
63 , fNumber(number)
64{
65 G4int i = 0;
66 G4double p = 0.0, qn = 0.0, sig = 0.0, un = 0.0;
67 const G4double maxDerivative = 0.99e30;
68 G4double* u = new G4double[fNumber - 1];
69
70 for(i = 0; i < fNumber; ++i)
71 {
72 fArgument[i] = pX[i];
73 fFunction[i] = pY[i];
74 }
75 if(pFirstDerStart > maxDerivative)
76 {
77 fSecondDerivative[0] = 0.0;
78 u[0] = 0.0;
79 }
80 else
81 {
82 fSecondDerivative[0] = -0.5;
83 u[0] = (3.0 / (fArgument[1] - fArgument[0])) *
84 ((fFunction[1] - fFunction[0]) / (fArgument[1] - fArgument[0]) -
85 pFirstDerStart);
86 }
87
88 // Decomposition loop for tridiagonal algorithm. fSecondDerivative[i]
89 // and u[i] are used for temporary storage of the decomposed factors.
90
91 for(i = 1; i < fNumber - 1; ++i)
92 {
93 sig =
94 (fArgument[i] - fArgument[i - 1]) / (fArgument[i + 1] - fArgument[i - 1]);
95 p = sig * fSecondDerivative[i - 1] + 2.0;
96 fSecondDerivative[i] = (sig - 1.0) / p;
97 u[i] =
98 (fFunction[i + 1] - fFunction[i]) / (fArgument[i + 1] - fArgument[i]) -
99 (fFunction[i] - fFunction[i - 1]) / (fArgument[i] - fArgument[i - 1]);
100 u[i] =
101 (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 =
112 (3.0 / (fArgument[fNumber - 1] - fArgument[fNumber - 2])) *
113 (pFirstDerFinish - (fFunction[fNumber - 1] - fFunction[fNumber - 2]) /
114 (fArgument[fNumber - 1] - fArgument[fNumber - 2]));
115 }
116 fSecondDerivative[fNumber - 1] =
117 (un - qn * u[fNumber - 2]) / (qn * fSecondDerivative[fNumber - 2] + 1.0);
118
119 // The backsubstitution loop for the triagonal algorithm of solving
120 // a linear system of equations.
121
122 for(G4int k = fNumber - 2; k >= 0; --k)
123 {
124 fSecondDerivative[k] =
125 fSecondDerivative[k] * fSecondDerivative[k + 1] + u[k];
126 }
127 delete[] u;
128}

◆ ~G4DataInterpolation()

G4DataInterpolation::~G4DataInterpolation ( )

Definition at line 135 of file G4DataInterpolation.cc.

136{
137 delete[] fArgument;
138 delete[] fFunction;
139 delete[] fSecondDerivative;
140}

◆ G4DataInterpolation() [3/3]

G4DataInterpolation::G4DataInterpolation ( const G4DataInterpolation &  )
delete

Member Function Documentation

◆ CorrelatedSearch()

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

Definition at line 412 of file G4DataInterpolation.cc.

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

◆ CubicSplineInterpolation()

G4double G4DataInterpolation::CubicSplineInterpolation ( G4double  pX) const

Definition at line 310 of file G4DataInterpolation.cc.

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()", "Error",
333 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]) *
343 deltaHL * deltaHL / 6.0;
344}
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:61
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35

◆ FastCubicSpline()

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

Definition at line 352 of file G4DataInterpolation.cc.

353{
354 G4double delta = fArgument[index + 1] - fArgument[index];
355 if(!(delta != 0.0))
356 {
357 G4Exception("G4DataInterpolation::FastCubicSpline()", "Error",
358 FatalException, "Bad fArgument input !");
359 }
360 G4double a = (fArgument[index + 1] - pX) / delta;
361 G4double b = (pX - fArgument[index]) / delta;
362
363 // Final evaluation of cubic spline polynomial for return
364
365 return a * fFunction[index] + b * fFunction[index + 1] +
366 ((a * a * a - a) * fSecondDerivative[index] +
367 (b * b * b - b) * fSecondDerivative[index + 1]) *
368 delta * delta / 6.0;
369}

◆ LocateArgument()

G4int G4DataInterpolation::LocateArgument ( G4double  pX) const

Definition at line 376 of file G4DataInterpolation.cc.

377{
378 G4int kLow = -1, kHigh = fNumber, k = 0;
379 G4bool ascend = (fArgument[fNumber - 1] >= fArgument[0]);
380 while((kHigh - kLow) > 1)
381 {
382 k = (kHigh + kLow) >> 1; // compute midpoint 'bisection'
383 if((pX >= fArgument[k]) == ascend)
384 {
385 kLow = k;
386 }
387 else
388 {
389 kHigh = k;
390 }
391 }
392 if(!(pX != fArgument[0]))
393 {
394 return 1;
395 }
396 else if(!(pX != fArgument[fNumber - 1]))
397 {
398 return fNumber - 2;
399 }
400 else
401 return kLow;
402}

◆ operator=()

G4DataInterpolation & G4DataInterpolation::operator= ( const G4DataInterpolation &  )
delete

◆ PolIntCoefficient()

void G4DataInterpolation::PolIntCoefficient ( G4double  cof[]) const

Definition at line 203 of file G4DataInterpolation.cc.

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

◆ PolynomInterpolation()

G4double G4DataInterpolation::PolynomInterpolation ( G4double  pX,
G4double &  deltaY 
) const

Definition at line 149 of file G4DataInterpolation.cc.

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

◆ 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()", "Error",
289 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: