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,v 1.10 2008/03/13 09:35:56 gcosmo Exp $ |
---|
28 | // GEANT4 tag $Name: geant4-09-03 $ |
---|
29 | // |
---|
30 | #include "G4DataInterpolation.hh" |
---|
31 | |
---|
32 | ////////////////////////////////////////////////////////////////////////////// |
---|
33 | // |
---|
34 | // Constructor for initializing of fArgument, fFunction and fNumber |
---|
35 | // data members |
---|
36 | |
---|
37 | G4DataInterpolation::G4DataInterpolation( G4double pX[], |
---|
38 | G4double pY[], |
---|
39 | G4int number ) |
---|
40 | : fArgument(new G4double[number]), |
---|
41 | fFunction(new G4double[number]), |
---|
42 | fSecondDerivative(0), |
---|
43 | fNumber(number) |
---|
44 | { |
---|
45 | for(G4int i=0;i<fNumber;i++) |
---|
46 | { |
---|
47 | fArgument[i] = pX[i] ; |
---|
48 | fFunction[i] = pY[i] ; |
---|
49 | } |
---|
50 | } |
---|
51 | |
---|
52 | //////////////////////////////////////////////////////////////////////////// |
---|
53 | // |
---|
54 | // Constructor for cubic spline interpolation. It creates the array |
---|
55 | // fSecondDerivative[0,...fNumber-1] which is used in this interpolation by |
---|
56 | // the function |
---|
57 | |
---|
58 | |
---|
59 | G4DataInterpolation::G4DataInterpolation( G4double pX[], |
---|
60 | G4double pY[], |
---|
61 | G4int number, |
---|
62 | G4double pFirstDerStart, |
---|
63 | G4double pFirstDerFinish ) |
---|
64 | : fArgument(new G4double[number]), |
---|
65 | fFunction(new G4double[number]), |
---|
66 | fSecondDerivative(new G4double[number]), |
---|
67 | fNumber(number) |
---|
68 | { |
---|
69 | G4int i=0 ; |
---|
70 | G4double p=0.0, qn=0.0, sig=0.0, un=0.0 ; |
---|
71 | const G4double maxDerivative = 0.99e30 ; |
---|
72 | G4double* u = new G4double[fNumber - 1] ; |
---|
73 | |
---|
74 | for(i=0;i<fNumber;i++) |
---|
75 | { |
---|
76 | fArgument[i] = pX[i] ; |
---|
77 | fFunction[i] = pY[i] ; |
---|
78 | } |
---|
79 | if(pFirstDerStart > maxDerivative) |
---|
80 | { |
---|
81 | fSecondDerivative[0] = 0.0 ; |
---|
82 | u[0] = 0.0 ; |
---|
83 | } |
---|
84 | else |
---|
85 | { |
---|
86 | fSecondDerivative[0] = -0.5 ; |
---|
87 | u[0] = (3.0/(fArgument[1]-fArgument[0])) |
---|
88 | * ((fFunction[1]-fFunction[0])/(fArgument[1]-fArgument[0]) |
---|
89 | - pFirstDerStart) ; |
---|
90 | } |
---|
91 | |
---|
92 | // Decomposition loop for tridiagonal algorithm. fSecondDerivative[i] |
---|
93 | // and u[i] are used for temporary storage of the decomposed factors. |
---|
94 | |
---|
95 | for(i=1;i<fNumber-1;i++) |
---|
96 | { |
---|
97 | sig = (fArgument[i]-fArgument[i-1])/(fArgument[i+1]-fArgument[i-1]) ; |
---|
98 | p = sig*fSecondDerivative[i-1] + 2.0 ; |
---|
99 | fSecondDerivative[i] = (sig - 1.0)/p ; |
---|
100 | u[i] = (fFunction[i+1]-fFunction[i])/(fArgument[i+1]-fArgument[i]) - |
---|
101 | (fFunction[i]-fFunction[i-1])/(fArgument[i]-fArgument[i-1]) ; |
---|
102 | u[i] =(6.0*u[i]/(fArgument[i+1]-fArgument[i-1]) - sig*u[i-1])/p ; |
---|
103 | } |
---|
104 | if(pFirstDerFinish > maxDerivative) |
---|
105 | { |
---|
106 | qn = 0.0 ; |
---|
107 | un = 0.0 ; |
---|
108 | } |
---|
109 | else |
---|
110 | { |
---|
111 | qn = 0.5 ; |
---|
112 | un = (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] = (un - qn*u[fNumber-2])/ |
---|
117 | (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] = fSecondDerivative[k]*fSecondDerivative[k+1] + u[k]; |
---|
125 | } |
---|
126 | delete[] u ; |
---|
127 | } |
---|
128 | |
---|
129 | ///////////////////////////////////////////////////////////////////////////// |
---|
130 | // |
---|
131 | // Destructor deletes dynamically created arrays for data members: fArgument, |
---|
132 | // fFunction and fSecondDerivative, all have dimension of fNumber |
---|
133 | |
---|
134 | G4DataInterpolation::~G4DataInterpolation() |
---|
135 | { |
---|
136 | delete [] fArgument ; |
---|
137 | delete [] fFunction ; |
---|
138 | if(fSecondDerivative) { delete [] fSecondDerivative; } |
---|
139 | } |
---|
140 | |
---|
141 | ///////////////////////////////////////////////////////////////////////////// |
---|
142 | // |
---|
143 | // This function returns the value P(pX), where P(x) is polynom of fNumber-1 |
---|
144 | // degree such that P(fArgument[i]) = fFunction[i], for i = 0, ..., fNumber-1. |
---|
145 | // This is Lagrange's form of interpolation and it is based on Neville's |
---|
146 | // algorithm |
---|
147 | |
---|
148 | G4double |
---|
149 | G4DataInterpolation::PolynomInterpolation(G4double pX, |
---|
150 | G4double& deltaY ) const |
---|
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, y=0.0 ; |
---|
154 | G4double* c = new G4double[fNumber] ; |
---|
155 | G4double* d = new G4double[fNumber] ; |
---|
156 | G4double diff = std::fabs(pX-fArgument[0]) ; |
---|
157 | for(i=0;i<fNumber;i++) |
---|
158 | { |
---|
159 | difi = std::fabs(pX-fArgument[i]) ; |
---|
160 | if(difi <diff) |
---|
161 | { |
---|
162 | k = i ; |
---|
163 | diff = difi ; |
---|
164 | } |
---|
165 | c[i] = fFunction[i] ; |
---|
166 | d[i] = fFunction[i] ; |
---|
167 | } |
---|
168 | y = fFunction[k--] ; |
---|
169 | for(j=1;j<fNumber;j++) |
---|
170 | { |
---|
171 | for(i=0;i<fNumber-j;i++) |
---|
172 | { |
---|
173 | deltaLow = fArgument[i] - pX ; |
---|
174 | deltaUp = fArgument[i+j] - pX ; |
---|
175 | cd = c[i+1] - d[i] ; |
---|
176 | mult = deltaLow - deltaUp ; |
---|
177 | if (!(mult != 0.0)) |
---|
178 | { |
---|
179 | G4Exception("G4DataInterpolation::PolynomInterpolation()", |
---|
180 | "Error", FatalException, "Coincident nodes !") ; |
---|
181 | } |
---|
182 | mult = cd/mult ; |
---|
183 | d[i] = deltaUp*mult ; |
---|
184 | c[i] = deltaLow*mult ; |
---|
185 | } |
---|
186 | y += (deltaY = (2*k < (fNumber - j -1) ? c[k+1] : d[k--] )) ; |
---|
187 | } |
---|
188 | delete[] c ; |
---|
189 | delete[] d ; |
---|
190 | |
---|
191 | return y ; |
---|
192 | } |
---|
193 | |
---|
194 | //////////////////////////////////////////////////////////////////////////// |
---|
195 | // |
---|
196 | // Given arrays fArgument[0,..,fNumber-1] and fFunction[0,..,fNumber-1], this |
---|
197 | // function calculates an array of coefficients. The coefficients don't provide |
---|
198 | // usually (fNumber>10) better accuracy for polynom interpolation, as compared |
---|
199 | // with PolynomInterpolation function. They could be used instead for derivate |
---|
200 | // calculations and some other applications. |
---|
201 | |
---|
202 | void |
---|
203 | G4DataInterpolation::PolIntCoefficient( G4double cof[]) const |
---|
204 | { |
---|
205 | G4int i=0, j=0 ; |
---|
206 | G4double factor=fNumber, 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 | |
---|
241 | ///////////////////////////////////////////////////////////////////////////// |
---|
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 |
---|
249 | G4DataInterpolation::RationalPolInterpolation(G4double pX, |
---|
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 | |
---|
303 | ///////////////////////////////////////////////////////////////////////////// |
---|
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 |
---|
311 | G4DataInterpolation::CubicSplineInterpolation(G4double pX) const |
---|
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 | |
---|
346 | /////////////////////////////////////////////////////////////////////////// |
---|
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 |
---|
353 | G4DataInterpolation::FastCubicSpline(G4double pX, |
---|
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 | |
---|
372 | //////////////////////////////////////////////////////////////////////////// |
---|
373 | // |
---|
374 | // Given argument pX, returns index k, so that pX bracketed by fArgument[k] |
---|
375 | // and fArgument[k+1] |
---|
376 | |
---|
377 | G4int |
---|
378 | G4DataInterpolation::LocateArgument(G4double pX) const |
---|
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 | |
---|
405 | ///////////////////////////////////////////////////////////////////////////// |
---|
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 |
---|
414 | G4DataInterpolation::CorrelatedSearch( G4double pX, |
---|
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 | // |
---|
499 | //////////////////////////////////////////////////////////////////////////// |
---|