Geant4-11
G4Integrator.icc
Go to the documentation of this file.
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
41template <class T, class F>
42G4double 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
69template <class T, class F>
70G4double 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
98template <class T, class F>
99G4double 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
129template <class T, class F>
130G4double 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
147template <class T, class F>
148G4double 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
157template <class T, class F>
158G4double 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
175template <class T, class F>
176void 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
204template <class T, class F>
205void 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//
215template <class T, class F>
216void 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
249template <class T, class F>
250G4double 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
264template <class T, class F>
265G4double 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
276template <class T, class F>
277G4double 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
309template <class T, class F>
310G4double 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
374template <class T, class F>
375G4double 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
385template <class T, class F>
386G4double 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
456template <class T, class F>
457G4double 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
486template <class T, class F>
487G4double 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
496template <class T, class F>
497G4double 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
532template <class T, class F>
533G4double 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
622template <class T, class F>
623G4double 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
632template <class T, class F>
633G4double 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
729template <class T, class F>
730G4double 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
766template <class T, class F>
767G4double 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
777template <class T, class F>
778G4double 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
823template <class T, class F>
824G4double 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
905template <class T, class F>
906G4double 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
916template <class T, class F>
917G4double 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
1002template <class T, class F>
1003G4double 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
1033template <class T, class F>
1034G4double 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
1121template <class T, class F>
1122G4double 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
1131template <class T, class F>
1132G4double 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
1224template <class T, class F>
1225G4double 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
1346template <class T, class F>
1347G4double 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
1357template <class T, class F>
1358G4double 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///////////////////////////////////////////////////////////////////
Generated by doxygen 1.9.3