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