| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4Integrator inline methods implementation
27 //
28 // Author: V.Grichine, 04.09.1999 - First impl
29 // G4SimpleIntegration class with H.P.
30 // E.TCherniaev advises
31 // -------------------------------------------
32
33 //////////////////////////////////////////////
34 //
35 // Sympson integration method
36 //
37 //////////////////////////////////////////////
38 //
39 // Integration of class member functions T::f
40
41 template <class T, class F>
42 G4double G4Integrator<T, F>::Simpson(T& typeT,
43 G4double
44 {
45 G4int i;
46 G4double step = (xFinal - xInitial) / itera
47 G4double x = xInitial;
48 G4double xPlus = xInitial + 0.5 * step;
49 G4double mean = ((typeT.*f)(xInitial) + (ty
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
67 // Convenient to use with 'this' pointer
68
69 template <class T, class F>
70 G4double G4Integrator<T, F>::Simpson(T* ptrT,
71 G4double
72 {
73 G4int i;
74 G4double step = (xFinal - xInitial) / itera
75 G4double x = xInitial;
76 G4double xPlus = xInitial + 0.5 * step;
77 G4double mean = ((ptrT->*f)(xInitial) + (pt
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
95 // Convenient to use, when function f is defin
96 // program
97
98 template <class T, class F>
99 G4double G4Integrator<T, F>::Simpson(G4double
100 G4double
101 {
102 G4int i;
103 G4double step = (xFinal - xInitial) / itera
104 G4double x = xInitial;
105 G4double xPlus = xInitial + 0.5 * step;
106 G4double mean = ((*f)(xInitial) + (*f)(xFin
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
131 G4double xF
132 {
133 static const G4double root = 1.0 / std::sqrt
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)
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
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 (*
159 G4double xF
160 {
161 static const G4double root = 1.0 / std::sqrt
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)
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,
177 G4double x
178 G4double&
179 {
180 if(depth > 100)
181 {
182 G4cout << "G4Integrator<T,F>::AdaptGauss:
183 G4cout << "Function varies too rapidly to
184 << G4endl;
185
186 return;
187 }
188 G4double xMean = (xInitial + xFinal) / 2
189 G4double leftHalf = Gauss(typeT, f, xInitia
190 G4double rightHalf = Gauss(typeT, f, xMean,
191 G4double full = Gauss(typeT, f, xInitia
192 if(std::fabs(leftHalf + rightHalf - full) <
193 {
194 sum += full;
195 }
196 else
197 {
198 ++depth;
199 AdaptGauss(typeT, f, xInitial, xMean, fTol
200 AdaptGauss(typeT, f, xMean, xFinal, fToler
201 }
202 }
203
204 template <class T, class F>
205 void G4Integrator<T, F>::AdaptGauss(T* ptrT, F
206 G4double x
207 G4double&
208 {
209 AdaptGauss(*ptrT, f, xInitial, xFinal, fTole
210 }
211
212 //////////////////////////////////////////////
213 //
214 //
215 template <class T, class F>
216 void G4Integrator<T, F>::AdaptGauss(G4double (
217 G4double x
218 G4double&
219 {
220 if(depth > 100)
221 {
222 G4cout << "G4SimpleIntegration::AdaptGauss
223 G4cout << "Function varies too rapidly to
224 << G4endl;
225
226 return;
227 }
228 G4double xMean = (xInitial + xFinal) / 2
229 G4double leftHalf = Gauss(f, xInitial, xMea
230 G4double rightHalf = Gauss(f, xMean, xFinal)
231 G4double full = Gauss(f, xInitial, xFin
232 if(std::fabs(leftHalf + rightHalf - full) <
233 {
234 sum += full;
235 }
236 else
237 {
238 ++depth;
239 AdaptGauss(f, xInitial, xMean, fTolerance,
240 AdaptGauss(f, xMean, xFinal, fTolerance, s
241 }
242 }
243
244 //////////////////////////////////////////////
245 //
246 // Adaptive Gauss integration with accuracy 'e
247 // Convenient for using with class object type
248
249 template <class T, class F>
250 G4double G4Integrator<T, F>::AdaptiveGauss(T&
251 G4d
252 {
253 G4int depth = 0;
254 G4double sum = 0.0;
255 AdaptGauss(typeT, f, xInitial, xFinal, e, su
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*
266 G4d
267 {
268 return AdaptiveGauss(*ptrT, f, xInitial, xFi
269 }
270
271 //////////////////////////////////////////////
272 //
273 // Adaptive Gauss integration with accuracy 'e
274 // Convenient for using with global scope func
275
276 template <class T, class F>
277 G4double G4Integrator<T, F>::AdaptiveGauss(G4d
278 G4d
279 G4d
280 {
281 G4int depth = 0;
282 G4double sum = 0.0;
283 AdaptGauss(f, xInitial, xFinal, e, sum, dept
284 return sum;
285 }
286
287 //////////////////////////////////////////////
288 // Gauss integration methods involving ortogon
289 //////////////////////////////////////////////
290 //
291 // Methods involving Legendre polynomials
292 //
293 //////////////////////////////////////////////
294 //
295 // The value nLegendre set the accuracy requir
296 // where the function pFunction will be evalua
297 // The function creates the arrays for absciss
298 // in Gauss-Legendre quadrature method.
299 // The values a and b are the limits of integr
300 // nLegendre MUST BE EVEN !!!
301 // Returns the integral of the function f betw
302 // Gauss-Legendre integration: the function is
303 // 2*fNumber times at interior points in the r
304 // Since the weights and abscissas are, in thi
305 // the midpoint of the range of integration, t
306 // fNumber distinct values of each.
307 // Convenient for using with some class object
308
309 template <class T, class F>
310 G4double G4Integrator<T, F>::Legendre(T& typeT
311 G4int nL
312 {
313 G4double nwt, nwt1, temp1, temp2, temp3, tem
314 G4double xDiff, xMean, dx, integral;
315
316 const G4double tolerance = 1.6e-10;
317 G4int i, j, k = nLegendre;
318 G4int fNumber = (nLegendre + 1) / 2;
319
320 if(2 * fNumber != k)
321 {
322 G4Exception("G4Integrator<T,F>::Legendre(T
323 FatalException, "Invalid (odd)
324 }
325
326 G4double* fAbscissa = new G4double[fNumber];
327 G4double* fWeight = new G4double[fNumber];
328
329 for(i = 1; i <= fNumber; ++i) // Loop over
330 {
331 nwt = std::cos(CLHEP::pi * (i - 0.25) /
332 (k + 0.5)); // Initial roo
333
334 do // loop of Newton's method
335 {
336 temp1 = 1.0;
337 temp2 = 0.0;
338 for(j = 1; j <= k; ++j)
339 {
340 temp3 = temp2;
341 temp2 = temp1;
342 temp1 = ((2.0 * j - 1.0) * nwt * temp2
343 }
344 temp = k * (nwt * temp1 - temp2) / (nwt
345 nwt1 = nwt;
346 nwt = nwt1 - temp1 / temp; // Newton's
347 } while(std::fabs(nwt - nwt1) > tolerance)
348
349 fAbscissa[fNumber - i] = nwt;
350 fWeight[fNumber - i] = 2.0 / ((1.0 - nwt
351 }
352
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
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,
376 G4int nL
377 {
378 return Legendre(*ptrT, f, a, b, nLegendre);
379 }
380
381 //////////////////////////////////////////////
382 //
383 // Convenient for using with global scope func
384
385 template <class T, class F>
386 G4double G4Integrator<T, F>::Legendre(G4double
387 G4double
388 {
389 G4double nwt, nwt1, temp1, temp2, temp3, tem
390 G4double xDiff, xMean, dx, integral;
391
392 const G4double tolerance = 1.6e-10;
393 G4int i, j, k = nLegendre;
394 G4int fNumber = (nLegendre + 1) / 2;
395
396 if(2 * fNumber != k)
397 {
398 G4Exception("G4Integrator<T,F>::Legendre(.
399 FatalException, "Invalid (odd)
400 }
401
402 G4double* fAbscissa = new G4double[fNumber];
403 G4double* fWeight = new G4double[fNumber];
404
405 for(i = 1; i <= fNumber; i++) // Loop over
406 {
407 nwt = std::cos(CLHEP::pi * (i - 0.25) /
408 (k + 0.5)); // Initial roo
409
410 do // loop of Newton's method
411 {
412 temp1 = 1.0;
413 temp2 = 0.0;
414 for(j = 1; j <= k; ++j)
415 {
416 temp3 = temp2;
417 temp2 = temp1;
418 temp1 = ((2.0 * j - 1.0) * nwt * temp2
419 }
420 temp = k * (nwt * temp1 - temp2) / (nwt
421 nwt1 = nwt;
422 nwt = nwt1 - temp1 / temp; // Newton's
423 } while(std::fabs(nwt - nwt1) > tolerance)
424
425 fAbscissa[fNumber - i] = nwt;
426 fWeight[fNumber - i] = 2.0 / ((1.0 - nwt
427 }
428
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)
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
450 // by ten point Gauss-Legendre integration: th
451 // ten times at interior points in the range o
452 // and abscissas are, in this case, symmetric
453 // range of integration, there are actually on
454 // Convenient for using with class object type
455
456 template <class T, class F>
457 G4double G4Integrator<T, F>::Legendre10(T& typ
458 {
459 G4int i;
460 G4double xDiff, xMean, dx, integral;
461
462 // From Abramowitz M., Stegan I.A. 1964 , Ha
463
464 static const G4double abscissa[] = { 0.14887
465 0.67940
466 0.97390
467
468 static const G4double weight[] = { 0.2955242
469 0.2190863
470 0.0666713
471 xMean = 0.5 * (a +
472 xDiff = 0.5 * (b -
473 integral = 0.0;
474 for(i = 0; i < 5; ++i)
475 {
476 dx = xDiff * abscissa[i];
477 integral += weight[i] * ((typeT.*f)(xMean
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* ptr
488 {
489 return Legendre10(*ptrT, f, a, b);
490 }
491
492 //////////////////////////////////////////////
493 //
494 // Convenient for using with global scope func
495
496 template <class T, class F>
497 G4double G4Integrator<T, F>::Legendre10(G4doub
498 G4doub
499 {
500 G4int i;
501 G4double xDiff, xMean, dx, integral;
502
503 // From Abramowitz M., Stegan I.A. 1964 , Ha
504
505 static const G4double abscissa[] = { 0.14887
506 0.67940
507 0.97390
508
509 static const G4double weight[] = { 0.2955242
510 0.2190863
511 0.0666713
512 xMean = 0.5 * (a +
513 xDiff = 0.5 * (b -
514 integral = 0.0;
515 for(i = 0; i < 5; ++i)
516 {
517 dx = xDiff * abscissa[i];
518 integral += weight[i] * ((*f)(xMean + dx)
519 }
520 return integral *= xDiff;
521 }
522
523 //////////////////////////////////////////////
524 //
525 // Returns the integral of the function to be
526 // by 96 point Gauss-Legendre integration: the
527 // ten Times at interior points in the range o
528 // and abscissas are, in this case, symmetric
529 // range of integration, there are actually on
530 // Convenient for using with some class object
531
532 template <class T, class F>
533 G4double G4Integrator<T, F>::Legendre96(T& typ
534 {
535 G4int i;
536 G4double xDiff, xMean, dx, integral;
537
538 // From Abramowitz M., Stegan I.A. 1964 , Ha
539
540 static const G4double abscissa[] = {
541 0.016276744849602969579, 0.048812985136049
542 0.081297495464425558994, 0.113695850110665
543 0.145973714654896941989, 0.178096882367618
544
545 0.210031310460567203603, 0.241743156163840
546 0.273198812591049141487, 0.304364944354496
547 0.335208522892625422616, 0.365696861472313
548
549 0.395797649828908603285, 0.425478988407300
550 0.454709422167743008636, 0.483457973920596
551 0.511694177154667673586, 0.539388108324357
552
553 0.566510418561397168404, 0.593032364777572
554 0.618925840125468570386, 0.644163403784967
555 0.668718310043916153953, 0.692564536642171
556
557 0.715676812348967626225, 0.738030643744400
558 0.759602341176647498703, 0.780369043867433
559 0.800308744139140817229, 0.819400310737931
560
561 0.837623511228187121494, 0.854959033434601
562 0.871388505909296502874, 0.886894517402420
563 0.901460635315852341319, 0.915071423120898
564
565 0.927712456722308690965, 0.939370339752755
566 0.950032717784437635756, 0.959688291448742
567 0.968326828463264212174, 0.975939174585136
568
569 0.982517263563014677447, 0.988054126329623
570 0.992543900323762624572, 0.995981842987209
571 0.998364375863181677724, 0.999689503883230
572 };
573
574 static const G4double weight[] = {
575 0.032550614492363166242, 0.032516118713868
576 0.032447163714064269364, 0.032343822568575
577 0.032206204794030250669, 0.032034456231992
578
579 0.031828758894411006535, 0.031589330770727
580 0.031316425596862355813, 0.031010332586313
581 0.030671376123669149014, 0.030299915420827
582
583 0.029896344136328385984, 0.029461089958167
584 0.028994614150555236543, 0.028497411065085
585 0.027970007616848334440, 0.027412962726029
586
587 0.026826866725591762198, 0.026212340735672
588 0.025570036005349361499, 0.024900633222483
589 0.024204841792364691282, 0.023483399085926
590
591 0.022737069658329374001, 0.021966644438744
592 0.021172939892191298988, 0.020356797154333
593 0.019519081140145022410, 0.018660679627411
594
595 0.017782502316045260838, 0.016885479864245
596 0.015970562902562291381, 0.015038721026994
597 0.014090941772314860916, 0.013128229566961
598
599 0.012151604671088319635, 0.011162102099838
600 0.010160770535008415758, 0.009148671230783
601 0.008126876925698759217, 0.007096470791153
602
603 0.006058545504235961683, 0.005014202742927
604 0.003964554338444686674, 0.002910731817934
605 0.001853960788946921732, 0.000796792065552
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
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* ptr
624 {
625 return Legendre96(*ptrT, f, a, b);
626 }
627
628 //////////////////////////////////////////////
629 //
630 // Convenient for using with global scope func
631
632 template <class T, class F>
633 G4double G4Integrator<T, F>::Legendre96(G4doub
634 G4doub
635 {
636 G4int i;
637 G4double xDiff, xMean, dx, integral;
638
639 // From Abramowitz M., Stegan I.A. 1964 , Ha
640
641 static const G4double abscissa[] = {
642 0.016276744849602969579, 0.048812985136049
643 0.081297495464425558994, 0.113695850110665
644 0.145973714654896941989, 0.178096882367618
645
646 0.210031310460567203603, 0.241743156163840
647 0.273198812591049141487, 0.304364944354496
648 0.335208522892625422616, 0.365696861472313
649
650 0.395797649828908603285, 0.425478988407300
651 0.454709422167743008636, 0.483457973920596
652 0.511694177154667673586, 0.539388108324357
653
654 0.566510418561397168404, 0.593032364777572
655 0.618925840125468570386, 0.644163403784967
656 0.668718310043916153953, 0.692564536642171
657
658 0.715676812348967626225, 0.738030643744400
659 0.759602341176647498703, 0.780369043867433
660 0.800308744139140817229, 0.819400310737931
661
662 0.837623511228187121494, 0.854959033434601
663 0.871388505909296502874, 0.886894517402420
664 0.901460635315852341319, 0.915071423120898
665
666 0.927712456722308690965, 0.939370339752755
667 0.950032717784437635756, 0.959688291448742
668 0.968326828463264212174, 0.975939174585136
669
670 0.982517263563014677447, 0.988054126329623
671 0.992543900323762624572, 0.995981842987209
672 0.998364375863181677724, 0.999689503883230
673 };
674
675 static const G4double weight[] = {
676 0.032550614492363166242, 0.032516118713868
677 0.032447163714064269364, 0.032343822568575
678 0.032206204794030250669, 0.032034456231992
679
680 0.031828758894411006535, 0.031589330770727
681 0.031316425596862355813, 0.031010332586313
682 0.030671376123669149014, 0.030299915420827
683
684 0.029896344136328385984, 0.029461089958167
685 0.028994614150555236543, 0.028497411065085
686 0.027970007616848334440, 0.027412962726029
687
688 0.026826866725591762198, 0.026212340735672
689 0.025570036005349361499, 0.024900633222483
690 0.024204841792364691282, 0.023483399085926
691
692 0.022737069658329374001, 0.021966644438744
693 0.021172939892191298988, 0.020356797154333
694 0.019519081140145022410, 0.018660679627411
695
696 0.017782502316045260838, 0.016885479864245
697 0.015970562902562291381, 0.015038721026994
698 0.014090941772314860916, 0.013128229566961
699
700 0.012151604671088319635, 0.011162102099838
701 0.010160770535008415758, 0.009148671230783
702 0.008126876925698759217, 0.007096470791153
703
704 0.006058545504235961683, 0.005014202742927
705 0.003964554338444686674, 0.002910731817934
706 0.001853960788946921732, 0.000796792065552
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)
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
726 // quadrature method.
727 // Convenient for using with class object type
728
729 template <class T, class F>
730 G4double G4Integrator<T, F>::Chebyshev(T& type
731 G4int n
732 {
733 G4int i;
734 G4double xDiff, xMean, dx, integral = 0.0;
735
736 G4int fNumber = nChebyshev; // Try to
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 - fAbscis
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
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
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(G4doubl
779 G4doubl
780 {
781 G4int i;
782 G4double xDiff, xMean, dx, integral = 0.0;
783
784 G4int fNumber = nChebyshev; // Try to
785 G4double cof = CLHEP::pi / fNumber;
786 G4double* fAbscissa = new G4double[fNumber];
787 G4double* fWeight = new G4double[fNumber];
788 for(i = 0; i < fNumber; ++i)
789 {
790 fAbscissa[i] = std::cos(cof * (i + 0.5));
791 fWeight[i] = cof * std::sqrt(1 - fAbscis
792 }
793
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(
817 // The value of nLaguerre sets the accuracy.
818 // The function creates arrays fAbscissa[0,..,
819 // fWeight[0,..,nLaguerre-1] .
820 // Convenient for using with class object 'typ
821 // (T::f)
822
823 template <class T, class F>
824 G4double G4Integrator<T, F>::Laguerre(T& typeT
825 G4int nL
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,
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
838 {
839 if(i == 1)
840 {
841 nwt = (1.0 + alpha) * (3.0 + 0.92 * alph
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.
847 }
848 else
849 {
850 cofi = i - 2;
851 nwt += ((1.0 + 2.55 * cofi) / (1.9 * cof
852 1.26 * cofi * alpha / (1.0 + 3.5
853 (nwt - fAbscissa[i - 3]) / (1.0 +
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 -
866 }
867 temp = (fNumber * temp1 - (fNumber + alp
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
879 FatalException, "Too many (>
880 }
881
882 fAbscissa[i - 1] = nwt;
883 fWeight[i - 1] = -std::exp(GammaLogarith
884 GammaLogarithm(
885 (temp * fNumber * temp2);
886 }
887
888 //
889 // Integral evaluation
890 //
891
892 for(i = 0; i < fNumber; ++i)
893 {
894 integral += fWeight[i] * (typeT.*f)(fAbsci
895 }
896 delete[] fAbscissa;
897 delete[] fWeight;
898 return integral;
899 }
900
901 //////////////////////////////////////////////
902 //
903 //
904
905 template <class T, class F>
906 G4double G4Integrator<T, F>::Laguerre(T* ptrT,
907 G4int nL
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
918 G4int nL
919 {
920 const G4double tolerance = 1.0e-10;
921 const G4int maxNumber = 12;
922 G4int i, j, k;
923 G4double nwt = 0., nwt1, temp1, temp2,
924 G4double integral = 0.0;
925
926 G4int fNumber = nLaguerre;
927 G4double* fAbscissa = new G4double[fNumber];
928 G4double* fWeight = new G4double[fNumber];
929
930 for(i = 1; i <= fNumber; ++i) // Loop over
931 {
932 if(i == 1)
933 {
934 nwt = (1.0 + alpha) * (3.0 + 0.92 * alph
935 (1.0 + 2.4 * fNumber + 1.8 * alpha
936 }
937 else if(i == 2)
938 {
939 nwt += (15.0 + 6.25 * alpha) / (1.0 + 0.
940 }
941 else
942 {
943 cofi = i - 2;
944 nwt += ((1.0 + 2.55 * cofi) / (1.9 * cof
945 1.26 * cofi * alpha / (1.0 + 3.5
946 (nwt - fAbscissa[i - 3]) / (1.0 +
947 }
948 for(k = 1; k <= maxNumber; ++k)
949 {
950 temp1 = 1.0;
951 temp2 = 0.0;
952
953 for(j = 1; j <= fNumber; ++j)
954 {
955 temp3 = temp2;
956 temp2 = temp1;
957 temp1 =
958 ((2 * j - 1 + alpha - nwt) * temp2 -
959 }
960 temp = (fNumber * temp1 - (fNumber + alp
961 nwt1 = nwt;
962 nwt = nwt1 - temp1 / temp;
963
964 if(std::fabs(nwt - nwt1) <= tolerance)
965 {
966 break;
967 }
968 }
969 if(k > maxNumber)
970 {
971 G4Exception("G4Integrator<T,F>::Laguerre
972 "Too many (>12) iterations."
973 }
974
975 fAbscissa[i - 1] = nwt;
976 fWeight[i - 1] = -std::exp(GammaLogarith
977 GammaLogarithm(
978 (temp * fNumber * temp2);
979 }
980
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
997 // Returns the value ln(Gamma(xx) for xx > 0.
998 // xx > 1. For 0 < xx < 1. the reflection form
999 // (Adapted from Numerical Recipes in C)
1000 //
1001
1002 template <class T, class F>
1003 G4double G4Integrator<T, F>::GammaLogarithm(G
1004 {
1005 static const G4double cof[6] = { 76.1800917
1006 24.0140982
1007 0.12086509
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 *
1020 }
1021
1022 /////////////////////////////////////////////
1023 //
1024 // Method involving Hermite polynomials
1025 //
1026 /////////////////////////////////////////////
1027 //
1028 //
1029 // Gauss-Hermite method for integration of st
1030 // from minus infinity to plus infinity .
1031 //
1032
1033 template <class T, class F>
1034 G4double G4Integrator<T, F>::Hermite(T& typeT
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,
1042
1043 G4double piInMinusQ =
1044 std::pow(CLHEP::pi, -0.25); // 1.0/std::
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
1055 1.85575001 * std::pow((G4double)(
1056 }
1057 else if(i == 2)
1058 {
1059 nwt -= 1.14001 * std::pow((G4double) nH
1060 }
1061 else if(i == 3)
1062 {
1063 nwt = 1.86002 * nwt - 0.86002 * fAbscis
1064 }
1065 else if(i == 4)
1066 {
1067 nwt = 1.91001 * nwt - 0.91001 * fAbscis
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) * te
1083 std::sqrt(((G4double)(j - 1))
1084 }
1085 temp = std::sqrt((G4double) 2 * nHermit
1086 nwt1 = nwt;
1087 nwt = nwt1 - temp1 / temp;
1088
1089 if(std::fabs(nwt - nwt1) <= tolerance)
1090 {
1091 break;
1092 }
1093 }
1094 if(k > maxNumber)
1095 {
1096 G4Exception("G4Integrator<T,F>::Hermite
1097 FatalException, "Too many (
1098 }
1099 fAbscissa[i - 1] = nwt;
1100 fWeight[i - 1] = 2.0 / (temp * temp);
1101 }
1102
1103 //
1104 // Integral calculation
1105 //
1106
1107 for(i = 0; i < fNumber; ++i)
1108 {
1109 integral +=
1110 fWeight[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,
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
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,
1140
1141 G4double piInMinusQ =
1142 std::pow(CLHEP::pi, -0.25); // 1.0/std::
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
1153 1.85575001 * std::pow((G4double)(
1154 }
1155 else if(i == 2)
1156 {
1157 nwt -= 1.14001 * std::pow((G4double) nH
1158 }
1159 else if(i == 3)
1160 {
1161 nwt = 1.86002 * nwt - 0.86002 * fAbscis
1162 }
1163 else if(i == 4)
1164 {
1165 nwt = 1.91001 * nwt - 0.91001 * fAbscis
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) * te
1181 std::sqrt(((G4double)(j - 1))
1182 }
1183 temp = std::sqrt((G4double) 2 * nHermit
1184 nwt1 = nwt;
1185 nwt = nwt1 - temp1 / temp;
1186
1187 if(std::fabs(nwt - nwt1) <= tolerance)
1188 {
1189 break;
1190 }
1191 }
1192 if(k > maxNumber)
1193 {
1194 G4Exception("G4Integrator<T,F>::Hermite
1195 "Too many (>12) iterations.
1196 }
1197 fAbscissa[i - 1] = nwt;
1198 fWeight[i - 1] = 2.0 / (temp * temp);
1199 }
1200
1201 //
1202 // Integral calculation
1203 //
1204
1205 for(i = 0; i < fNumber; ++i)
1206 {
1207 integral += fWeight[i] * ((*f)(fAbscissa[
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
1221 // from minus unit to plus unit .
1222 //
1223
1224 template <class T, class F>
1225 G4double G4Integrator<T, F>::Jacobi(T& typeT,
1226 G4double
1227 {
1228 const G4double tolerance = 1.0e-12;
1229 const G4double maxNumber = 12;
1230 G4int i, k, j;
1231 G4double alphaBeta, alphaReduced, betaReduc
1232
1233 G4double a, b, c, nwt1, nwt2, nwt3, nwt, te
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
1246 0.767999 * alp
1247 root2 = 1.0 + 1.48 * alphaReduce
1248 0.451998 * alphaReduced * alpha
1249 0.83001 * alphaReduced * betaRe
1250 root = 1.0 - root1 / root2;
1251 }
1252 else if(i == 2)
1253 {
1254 root1 = (4.1002 + alpha) / ((1.0 + alph
1255 root2 = 1.0 + 0.06 * (nJacobi - 8.0) *
1256 root3 =
1257 1.0 + 0.012002 * beta * (1.0 + 0.2499
1258 root -= (1.0 - root) * root1 * root2 *
1259 }
1260 else if(i == 3)
1261 {
1262 root1 = (1.67001 + 0.27998 * alpha) / (
1263 root2 = 1.0 + 0.22 * (nJacobi - 8.0) /
1264 root3 = 1.0 + 8.0 * beta / ((6.28001 +
1265 root -= (fAbscissa[0] - root) * root1 *
1266 }
1267 else if(i == nJacobi - 1)
1268 {
1269 root1 = (1.0 + 0.235002 * beta) / (0.76
1270 root2 = 1.0 / (1.0 + 0.639002 * (nJacob
1271 (1.0 + 0.71001 *
1272 root3 = 1.0 / (1.0 + 20.0 * alpha / ((7
1273 root += (root - fAbscissa[nJacobi - 4])
1274 }
1275 else if(i == nJacobi)
1276 {
1277 root1 = (1.0 + 0.37002 * beta) / (1.670
1278 root2 = 1.0 / (1.0 + 0.22 * (nJacobi -
1279 root3 =
1280 1.0 / (1.0 + 8.0 * alpha / ((6.28002
1281 root += (root - fAbscissa[nJacobi - 3])
1282 }
1283 else
1284 {
1285 root = 3.0 * fAbscissa[i - 2] - 3.0 * f
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
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) * (tem
1299 b = (temp - 1.0) *
1300 (alpha * alpha - beta * beta + te
1301 c = 2.0 * (j - 1 + alpha) * (j - 1
1302 nwt1 = (b * nwt2 - c * nwt3) / a;
1303 }
1304 nwt = (nJacobi * (alpha - beta - temp *
1305 2.0 * (nJacobi + alpha) * (nJaco
1306 (temp * (1.0 - root * root));
1307 rootTemp = root;
1308 root = rootTemp - nwt1 / nwt;
1309 if(std::fabs(root - rootTemp) <= tolera
1310 {
1311 break;
1312 }
1313 }
1314 if(k > maxNumber)
1315 {
1316 G4Exception("G4Integrator<T,F>::Jacobi(
1317 FatalException, "Too many (
1318 }
1319 fAbscissa[i - 1] = root;
1320 fWeight[i - 1] =
1321 std::exp(GammaLogarithm((G4double)(alph
1322 GammaLogarithm((G4double)(beta
1323 GammaLogarithm((G4double)(nJac
1324 GammaLogarithm((G4double)(nJac
1325 temp * std::pow(2.0, alphaBeta) / (nwt
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)(fAbsc
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,
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
1359 G4double
1360 {
1361 const G4double tolerance = 1.0e-12;
1362 const G4double maxNumber = 12;
1363 G4int i, k, j;
1364 G4double alphaBeta, alphaReduced, betaReduc
1365
1366 G4double a, b, c, nwt1, nwt2, nwt3, nwt, te
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
1379 0.767999 * alp
1380 root2 = 1.0 + 1.48 * alphaReduce
1381 0.451998 * alphaReduced * alpha
1382 0.83001 * alphaReduced * betaRe
1383 root = 1.0 - root1 / root2;
1384 }
1385 else if(i == 2)
1386 {
1387 root1 = (4.1002 + alpha) / ((1.0 + alph
1388 root2 = 1.0 + 0.06 * (nJacobi - 8.0) *
1389 root3 =
1390 1.0 + 0.012002 * beta * (1.0 + 0.2499
1391 root -= (1.0 - root) * root1 * root2 *
1392 }
1393 else if(i == 3)
1394 {
1395 root1 = (1.67001 + 0.27998 * alpha) / (
1396 root2 = 1.0 + 0.22 * (nJacobi - 8.0) /
1397 root3 = 1.0 + 8.0 * beta / ((6.28001 +
1398 root -= (fAbscissa[0] - root) * root1 *
1399 }
1400 else if(i == nJacobi - 1)
1401 {
1402 root1 = (1.0 + 0.235002 * beta) / (0.76
1403 root2 = 1.0 / (1.0 + 0.639002 * (nJacob
1404 (1.0 + 0.71001 *
1405 root3 = 1.0 / (1.0 + 20.0 * alpha / ((7
1406 root += (root - fAbscissa[nJacobi - 4])
1407 }
1408 else if(i == nJacobi)
1409 {
1410 root1 = (1.0 + 0.37002 * beta) / (1.670
1411 root2 = 1.0 / (1.0 + 0.22 * (nJacobi -
1412 root3 =
1413 1.0 / (1.0 + 8.0 * alpha / ((6.28002
1414 root += (root - fAbscissa[nJacobi - 3])
1415 }
1416 else
1417 {
1418 root = 3.0 * fAbscissa[i - 2] - 3.0 * f
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
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) * (tem
1432 b = (temp - 1.0) *
1433 (alpha * alpha - beta * beta + te
1434 c = 2.0 * (j - 1 + alpha) * (j - 1
1435 nwt1 = (b * nwt2 - c * nwt3) / a;
1436 }
1437 nwt = (nJacobi * (alpha - beta - temp *
1438 2.0 * (nJacobi + alpha) * (nJaco
1439 (temp * (1.0 - root * root));
1440 rootTemp = root;
1441 root = rootTemp - nwt1 / nwt;
1442 if(std::fabs(root - rootTemp) <= tolera
1443 {
1444 break;
1445 }
1446 }
1447 if(k > maxNumber)
1448 {
1449 G4Exception("G4Integrator<T,F>::Jacobi(
1450 "Too many (>12) iterations.
1451 }
1452 fAbscissa[i - 1] = root;
1453 fWeight[i - 1] =
1454 std::exp(GammaLogarithm((G4double)(alph
1455 GammaLogarithm((G4double)(beta
1456 GammaLogarithm((G4double)(nJac
1457 GammaLogarithm((G4double)(nJac
1458 temp * std::pow(2.0, alphaBeta) / (nwt
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