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