| Geant4 Cross Reference |
>> 1 // This code implementation is the intellectual property of
>> 2 // the GEANT4 collaboration.
1 // 3 //
2 // ******************************************* << 4 // By copying, distributing or modifying the Program (or any work
3 // * License and Disclaimer << 5 // based on the Program) you indicate your acceptance of this statement,
4 // * << 6 // and all its terms.
5 // * The Geant4 software is copyright of th << 7 //
6 // * the Geant4 Collaboration. It is provided << 8 // $Id: G4Integrator.icc,v 1.4 1999/11/23 14:59:59 gcosmo Exp $
7 // * conditions of the Geant4 Software License << 9 // GEANT4 tag $Name: geant4-01-01 $
8 // * LICENSE and available at http://cern.ch/ << 10 //
9 // * include a list of copyright holders. << 11 // Implementation of G4Integrator methods.
10 // * << 12 //
11 // * Neither the authors of this software syst << 13 //
12 // * institutes,nor the agencies providing fin <<
13 // * work make any representation or warran <<
14 // * regarding this software system or assum <<
15 // * use. Please see the license in the file <<
16 // * for the full disclaimer and the limitatio <<
17 // * <<
18 // * This code implementation is the result <<
19 // * technical work of the GEANT4 collaboratio <<
20 // * By using, copying, modifying or distri <<
21 // * any work based on the software) you ag <<
22 // * use in resulting scientific publicati <<
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* <<
25 // <<
26 // G4Integrator inline methods implementation <<
27 // <<
28 // Author: V.Grichine, 04.09.1999 - First impl <<
29 // G4SimpleIntegration class with H.P. <<
30 // E.TCherniaev advises <<
31 // ------------------------------------------- <<
32 14
33 ////////////////////////////////////////////// 15 /////////////////////////////////////////////////////////////////////
34 // 16 //
35 // Sympson integration method 17 // Sympson integration method
36 // 18 //
37 ////////////////////////////////////////////// 19 /////////////////////////////////////////////////////////////////////
38 // 20 //
39 // Integration of class member functions T::f << 21 // Integration of class member functions T::f by Simpson method.
40 22
41 template <class T, class F> << 23 template <class T, class F>
42 G4double G4Integrator<T, F>::Simpson(T& typeT, << 24 G4double G4Integrator::Simpson( T& typeT,
43 G4double << 25 F f,
44 { << 26 G4double xInitial,
45 G4int i; << 27 G4double xFinal,
46 G4double step = (xFinal - xInitial) / itera << 28 G4int iterationNumber )
47 G4double x = xInitial; << 29 {
48 G4double xPlus = xInitial + 0.5 * step; << 30 G4int i ;
49 G4double mean = ((typeT.*f)(xInitial) + (ty << 31 G4double step = (xFinal - xInitial)/iterationNumber ;
50 G4double sum = (typeT.*f)(xPlus); << 32 G4double x = xInitial ;
51 << 33 G4double xPlus = xInitial + 0.5*step ;
52 for(i = 1; i < iterationNumber; ++i) << 34 G4double mean = ( (typeT.*f)(xInitial) + (typeT.*f)(xFinal) )*0.5 ;
53 { << 35 G4double sum = (typeT.*f)(xPlus) ;
54 x += step; << 36
55 xPlus += step; << 37 for(i=1;i<iterationNumber;i++)
56 mean += (typeT.*f)(x); << 38 {
57 sum += (typeT.*f)(xPlus); << 39 x += step ;
58 } << 40 xPlus += step ;
59 mean += 2.0 * sum; << 41 mean += (typeT.*f)(x) ;
>> 42 sum += (typeT.*f)(xPlus) ;
>> 43 }
>> 44 mean += 2.0*sum ;
60 45
61 return mean * step / 3.0; << 46 return mean*step/3.0 ;
62 } 47 }
63 48
64 ////////////////////////////////////////////// 49 /////////////////////////////////////////////////////////////////////
65 // 50 //
66 // Integration of class member functions T::f 51 // Integration of class member functions T::f by Simpson method.
67 // Convenient to use with 'this' pointer 52 // Convenient to use with 'this' pointer
68 53
69 template <class T, class F> << 54 template <class T, class F>
70 G4double G4Integrator<T, F>::Simpson(T* ptrT, << 55 G4double G4Integrator::Simpson( T* ptrT,
71 G4double << 56 F f,
72 { << 57 G4double xInitial,
73 G4int i; << 58 G4double xFinal,
74 G4double step = (xFinal - xInitial) / itera << 59 G4int iterationNumber )
75 G4double x = xInitial; << 60 {
76 G4double xPlus = xInitial + 0.5 * step; << 61 G4int i ;
77 G4double mean = ((ptrT->*f)(xInitial) + (pt << 62 G4double step = (xFinal - xInitial)/iterationNumber ;
78 G4double sum = (ptrT->*f)(xPlus); << 63 G4double x = xInitial ;
>> 64 G4double xPlus = xInitial + 0.5*step ;
>> 65 G4double mean = ( (ptrT->*f)(xInitial) + (ptrT->*f)(xFinal) )*0.5 ;
>> 66 G4double sum = (ptrT->*f)(xPlus) ;
>> 67
>> 68 for(i=1;i<iterationNumber;i++)
>> 69 {
>> 70 x += step ;
>> 71 xPlus += step ;
>> 72 mean += (ptrT->*f)(x) ;
>> 73 sum += (ptrT->*f)(xPlus) ;
>> 74 }
>> 75 mean += 2.0*sum ;
79 76
80 for(i = 1; i < iterationNumber; ++i) << 77 return mean*step/3.0 ;
81 { <<
82 x += step; <<
83 xPlus += step; <<
84 mean += (ptrT->*f)(x); <<
85 sum += (ptrT->*f)(xPlus); <<
86 } <<
87 mean += 2.0 * sum; <<
88 <<
89 return mean * step / 3.0; <<
90 } 78 }
91 79
92 ////////////////////////////////////////////// 80 /////////////////////////////////////////////////////////////////////
93 // 81 //
94 // Integration of class member functions T::f 82 // Integration of class member functions T::f by Simpson method.
95 // Convenient to use, when function f is defin 83 // Convenient to use, when function f is defined in global scope, i.e. in main()
96 // program 84 // program
97 85
98 template <class T, class F> << 86 G4double G4Integrator::Simpson( G4double (*f)(G4double),
99 G4double G4Integrator<T, F>::Simpson(G4double << 87 G4double xInitial,
100 G4double << 88 G4double xFinal,
101 { << 89 G4int iterationNumber )
102 G4int i; << 90 {
103 G4double step = (xFinal - xInitial) / itera << 91 G4int i ;
104 G4double x = xInitial; << 92 G4double step = (xFinal - xInitial)/iterationNumber ;
105 G4double xPlus = xInitial + 0.5 * step; << 93 G4double x = xInitial ;
106 G4double mean = ((*f)(xInitial) + (*f)(xFin << 94 G4double xPlus = xInitial + 0.5*step ;
107 G4double sum = (*f)(xPlus); << 95 G4double mean = ( (*f)(xInitial) + (*f)(xFinal) )*0.5 ;
108 << 96 G4double sum = (*f)(xPlus) ;
109 for(i = 1; i < iterationNumber; ++i) << 97
110 { << 98 for(i=1;i<iterationNumber;i++)
111 x += step; << 99 {
112 xPlus += step; << 100 x += step ;
113 mean += (*f)(x); << 101 xPlus += step ;
114 sum += (*f)(xPlus); << 102 mean += (*f)(x) ;
115 } << 103 sum += (*f)(xPlus) ;
116 mean += 2.0 * sum; << 104 }
>> 105 mean += 2.0*sum ;
117 106
118 return mean * step / 3.0; << 107 return mean*step/3.0 ;
119 } 108 }
120 109
121 ////////////////////////////////////////////// 110 //////////////////////////////////////////////////////////////////////////
122 // 111 //
123 // Adaptive Gauss method 112 // Adaptive Gauss method
124 // 113 //
125 ////////////////////////////////////////////// 114 //////////////////////////////////////////////////////////////////////////
126 // 115 //
127 // 116 //
128 117
129 template <class T, class F> << 118 template <class T, class F>
130 G4double G4Integrator<T, F>::Gauss(T& typeT, F << 119 G4double G4Integrator::Gauss( T& typeT, F f,
131 G4double xF << 120 G4double xInitial, G4double xFinal )
132 { << 121 {
133 static const G4double root = 1.0 / std::sqrt << 122 static G4double root = 1.0/sqrt(3.0) ;
134 << 123
135 G4double xMean = (xInitial + xFinal) / 2.0; << 124 G4double xMean = (xInitial + xFinal)/2.0 ;
136 G4double Step = (xFinal - xInitial) / 2.0; << 125 G4double Step = (xFinal - xInitial)/2.0 ;
137 G4double delta = Step * root; << 126 G4double delta = Step*root ;
138 G4double sum = ((typeT.*f)(xMean + delta) << 127 G4double sum = ((typeT.*f)(xMean + delta) +
139 << 128 (typeT.*f)(xMean - delta)) ;
140 return sum * Step; << 129
>> 130 return sum*Step ;
141 } 131 }
142 132
143 ////////////////////////////////////////////// 133 //////////////////////////////////////////////////////////////////////
144 // 134 //
145 // 135 //
146 136
147 template <class T, class F> << 137 template <class T, class F> G4double
148 G4double G4Integrator<T, F>::Gauss(T* ptrT, F << 138 G4Integrator::Gauss( T* ptrT, F f, G4double a, G4double b )
149 { 139 {
150 return Gauss(*ptrT, f, a, b); << 140 return Gauss(*ptrT,f,a,b) ;
151 } 141 }
152 142
153 ////////////////////////////////////////////// 143 ///////////////////////////////////////////////////////////////////////
154 // 144 //
155 // 145 //
156 146
157 template <class T, class F> << 147 G4double G4Integrator::Gauss( G4double (*f)(G4double),
158 G4double G4Integrator<T, F>::Gauss(G4double (* << 148 G4double xInitial, G4double xFinal)
159 G4double xF << 149 {
160 { << 150 static G4double root = 1.0/sqrt(3.0) ;
161 static const G4double root = 1.0 / std::sqrt << 151
162 << 152 G4double xMean = (xInitial + xFinal)/2.0 ;
163 G4double xMean = (xInitial + xFinal) / 2.0; << 153 G4double Step = (xFinal - xInitial)/2.0 ;
164 G4double Step = (xFinal - xInitial) / 2.0; << 154 G4double delta = Step*root ;
165 G4double delta = Step * root; << 155 G4double sum = ( (*f)(xMean + delta) + (*f)(xMean - delta) ) ;
166 G4double sum = ((*f)(xMean + delta) + (*f) << 156
167 << 157 return sum*Step ;
168 return sum * Step; <<
169 } 158 }
170 159
171 ////////////////////////////////////////////// 160 ///////////////////////////////////////////////////////////////////////////
172 // 161 //
173 // 162 //
174 163
175 template <class T, class F> << 164 template <class T, class F>
176 void G4Integrator<T, F>::AdaptGauss(T& typeT, << 165 void G4Integrator::AdaptGauss( T& typeT, F f, G4double xInitial,
177 G4double x << 166 G4double xFinal, G4double fTolerance,
178 G4double& << 167 G4double& sum,
>> 168 G4int& depth )
>> 169 {
>> 170 if(depth > 100)
>> 171 {
>> 172 G4cout<<"G4Integrator::AdaptGauss: WARNING !!!"<<G4endl ;
>> 173 G4cout
>> 174 <<"Function varies too rapidly to get stated accuracy in 100 steps "<<G4endl ;
>> 175
>> 176 return ;
>> 177 }
>> 178 G4double xMean = (xInitial + xFinal)/2.0 ;
>> 179 G4double leftHalf = Gauss(typeT,f,xInitial,xMean) ;
>> 180 G4double rightHalf = Gauss(typeT,f,xMean,xFinal) ;
>> 181 G4double full = Gauss(typeT,f,xInitial,xFinal) ;
>> 182 if(fabs(leftHalf+rightHalf-full) < fTolerance)
>> 183 {
>> 184 sum += full ;
>> 185 }
>> 186 else
>> 187 {
>> 188 depth++ ;
>> 189 AdaptGauss(typeT,f,xInitial,xMean,fTolerance,sum,depth) ;
>> 190 AdaptGauss(typeT,f,xMean,xFinal,fTolerance,sum,depth) ;
>> 191 }
>> 192 }
>> 193
>> 194 template <class T, class F>
>> 195 void G4Integrator::AdaptGauss( T* ptrT, F f, G4double xInitial,
>> 196 G4double xFinal, G4double fTolerance,
>> 197 G4double& sum,
>> 198 G4int& depth )
179 { 199 {
180 if(depth > 100) << 200 return AdaptGauss(*ptrT,f,xInitial,xFinal,fTolerance,sum,depth) ;
181 { <<
182 G4cout << "G4Integrator<T,F>::AdaptGauss: <<
183 G4cout << "Function varies too rapidly to <<
184 << G4endl; <<
185 <<
186 return; <<
187 } <<
188 G4double xMean = (xInitial + xFinal) / 2 <<
189 G4double leftHalf = Gauss(typeT, f, xInitia <<
190 G4double rightHalf = Gauss(typeT, f, xMean, <<
191 G4double full = Gauss(typeT, f, xInitia <<
192 if(std::fabs(leftHalf + rightHalf - full) < <<
193 { <<
194 sum += full; <<
195 } <<
196 else <<
197 { <<
198 ++depth; <<
199 AdaptGauss(typeT, f, xInitial, xMean, fTol <<
200 AdaptGauss(typeT, f, xMean, xFinal, fToler <<
201 } <<
202 } <<
203 <<
204 template <class T, class F> <<
205 void G4Integrator<T, F>::AdaptGauss(T* ptrT, F <<
206 G4double x <<
207 G4double& <<
208 { <<
209 AdaptGauss(*ptrT, f, xInitial, xFinal, fTole <<
210 } 201 }
211 202
212 ////////////////////////////////////////////// 203 /////////////////////////////////////////////////////////////////////////
213 // 204 //
214 // 205 //
215 template <class T, class F> <<
216 void G4Integrator<T, F>::AdaptGauss(G4double ( <<
217 G4double x <<
218 G4double& <<
219 { <<
220 if(depth > 100) <<
221 { <<
222 G4cout << "G4SimpleIntegration::AdaptGauss <<
223 G4cout << "Function varies too rapidly to <<
224 << G4endl; <<
225 206
226 return; << 207 void G4Integrator::AdaptGauss( G4double (*f)(G4double),
227 } << 208 G4double xInitial, G4double xFinal,
228 G4double xMean = (xInitial + xFinal) / 2 << 209 G4double fTolerance, G4double& sum,
229 G4double leftHalf = Gauss(f, xInitial, xMea << 210 G4int& depth )
230 G4double rightHalf = Gauss(f, xMean, xFinal) << 211 {
231 G4double full = Gauss(f, xInitial, xFin << 212 if(depth > 100)
232 if(std::fabs(leftHalf + rightHalf - full) < << 213 {
233 { << 214 G4cout<<"G4SimpleIntegration::AdaptGauss: WARNING !!!"<<G4endl ;
234 sum += full; << 215 G4cout<<"Function varies too rapidly to get stated accuracy in 100 steps "
235 } << 216 <<G4endl ;
236 else << 217
237 { << 218 return ;
238 ++depth; << 219 }
239 AdaptGauss(f, xInitial, xMean, fTolerance, << 220 G4double xMean = (xInitial + xFinal)/2.0 ;
240 AdaptGauss(f, xMean, xFinal, fTolerance, s << 221 G4double leftHalf = Gauss(f,xInitial,xMean) ;
241 } << 222 G4double rightHalf = Gauss(f,xMean,xFinal) ;
>> 223 G4double full = Gauss(f,xInitial,xFinal) ;
>> 224 if(fabs(leftHalf+rightHalf-full) < fTolerance)
>> 225 {
>> 226 sum += full ;
>> 227 }
>> 228 else
>> 229 {
>> 230 depth++ ;
>> 231 AdaptGauss(f,xInitial,xMean,fTolerance,sum,depth) ;
>> 232 AdaptGauss(f,xMean,xFinal,fTolerance,sum,depth) ;
>> 233 }
242 } 234 }
243 235
>> 236
>> 237
>> 238
244 ////////////////////////////////////////////// 239 ////////////////////////////////////////////////////////////////////////
245 // 240 //
246 // Adaptive Gauss integration with accuracy 'e 241 // Adaptive Gauss integration with accuracy 'e'
247 // Convenient for using with class object type 242 // Convenient for using with class object typeT
248 << 243
249 template <class T, class F> << 244 template<class T, class F> G4double
250 G4double G4Integrator<T, F>::AdaptiveGauss(T& << 245 G4Integrator::AdaptiveGauss( T& typeT, F f, G4double xInitial,
251 G4d << 246 G4double xFinal, G4double e )
252 { << 247 {
253 G4int depth = 0; << 248 G4int depth = 0 ;
254 G4double sum = 0.0; << 249 G4double sum = 0.0 ;
255 AdaptGauss(typeT, f, xInitial, xFinal, e, su << 250 AdaptGauss(typeT,f,xInitial,xFinal,e,sum,depth) ;
256 return sum; << 251 return sum ;
257 } 252 }
258 253
259 ////////////////////////////////////////////// 254 ////////////////////////////////////////////////////////////////////////
260 // 255 //
261 // Adaptive Gauss integration with accuracy 'e 256 // Adaptive Gauss integration with accuracy 'e'
262 // Convenient for using with 'this' pointer 257 // Convenient for using with 'this' pointer
263 << 258
264 template <class T, class F> << 259 template<class T, class F> G4double
265 G4double G4Integrator<T, F>::AdaptiveGauss(T* << 260 G4Integrator::AdaptiveGauss( T* ptrT, F f, G4double xInitial,
266 G4d << 261 G4double xFinal, G4double e )
267 { 262 {
268 return AdaptiveGauss(*ptrT, f, xInitial, xFi << 263 return AdaptiveGauss(*ptrT,f,xInitial,xFinal,e) ;
269 } 264 }
270 265
271 ////////////////////////////////////////////// 266 ////////////////////////////////////////////////////////////////////////
272 // 267 //
273 // Adaptive Gauss integration with accuracy 'e 268 // Adaptive Gauss integration with accuracy 'e'
274 // Convenient for using with global scope func 269 // Convenient for using with global scope function f
275 << 270
276 template <class T, class F> << 271 G4double
277 G4double G4Integrator<T, F>::AdaptiveGauss(G4d << 272 G4Integrator::AdaptiveGauss( G4double (*f)(G4double),
278 G4d << 273 G4double xInitial, G4double xFinal, G4double e )
279 G4d << 274 {
280 { << 275 G4int depth = 0 ;
281 G4int depth = 0; << 276 G4double sum = 0.0 ;
282 G4double sum = 0.0; << 277 AdaptGauss(f,xInitial,xFinal,e,sum,depth) ;
283 AdaptGauss(f, xInitial, xFinal, e, sum, dept << 278 return sum ;
284 return sum; <<
285 } 279 }
286 280
287 ////////////////////////////////////////////// 281 ////////////////////////////////////////////////////////////////////////////
288 // Gauss integration methods involving ortogon 282 // Gauss integration methods involving ortogonal polynomials
289 ////////////////////////////////////////////// 283 ////////////////////////////////////////////////////////////////////////////
290 // 284 //
291 // Methods involving Legendre polynomials << 285 // Methods involving Legendre polynomials
292 // 286 //
293 ////////////////////////////////////////////// 287 /////////////////////////////////////////////////////////////////////////
294 // 288 //
295 // The value nLegendre set the accuracy requir 289 // The value nLegendre set the accuracy required, i.e the number of points
296 // where the function pFunction will be evalua 290 // where the function pFunction will be evaluated during integration.
297 // The function creates the arrays for absciss << 291 // The function creates the arrays for abscissas and weights that used
298 // in Gauss-Legendre quadrature method. << 292 // in Gauss-Legendre quadrature method.
299 // The values a and b are the limits of integr 293 // The values a and b are the limits of integration of the function f .
300 // nLegendre MUST BE EVEN !!! 294 // nLegendre MUST BE EVEN !!!
301 // Returns the integral of the function f betw << 295 // Returns the integral of the function f between a and b, by 2*fNumber point
302 // Gauss-Legendre integration: the function is 296 // Gauss-Legendre integration: the function is evaluated exactly
303 // 2*fNumber times at interior points in the r << 297 // 2*fNumber times at interior points in the range of integration.
304 // Since the weights and abscissas are, in thi << 298 // Since the weights and abscissas are, in this case, symmetric around
305 // the midpoint of the range of integration, t << 299 // the midpoint of the range of integration, there are actually only
306 // fNumber distinct values of each. 300 // fNumber distinct values of each.
307 // Convenient for using with some class object 301 // Convenient for using with some class object dataT
308 302
309 template <class T, class F> << 303 template <class T, class F> G4double
310 G4double G4Integrator<T, F>::Legendre(T& typeT << 304 G4Integrator::Legendre( T& typeT, F f, G4double a, G4double b, G4int nLegendre)
311 G4int nL << 305 {
312 { << 306 G4double newton, newton1, temp1, temp2, temp3, temp ;
313 G4double nwt, nwt1, temp1, temp2, temp3, tem << 307 G4double xDiff, xMean, dx, integral ;
314 G4double xDiff, xMean, dx, integral; <<
315 <<
316 const G4double tolerance = 1.6e-10; <<
317 G4int i, j, k = nLegendre; <<
318 G4int fNumber = (nLegendre + 1) / 2; <<
319 <<
320 if(2 * fNumber != k) <<
321 { <<
322 G4Exception("G4Integrator<T,F>::Legendre(T <<
323 FatalException, "Invalid (odd) <<
324 } <<
325 <<
326 G4double* fAbscissa = new G4double[fNumber]; <<
327 G4double* fWeight = new G4double[fNumber]; <<
328 <<
329 for(i = 1; i <= fNumber; ++i) // Loop over <<
330 { <<
331 nwt = std::cos(CLHEP::pi * (i - 0.25) / <<
332 (k + 0.5)); // Initial roo <<
333 <<
334 do // loop of Newton's method <<
335 { <<
336 temp1 = 1.0; <<
337 temp2 = 0.0; <<
338 for(j = 1; j <= k; ++j) <<
339 { <<
340 temp3 = temp2; <<
341 temp2 = temp1; <<
342 temp1 = ((2.0 * j - 1.0) * nwt * temp2 <<
343 } <<
344 temp = k * (nwt * temp1 - temp2) / (nwt <<
345 nwt1 = nwt; <<
346 nwt = nwt1 - temp1 / temp; // Newton's <<
347 } while(std::fabs(nwt - nwt1) > tolerance) <<
348 <<
349 fAbscissa[fNumber - i] = nwt; <<
350 fWeight[fNumber - i] = 2.0 / ((1.0 - nwt <<
351 } <<
352 308
353 // << 309 const G4double tolerance = 1.6e-10 ;
354 // Now we ready to get integral << 310 G4int i, j, k = nLegendre ;
355 // << 311 G4int fNumber = (nLegendre + 1)/2 ;
356 << 312
357 xMean = 0.5 * (a + b); << 313 if(2*fNumber != k)
358 xDiff = 0.5 * (b - a); << 314 {
359 integral = 0.0; << 315 G4Exception("Invalid (odd) n Legendre in G4Integrator::Legendre") ;
360 for(i = 0; i < fNumber; ++i) << 316 }
361 { << 317
362 dx = xDiff * fAbscissa[i]; << 318 G4double* fAbscissa = new G4double[fNumber] ;
363 integral += fWeight[i] * ((typeT.*f)(xMean << 319 G4double* fWeight = new G4double[fNumber] ;
364 } << 320
365 delete[] fAbscissa; << 321 for(i=1;i<=fNumber;i++) // Loop over the desired roots
366 delete[] fWeight; << 322 {
367 return integral *= xDiff; << 323 newton = cos(pi*(i - 0.25)/(k + 0.5)) ; // Initial root approximation
368 } << 324
>> 325 do // loop of Newton's method
>> 326 {
>> 327 temp1 = 1.0 ;
>> 328 temp2 = 0.0 ;
>> 329 for(j=1;j<=k;j++)
>> 330 {
>> 331 temp3 = temp2 ;
>> 332 temp2 = temp1 ;
>> 333 temp1 = ((2.0*j - 1.0)*newton*temp2 - (j - 1.0)*temp3)/j ;
>> 334 }
>> 335 temp = k*(newton*temp1 - temp2)/(newton*newton - 1.0) ;
>> 336 newton1 = newton ;
>> 337 newton = newton1 - temp1/temp ; // Newton's method
>> 338 }
>> 339 while(fabs(newton - newton1) > tolerance) ;
>> 340
>> 341 fAbscissa[fNumber-i] = newton ;
>> 342 fWeight[fNumber-i] = 2.0/((1.0 - newton*newton)*temp*temp) ;
>> 343 }
>> 344 //
>> 345 // Now we ready to get integral
>> 346 //
>> 347
>> 348 xMean = 0.5*(a + b) ;
>> 349 xDiff = 0.5*(b - a) ;
>> 350 integral = 0.0 ;
>> 351 for(i=0;i<fNumber;i++)
>> 352 {
>> 353 dx = xDiff*fAbscissa[i] ;
>> 354 integral += fWeight[i]*( (typeT.*f)(xMean + dx) +
>> 355 (typeT.*f)(xMean - dx) ) ;
>> 356 }
>> 357 return integral *= xDiff ;
>> 358 }
369 359
370 ////////////////////////////////////////////// 360 ///////////////////////////////////////////////////////////////////////
371 // 361 //
372 // Convenient for using with the pointer 'this 362 // Convenient for using with the pointer 'this'
373 363
374 template <class T, class F> << 364 template <class T, class F> G4double
375 G4double G4Integrator<T, F>::Legendre(T* ptrT, << 365 G4Integrator::Legendre( T* ptrT, F f, G4double a, G4double b, G4int nLegendre)
376 G4int nL <<
377 { 366 {
378 return Legendre(*ptrT, f, a, b, nLegendre); << 367 return Legendre(*ptrT,f,a,b,nLegendre) ;
379 } 368 }
380 369
381 ////////////////////////////////////////////// 370 ///////////////////////////////////////////////////////////////////////
382 // 371 //
383 // Convenient for using with global scope func 372 // Convenient for using with global scope function f
384 373
385 template <class T, class F> <<
386 G4double G4Integrator<T, F>::Legendre(G4double <<
387 G4double <<
388 { <<
389 G4double nwt, nwt1, temp1, temp2, temp3, tem <<
390 G4double xDiff, xMean, dx, integral; <<
391 <<
392 const G4double tolerance = 1.6e-10; <<
393 G4int i, j, k = nLegendre; <<
394 G4int fNumber = (nLegendre + 1) / 2; <<
395 <<
396 if(2 * fNumber != k) <<
397 { <<
398 G4Exception("G4Integrator<T,F>::Legendre(. <<
399 FatalException, "Invalid (odd) <<
400 } <<
401 <<
402 G4double* fAbscissa = new G4double[fNumber]; <<
403 G4double* fWeight = new G4double[fNumber]; <<
404 <<
405 for(i = 1; i <= fNumber; i++) // Loop over <<
406 { <<
407 nwt = std::cos(CLHEP::pi * (i - 0.25) / <<
408 (k + 0.5)); // Initial roo <<
409 <<
410 do // loop of Newton's method <<
411 { <<
412 temp1 = 1.0; <<
413 temp2 = 0.0; <<
414 for(j = 1; j <= k; ++j) <<
415 { <<
416 temp3 = temp2; <<
417 temp2 = temp1; <<
418 temp1 = ((2.0 * j - 1.0) * nwt * temp2 <<
419 } <<
420 temp = k * (nwt * temp1 - temp2) / (nwt <<
421 nwt1 = nwt; <<
422 nwt = nwt1 - temp1 / temp; // Newton's <<
423 } while(std::fabs(nwt - nwt1) > tolerance) <<
424 <<
425 fAbscissa[fNumber - i] = nwt; <<
426 fWeight[fNumber - i] = 2.0 / ((1.0 - nwt <<
427 } <<
428 374
429 // << 375 G4double G4Integrator::
430 // Now we ready to get integral << 376 Legendre( G4double (*f)(G4double), G4double a, G4double b, G4int nLegendre)
431 // << 377 {
432 << 378 G4double newton, newton1, temp1, temp2, temp3, temp ;
433 xMean = 0.5 * (a + b); << 379 G4double xDiff, xMean, dx, integral ;
434 xDiff = 0.5 * (b - a); <<
435 integral = 0.0; <<
436 for(i = 0; i < fNumber; ++i) <<
437 { <<
438 dx = xDiff * fAbscissa[i]; <<
439 integral += fWeight[i] * ((*f)(xMean + dx) <<
440 } <<
441 delete[] fAbscissa; <<
442 delete[] fWeight; <<
443 380
444 return integral *= xDiff; << 381 const G4double tolerance = 1.6e-10 ;
445 } << 382 G4int i, j, k = nLegendre ;
>> 383 G4int fNumber = (nLegendre + 1)/2 ;
>> 384
>> 385 if(2*fNumber != k)
>> 386 {
>> 387 G4Exception("Invalid (odd) n Legendre in G4Integrator::Legendre") ;
>> 388 }
>> 389
>> 390 G4double* fAbscissa = new G4double[fNumber] ;
>> 391 G4double* fWeight = new G4double[fNumber] ;
>> 392
>> 393 for(i=1;i<=fNumber;i++) // Loop over the desired roots
>> 394 {
>> 395 newton = cos(pi*(i - 0.25)/(k + 0.5)) ; // Initial root approximation
>> 396
>> 397 do // loop of Newton's method
>> 398 {
>> 399 temp1 = 1.0 ;
>> 400 temp2 = 0.0 ;
>> 401 for(j=1;j<=k;j++)
>> 402 {
>> 403 temp3 = temp2 ;
>> 404 temp2 = temp1 ;
>> 405 temp1 = ((2.0*j - 1.0)*newton*temp2 - (j - 1.0)*temp3)/j ;
>> 406 }
>> 407 temp = k*(newton*temp1 - temp2)/(newton*newton - 1.0) ;
>> 408 newton1 = newton ;
>> 409 newton = newton1 - temp1/temp ; // Newton's method
>> 410 }
>> 411 while(fabs(newton - newton1) > tolerance) ;
>> 412
>> 413 fAbscissa[fNumber-i] = newton ;
>> 414 fWeight[fNumber-i] = 2.0/((1.0 - newton*newton)*temp*temp) ;
>> 415 }
>> 416 //
>> 417 // Now we ready to get integral
>> 418 //
>> 419
>> 420 xMean = 0.5*(a + b) ;
>> 421 xDiff = 0.5*(b - a) ;
>> 422 integral = 0.0 ;
>> 423 for(i=0;i<fNumber;i++)
>> 424 {
>> 425 dx = xDiff*fAbscissa[i] ;
>> 426 integral += fWeight[i]*( (*f)(xMean + dx) + (*f)(xMean - dx) ) ;
>> 427 }
>> 428 return integral *= xDiff ;
>> 429 }
446 430
447 ////////////////////////////////////////////// 431 ////////////////////////////////////////////////////////////////////////////
448 // 432 //
449 // Returns the integral of the function to be 433 // Returns the integral of the function to be pointed by T::f between a and b,
450 // by ten point Gauss-Legendre integration: th 434 // by ten point Gauss-Legendre integration: the function is evaluated exactly
451 // ten times at interior points in the range o 435 // ten times at interior points in the range of integration. Since the weights
452 // and abscissas are, in this case, symmetric << 436 // and abscissas are, in this case, symmetric around the midpoint of the
453 // range of integration, there are actually on 437 // range of integration, there are actually only five distinct values of each
454 // Convenient for using with class object type 438 // Convenient for using with class object typeT
455 439
456 template <class T, class F> << 440 template <class T, class F>
457 G4double G4Integrator<T, F>::Legendre10(T& typ << 441 G4double G4Integrator::Legendre10( T& typeT, F f,G4double a, G4double b)
458 { 442 {
459 G4int i; << 443 G4int i ;
460 G4double xDiff, xMean, dx, integral; << 444 G4double xDiff, xMean, dx, integral ;
461 << 445
462 // From Abramowitz M., Stegan I.A. 1964 , Ha << 446 // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 916
463 << 447
464 static const G4double abscissa[] = { 0.14887 << 448 static G4double abscissa[] = { 0.148874338981631, 0.433395394129247,
465 0.67940 << 449 0.679409568299024, 0.865063366688985,
466 0.97390 << 450 0.973906528517172 } ;
467 << 451
468 static const G4double weight[] = { 0.2955242 << 452 static G4double weight[] = { 0.295524224714753, 0.269266719309996,
469 0.2190863 << 453 0.219086362515982, 0.149451349150581,
470 0.0666713 << 454 0.066671344308688 } ;
471 xMean = 0.5 * (a + << 455 xMean = 0.5*(a + b) ;
472 xDiff = 0.5 * (b - << 456 xDiff = 0.5*(b - a) ;
473 integral = 0.0; << 457 integral = 0.0 ;
474 for(i = 0; i < 5; ++i) << 458 for(i=0;i<5;i++)
475 { << 459 {
476 dx = xDiff * abscissa[i]; << 460 dx = xDiff*abscissa[i] ;
477 integral += weight[i] * ((typeT.*f)(xMean << 461 integral += weight[i]*( (typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx)) ;
478 } << 462 }
479 return integral *= xDiff; << 463 return integral *= xDiff ;
480 } 464 }
481 465
482 ////////////////////////////////////////////// 466 ///////////////////////////////////////////////////////////////////////////
483 // 467 //
484 // Convenient for using with the pointer 'this 468 // Convenient for using with the pointer 'this'
485 469
486 template <class T, class F> << 470 template <class T, class F>
487 G4double G4Integrator<T, F>::Legendre10(T* ptr << 471 G4double G4Integrator::Legendre10( T* ptrT, F f,G4double a, G4double b)
488 { 472 {
489 return Legendre10(*ptrT, f, a, b); << 473 return Legendre10(*ptrT,f,a,b) ;
490 } << 474 }
491 475
492 ////////////////////////////////////////////// 476 //////////////////////////////////////////////////////////////////////////
493 // 477 //
494 // Convenient for using with global scope func 478 // Convenient for using with global scope functions
495 479
496 template <class T, class F> << 480 G4double
497 G4double G4Integrator<T, F>::Legendre10(G4doub << 481 G4Integrator::Legendre10( G4double (*f)(G4double), G4double a, G4double b)
498 G4doub << 482 {
499 { << 483 G4int i ;
500 G4int i; << 484 G4double xDiff, xMean, dx, integral ;
501 G4double xDiff, xMean, dx, integral; << 485
502 << 486 // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 916
503 // From Abramowitz M., Stegan I.A. 1964 , Ha << 487
504 << 488 static G4double abscissa[] = { 0.148874338981631, 0.433395394129247,
505 static const G4double abscissa[] = { 0.14887 << 489 0.679409568299024, 0.865063366688985,
506 0.67940 << 490 0.973906528517172 } ;
507 0.97390 << 491
508 << 492 static G4double weight[] = { 0.295524224714753, 0.269266719309996,
509 static const G4double weight[] = { 0.2955242 << 493 0.219086362515982, 0.149451349150581,
510 0.2190863 << 494 0.066671344308688 } ;
511 0.0666713 << 495 xMean = 0.5*(a + b) ;
512 xMean = 0.5 * (a + << 496 xDiff = 0.5*(b - a) ;
513 xDiff = 0.5 * (b - << 497 integral = 0.0 ;
514 integral = 0.0; << 498 for(i=0;i<5;i++)
515 for(i = 0; i < 5; ++i) << 499 {
516 { << 500 dx = xDiff*abscissa[i] ;
517 dx = xDiff * abscissa[i]; << 501 integral += weight[i]*( (*f)(xMean + dx) + (*f)(xMean - dx)) ;
518 integral += weight[i] * ((*f)(xMean + dx) << 502 }
519 } << 503 return integral *= xDiff ;
520 return integral *= xDiff; <<
521 } 504 }
522 505
523 ////////////////////////////////////////////// 506 ///////////////////////////////////////////////////////////////////////
524 // 507 //
525 // Returns the integral of the function to be 508 // Returns the integral of the function to be pointed by T::f between a and b,
526 // by 96 point Gauss-Legendre integration: the 509 // by 96 point Gauss-Legendre integration: the function is evaluated exactly
527 // ten Times at interior points in the range o 510 // ten Times at interior points in the range of integration. Since the weights
528 // and abscissas are, in this case, symmetric << 511 // and abscissas are, in this case, symmetric around the midpoint of the
529 // range of integration, there are actually on 512 // range of integration, there are actually only five distinct values of each
530 // Convenient for using with some class object 513 // Convenient for using with some class object typeT
531 514
532 template <class T, class F> << 515 template <class T, class F>
533 G4double G4Integrator<T, F>::Legendre96(T& typ << 516 G4double G4Integrator::Legendre96( T& typeT, F f,G4double a, G4double b)
534 { 517 {
535 G4int i; << 518 G4int i ;
536 G4double xDiff, xMean, dx, integral; << 519 G4double xDiff, xMean, dx, integral ;
537 << 520
538 // From Abramowitz M., Stegan I.A. 1964 , Ha << 521 // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919
539 << 522
540 static const G4double abscissa[] = { << 523 static G4double
541 0.016276744849602969579, 0.048812985136049 << 524 abscissa[] = {
542 0.081297495464425558994, 0.113695850110665 << 525 0.016276744849602969579, 0.048812985136049731112,
543 0.145973714654896941989, 0.178096882367618 << 526 0.081297495464425558994, 0.113695850110665920911,
544 << 527 0.145973714654896941989, 0.178096882367618602759, // 6
545 0.210031310460567203603, 0.241743156163840 << 528
546 0.273198812591049141487, 0.304364944354496 << 529 0.210031310460567203603, 0.241743156163840012328,
547 0.335208522892625422616, 0.365696861472313 << 530 0.273198812591049141487, 0.304364944354496353024,
548 << 531 0.335208522892625422616, 0.365696861472313635031, // 12
549 0.395797649828908603285, 0.425478988407300 << 532
550 0.454709422167743008636, 0.483457973920596 << 533 0.395797649828908603285, 0.425478988407300545365,
551 0.511694177154667673586, 0.539388108324357 << 534 0.454709422167743008636, 0.483457973920596359768,
552 << 535 0.511694177154667673586, 0.539388108324357436227, // 18
553 0.566510418561397168404, 0.593032364777572 << 536
554 0.618925840125468570386, 0.644163403784967 << 537 0.566510418561397168404, 0.593032364777572080684,
555 0.668718310043916153953, 0.692564536642171 << 538 0.618925840125468570386, 0.644163403784967106798,
556 << 539 0.668718310043916153953, 0.692564536642171561344, // 24
557 0.715676812348967626225, 0.738030643744400 << 540
558 0.759602341176647498703, 0.780369043867433 << 541 0.715676812348967626225, 0.738030643744400132851,
559 0.800308744139140817229, 0.819400310737931 << 542 0.759602341176647498703, 0.780369043867433217604,
560 << 543 0.800308744139140817229, 0.819400310737931675539, // 30
561 0.837623511228187121494, 0.854959033434601 << 544
562 0.871388505909296502874, 0.886894517402420 << 545 0.837623511228187121494, 0.854959033434601455463,
563 0.901460635315852341319, 0.915071423120898 << 546 0.871388505909296502874, 0.886894517402420416057,
564 << 547 0.901460635315852341319, 0.915071423120898074206, // 36
565 0.927712456722308690965, 0.939370339752755 << 548
566 0.950032717784437635756, 0.959688291448742 << 549 0.927712456722308690965, 0.939370339752755216932,
567 0.968326828463264212174, 0.975939174585136 << 550 0.950032717784437635756, 0.959688291448742539300,
568 << 551 0.968326828463264212174, 0.975939174585136466453, // 42
569 0.982517263563014677447, 0.988054126329623 << 552
570 0.992543900323762624572, 0.995981842987209 << 553 0.982517263563014677447, 0.988054126329623799481,
571 0.998364375863181677724, 0.999689503883230 << 554 0.992543900323762624572, 0.995981842987209290650,
572 }; << 555 0.998364375863181677724, 0.999689503883230766828 // 48
573 << 556 } ;
574 static const G4double weight[] = { << 557
575 0.032550614492363166242, 0.032516118713868 << 558 static G4double
576 0.032447163714064269364, 0.032343822568575 << 559 weight[] = {
577 0.032206204794030250669, 0.032034456231992 << 560 0.032550614492363166242, 0.032516118713868835987,
578 << 561 0.032447163714064269364, 0.032343822568575928429,
579 0.031828758894411006535, 0.031589330770727 << 562 0.032206204794030250669, 0.032034456231992663218, // 6
580 0.031316425596862355813, 0.031010332586313 << 563
581 0.030671376123669149014, 0.030299915420827 << 564 0.031828758894411006535, 0.031589330770727168558,
582 << 565 0.031316425596862355813, 0.031010332586313837423,
583 0.029896344136328385984, 0.029461089958167 << 566 0.030671376123669149014, 0.030299915420827593794, // 12
584 0.028994614150555236543, 0.028497411065085 << 567
585 0.027970007616848334440, 0.027412962726029 << 568 0.029896344136328385984, 0.029461089958167905970,
586 << 569 0.028994614150555236543, 0.028497411065085385646,
587 0.026826866725591762198, 0.026212340735672 << 570 0.027970007616848334440, 0.027412962726029242823, // 18
588 0.025570036005349361499, 0.024900633222483 << 571
589 0.024204841792364691282, 0.023483399085926 << 572 0.026826866725591762198, 0.026212340735672413913,
590 << 573 0.025570036005349361499, 0.024900633222483610288,
591 0.022737069658329374001, 0.021966644438744 << 574 0.024204841792364691282, 0.023483399085926219842, // 24
592 0.021172939892191298988, 0.020356797154333 << 575
593 0.019519081140145022410, 0.018660679627411 << 576 0.022737069658329374001, 0.021966644438744349195,
594 << 577 0.021172939892191298988, 0.020356797154333324595,
595 0.017782502316045260838, 0.016885479864245 << 578 0.019519081140145022410, 0.018660679627411467385, // 30
596 0.015970562902562291381, 0.015038721026994 << 579
597 0.014090941772314860916, 0.013128229566961 << 580 0.017782502316045260838, 0.016885479864245172450,
598 << 581 0.015970562902562291381, 0.015038721026994938006,
599 0.012151604671088319635, 0.011162102099838 << 582 0.014090941772314860916, 0.013128229566961572637, // 36
600 0.010160770535008415758, 0.009148671230783 << 583
601 0.008126876925698759217, 0.007096470791153 << 584 0.012151604671088319635, 0.011162102099838498591,
602 << 585 0.010160770535008415758, 0.009148671230783386633,
603 0.006058545504235961683, 0.005014202742927 << 586 0.008126876925698759217, 0.007096470791153865269, // 42
604 0.003964554338444686674, 0.002910731817934 << 587
605 0.001853960788946921732, 0.000796792065552 << 588 0.006058545504235961683, 0.005014202742927517693,
606 }; << 589 0.003964554338444686674, 0.002910731817934946408,
607 xMean = 0.5 * (a + b); << 590 0.001853960788946921732, 0.000796792065552012429 // 48
608 xDiff = 0.5 * (b - a); << 591 } ;
609 integral = 0.0; << 592 xMean = 0.5*(a + b) ;
610 for(i = 0; i < 48; ++i) << 593 xDiff = 0.5*(b - a) ;
611 { << 594 integral = 0.0 ;
612 dx = xDiff * abscissa[i]; << 595 for(i=0;i<48;i++)
613 integral += weight[i] * ((typeT.*f)(xMean << 596 {
614 } << 597 dx = xDiff*abscissa[i] ;
615 return integral *= xDiff; << 598 integral += weight[i]*((typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx)) ;
>> 599 }
>> 600 return integral *= xDiff ;
616 } 601 }
617 602
618 ////////////////////////////////////////////// 603 ///////////////////////////////////////////////////////////////////////
619 // 604 //
620 // Convenient for using with the pointer 'this 605 // Convenient for using with the pointer 'this'
621 606
622 template <class T, class F> << 607 template <class T, class F>
623 G4double G4Integrator<T, F>::Legendre96(T* ptr << 608 G4double G4Integrator::Legendre96( T* ptrT, F f,G4double a, G4double b)
624 { 609 {
625 return Legendre96(*ptrT, f, a, b); << 610 return Legendre96(*ptrT,f,a,b) ;
626 } << 611 }
627 612
628 ////////////////////////////////////////////// 613 ///////////////////////////////////////////////////////////////////////
629 // 614 //
630 // Convenient for using with global scope func << 615 // Convenient for using with global scope function f
631 616
632 template <class T, class F> << 617 G4double
633 G4double G4Integrator<T, F>::Legendre96(G4doub << 618 G4Integrator::Legendre96( G4double (*f)(G4double), G4double a, G4double b)
634 G4doub << 619 {
635 { << 620 G4int i ;
636 G4int i; << 621 G4double xDiff, xMean, dx, integral ;
637 G4double xDiff, xMean, dx, integral; << 622
638 << 623 // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919
639 // From Abramowitz M., Stegan I.A. 1964 , Ha << 624
640 << 625 static G4double
641 static const G4double abscissa[] = { << 626 abscissa[] = {
642 0.016276744849602969579, 0.048812985136049 << 627 0.016276744849602969579, 0.048812985136049731112,
643 0.081297495464425558994, 0.113695850110665 << 628 0.081297495464425558994, 0.113695850110665920911,
644 0.145973714654896941989, 0.178096882367618 << 629 0.145973714654896941989, 0.178096882367618602759, // 6
645 << 630
646 0.210031310460567203603, 0.241743156163840 << 631 0.210031310460567203603, 0.241743156163840012328,
647 0.273198812591049141487, 0.304364944354496 << 632 0.273198812591049141487, 0.304364944354496353024,
648 0.335208522892625422616, 0.365696861472313 << 633 0.335208522892625422616, 0.365696861472313635031, // 12
649 << 634
650 0.395797649828908603285, 0.425478988407300 << 635 0.395797649828908603285, 0.425478988407300545365,
651 0.454709422167743008636, 0.483457973920596 << 636 0.454709422167743008636, 0.483457973920596359768,
652 0.511694177154667673586, 0.539388108324357 << 637 0.511694177154667673586, 0.539388108324357436227, // 18
653 << 638
654 0.566510418561397168404, 0.593032364777572 << 639 0.566510418561397168404, 0.593032364777572080684,
655 0.618925840125468570386, 0.644163403784967 << 640 0.618925840125468570386, 0.644163403784967106798,
656 0.668718310043916153953, 0.692564536642171 << 641 0.668718310043916153953, 0.692564536642171561344, // 24
657 << 642
658 0.715676812348967626225, 0.738030643744400 << 643 0.715676812348967626225, 0.738030643744400132851,
659 0.759602341176647498703, 0.780369043867433 << 644 0.759602341176647498703, 0.780369043867433217604,
660 0.800308744139140817229, 0.819400310737931 << 645 0.800308744139140817229, 0.819400310737931675539, // 30
661 << 646
662 0.837623511228187121494, 0.854959033434601 << 647 0.837623511228187121494, 0.854959033434601455463,
663 0.871388505909296502874, 0.886894517402420 << 648 0.871388505909296502874, 0.886894517402420416057,
664 0.901460635315852341319, 0.915071423120898 << 649 0.901460635315852341319, 0.915071423120898074206, // 36
665 << 650
666 0.927712456722308690965, 0.939370339752755 << 651 0.927712456722308690965, 0.939370339752755216932,
667 0.950032717784437635756, 0.959688291448742 << 652 0.950032717784437635756, 0.959688291448742539300,
668 0.968326828463264212174, 0.975939174585136 << 653 0.968326828463264212174, 0.975939174585136466453, // 42
669 << 654
670 0.982517263563014677447, 0.988054126329623 << 655 0.982517263563014677447, 0.988054126329623799481,
671 0.992543900323762624572, 0.995981842987209 << 656 0.992543900323762624572, 0.995981842987209290650,
672 0.998364375863181677724, 0.999689503883230 << 657 0.998364375863181677724, 0.999689503883230766828 // 48
673 }; << 658 } ;
674 << 659
675 static const G4double weight[] = { << 660 static G4double
676 0.032550614492363166242, 0.032516118713868 << 661 weight[] = {
677 0.032447163714064269364, 0.032343822568575 << 662 0.032550614492363166242, 0.032516118713868835987,
678 0.032206204794030250669, 0.032034456231992 << 663 0.032447163714064269364, 0.032343822568575928429,
679 << 664 0.032206204794030250669, 0.032034456231992663218, // 6
680 0.031828758894411006535, 0.031589330770727 << 665
681 0.031316425596862355813, 0.031010332586313 << 666 0.031828758894411006535, 0.031589330770727168558,
682 0.030671376123669149014, 0.030299915420827 << 667 0.031316425596862355813, 0.031010332586313837423,
683 << 668 0.030671376123669149014, 0.030299915420827593794, // 12
684 0.029896344136328385984, 0.029461089958167 << 669
685 0.028994614150555236543, 0.028497411065085 << 670 0.029896344136328385984, 0.029461089958167905970,
686 0.027970007616848334440, 0.027412962726029 << 671 0.028994614150555236543, 0.028497411065085385646,
687 << 672 0.027970007616848334440, 0.027412962726029242823, // 18
688 0.026826866725591762198, 0.026212340735672 << 673
689 0.025570036005349361499, 0.024900633222483 << 674 0.026826866725591762198, 0.026212340735672413913,
690 0.024204841792364691282, 0.023483399085926 << 675 0.025570036005349361499, 0.024900633222483610288,
691 << 676 0.024204841792364691282, 0.023483399085926219842, // 24
692 0.022737069658329374001, 0.021966644438744 << 677
693 0.021172939892191298988, 0.020356797154333 << 678 0.022737069658329374001, 0.021966644438744349195,
694 0.019519081140145022410, 0.018660679627411 << 679 0.021172939892191298988, 0.020356797154333324595,
695 << 680 0.019519081140145022410, 0.018660679627411467385, // 30
696 0.017782502316045260838, 0.016885479864245 << 681
697 0.015970562902562291381, 0.015038721026994 << 682 0.017782502316045260838, 0.016885479864245172450,
698 0.014090941772314860916, 0.013128229566961 << 683 0.015970562902562291381, 0.015038721026994938006,
699 << 684 0.014090941772314860916, 0.013128229566961572637, // 36
700 0.012151604671088319635, 0.011162102099838 << 685
701 0.010160770535008415758, 0.009148671230783 << 686 0.012151604671088319635, 0.011162102099838498591,
702 0.008126876925698759217, 0.007096470791153 << 687 0.010160770535008415758, 0.009148671230783386633,
703 << 688 0.008126876925698759217, 0.007096470791153865269, // 42
704 0.006058545504235961683, 0.005014202742927 << 689
705 0.003964554338444686674, 0.002910731817934 << 690 0.006058545504235961683, 0.005014202742927517693,
706 0.001853960788946921732, 0.000796792065552 << 691 0.003964554338444686674, 0.002910731817934946408,
707 }; << 692 0.001853960788946921732, 0.000796792065552012429 // 48
708 xMean = 0.5 * (a + b); << 693 } ;
709 xDiff = 0.5 * (b - a); << 694 xMean = 0.5*(a + b) ;
710 integral = 0.0; << 695 xDiff = 0.5*(b - a) ;
711 for(i = 0; i < 48; ++i) << 696 integral = 0.0 ;
712 { << 697 for(i=0;i<48;i++)
713 dx = xDiff * abscissa[i]; << 698 {
714 integral += weight[i] * ((*f)(xMean + dx) << 699 dx = xDiff*abscissa[i] ;
715 } << 700 integral += weight[i]*((*f)(xMean + dx) + (*f)(xMean - dx)) ;
716 return integral *= xDiff; << 701 }
>> 702 return integral *= xDiff ;
717 } 703 }
718 704
719 ////////////////////////////////////////////// 705 //////////////////////////////////////////////////////////////////////////////
720 // 706 //
721 // Methods involving Chebyshev polynomials << 707 // Methods involving Chebyshev polynomials
722 // 708 //
723 ////////////////////////////////////////////// 709 ///////////////////////////////////////////////////////////////////////////
724 // 710 //
725 // Integrates function pointed by T::f from a << 711 // Integrates function pointed by T::f from a to b by Gauss-Chebyshev
726 // quadrature method. 712 // quadrature method.
727 // Convenient for using with class object type 713 // Convenient for using with class object typeT
728 714
729 template <class T, class F> << 715 template <class T, class F> G4double
730 G4double G4Integrator<T, F>::Chebyshev(T& type << 716 G4Integrator::Chebyshev( T& typeT, F f, G4double a,
731 G4int n << 717 G4double b, G4int nChebyshev )
732 { << 718 {
733 G4int i; << 719 G4int i ;
734 G4double xDiff, xMean, dx, integral = 0.0; << 720 G4double xDiff, xMean, dx, integral = 0.0 ;
735 << 721
736 G4int fNumber = nChebyshev; // Try to << 722 G4int fNumber = nChebyshev ; // Try to reduce fNumber twice ??
737 G4double cof = CLHEP::pi / fNumber; << 723 G4double cof = pi/fNumber ;
738 G4double* fAbscissa = new G4double[fNumber]; << 724 G4double* fAbscissa = new G4double[fNumber] ;
739 G4double* fWeight = new G4double[fNumber]; << 725 G4double* fWeight = new G4double[fNumber] ;
740 for(i = 0; i < fNumber; ++i) << 726 for(i=0;i<fNumber;i++)
741 { << 727 {
742 fAbscissa[i] = std::cos(cof * (i + 0.5)); << 728 fAbscissa[i] = cos(cof*(i + 0.5)) ;
743 fWeight[i] = cof * std::sqrt(1 - fAbscis << 729 fWeight[i] = cof*sqrt(1 - fAbscissa[i]*fAbscissa[i]) ;
744 } << 730 }
745 << 731 //
746 // << 732 // Now we ready to estimate the integral
747 // Now we ready to estimate the integral << 733 //
748 // << 734 xMean = 0.5*(a + b) ;
749 << 735 xDiff = 0.5*(b - a) ;
750 xMean = 0.5 * (a + b); << 736 for(i=0;i<fNumber;i++)
751 xDiff = 0.5 * (b - a); << 737 {
752 for(i = 0; i < fNumber; ++i) << 738 dx = xDiff*fAbscissa[i] ;
753 { << 739 integral += fWeight[i]*(typeT.*f)(xMean + dx) ;
754 dx = xDiff * fAbscissa[i]; << 740 }
755 integral += fWeight[i] * (typeT.*f)(xMean << 741 return integral *= xDiff ;
756 } <<
757 delete[] fAbscissa; <<
758 delete[] fWeight; <<
759 return integral *= xDiff; <<
760 } 742 }
761 743
762 ////////////////////////////////////////////// 744 ///////////////////////////////////////////////////////////////////////
763 // 745 //
764 // Convenient for using with 'this' pointer 746 // Convenient for using with 'this' pointer
765 747
766 template <class T, class F> << 748 template <class T, class F> G4double
767 G4double G4Integrator<T, F>::Chebyshev(T* ptrT << 749 G4Integrator::Chebyshev( T* ptrT, F f, G4double a, G4double b, G4int n)
768 G4int n <<
769 { 750 {
770 return Chebyshev(*ptrT, f, a, b, n); << 751 return Chebyshev(*ptrT,f,a,b,n) ;
771 } << 752 }
772 753
773 ////////////////////////////////////////////// 754 ////////////////////////////////////////////////////////////////////////
774 // 755 //
775 // For use with global scope functions f << 756 // For use with global scope functions f
776 <<
777 template <class T, class F> <<
778 G4double G4Integrator<T, F>::Chebyshev(G4doubl <<
779 G4doubl <<
780 { <<
781 G4int i; <<
782 G4double xDiff, xMean, dx, integral = 0.0; <<
783 <<
784 G4int fNumber = nChebyshev; // Try to <<
785 G4double cof = CLHEP::pi / fNumber; <<
786 G4double* fAbscissa = new G4double[fNumber]; <<
787 G4double* fWeight = new G4double[fNumber]; <<
788 for(i = 0; i < fNumber; ++i) <<
789 { <<
790 fAbscissa[i] = std::cos(cof * (i + 0.5)); <<
791 fWeight[i] = cof * std::sqrt(1 - fAbscis <<
792 } <<
793 757
794 // << 758 G4double
795 // Now we ready to estimate the integral << 759 G4Integrator::Chebyshev( G4double (*f)(G4double),
796 // << 760 G4double a, G4double b, G4int nChebyshev)
797 << 761 {
798 xMean = 0.5 * (a + b); << 762 G4int i ;
799 xDiff = 0.5 * (b - a); << 763 G4double xDiff, xMean, dx, integral = 0.0 ;
800 for(i = 0; i < fNumber; ++i) << 764
801 { << 765 G4int fNumber = nChebyshev ; // Try to reduce fNumber twice ??
802 dx = xDiff * fAbscissa[i]; << 766 G4double cof = pi/fNumber ;
803 integral += fWeight[i] * (*f)(xMean + dx); << 767 G4double* fAbscissa = new G4double[fNumber] ;
804 } << 768 G4double* fWeight = new G4double[fNumber] ;
805 delete[] fAbscissa; << 769 for(i=0;i<fNumber;i++)
806 delete[] fWeight; << 770 {
807 return integral *= xDiff; << 771 fAbscissa[i] = cos(cof*(i + 0.5)) ;
>> 772 fWeight[i] = cof*sqrt(1 - fAbscissa[i]*fAbscissa[i]) ;
>> 773 }
>> 774 //
>> 775 // Now we ready to estimate the integral
>> 776 //
>> 777 xMean = 0.5*(a + b) ;
>> 778 xDiff = 0.5*(b - a) ;
>> 779 for(i=0;i<fNumber;i++)
>> 780 {
>> 781 dx = xDiff*fAbscissa[i] ;
>> 782 integral += fWeight[i]*(*f)(xMean + dx) ;
>> 783 }
>> 784 return integral *= xDiff ;
808 } 785 }
809 786
810 ////////////////////////////////////////////// 787 //////////////////////////////////////////////////////////////////////
811 // 788 //
812 // Method involving Laguerre polynomials 789 // Method involving Laguerre polynomials
813 // 790 //
814 ////////////////////////////////////////////// 791 //////////////////////////////////////////////////////////////////////
815 // 792 //
816 // Integral from zero to infinity of std::pow( << 793 // Integral from zero to infinity of pow(x,alpha)*exp(-x)*f(x).
817 // The value of nLaguerre sets the accuracy. 794 // The value of nLaguerre sets the accuracy.
818 // The function creates arrays fAbscissa[0,.., << 795 // The function creates arrays fAbscissa[0,..,nLaguerre-1] and
819 // fWeight[0,..,nLaguerre-1] . << 796 // fWeight[0,..,nLaguerre-1] .
820 // Convenient for using with class object 'typ 797 // Convenient for using with class object 'typeT' and (typeT.*f) function
821 // (T::f) 798 // (T::f)
822 799
823 template <class T, class F> << 800 template <class T, class F> G4double
824 G4double G4Integrator<T, F>::Laguerre(T& typeT << 801 G4Integrator::Laguerre( T& typeT, F f, G4double alpha, G4int nLaguerre )
825 G4int nL << 802 {
826 { << 803 const G4double tolerance = 1.0e-10 ;
827 const G4double tolerance = 1.0e-10; << 804 const G4int maxNumber = 12 ;
828 const G4int maxNumber = 12; << 805 G4int i, j, k ;
829 G4int i, j, k; << 806 G4double newton, newton1, temp1, temp2, temp3, temp, cofi ;
830 G4double nwt = 0., nwt1, temp1, temp2, << 807 G4double integral = 0.0 ;
831 G4double integral = 0.0; << 808
832 << 809 G4int fNumber = nLaguerre ;
833 G4int fNumber = nLaguerre; << 810 G4double* fAbscissa = new G4double[fNumber] ;
834 G4double* fAbscissa = new G4double[fNumber]; << 811 G4double* fWeight = new G4double[fNumber] ;
835 G4double* fWeight = new G4double[fNumber]; << 812
>> 813 for(i=1;i<=fNumber;i++) // Loop over the desired roots
>> 814 {
>> 815 if(i == 1)
>> 816 {
>> 817 newton = (1.0 + alpha)*(3.0 + 0.92*alpha)/(1.0 + 2.4*fNumber + 1.8*alpha) ;
>> 818 }
>> 819 else if(i == 2)
>> 820 {
>> 821 newton += (15.0 + 6.25*alpha)/(1.0 + 0.9*alpha + 2.5*fNumber) ;
>> 822 }
>> 823 else
>> 824 {
>> 825 cofi = i - 2 ;
>> 826 newton += ((1.0+2.55*cofi)/(1.9*cofi) + 1.26*cofi*alpha/(1.0+3.5*cofi))*
>> 827 (newton - fAbscissa[i-3])/(1.0 + 0.3*alpha) ;
>> 828 }
>> 829 for(k=1;k<=maxNumber;k++)
>> 830 {
>> 831 temp1 = 1.0 ;
>> 832 temp2 = 0.0 ;
836 833
837 for(i = 1; i <= fNumber; ++i) // Loop over << 834 for(j=1;j<=fNumber;j++)
838 { << 835 {
839 if(i == 1) << 836 temp3 = temp2 ;
840 { << 837 temp2 = temp1 ;
841 nwt = (1.0 + alpha) * (3.0 + 0.92 * alph << 838 temp1 = ((2*j - 1 + alpha - newton)*temp2 - (j - 1 + alpha)*temp3)/j ;
842 (1.0 + 2.4 * fNumber + 1.8 * alpha << 839 }
843 } << 840 temp = (fNumber*temp1 - (fNumber +alpha)*temp2)/newton ;
844 else if(i == 2) << 841 newton1 = newton ;
845 { << 842 newton = newton1 - temp1/temp ;
846 nwt += (15.0 + 6.25 * alpha) / (1.0 + 0. << 843
847 } << 844 if(fabs(newton - newton1) <= tolerance)
848 else << 845 {
849 { << 846 break ;
850 cofi = i - 2; << 847 }
851 nwt += ((1.0 + 2.55 * cofi) / (1.9 * cof << 848 }
852 1.26 * cofi * alpha / (1.0 + 3.5 << 849 if(k > maxNumber)
853 (nwt - fAbscissa[i - 3]) / (1.0 + << 850 {
854 } << 851 G4Exception("Too many (>12) iterations in G4Integration::Laguerre") ;
855 for(k = 1; k <= maxNumber; ++k) << 852 }
856 { << 853
857 temp1 = 1.0; << 854 fAbscissa[i-1] = newton ;
858 temp2 = 0.0; << 855 fWeight[i-1] = -exp(GammaLogarithm(alpha + fNumber) -
859 << 856 GammaLogarithm((G4double)fNumber))/(temp*fNumber*temp2) ;
860 for(j = 1; j <= fNumber; ++j) << 857 }
861 { << 858 //
862 temp3 = temp2; << 859 // Integral evaluation
863 temp2 = temp1; << 860 //
864 temp1 = << 861 for(i=0;i<fNumber;i++)
865 ((2 * j - 1 + alpha - nwt) * temp2 - << 862 {
866 } << 863 integral += fWeight[i]*(typeT.*f)(fAbscissa[i]) ;
867 temp = (fNumber * temp1 - (fNumber + alp << 864 }
868 nwt1 = nwt; << 865 return integral ;
869 nwt = nwt1 - temp1 / temp; << 866 }
870 <<
871 if(std::fabs(nwt - nwt1) <= tolerance) <<
872 { <<
873 break; <<
874 } <<
875 } <<
876 if(k > maxNumber) <<
877 { <<
878 G4Exception("G4Integrator<T,F>::Laguerre <<
879 FatalException, "Too many (> <<
880 } <<
881 <<
882 fAbscissa[i - 1] = nwt; <<
883 fWeight[i - 1] = -std::exp(GammaLogarith <<
884 GammaLogarithm( <<
885 (temp * fNumber * temp2); <<
886 } <<
887 867
888 // <<
889 // Integral evaluation <<
890 // <<
891 868
892 for(i = 0; i < fNumber; ++i) <<
893 { <<
894 integral += fWeight[i] * (typeT.*f)(fAbsci <<
895 } <<
896 delete[] fAbscissa; <<
897 delete[] fWeight; <<
898 return integral; <<
899 } <<
900 869
901 ////////////////////////////////////////////// 870 //////////////////////////////////////////////////////////////////////
902 // 871 //
903 // 872 //
904 873
905 template <class T, class F> << 874 template <class T, class F> G4double
906 G4double G4Integrator<T, F>::Laguerre(T* ptrT, << 875 G4Integrator::Laguerre( T* ptrT, F f, G4double alpha, G4int nLaguerre )
907 G4int nL <<
908 { 876 {
909 return Laguerre(*ptrT, f, alpha, nLaguerre); << 877 return Laguerre(*ptrT,f,alpha,nLaguerre) ;
910 } 878 }
911 879
912 ////////////////////////////////////////////// 880 ////////////////////////////////////////////////////////////////////////
913 // 881 //
914 // For use with global scope functions f << 882 // For use with global scope functions f
915 <<
916 template <class T, class F> <<
917 G4double G4Integrator<T, F>::Laguerre(G4double <<
918 G4int nL <<
919 { <<
920 const G4double tolerance = 1.0e-10; <<
921 const G4int maxNumber = 12; <<
922 G4int i, j, k; <<
923 G4double nwt = 0., nwt1, temp1, temp2, <<
924 G4double integral = 0.0; <<
925 <<
926 G4int fNumber = nLaguerre; <<
927 G4double* fAbscissa = new G4double[fNumber]; <<
928 G4double* fWeight = new G4double[fNumber]; <<
929 <<
930 for(i = 1; i <= fNumber; ++i) // Loop over <<
931 { <<
932 if(i == 1) <<
933 { <<
934 nwt = (1.0 + alpha) * (3.0 + 0.92 * alph <<
935 (1.0 + 2.4 * fNumber + 1.8 * alpha <<
936 } <<
937 else if(i == 2) <<
938 { <<
939 nwt += (15.0 + 6.25 * alpha) / (1.0 + 0. <<
940 } <<
941 else <<
942 { <<
943 cofi = i - 2; <<
944 nwt += ((1.0 + 2.55 * cofi) / (1.9 * cof <<
945 1.26 * cofi * alpha / (1.0 + 3.5 <<
946 (nwt - fAbscissa[i - 3]) / (1.0 + <<
947 } <<
948 for(k = 1; k <= maxNumber; ++k) <<
949 { <<
950 temp1 = 1.0; <<
951 temp2 = 0.0; <<
952 <<
953 for(j = 1; j <= fNumber; ++j) <<
954 { <<
955 temp3 = temp2; <<
956 temp2 = temp1; <<
957 temp1 = <<
958 ((2 * j - 1 + alpha - nwt) * temp2 - <<
959 } <<
960 temp = (fNumber * temp1 - (fNumber + alp <<
961 nwt1 = nwt; <<
962 nwt = nwt1 - temp1 / temp; <<
963 <<
964 if(std::fabs(nwt - nwt1) <= tolerance) <<
965 { <<
966 break; <<
967 } <<
968 } <<
969 if(k > maxNumber) <<
970 { <<
971 G4Exception("G4Integrator<T,F>::Laguerre <<
972 "Too many (>12) iterations." <<
973 } <<
974 <<
975 fAbscissa[i - 1] = nwt; <<
976 fWeight[i - 1] = -std::exp(GammaLogarith <<
977 GammaLogarithm( <<
978 (temp * fNumber * temp2); <<
979 } <<
980 883
981 // << 884 G4double
982 // Integral evaluation << 885 G4Integrator::Laguerre( G4double (*f)(G4double),
983 // << 886 G4double alpha, G4int nLaguerre)
>> 887 {
>> 888 const G4double tolerance = 1.0e-10 ;
>> 889 const G4int maxNumber = 12 ;
>> 890 G4int i, j, k ;
>> 891 G4double newton, newton1, temp1, temp2, temp3, temp, cofi ;
>> 892 G4double integral = 0.0 ;
>> 893
>> 894 G4int fNumber = nLaguerre ;
>> 895 G4double* fAbscissa = new G4double[fNumber] ;
>> 896 G4double* fWeight = new G4double[fNumber] ;
>> 897
>> 898 for(i=1;i<=fNumber;i++) // Loop over the desired roots
>> 899 {
>> 900 if(i == 1)
>> 901 {
>> 902 newton = (1.0 + alpha)*(3.0 + 0.92*alpha)/(1.0 + 2.4*fNumber + 1.8*alpha) ;
>> 903 }
>> 904 else if(i == 2)
>> 905 {
>> 906 newton += (15.0 + 6.25*alpha)/(1.0 + 0.9*alpha + 2.5*fNumber) ;
>> 907 }
>> 908 else
>> 909 {
>> 910 cofi = i - 2 ;
>> 911 newton += ((1.0+2.55*cofi)/(1.9*cofi) + 1.26*cofi*alpha/(1.0+3.5*cofi))*
>> 912 (newton - fAbscissa[i-3])/(1.0 + 0.3*alpha) ;
>> 913 }
>> 914 for(k=1;k<=maxNumber;k++)
>> 915 {
>> 916 temp1 = 1.0 ;
>> 917 temp2 = 0.0 ;
984 918
985 for(i = 0; i < fNumber; i++) << 919 for(j=1;j<=fNumber;j++)
986 { << 920 {
987 integral += fWeight[i] * (*f)(fAbscissa[i] << 921 temp3 = temp2 ;
988 } << 922 temp2 = temp1 ;
989 delete[] fAbscissa; << 923 temp1 = ((2*j - 1 + alpha - newton)*temp2 - (j - 1 + alpha)*temp3)/j ;
990 delete[] fWeight; << 924 }
991 return integral; << 925 temp = (fNumber*temp1 - (fNumber +alpha)*temp2)/newton ;
>> 926 newton1 = newton ;
>> 927 newton = newton1 - temp1/temp ;
>> 928
>> 929 if(fabs(newton - newton1) <= tolerance)
>> 930 {
>> 931 break ;
>> 932 }
>> 933 }
>> 934 if(k > maxNumber)
>> 935 {
>> 936 G4Exception("Too many (>12) iterations in G4Integration::Laguerre") ;
>> 937 }
>> 938
>> 939 fAbscissa[i-1] = newton ;
>> 940 fWeight[i-1] = -exp(GammaLogarithm(alpha + fNumber) -
>> 941 GammaLogarithm((G4double)fNumber))/(temp*fNumber*temp2) ;
>> 942 }
>> 943 //
>> 944 // Integral evaluation
>> 945 //
>> 946 for(i=0;i<fNumber;i++)
>> 947 {
>> 948 integral += fWeight[i]*(*f)(fAbscissa[i]) ;
>> 949 }
>> 950 return integral ;
992 } 951 }
993 952
994 ////////////////////////////////////////////// 953 ///////////////////////////////////////////////////////////////////////
995 // 954 //
996 // Auxiliary function which returns the value << 955 // Auxiliary function which returns the value of log(gamma-function(x))
997 // Returns the value ln(Gamma(xx) for xx > 0. << 956 // Returns the value ln(Gamma(xx) for xx > 0. Full accuracy is obtained for
998 // xx > 1. For 0 < xx < 1. the reflection form 957 // xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first.
999 // (Adapted from Numerical Recipes in C) 958 // (Adapted from Numerical Recipes in C)
1000 // 959 //
1001 960
1002 template <class T, class F> << 961 G4double G4Integrator::GammaLogarithm(G4double xx)
1003 G4double G4Integrator<T, F>::GammaLogarithm(G <<
1004 { 962 {
1005 static const G4double cof[6] = { 76.1800917 << 963 static G4double cof[6] = { 76.18009172947146, -86.50532032941677,
1006 24.0140982 << 964 24.01409824083091, -1.231739572450155,
1007 0.12086509 << 965 0.1208650973866179e-2, -0.5395239384953e-5 } ;
1008 G4int j; << 966 register HepInt j;
1009 G4double x = xx - 1.0; << 967 G4double x = xx - 1.0 ;
1010 G4double tmp = x + 5.5; << 968 G4double tmp = x + 5.5 ;
1011 tmp -= (x + 0.5) * std::log(tmp); << 969 tmp -= (x + 0.5) * log(tmp) ;
1012 G4double ser = 1.000000000190015; << 970 G4double ser = 1.000000000190015 ;
1013 971
1014 for(j = 0; j <= 5; ++j) << 972 for ( j = 0; j <= 5; j++ )
1015 { 973 {
1016 x += 1.0; << 974 x += 1.0 ;
1017 ser += cof[j] / x; << 975 ser += cof[j]/x ;
1018 } 976 }
1019 return -tmp + std::log(2.5066282746310005 * << 977 return -tmp + log(2.5066282746310005*ser) ;
1020 } 978 }
1021 979
1022 ///////////////////////////////////////////// 980 ///////////////////////////////////////////////////////////////////////
1023 // 981 //
1024 // Method involving Hermite polynomials 982 // Method involving Hermite polynomials
1025 // 983 //
1026 ///////////////////////////////////////////// 984 ///////////////////////////////////////////////////////////////////////
1027 // 985 //
1028 // 986 //
1029 // Gauss-Hermite method for integration of st << 987 // Gauss-Hermite method for integration of exp(-x*x)*f(x)
1030 // from minus infinity to plus infinity . << 988 // from minus infinity to plus infinity .
1031 // 989 //
1032 990
1033 template <class T, class F> << 991 template <class T, class F>
1034 G4double G4Integrator<T, F>::Hermite(T& typeT << 992 G4double G4Integrator::Hermite( T& typeT, F f, G4int nHermite)
1035 { 993 {
1036 const G4double tolerance = 1.0e-12; << 994 const G4double tolerance = 1.0e-12 ;
1037 const G4int maxNumber = 12; << 995 const G4int maxNumber = 12 ;
1038 << 996
1039 G4int i, j, k; << 997 G4int i, j, k ;
1040 G4double integral = 0.0; << 998 G4double integral = 0.0 ;
1041 G4double nwt = 0., nwt1, temp1, temp2, << 999 G4double newton, newton1, temp1, temp2, temp3, temp ;
1042 1000
1043 G4double piInMinusQ = << 1001 G4double piInMinusQ = pow(pi,-0.25) ; // 1.0/sqrt(sqrt(pi)) ??
1044 std::pow(CLHEP::pi, -0.25); // 1.0/std:: <<
1045 1002
1046 G4int fNumber = (nHermite + 1) / 2; << 1003 G4int fNumber = (nHermite +1)/2 ;
1047 G4double* fAbscissa = new G4double[fNumber] << 1004 G4double* fAbscissa = new G4double[fNumber] ;
1048 G4double* fWeight = new G4double[fNumber] << 1005 G4double* fWeight = new G4double[fNumber] ;
1049 1006
1050 for(i = 1; i <= fNumber; ++i) << 1007 for(i=1;i<=fNumber;i++)
1051 { << 1008 {
1052 if(i == 1) << 1009 if(i == 1)
1053 { << 1010 {
1054 nwt = std::sqrt((G4double)(2 * nHermite << 1011 newton = sqrt((G4double)(2*nHermite + 1)) -
1055 1.85575001 * std::pow((G4double)( << 1012 1.85575001*pow((G4double)(2*nHermite + 1),-0.16666999) ;
1056 } << 1013 }
1057 else if(i == 2) << 1014 else if(i == 2)
1058 { << 1015 {
1059 nwt -= 1.14001 * std::pow((G4double) nH << 1016 newton -= 1.14001*pow((G4double)nHermite,0.425999)/newton ;
1060 } << 1017 }
1061 else if(i == 3) << 1018 else if(i == 3)
1062 { << 1019 {
1063 nwt = 1.86002 * nwt - 0.86002 * fAbscis << 1020 newton = 1.86002*newton - 0.86002*fAbscissa[0] ;
1064 } << 1021 }
1065 else if(i == 4) << 1022 else if(i == 4)
1066 { << 1023 {
1067 nwt = 1.91001 * nwt - 0.91001 * fAbscis << 1024 newton = 1.91001*newton - 0.91001*fAbscissa[1] ;
1068 } << 1025 }
1069 else << 1026 else
1070 { << 1027 {
1071 nwt = 2.0 * nwt - fAbscissa[i - 3]; << 1028 newton = 2.0*newton - fAbscissa[i - 3] ;
1072 } << 1029 }
1073 for(k = 1; k <= maxNumber; ++k) << 1030 for(k=1;k<=maxNumber;k++)
1074 { << 1031 {
1075 temp1 = piInMinusQ; << 1032 temp1 = piInMinusQ ;
1076 temp2 = 0.0; << 1033 temp2 = 0.0 ;
1077 <<
1078 for(j = 1; j <= nHermite; ++j) <<
1079 { <<
1080 temp3 = temp2; <<
1081 temp2 = temp1; <<
1082 temp1 = nwt * std::sqrt(2.0 / j) * te <<
1083 std::sqrt(((G4double)(j - 1)) <<
1084 } <<
1085 temp = std::sqrt((G4double) 2 * nHermit <<
1086 nwt1 = nwt; <<
1087 nwt = nwt1 - temp1 / temp; <<
1088 <<
1089 if(std::fabs(nwt - nwt1) <= tolerance) <<
1090 { <<
1091 break; <<
1092 } <<
1093 } <<
1094 if(k > maxNumber) <<
1095 { <<
1096 G4Exception("G4Integrator<T,F>::Hermite <<
1097 FatalException, "Too many ( <<
1098 } <<
1099 fAbscissa[i - 1] = nwt; <<
1100 fWeight[i - 1] = 2.0 / (temp * temp); <<
1101 } <<
1102 <<
1103 // <<
1104 // Integral calculation <<
1105 // <<
1106 1034
1107 for(i = 0; i < fNumber; ++i) << 1035 for(j=1;j<=nHermite;j++)
1108 { << 1036 {
1109 integral += << 1037 temp3 = temp2 ;
1110 fWeight[i] * ((typeT.*f)(fAbscissa[i]) << 1038 temp2 = temp1 ;
1111 } << 1039 temp1 = newton*sqrt(2.0/j)*temp2 -
1112 delete[] fAbscissa; << 1040 sqrt(((G4double)(j - 1))/j)*temp3 ;
1113 delete[] fWeight; << 1041 }
1114 return integral; << 1042 temp = sqrt((G4double)2*nHermite)*temp2 ;
>> 1043 newton1 = newton ;
>> 1044 newton = newton1 - temp1/temp ;
>> 1045
>> 1046 if(fabs(newton - newton1) <= tolerance)
>> 1047 {
>> 1048 break ;
>> 1049 }
>> 1050 }
>> 1051 if(k > maxNumber)
>> 1052 {
>> 1053 G4Exception("Too many (>12) iterations in G4Integrator::Hermite") ;
>> 1054 }
>> 1055 fAbscissa[i-1] = newton ;
>> 1056 fWeight[i-1] = 2.0/(temp*temp) ;
>> 1057 }
>> 1058 //
>> 1059 // Integral calculation
>> 1060 //
>> 1061 for(i=0;i<fNumber;i++)
>> 1062 {
>> 1063 integral += fWeight[i]*( (typeT.*f)(fAbscissa[i]) +
>> 1064 (typeT.*f)(-fAbscissa[i]) ) ;
>> 1065 }
>> 1066 return integral ;
1115 } 1067 }
1116 1068
>> 1069
1117 ///////////////////////////////////////////// 1070 ////////////////////////////////////////////////////////////////////////
1118 // 1071 //
1119 // For use with 'this' pointer 1072 // For use with 'this' pointer
1120 1073
1121 template <class T, class F> << 1074 template <class T, class F>
1122 G4double G4Integrator<T, F>::Hermite(T* ptrT, << 1075 G4double G4Integrator::Hermite( T* ptrT, F f, G4int n)
1123 { 1076 {
1124 return Hermite(*ptrT, f, n); << 1077 return Hermite(*ptrT,f,n) ;
1125 } << 1078 }
1126 1079
1127 ///////////////////////////////////////////// 1080 ////////////////////////////////////////////////////////////////////////
1128 // 1081 //
1129 // For use with global scope f 1082 // For use with global scope f
1130 1083
1131 template <class T, class F> << 1084 G4double G4Integrator::Hermite( G4double (*f)(G4double), G4int nHermite)
1132 G4double G4Integrator<T, F>::Hermite(G4double <<
1133 { 1085 {
1134 const G4double tolerance = 1.0e-12; << 1086 const G4double tolerance = 1.0e-12 ;
1135 const G4int maxNumber = 12; << 1087 const G4int maxNumber = 12 ;
1136 << 1088
1137 G4int i, j, k; << 1089 G4int i, j, k ;
1138 G4double integral = 0.0; << 1090 G4double integral = 0.0 ;
1139 G4double nwt = 0., nwt1, temp1, temp2, << 1091 G4double newton, newton1, temp1, temp2, temp3, temp ;
1140 << 1092
1141 G4double piInMinusQ = << 1093 G4double piInMinusQ = pow(pi,-0.25) ; // 1.0/sqrt(sqrt(pi)) ??
1142 std::pow(CLHEP::pi, -0.25); // 1.0/std:: << 1094
1143 << 1095 G4int fNumber = (nHermite +1)/2 ;
1144 G4int fNumber = (nHermite + 1) / 2; << 1096 G4double* fAbscissa = new G4double[fNumber] ;
1145 G4double* fAbscissa = new G4double[fNumber] << 1097 G4double* fWeight = new G4double[fNumber] ;
1146 G4double* fWeight = new G4double[fNumber] << 1098
1147 << 1099 for(i=1;i<=fNumber;i++)
1148 for(i = 1; i <= fNumber; ++i) << 1100 {
1149 { << 1101 if(i == 1)
1150 if(i == 1) << 1102 {
1151 { << 1103 newton = sqrt((G4double)(2*nHermite + 1)) -
1152 nwt = std::sqrt((G4double)(2 * nHermite << 1104 1.85575001*pow((G4double)(2*nHermite + 1),-0.16666999) ;
1153 1.85575001 * std::pow((G4double)( << 1105 }
1154 } << 1106 else if(i == 2)
1155 else if(i == 2) << 1107 {
1156 { << 1108 newton -= 1.14001*pow((G4double)nHermite,0.425999)/newton ;
1157 nwt -= 1.14001 * std::pow((G4double) nH << 1109 }
1158 } << 1110 else if(i == 3)
1159 else if(i == 3) << 1111 {
1160 { << 1112 newton = 1.86002*newton - 0.86002*fAbscissa[0] ;
1161 nwt = 1.86002 * nwt - 0.86002 * fAbscis << 1113 }
1162 } << 1114 else if(i == 4)
1163 else if(i == 4) << 1115 {
1164 { << 1116 newton = 1.91001*newton - 0.91001*fAbscissa[1] ;
1165 nwt = 1.91001 * nwt - 0.91001 * fAbscis << 1117 }
1166 } << 1118 else
1167 else << 1119 {
1168 { << 1120 newton = 2.0*newton - fAbscissa[i - 3] ;
1169 nwt = 2.0 * nwt - fAbscissa[i - 3]; << 1121 }
1170 } << 1122 for(k=1;k<=maxNumber;k++)
1171 for(k = 1; k <= maxNumber; ++k) << 1123 {
1172 { << 1124 temp1 = piInMinusQ ;
1173 temp1 = piInMinusQ; << 1125 temp2 = 0.0 ;
1174 temp2 = 0.0; <<
1175 <<
1176 for(j = 1; j <= nHermite; ++j) <<
1177 { <<
1178 temp3 = temp2; <<
1179 temp2 = temp1; <<
1180 temp1 = nwt * std::sqrt(2.0 / j) * te <<
1181 std::sqrt(((G4double)(j - 1)) <<
1182 } <<
1183 temp = std::sqrt((G4double) 2 * nHermit <<
1184 nwt1 = nwt; <<
1185 nwt = nwt1 - temp1 / temp; <<
1186 <<
1187 if(std::fabs(nwt - nwt1) <= tolerance) <<
1188 { <<
1189 break; <<
1190 } <<
1191 } <<
1192 if(k > maxNumber) <<
1193 { <<
1194 G4Exception("G4Integrator<T,F>::Hermite <<
1195 "Too many (>12) iterations. <<
1196 } <<
1197 fAbscissa[i - 1] = nwt; <<
1198 fWeight[i - 1] = 2.0 / (temp * temp); <<
1199 } <<
1200 <<
1201 // <<
1202 // Integral calculation <<
1203 // <<
1204 1126
1205 for(i = 0; i < fNumber; ++i) << 1127 for(j=1;j<=nHermite;j++)
1206 { << 1128 {
1207 integral += fWeight[i] * ((*f)(fAbscissa[ << 1129 temp3 = temp2 ;
1208 } << 1130 temp2 = temp1 ;
1209 delete[] fAbscissa; << 1131 temp1 = newton*sqrt(2.0/j)*temp2 -
1210 delete[] fWeight; << 1132 sqrt(((G4double)(j - 1))/j)*temp3 ;
1211 return integral; << 1133 }
>> 1134 temp = sqrt((G4double)2*nHermite)*temp2 ;
>> 1135 newton1 = newton ;
>> 1136 newton = newton1 - temp1/temp ;
>> 1137
>> 1138 if(fabs(newton - newton1) <= tolerance)
>> 1139 {
>> 1140 break ;
>> 1141 }
>> 1142 }
>> 1143 if(k > maxNumber)
>> 1144 {
>> 1145 G4Exception("Too many (>12) iterations in G4Integrator::Hermite") ;
>> 1146 }
>> 1147 fAbscissa[i-1] = newton ;
>> 1148 fWeight[i-1] = 2.0/(temp*temp) ;
>> 1149 }
>> 1150 //
>> 1151 // Integral calculation
>> 1152 //
>> 1153 for(i=0;i<fNumber;i++)
>> 1154 {
>> 1155 integral += fWeight[i]*( (*f)(fAbscissa[i]) + (*f)(-fAbscissa[i]) ) ;
>> 1156 }
>> 1157 return integral ;
1212 } 1158 }
1213 1159
1214 ///////////////////////////////////////////// 1160 ////////////////////////////////////////////////////////////////////////////
1215 // 1161 //
1216 // Method involving Jacobi polynomials 1162 // Method involving Jacobi polynomials
1217 // 1163 //
1218 ///////////////////////////////////////////// 1164 ////////////////////////////////////////////////////////////////////////////
1219 // 1165 //
1220 // Gauss-Jacobi method for integration of ((1 1166 // Gauss-Jacobi method for integration of ((1-x)^alpha)*((1+x)^beta)*f(x)
1221 // from minus unit to plus unit . 1167 // from minus unit to plus unit .
1222 // 1168 //
1223 1169
1224 template <class T, class F> << 1170 template <class T, class F>
1225 G4double G4Integrator<T, F>::Jacobi(T& typeT, << 1171 G4double G4Integrator::Jacobi( T& typeT, F f, G4double alpha,
1226 G4double << 1172 G4double beta, G4int nJacobi)
1227 { << 1173 {
1228 const G4double tolerance = 1.0e-12; << 1174 const G4double tolerance = 1.0e-12 ;
1229 const G4double maxNumber = 12; << 1175 const G4double maxNumber = 12 ;
1230 G4int i, k, j; << 1176 G4int i, k, j ;
1231 G4double alphaBeta, alphaReduced, betaReduc << 1177 G4double alphaBeta, alphaReduced, betaReduced, root1, root2, root3 ;
1232 << 1178 G4double a, b, c, newton1, newton2, newton3, newton, temp, root, rootTemp ;
1233 G4double a, b, c, nwt1, nwt2, nwt3, nwt, te << 1179
1234 << 1180 G4int fNumber = nJacobi ;
1235 G4int fNumber = nJacobi; << 1181 G4double* fAbscissa = new G4double[fNumber] ;
1236 G4double* fAbscissa = new G4double[fNumber] << 1182 G4double* fWeight = new G4double[fNumber] ;
1237 G4double* fWeight = new G4double[fNumber] << 1183
1238 << 1184 for (i=1;i<=nJacobi;i++)
1239 for(i = 1; i <= nJacobi; ++i) << 1185 {
1240 { << 1186 if (i == 1)
1241 if(i == 1) << 1187 {
1242 { << 1188 alphaReduced = alpha/nJacobi ;
1243 alphaReduced = alpha / nJacobi; << 1189 betaReduced = beta/nJacobi ;
1244 betaReduced = beta / nJacobi; << 1190 root1 = (1.0+alpha)*(2.78002/(4.0+nJacobi*nJacobi)+
1245 root1 = (1.0 + alpha) * (2.78002 << 1191 0.767999*alphaReduced/nJacobi) ;
1246 0.767999 * alp << 1192 root2 = 1.0+1.48*alphaReduced+0.96002*betaReduced +
1247 root2 = 1.0 + 1.48 * alphaReduce << 1193 0.451998*alphaReduced*alphaReduced +
1248 0.451998 * alphaReduced * alpha << 1194 0.83001*alphaReduced*betaReduced ;
1249 0.83001 * alphaReduced * betaRe << 1195 root = 1.0-root1/root2 ;
1250 root = 1.0 - root1 / root2; << 1196 }
1251 } << 1197 else if (i == 2)
1252 else if(i == 2) << 1198 {
1253 { << 1199 root1=(4.1002+alpha)/((1.0+alpha)*(1.0+0.155998*alpha)) ;
1254 root1 = (4.1002 + alpha) / ((1.0 + alph << 1200 root2=1.0+0.06*(nJacobi-8.0)*(1.0+0.12*alpha)/nJacobi ;
1255 root2 = 1.0 + 0.06 * (nJacobi - 8.0) * << 1201 root3=1.0+0.012002*beta*(1.0+0.24997*fabs(alpha))/nJacobi ;
1256 root3 = << 1202 root -= (1.0-root)*root1*root2*root3 ;
1257 1.0 + 0.012002 * beta * (1.0 + 0.2499 << 1203 }
1258 root -= (1.0 - root) * root1 * root2 * << 1204 else if (i == 3)
1259 } << 1205 {
1260 else if(i == 3) << 1206 root1=(1.67001+0.27998*alpha)/(1.0+0.37002*alpha) ;
1261 { << 1207 root2=1.0+0.22*(nJacobi-8.0)/nJacobi ;
1262 root1 = (1.67001 + 0.27998 * alpha) / ( << 1208 root3=1.0+8.0*beta/((6.28001+beta)*nJacobi*nJacobi) ;
1263 root2 = 1.0 + 0.22 * (nJacobi - 8.0) / << 1209 root -= (fAbscissa[0]-root)*root1*root2*root3 ;
1264 root3 = 1.0 + 8.0 * beta / ((6.28001 + << 1210 }
1265 root -= (fAbscissa[0] - root) * root1 * << 1211 else if (i == nJacobi-1)
1266 } << 1212 {
1267 else if(i == nJacobi - 1) << 1213 root1=(1.0+0.235002*beta)/(0.766001+0.118998*beta) ;
1268 { << 1214 root2=1.0/(1.0+0.639002*(nJacobi-4.0)/(1.0+0.71001*(nJacobi-4.0))) ;
1269 root1 = (1.0 + 0.235002 * beta) / (0.76 << 1215 root3=1.0/(1.0+20.0*alpha/((7.5+alpha)*nJacobi*nJacobi)) ;
1270 root2 = 1.0 / (1.0 + 0.639002 * (nJacob << 1216 root += (root-fAbscissa[nJacobi-4])*root1*root2*root3 ;
1271 (1.0 + 0.71001 * << 1217 }
1272 root3 = 1.0 / (1.0 + 20.0 * alpha / ((7 << 1218 else if (i == nJacobi)
1273 root += (root - fAbscissa[nJacobi - 4]) << 1219 {
1274 } << 1220 root1 = (1.0+0.37002*beta)/(1.67001+0.27998*beta) ;
1275 else if(i == nJacobi) << 1221 root2 = 1.0/(1.0+0.22*(nJacobi-8.0)/nJacobi) ;
1276 { << 1222 root3 = 1.0/(1.0+8.0*alpha/((6.28002+alpha)*nJacobi*nJacobi)) ;
1277 root1 = (1.0 + 0.37002 * beta) / (1.670 << 1223 root += (root-fAbscissa[nJacobi-3])*root1*root2*root3 ;
1278 root2 = 1.0 / (1.0 + 0.22 * (nJacobi - << 1224 }
1279 root3 = << 1225 else
1280 1.0 / (1.0 + 8.0 * alpha / ((6.28002 << 1226 {
1281 root += (root - fAbscissa[nJacobi - 3]) << 1227 root = 3.0*fAbscissa[i-2]-3.0*fAbscissa[i-3]+fAbscissa[i-4] ;
1282 } << 1228 }
1283 else << 1229 alphaBeta = alpha + beta ;
1284 { << 1230 for (k=1;k<=maxNumber;k++)
1285 root = 3.0 * fAbscissa[i - 2] - 3.0 * f << 1231 {
1286 } << 1232 temp = 2.0 + alphaBeta ;
1287 alphaBeta = alpha + beta; << 1233 newton1 = (alpha-beta+temp*root)/2.0 ;
1288 for(k = 1; k <= maxNumber; ++k) << 1234 newton2 = 1.0 ;
1289 { << 1235 for (j=2;j<=nJacobi;j++)
1290 temp = 2.0 + alphaBeta; << 1236 {
1291 nwt1 = (alpha - beta + temp * root) / 2 << 1237 newton3 = newton2 ;
1292 nwt2 = 1.0; << 1238 newton2 = newton1 ;
1293 for(j = 2; j <= nJacobi; ++j) << 1239 temp = 2*j+alphaBeta ;
1294 { << 1240 a = 2*j*(j+alphaBeta)*(temp-2.0) ;
1295 nwt3 = nwt2; << 1241 b = (temp-1.0)*(alpha*alpha-beta*beta+temp*(temp-2.0)*root) ;
1296 nwt2 = nwt1; << 1242 c = 2.0*(j-1+alpha)*(j-1+beta)*temp ;
1297 temp = 2 * j + alphaBeta; << 1243 newton1 = (b*newton2-c*newton3)/a ;
1298 a = 2 * j * (j + alphaBeta) * (tem << 1244 }
1299 b = (temp - 1.0) * << 1245 newton = (nJacobi*(alpha - beta - temp*root)*newton1 +
1300 (alpha * alpha - beta * beta + te << 1246 2.0*(nJacobi + alpha)*(nJacobi + beta)*newton2)/
1301 c = 2.0 * (j - 1 + alpha) * (j - 1 << 1247 (temp*(1.0 - root*root)) ;
1302 nwt1 = (b * nwt2 - c * nwt3) / a; << 1248 rootTemp = root ;
1303 } << 1249 root = rootTemp - newton1/newton ;
1304 nwt = (nJacobi * (alpha - beta - temp * << 1250 if (fabs(root-rootTemp) <= tolerance)
1305 2.0 * (nJacobi + alpha) * (nJaco << 1251 {
1306 (temp * (1.0 - root * root)); << 1252 break ;
1307 rootTemp = root; << 1253 }
1308 root = rootTemp - nwt1 / nwt; << 1254 }
1309 if(std::fabs(root - rootTemp) <= tolera << 1255 if (k > maxNumber)
1310 { << 1256 {
1311 break; << 1257 G4Exception("Too many iterations (>12) in G4Integrator::Jacobi") ;
1312 } << 1258 }
1313 } << 1259 fAbscissa[i-1] = root ;
1314 if(k > maxNumber) << 1260 fWeight[i-1] = exp(GammaLogarithm((G4double)(alpha+nJacobi)) +
1315 { << 1261 GammaLogarithm((G4double)(beta+nJacobi)) -
1316 G4Exception("G4Integrator<T,F>::Jacobi( << 1262 GammaLogarithm((G4double)(nJacobi+1.0)) -
1317 FatalException, "Too many ( << 1263 GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0)))
1318 } << 1264 *temp*pow(2.0,alphaBeta)/(newton*newton2) ;
1319 fAbscissa[i - 1] = root; << 1265 }
1320 fWeight[i - 1] = << 1266 //
1321 std::exp(GammaLogarithm((G4double)(alph << 1267 // Calculation of the integral
1322 GammaLogarithm((G4double)(beta << 1268 //
1323 GammaLogarithm((G4double)(nJac << 1269 G4double integral = 0.0 ;
1324 GammaLogarithm((G4double)(nJac << 1270 for(i=0;i<fNumber;i++)
1325 temp * std::pow(2.0, alphaBeta) / (nwt << 1271 {
1326 } << 1272 integral += fWeight[i]*(typeT.*f)(fAbscissa[i]) ;
1327 << 1273 }
1328 // << 1274 return integral ;
1329 // Calculation of the integral <<
1330 // <<
1331 <<
1332 G4double integral = 0.0; <<
1333 for(i = 0; i < fNumber; ++i) <<
1334 { <<
1335 integral += fWeight[i] * (typeT.*f)(fAbsc <<
1336 } <<
1337 delete[] fAbscissa; <<
1338 delete[] fWeight; <<
1339 return integral; <<
1340 } 1275 }
1341 1276
>> 1277
1342 ///////////////////////////////////////////// 1278 /////////////////////////////////////////////////////////////////////////
1343 // 1279 //
1344 // For use with 'this' pointer 1280 // For use with 'this' pointer
1345 1281
1346 template <class T, class F> << 1282 template <class T, class F>
1347 G4double G4Integrator<T, F>::Jacobi(T* ptrT, << 1283 G4double G4Integrator::Jacobi( T* ptrT, F f, G4double alpha,
1348 G4int n) << 1284 G4double beta, G4int n)
1349 { 1285 {
1350 return Jacobi(*ptrT, f, alpha, beta, n); << 1286 return Jacobi(*ptrT,f,alpha,beta,n) ;
1351 } << 1287 }
1352 1288
1353 ///////////////////////////////////////////// 1289 /////////////////////////////////////////////////////////////////////////
1354 // 1290 //
1355 // For use with global scope f << 1291 // For use with global scope f
1356 1292
1357 template <class T, class F> << 1293 G4double G4Integrator::Jacobi( G4double (*f)(G4double), G4double alpha,
1358 G4double G4Integrator<T, F>::Jacobi(G4double << 1294 G4double beta, G4int nJacobi)
1359 G4double << 1295 {
1360 { << 1296 const G4double tolerance = 1.0e-12 ;
1361 const G4double tolerance = 1.0e-12; << 1297 const G4double maxNumber = 12 ;
1362 const G4double maxNumber = 12; << 1298 G4int i, k, j ;
1363 G4int i, k, j; << 1299 G4double alphaBeta, alphaReduced, betaReduced, root1, root2, root3 ;
1364 G4double alphaBeta, alphaReduced, betaReduc << 1300 G4double a, b, c, newton1, newton2, newton3, newton, temp, root, rootTemp ;
1365 << 1301
1366 G4double a, b, c, nwt1, nwt2, nwt3, nwt, te << 1302 G4int fNumber = nJacobi ;
1367 << 1303 G4double* fAbscissa = new G4double[fNumber] ;
1368 G4int fNumber = nJacobi; << 1304 G4double* fWeight = new G4double[fNumber] ;
1369 G4double* fAbscissa = new G4double[fNumber] << 1305
1370 G4double* fWeight = new G4double[fNumber] << 1306 for (i=1;i<=nJacobi;i++)
1371 << 1307 {
1372 for(i = 1; i <= nJacobi; ++i) << 1308 if (i == 1)
1373 { << 1309 {
1374 if(i == 1) << 1310 alphaReduced = alpha/nJacobi ;
1375 { << 1311 betaReduced = beta/nJacobi ;
1376 alphaReduced = alpha / nJacobi; << 1312 root1 = (1.0+alpha)*(2.78002/(4.0+nJacobi*nJacobi)+
1377 betaReduced = beta / nJacobi; << 1313 0.767999*alphaReduced/nJacobi) ;
1378 root1 = (1.0 + alpha) * (2.78002 << 1314 root2 = 1.0+1.48*alphaReduced+0.96002*betaReduced +
1379 0.767999 * alp << 1315 0.451998*alphaReduced*alphaReduced +
1380 root2 = 1.0 + 1.48 * alphaReduce << 1316 0.83001*alphaReduced*betaReduced ;
1381 0.451998 * alphaReduced * alpha << 1317 root = 1.0-root1/root2 ;
1382 0.83001 * alphaReduced * betaRe << 1318 }
1383 root = 1.0 - root1 / root2; << 1319 else if (i == 2)
1384 } << 1320 {
1385 else if(i == 2) << 1321 root1=(4.1002+alpha)/((1.0+alpha)*(1.0+0.155998*alpha)) ;
1386 { << 1322 root2=1.0+0.06*(nJacobi-8.0)*(1.0+0.12*alpha)/nJacobi ;
1387 root1 = (4.1002 + alpha) / ((1.0 + alph << 1323 root3=1.0+0.012002*beta*(1.0+0.24997*fabs(alpha))/nJacobi ;
1388 root2 = 1.0 + 0.06 * (nJacobi - 8.0) * << 1324 root -= (1.0-root)*root1*root2*root3 ;
1389 root3 = << 1325 }
1390 1.0 + 0.012002 * beta * (1.0 + 0.2499 << 1326 else if (i == 3)
1391 root -= (1.0 - root) * root1 * root2 * << 1327 {
1392 } << 1328 root1=(1.67001+0.27998*alpha)/(1.0+0.37002*alpha) ;
1393 else if(i == 3) << 1329 root2=1.0+0.22*(nJacobi-8.0)/nJacobi ;
1394 { << 1330 root3=1.0+8.0*beta/((6.28001+beta)*nJacobi*nJacobi) ;
1395 root1 = (1.67001 + 0.27998 * alpha) / ( << 1331 root -= (fAbscissa[0]-root)*root1*root2*root3 ;
1396 root2 = 1.0 + 0.22 * (nJacobi - 8.0) / << 1332 }
1397 root3 = 1.0 + 8.0 * beta / ((6.28001 + << 1333 else if (i == nJacobi-1)
1398 root -= (fAbscissa[0] - root) * root1 * << 1334 {
1399 } << 1335 root1=(1.0+0.235002*beta)/(0.766001+0.118998*beta) ;
1400 else if(i == nJacobi - 1) << 1336 root2=1.0/(1.0+0.639002*(nJacobi-4.0)/(1.0+0.71001*(nJacobi-4.0))) ;
1401 { << 1337 root3=1.0/(1.0+20.0*alpha/((7.5+alpha)*nJacobi*nJacobi)) ;
1402 root1 = (1.0 + 0.235002 * beta) / (0.76 << 1338 root += (root-fAbscissa[nJacobi-4])*root1*root2*root3 ;
1403 root2 = 1.0 / (1.0 + 0.639002 * (nJacob << 1339 }
1404 (1.0 + 0.71001 * << 1340 else if (i == nJacobi)
1405 root3 = 1.0 / (1.0 + 20.0 * alpha / ((7 << 1341 {
1406 root += (root - fAbscissa[nJacobi - 4]) << 1342 root1 = (1.0+0.37002*beta)/(1.67001+0.27998*beta) ;
1407 } << 1343 root2 = 1.0/(1.0+0.22*(nJacobi-8.0)/nJacobi) ;
1408 else if(i == nJacobi) << 1344 root3 = 1.0/(1.0+8.0*alpha/((6.28002+alpha)*nJacobi*nJacobi)) ;
1409 { << 1345 root += (root-fAbscissa[nJacobi-3])*root1*root2*root3 ;
1410 root1 = (1.0 + 0.37002 * beta) / (1.670 << 1346 }
1411 root2 = 1.0 / (1.0 + 0.22 * (nJacobi - << 1347 else
1412 root3 = << 1348 {
1413 1.0 / (1.0 + 8.0 * alpha / ((6.28002 << 1349 root = 3.0*fAbscissa[i-2]-3.0*fAbscissa[i-3]+fAbscissa[i-4] ;
1414 root += (root - fAbscissa[nJacobi - 3]) << 1350 }
1415 } << 1351 alphaBeta = alpha + beta ;
1416 else << 1352 for (k=1;k<=maxNumber;k++)
1417 { << 1353 {
1418 root = 3.0 * fAbscissa[i - 2] - 3.0 * f << 1354 temp = 2.0 + alphaBeta ;
1419 } << 1355 newton1 = (alpha-beta+temp*root)/2.0 ;
1420 alphaBeta = alpha + beta; << 1356 newton2 = 1.0 ;
1421 for(k = 1; k <= maxNumber; ++k) << 1357 for (j=2;j<=nJacobi;j++)
1422 { << 1358 {
1423 temp = 2.0 + alphaBeta; << 1359 newton3 = newton2 ;
1424 nwt1 = (alpha - beta + temp * root) / 2 << 1360 newton2 = newton1 ;
1425 nwt2 = 1.0; << 1361 temp = 2*j+alphaBeta ;
1426 for(j = 2; j <= nJacobi; ++j) << 1362 a = 2*j*(j+alphaBeta)*(temp-2.0) ;
1427 { << 1363 b = (temp-1.0)*(alpha*alpha-beta*beta+temp*(temp-2.0)*root) ;
1428 nwt3 = nwt2; << 1364 c = 2.0*(j-1+alpha)*(j-1+beta)*temp ;
1429 nwt2 = nwt1; << 1365 newton1 = (b*newton2-c*newton3)/a ;
1430 temp = 2 * j + alphaBeta; << 1366 }
1431 a = 2 * j * (j + alphaBeta) * (tem << 1367 newton = (nJacobi*(alpha - beta - temp*root)*newton1 +
1432 b = (temp - 1.0) * << 1368 2.0*(nJacobi + alpha)*(nJacobi + beta)*newton2)/
1433 (alpha * alpha - beta * beta + te << 1369 (temp*(1.0 - root*root)) ;
1434 c = 2.0 * (j - 1 + alpha) * (j - 1 << 1370 rootTemp = root ;
1435 nwt1 = (b * nwt2 - c * nwt3) / a; << 1371 root = rootTemp - newton1/newton ;
1436 } << 1372 if (fabs(root-rootTemp) <= tolerance)
1437 nwt = (nJacobi * (alpha - beta - temp * << 1373 {
1438 2.0 * (nJacobi + alpha) * (nJaco << 1374 break ;
1439 (temp * (1.0 - root * root)); << 1375 }
1440 rootTemp = root; << 1376 }
1441 root = rootTemp - nwt1 / nwt; << 1377 if (k > maxNumber)
1442 if(std::fabs(root - rootTemp) <= tolera << 1378 {
1443 { << 1379 G4Exception("Too many iterations (>12) in G4Integrator::Jacobi") ;
1444 break; << 1380 }
1445 } << 1381 fAbscissa[i-1] = root ;
1446 } << 1382 fWeight[i-1] = exp(GammaLogarithm((G4double)(alpha+nJacobi)) +
1447 if(k > maxNumber) << 1383 GammaLogarithm((G4double)(beta+nJacobi)) -
1448 { << 1384 GammaLogarithm((G4double)(nJacobi+1.0)) -
1449 G4Exception("G4Integrator<T,F>::Jacobi( << 1385 GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0)))
1450 "Too many (>12) iterations. << 1386 *temp*pow(2.0,alphaBeta)/(newton*newton2) ;
1451 } << 1387 }
1452 fAbscissa[i - 1] = root; << 1388 //
1453 fWeight[i - 1] = << 1389 // Calculation of the integral
1454 std::exp(GammaLogarithm((G4double)(alph << 1390 //
1455 GammaLogarithm((G4double)(beta << 1391 G4double integral = 0.0 ;
1456 GammaLogarithm((G4double)(nJac << 1392 for(i=0;i<fNumber;i++)
1457 GammaLogarithm((G4double)(nJac << 1393 {
1458 temp * std::pow(2.0, alphaBeta) / (nwt << 1394 integral += fWeight[i]*(*f)(fAbscissa[i]) ;
1459 } << 1395 }
>> 1396 return integral ;
>> 1397 }
1460 1398
1461 // <<
1462 // Calculation of the integral <<
1463 // <<
1464 1399
1465 G4double integral = 0.0; <<
1466 for(i = 0; i < fNumber; ++i) <<
1467 { <<
1468 integral += fWeight[i] * (*f)(fAbscissa[i <<
1469 } <<
1470 delete[] fAbscissa; <<
1471 delete[] fWeight; <<
1472 return integral; <<
1473 } <<
1474 1400
1475 // 1401 //
1476 // 1402 //
1477 ///////////////////////////////////////////// 1403 ///////////////////////////////////////////////////////////////////
>> 1404
>> 1405
>> 1406
>> 1407
1478 1408