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