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