| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4DataInterpolation class implementation 26 // G4DataInterpolation class implementation
27 // 27 //
28 // Author: V.Grichine, 03.04.1997 28 // Author: V.Grichine, 03.04.1997
29 // ------------------------------------------- 29 // --------------------------------------------------------------------
30 30
31 #include "G4DataInterpolation.hh" 31 #include "G4DataInterpolation.hh"
32 32
33 ////////////////////////////////////////////// 33 //////////////////////////////////////////////////////////////////////////////
34 // 34 //
35 // Constructor for initializing of fArgument, 35 // Constructor for initializing of fArgument, fFunction and fNumber
36 // data members 36 // data members
37 37
38 G4DataInterpolation::G4DataInterpolation(G4dou 38 G4DataInterpolation::G4DataInterpolation(G4double pX[], G4double pY[],
39 G4int 39 G4int number)
40 : fArgument(new G4double[number]) 40 : fArgument(new G4double[number])
41 , fFunction(new G4double[number]) 41 , fFunction(new G4double[number])
42 , fNumber(number) 42 , fNumber(number)
43 { 43 {
44 for(G4int i = 0; i < fNumber; ++i) 44 for(G4int i = 0; i < fNumber; ++i)
45 { 45 {
46 fArgument[i] = pX[i]; 46 fArgument[i] = pX[i];
47 fFunction[i] = pY[i]; 47 fFunction[i] = pY[i];
48 } 48 }
49 } 49 }
50 50
51 ////////////////////////////////////////////// 51 ////////////////////////////////////////////////////////////////////////////
52 // 52 //
53 // Constructor for cubic spline interpolation. 53 // Constructor for cubic spline interpolation. It creates the array
54 // fSecondDerivative[0,...fNumber-1] which is 54 // fSecondDerivative[0,...fNumber-1] which is used in this interpolation by
55 // the function 55 // the function
56 56
57 G4DataInterpolation::G4DataInterpolation(G4dou 57 G4DataInterpolation::G4DataInterpolation(G4double pX[], G4double pY[],
58 G4int 58 G4int number, G4double pFirstDerStart,
59 G4dou 59 G4double pFirstDerFinish)
60 : fArgument(new G4double[number]) 60 : fArgument(new G4double[number])
61 , fFunction(new G4double[number]) 61 , fFunction(new G4double[number])
62 , fSecondDerivative(new G4double[number]) 62 , fSecondDerivative(new G4double[number])
63 , fNumber(number) 63 , fNumber(number)
64 { 64 {
65 G4int i = 0; 65 G4int i = 0;
66 G4double p = 0.0, qn = 0.0, sig = 0.0, un = 66 G4double p = 0.0, qn = 0.0, sig = 0.0, un = 0.0;
67 const G4double maxDerivative = 0.99e30; 67 const G4double maxDerivative = 0.99e30;
68 auto u = new G4double[fNumb << 68 G4double* u = new G4double[fNumber - 1];
69 69
70 for(i = 0; i < fNumber; ++i) 70 for(i = 0; i < fNumber; ++i)
71 { 71 {
72 fArgument[i] = pX[i]; 72 fArgument[i] = pX[i];
73 fFunction[i] = pY[i]; 73 fFunction[i] = pY[i];
74 } 74 }
75 if(pFirstDerStart > maxDerivative) 75 if(pFirstDerStart > maxDerivative)
76 { 76 {
77 fSecondDerivative[0] = 0.0; 77 fSecondDerivative[0] = 0.0;
78 u[0] = 0.0; 78 u[0] = 0.0;
79 } 79 }
80 else 80 else
81 { 81 {
82 fSecondDerivative[0] = -0.5; 82 fSecondDerivative[0] = -0.5;
83 u[0] = (3.0 / (fArgument[1 83 u[0] = (3.0 / (fArgument[1] - fArgument[0])) *
84 ((fFunction[1] - fFunction[0]) / (f 84 ((fFunction[1] - fFunction[0]) / (fArgument[1] - fArgument[0]) -
85 pFirstDerStart); 85 pFirstDerStart);
86 } 86 }
87 87
88 // Decomposition loop for tridiagonal algori 88 // Decomposition loop for tridiagonal algorithm. fSecondDerivative[i]
89 // and u[i] are used for temporary storage o 89 // and u[i] are used for temporary storage of the decomposed factors.
90 90
91 for(i = 1; i < fNumber - 1; ++i) 91 for(i = 1; i < fNumber - 1; ++i)
92 { 92 {
93 sig = 93 sig =
94 (fArgument[i] - fArgument[i - 1]) / (fAr 94 (fArgument[i] - fArgument[i - 1]) / (fArgument[i + 1] - fArgument[i - 1]);
95 p = sig * fSecondDeriva 95 p = sig * fSecondDerivative[i - 1] + 2.0;
96 fSecondDerivative[i] = (sig - 1.0) / p; 96 fSecondDerivative[i] = (sig - 1.0) / p;
97 u[i] = 97 u[i] =
98 (fFunction[i + 1] - fFunction[i]) / (fAr 98 (fFunction[i + 1] - fFunction[i]) / (fArgument[i + 1] - fArgument[i]) -
99 (fFunction[i] - fFunction[i - 1]) / (fAr 99 (fFunction[i] - fFunction[i - 1]) / (fArgument[i] - fArgument[i - 1]);
100 u[i] = 100 u[i] =
101 (6.0 * u[i] / (fArgument[i + 1] - fArgum 101 (6.0 * u[i] / (fArgument[i + 1] - fArgument[i - 1]) - sig * u[i - 1]) / p;
102 } 102 }
103 if(pFirstDerFinish > maxDerivative) 103 if(pFirstDerFinish > maxDerivative)
104 { 104 {
105 qn = 0.0; 105 qn = 0.0;
106 un = 0.0; 106 un = 0.0;
107 } 107 }
108 else 108 else
109 { 109 {
110 qn = 0.5; 110 qn = 0.5;
111 un = 111 un =
112 (3.0 / (fArgument[fNumber - 1] - fArgume 112 (3.0 / (fArgument[fNumber - 1] - fArgument[fNumber - 2])) *
113 (pFirstDerFinish - (fFunction[fNumber - 113 (pFirstDerFinish - (fFunction[fNumber - 1] - fFunction[fNumber - 2]) /
114 (fArgument[fNumber 114 (fArgument[fNumber - 1] - fArgument[fNumber - 2]));
115 } 115 }
116 fSecondDerivative[fNumber - 1] = 116 fSecondDerivative[fNumber - 1] =
117 (un - qn * u[fNumber - 2]) / (qn * fSecond 117 (un - qn * u[fNumber - 2]) / (qn * fSecondDerivative[fNumber - 2] + 1.0);
118 118
119 // The backsubstitution loop for the triagon 119 // The backsubstitution loop for the triagonal algorithm of solving
120 // a linear system of equations. 120 // a linear system of equations.
121 121
122 for(G4int k = fNumber - 2; k >= 0; --k) 122 for(G4int k = fNumber - 2; k >= 0; --k)
123 { 123 {
124 fSecondDerivative[k] = 124 fSecondDerivative[k] =
125 fSecondDerivative[k] * fSecondDerivative 125 fSecondDerivative[k] * fSecondDerivative[k + 1] + u[k];
126 } 126 }
127 delete[] u; 127 delete[] u;
128 } 128 }
129 129
130 ////////////////////////////////////////////// 130 /////////////////////////////////////////////////////////////////////////////
131 // 131 //
132 // Destructor deletes dynamically created arra 132 // Destructor deletes dynamically created arrays for data members: fArgument,
133 // fFunction and fSecondDerivative, all have d 133 // fFunction and fSecondDerivative, all have dimension of fNumber
134 134
135 G4DataInterpolation::~G4DataInterpolation() 135 G4DataInterpolation::~G4DataInterpolation()
136 { 136 {
137 delete[] fArgument; 137 delete[] fArgument;
138 delete[] fFunction; 138 delete[] fFunction;
139 delete[] fSecondDerivative; 139 delete[] fSecondDerivative;
140 } 140 }
141 141
142 ////////////////////////////////////////////// 142 /////////////////////////////////////////////////////////////////////////////
143 // 143 //
144 // This function returns the value P(pX), wher 144 // This function returns the value P(pX), where P(x) is polynom of fNumber-1
145 // degree such that P(fArgument[i]) = fFunctio 145 // degree such that P(fArgument[i]) = fFunction[i], for i = 0, ..., fNumber-1.
146 // This is Lagrange's form of interpolation an 146 // This is Lagrange's form of interpolation and it is based on Neville's
147 // algorithm 147 // algorithm
148 148
149 G4double G4DataInterpolation::PolynomInterpola 149 G4double G4DataInterpolation::PolynomInterpolation(G4double pX,
150 150 G4double& deltaY) const
151 { 151 {
152 G4int i = 0, j = 1, k = 0; 152 G4int i = 0, j = 1, k = 0;
153 G4double mult = 0.0, difi = 0.0, deltaLow = 153 G4double mult = 0.0, difi = 0.0, deltaLow = 0.0, deltaUp = 0.0, cd = 0.0,
154 y = 0.0; 154 y = 0.0;
155 auto c = new G4double[fNumber]; << 155 G4double* c = new G4double[fNumber];
156 auto d = new G4double[fNumber]; << 156 G4double* d = new G4double[fNumber];
157 G4double diff = std::fabs(pX - fArgument[0]) 157 G4double diff = std::fabs(pX - fArgument[0]);
158 for(i = 0; i < fNumber; ++i) 158 for(i = 0; i < fNumber; ++i)
159 { 159 {
160 difi = std::fabs(pX - fArgument[i]); 160 difi = std::fabs(pX - fArgument[i]);
161 if(difi < diff) 161 if(difi < diff)
162 { 162 {
163 k = i; 163 k = i;
164 diff = difi; 164 diff = difi;
165 } 165 }
166 c[i] = fFunction[i]; 166 c[i] = fFunction[i];
167 d[i] = fFunction[i]; 167 d[i] = fFunction[i];
168 } 168 }
169 y = fFunction[k--]; 169 y = fFunction[k--];
170 for(j = 1; j < fNumber; ++j) 170 for(j = 1; j < fNumber; ++j)
171 { 171 {
172 for(i = 0; i < fNumber - j; ++i) 172 for(i = 0; i < fNumber - j; ++i)
173 { 173 {
174 deltaLow = fArgument[i] - pX; 174 deltaLow = fArgument[i] - pX;
175 deltaUp = fArgument[i + j] - pX; 175 deltaUp = fArgument[i + j] - pX;
176 cd = c[i + 1] - d[i]; 176 cd = c[i + 1] - d[i];
177 mult = deltaLow - deltaUp; 177 mult = deltaLow - deltaUp;
178 if(!(mult != 0.0)) 178 if(!(mult != 0.0))
179 { 179 {
180 G4Exception("G4DataInterpolation::Poly 180 G4Exception("G4DataInterpolation::PolynomInterpolation()", "Error",
181 FatalException, "Coinciden 181 FatalException, "Coincident nodes !");
182 } 182 }
183 mult = cd / mult; 183 mult = cd / mult;
184 d[i] = deltaUp * mult; 184 d[i] = deltaUp * mult;
185 c[i] = deltaLow * mult; 185 c[i] = deltaLow * mult;
186 } 186 }
187 y += (deltaY = (2 * k < (fNumber - j - 1) 187 y += (deltaY = (2 * k < (fNumber - j - 1) ? c[k + 1] : d[k--]));
188 } 188 }
189 delete[] c; 189 delete[] c;
190 delete[] d; 190 delete[] d;
191 191
192 return y; 192 return y;
193 } 193 }
194 194
195 ////////////////////////////////////////////// 195 ////////////////////////////////////////////////////////////////////////////
196 // 196 //
197 // Given arrays fArgument[0,..,fNumber-1] and 197 // Given arrays fArgument[0,..,fNumber-1] and fFunction[0,..,fNumber-1], this
198 // function calculates an array of coefficient 198 // function calculates an array of coefficients. The coefficients don't provide
199 // usually (fNumber>10) better accuracy for po 199 // usually (fNumber>10) better accuracy for polynom interpolation, as compared
200 // with PolynomInterpolation function. They co 200 // with PolynomInterpolation function. They could be used instead for derivate
201 // calculations and some other applications. 201 // calculations and some other applications.
202 202
203 void G4DataInterpolation::PolIntCoefficient(G4 203 void G4DataInterpolation::PolIntCoefficient(G4double cof[]) const
204 { 204 {
205 G4int i = 0, j = 0; 205 G4int i = 0, j = 0;
206 G4double factor; 206 G4double factor;
207 G4double reducedY = 0.0, mult = 1.0; 207 G4double reducedY = 0.0, mult = 1.0;
208 auto tempArgument = new G4double[fNumber]; << 208 G4double* tempArgument = new G4double[fNumber];
209 209
210 for(i = 0; i < fNumber; ++i) 210 for(i = 0; i < fNumber; ++i)
211 { 211 {
212 tempArgument[i] = cof[i] = 0.0; 212 tempArgument[i] = cof[i] = 0.0;
213 } 213 }
214 tempArgument[fNumber - 1] = -fArgument[0]; 214 tempArgument[fNumber - 1] = -fArgument[0];
215 215
216 for(i = 1; i < fNumber; ++i) 216 for(i = 1; i < fNumber; ++i)
217 { 217 {
218 for(j = fNumber - 1 - i; j < fNumber - 1; 218 for(j = fNumber - 1 - i; j < fNumber - 1; ++j)
219 { 219 {
220 tempArgument[j] -= fArgument[i] * tempAr 220 tempArgument[j] -= fArgument[i] * tempArgument[j + 1];
221 } 221 }
222 tempArgument[fNumber - 1] -= fArgument[i]; 222 tempArgument[fNumber - 1] -= fArgument[i];
223 } 223 }
224 for(i = 0; i < fNumber; ++i) 224 for(i = 0; i < fNumber; ++i)
225 { 225 {
226 factor = fNumber; 226 factor = fNumber;
227 for(j = fNumber - 1; j >= 1; --j) 227 for(j = fNumber - 1; j >= 1; --j)
228 { 228 {
229 factor = j * tempArgument[j] + factor * 229 factor = j * tempArgument[j] + factor * fArgument[i];
230 } 230 }
231 reducedY = fFunction[i] / factor; 231 reducedY = fFunction[i] / factor;
232 mult = 1.0; 232 mult = 1.0;
233 for(j = fNumber - 1; j >= 0; --j) 233 for(j = fNumber - 1; j >= 0; --j)
234 { 234 {
235 cof[j] += mult * reducedY; 235 cof[j] += mult * reducedY;
236 mult = tempArgument[j] + mult * fArgumen 236 mult = tempArgument[j] + mult * fArgument[i];
237 } 237 }
238 } 238 }
239 delete[] tempArgument; 239 delete[] tempArgument;
240 } 240 }
241 241
242 ////////////////////////////////////////////// 242 /////////////////////////////////////////////////////////////////////////////
243 // 243 //
244 // The function returns diagonal rational func 244 // The function returns diagonal rational function (Bulirsch and Stoer
245 // algorithm of Neville type) Pn(x)/Qm(x) wher 245 // algorithm of Neville type) Pn(x)/Qm(x) where P and Q are polynoms.
246 // Tests showed the method is not stable and h 246 // Tests showed the method is not stable and hasn't advantage if compared
247 // with polynomial interpolation ?! 247 // with polynomial interpolation ?!
248 248
249 G4double G4DataInterpolation::RationalPolInter 249 G4double G4DataInterpolation::RationalPolInterpolation(G4double pX,
250 250 G4double& deltaY) const
251 { 251 {
252 G4int i = 0, j = 1, k = 0; 252 G4int i = 0, j = 1, k = 0;
253 const G4double tolerance = 1.6e-24; 253 const G4double tolerance = 1.6e-24;
254 G4double mult = 0.0, difi = 0.0, cd = 0.0, y 254 G4double mult = 0.0, difi = 0.0, cd = 0.0, y = 0.0, cof = 0.0;
255 auto c = new G4double[fNumber]; << 255 G4double* c = new G4double[fNumber];
256 auto d = new G4double[fNumber]; << 256 G4double* d = new G4double[fNumber];
257 G4double diff = std::fabs(pX - fArgument[0]) 257 G4double diff = std::fabs(pX - fArgument[0]);
258 for(i = 0; i < fNumber; ++i) 258 for(i = 0; i < fNumber; ++i)
259 { 259 {
260 difi = std::fabs(pX - fArgument[i]); 260 difi = std::fabs(pX - fArgument[i]);
261 if(!(difi != 0.0)) 261 if(!(difi != 0.0))
262 { 262 {
263 y = fFunction[i]; 263 y = fFunction[i];
264 deltaY = 0.0; 264 deltaY = 0.0;
265 delete[] c; 265 delete[] c;
266 delete[] d; 266 delete[] d;
267 return y; 267 return y;
268 } 268 }
269 if(difi < diff) << 269 else if(difi < diff)
270 { 270 {
271 k = i; 271 k = i;
272 diff = difi; 272 diff = difi;
273 } 273 }
274 c[i] = fFunction[i]; 274 c[i] = fFunction[i];
275 d[i] = fFunction[i] + tolerance; // to pr 275 d[i] = fFunction[i] + tolerance; // to prevent rare zero/zero cases
276 } 276 }
277 y = fFunction[k--]; 277 y = fFunction[k--];
278 for(j = 1; j < fNumber; ++j) 278 for(j = 1; j < fNumber; ++j)
279 { 279 {
280 for(i = 0; i < fNumber - j; ++i) 280 for(i = 0; i < fNumber - j; ++i)
281 { 281 {
282 cd = c[i + 1] - d[i]; 282 cd = c[i + 1] - d[i];
283 difi = fArgument[i + j] - pX; 283 difi = fArgument[i + j] - pX;
284 cof = (fArgument[i] - pX) * d[i] / difi 284 cof = (fArgument[i] - pX) * d[i] / difi;
285 mult = cof - c[i + 1]; 285 mult = cof - c[i + 1];
286 if(!(mult != 0.0)) // function to be in 286 if(!(mult != 0.0)) // function to be interpolated has pole at pX
287 { 287 {
288 G4Exception("G4DataInterpolation::Rati 288 G4Exception("G4DataInterpolation::RationalPolInterpolation()", "Error",
289 FatalException, "Coinciden 289 FatalException, "Coincident nodes !");
290 } 290 }
291 mult = cd / mult; 291 mult = cd / mult;
292 d[i] = c[i + 1] * mult; 292 d[i] = c[i + 1] * mult;
293 c[i] = cof * mult; 293 c[i] = cof * mult;
294 } 294 }
295 y += (deltaY = (2 * k < (fNumber - j - 1) 295 y += (deltaY = (2 * k < (fNumber - j - 1) ? c[k + 1] : d[k--]));
296 } 296 }
297 delete[] c; 297 delete[] c;
298 delete[] d; 298 delete[] d;
299 299
300 return y; 300 return y;
301 } 301 }
302 302
303 ////////////////////////////////////////////// 303 /////////////////////////////////////////////////////////////////////////////
304 // 304 //
305 // Cubic spline interpolation in point pX for 305 // Cubic spline interpolation in point pX for function given by the table:
306 // fArgument, fFunction. The constructor, whic 306 // fArgument, fFunction. The constructor, which creates fSecondDerivative,
307 // must be called before. The function works o 307 // must be called before. The function works optimal, if sequential calls
308 // are in random values of pX. 308 // are in random values of pX.
309 309
310 G4double G4DataInterpolation::CubicSplineInter 310 G4double G4DataInterpolation::CubicSplineInterpolation(G4double pX) const
311 { 311 {
312 G4int kLow = 0, kHigh = fNumber - 1, k = 0; 312 G4int kLow = 0, kHigh = fNumber - 1, k = 0;
313 313
314 // Searching in the table by means of bisect 314 // Searching in the table by means of bisection method.
315 // fArgument must be monotonic, either incre 315 // fArgument must be monotonic, either increasing or decreasing
316 316
317 while((kHigh - kLow) > 1) 317 while((kHigh - kLow) > 1)
318 { 318 {
319 k = (kHigh + kLow) >> 1; // compute midpo 319 k = (kHigh + kLow) >> 1; // compute midpoint 'bisection'
320 if(fArgument[k] > pX) 320 if(fArgument[k] > pX)
321 { 321 {
322 kHigh = k; 322 kHigh = k;
323 } 323 }
324 else 324 else
325 { 325 {
326 kLow = k; 326 kLow = k;
327 } 327 }
328 } // kLow and kHigh now bracket the input v 328 } // kLow and kHigh now bracket the input value of pX
329 G4double deltaHL = fArgument[kHigh] - fArgum 329 G4double deltaHL = fArgument[kHigh] - fArgument[kLow];
330 if(!(deltaHL != 0.0)) 330 if(!(deltaHL != 0.0))
331 { 331 {
332 G4Exception("G4DataInterpolation::CubicSpl 332 G4Exception("G4DataInterpolation::CubicSplineInterpolation()", "Error",
333 FatalException, "Bad fArgument 333 FatalException, "Bad fArgument input !");
334 } 334 }
335 G4double a = (fArgument[kHigh] - pX) / delta 335 G4double a = (fArgument[kHigh] - pX) / deltaHL;
336 G4double b = (pX - fArgument[kLow]) / deltaH 336 G4double b = (pX - fArgument[kLow]) / deltaHL;
337 337
338 // Final evaluation of cubic spline polynomi 338 // Final evaluation of cubic spline polynomial for return
339 339
340 return a * fFunction[kLow] + b * fFunction[k 340 return a * fFunction[kLow] + b * fFunction[kHigh] +
341 ((a * a * a - a) * fSecondDerivative[ 341 ((a * a * a - a) * fSecondDerivative[kLow] +
342 (b * b * b - b) * fSecondDerivative[ 342 (b * b * b - b) * fSecondDerivative[kHigh]) *
343 deltaHL * deltaHL / 6.0; 343 deltaHL * deltaHL / 6.0;
344 } 344 }
345 345
346 ////////////////////////////////////////////// 346 ///////////////////////////////////////////////////////////////////////////
347 // 347 //
348 // Return cubic spline interpolation in the po 348 // Return cubic spline interpolation in the point pX which is located between
349 // fArgument[index] and fArgument[index+1]. It 349 // fArgument[index] and fArgument[index+1]. It is usually called in sequence
350 // of known from external analysis values of i 350 // of known from external analysis values of index.
351 351
352 G4double G4DataInterpolation::FastCubicSpline( 352 G4double G4DataInterpolation::FastCubicSpline(G4double pX, G4int index) const
353 { 353 {
354 G4double delta = fArgument[index + 1] - fArg 354 G4double delta = fArgument[index + 1] - fArgument[index];
355 if(!(delta != 0.0)) 355 if(!(delta != 0.0))
356 { 356 {
357 G4Exception("G4DataInterpolation::FastCubi 357 G4Exception("G4DataInterpolation::FastCubicSpline()", "Error",
358 FatalException, "Bad fArgument 358 FatalException, "Bad fArgument input !");
359 } 359 }
360 G4double a = (fArgument[index + 1] - pX) / d 360 G4double a = (fArgument[index + 1] - pX) / delta;
361 G4double b = (pX - fArgument[index]) / delta 361 G4double b = (pX - fArgument[index]) / delta;
362 362
363 // Final evaluation of cubic spline polynomi 363 // Final evaluation of cubic spline polynomial for return
364 364
365 return a * fFunction[index] + b * fFunction[ 365 return a * fFunction[index] + b * fFunction[index + 1] +
366 ((a * a * a - a) * fSecondDerivative[ 366 ((a * a * a - a) * fSecondDerivative[index] +
367 (b * b * b - b) * fSecondDerivative[ 367 (b * b * b - b) * fSecondDerivative[index + 1]) *
368 delta * delta / 6.0; 368 delta * delta / 6.0;
369 } 369 }
370 370
371 ////////////////////////////////////////////// 371 ////////////////////////////////////////////////////////////////////////////
372 // 372 //
373 // Given argument pX, returns index k, so that 373 // Given argument pX, returns index k, so that pX bracketed by fArgument[k]
374 // and fArgument[k+1] 374 // and fArgument[k+1]
375 375
376 G4int G4DataInterpolation::LocateArgument(G4do 376 G4int G4DataInterpolation::LocateArgument(G4double pX) const
377 { 377 {
378 G4int kLow = -1, kHigh = fNumber, k = 0; 378 G4int kLow = -1, kHigh = fNumber, k = 0;
379 G4bool ascend = (fArgument[fNumber - 1] >= f 379 G4bool ascend = (fArgument[fNumber - 1] >= fArgument[0]);
380 while((kHigh - kLow) > 1) 380 while((kHigh - kLow) > 1)
381 { 381 {
382 k = (kHigh + kLow) >> 1; // compute midpo 382 k = (kHigh + kLow) >> 1; // compute midpoint 'bisection'
383 if((pX >= fArgument[k]) == ascend) 383 if((pX >= fArgument[k]) == ascend)
384 { 384 {
385 kLow = k; 385 kLow = k;
386 } 386 }
387 else 387 else
388 { 388 {
389 kHigh = k; 389 kHigh = k;
390 } 390 }
391 } 391 }
392 if(!(pX != fArgument[0])) 392 if(!(pX != fArgument[0]))
393 { 393 {
394 return 1; 394 return 1;
395 } 395 }
396 if(!(pX != fArgument[fNumber - 1])) << 396 else if(!(pX != fArgument[fNumber - 1]))
397 { 397 {
398 return fNumber - 2; 398 return fNumber - 2;
399 } 399 }
400 << 400 else
401 return kLow; << 401 return kLow;
402 <<
403 } 402 }
404 403
405 ////////////////////////////////////////////// 404 /////////////////////////////////////////////////////////////////////////////
406 // 405 //
407 // Given a value pX, returns a value 'index' s 406 // Given a value pX, returns a value 'index' such that pX is between
408 // fArgument[index] and fArgument[index+1]. fA 407 // fArgument[index] and fArgument[index+1]. fArgument MUST BE MONOTONIC,
409 // either increasing or decreasing. If index = 408 // either increasing or decreasing. If index = -1 or fNumber, this indicates
410 // that pX is out of range. The value index on 409 // that pX is out of range. The value index on input is taken as the initial
411 // approximation for index on output. 410 // approximation for index on output.
412 411
413 void G4DataInterpolation::CorrelatedSearch(G4d 412 void G4DataInterpolation::CorrelatedSearch(G4double pX, G4int& index) const
414 { 413 {
415 G4int kHigh = 0, k = 0, Increment = 0; 414 G4int kHigh = 0, k = 0, Increment = 0;
416 // ascend = true for ascending order of tabl 415 // ascend = true for ascending order of table, false otherwise
417 G4bool ascend = (fArgument[fNumber - 1] >= f 416 G4bool ascend = (fArgument[fNumber - 1] >= fArgument[0]);
418 if(index < 0 || index > fNumber - 1) 417 if(index < 0 || index > fNumber - 1)
419 { 418 {
420 index = -1; 419 index = -1;
421 kHigh = fNumber; 420 kHigh = fNumber;
422 } 421 }
423 else 422 else
424 { 423 {
425 Increment = 1; // What value would be t 424 Increment = 1; // What value would be the best ?
426 if((pX >= fArgument[index]) == ascend) 425 if((pX >= fArgument[index]) == ascend)
427 { 426 {
428 if(index == fNumber - 1) 427 if(index == fNumber - 1)
429 { 428 {
430 index = fNumber; 429 index = fNumber;
431 return; 430 return;
432 } 431 }
433 kHigh = index + 1; 432 kHigh = index + 1;
434 while((pX >= fArgument[kHigh]) == ascend 433 while((pX >= fArgument[kHigh]) == ascend)
435 { 434 {
436 index = kHigh; 435 index = kHigh;
437 Increment += Increment; // double the 436 Increment += Increment; // double the Increment
438 kHigh = index + Increment; 437 kHigh = index + Increment;
439 if(kHigh > (fNumber - 1)) 438 if(kHigh > (fNumber - 1))
440 { 439 {
441 kHigh = fNumber; 440 kHigh = fNumber;
442 break; 441 break;
443 } 442 }
444 } 443 }
445 } 444 }
446 else 445 else
447 { 446 {
448 if(index == 0) 447 if(index == 0)
449 { 448 {
450 index = -1; 449 index = -1;
451 return; 450 return;
452 } 451 }
453 kHigh = --index; 452 kHigh = --index;
454 while((pX < fArgument[index]) == ascend) 453 while((pX < fArgument[index]) == ascend)
455 { 454 {
456 kHigh = index; 455 kHigh = index;
457 Increment <<= 1; // double the Increm 456 Increment <<= 1; // double the Increment
458 if(Increment >= kHigh) 457 if(Increment >= kHigh)
459 { 458 {
460 index = -1; 459 index = -1;
461 break; 460 break;
462 } 461 }
463 << 462 else
464 index = kHigh - Increment; << 463 {
465 << 464 index = kHigh - Increment;
>> 465 }
466 } 466 }
467 } // Value bracketed 467 } // Value bracketed
468 } 468 }
469 // final bisection searching 469 // final bisection searching
470 470
471 while((kHigh - index) != 1) 471 while((kHigh - index) != 1)
472 { 472 {
473 k = (kHigh + index) >> 1; 473 k = (kHigh + index) >> 1;
474 if((pX >= fArgument[k]) == ascend) 474 if((pX >= fArgument[k]) == ascend)
475 { 475 {
476 index = k; 476 index = k;
477 } 477 }
478 else 478 else
479 { 479 {
480 kHigh = k; 480 kHigh = k;
481 } 481 }
482 } 482 }
483 if(!(pX != fArgument[fNumber - 1])) 483 if(!(pX != fArgument[fNumber - 1]))
484 { 484 {
485 index = fNumber - 2; 485 index = fNumber - 2;
486 } 486 }
487 if(!(pX != fArgument[0])) 487 if(!(pX != fArgument[0]))
488 { 488 {
489 index = 0; 489 index = 0;
490 } 490 }
491 return; 491 return;
492 } 492 }
493 493
494 // 494 //
495 // 495 //
496 ////////////////////////////////////////////// 496 ////////////////////////////////////////////////////////////////////////////
497 497