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