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