| 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 // G4GaussLaguerreQ class implementation
27 //
28 // Author: V.Grichine, 13.05.1997
29 // -------------------------------------------
30
31 #include "G4GaussLaguerreQ.hh"
32
33 // -------------------------------------------
34 //
35 // Constructor for Gauss-Laguerre quadrature m
36 // infinity of std::pow(x,alpha)*std::exp(-x)*
37 // The value of nLaguerre sets the accuracy.
38 // The constructor creates arrays fAbscissa[0,
39 // fWeight[0,..,nLaguerre-1] .
40 //
41
42 G4GaussLaguerreQ::G4GaussLaguerreQ(function pF
43 G4int nLagu
44 : G4VGaussianQuadrature(pFunction)
45 {
46 const G4double tolerance = 1.0e-10;
47 const G4int maxNumber = 12;
48 G4int i = 1, k = 1;
49 G4double newton0 = 0.0, newton1 = 0.0, temp1
50 temp = 0.0, cofi = 0.0;
51
52 fNumber = nLaguerre;
53 fAbscissa = new G4double[fNumber];
54 fWeight = new G4double[fNumber];
55
56 for(i = 1; i <= fNumber; ++i) // Loop over
57 {
58 if(i == 1)
59 {
60 newton0 = (1.0 + alpha) * (3.0 + 0.92 *
61 (1.0 + 2.4 * fNumber + 1.8 * a
62 }
63 else if(i == 2)
64 {
65 newton0 += (15.0 + 6.25 * alpha) / (1.0
66 }
67 else
68 {
69 cofi = i - 2;
70 newton0 += ((1.0 + 2.55 * cofi) / (1.9 *
71 1.26 * cofi * alpha / (1.0 +
72 (newton0 - fAbscissa[i - 3])
73 }
74 for(k = 1; k <= maxNumber; ++k)
75 {
76 temp1 = 1.0;
77 temp2 = 0.0;
78 for(G4int j = 1; j <= fNumber; ++j)
79 {
80 temp3 = temp2;
81 temp2 = temp1;
82 temp1 =
83 ((2 * j - 1 + alpha - newton0) * tem
84 }
85 temp = (fNumber * temp1 - (fNumber +
86 newton1 = newton0;
87 newton0 = newton1 - temp1 / temp;
88 if(std::fabs(newton0 - newton1) <= toler
89 {
90 break;
91 }
92 }
93 if(k > maxNumber)
94 {
95 G4Exception("G4GaussLaguerreQ::G4GaussLa
96 FatalException,
97 "Too many iterations in Gaus
98 }
99
100 fAbscissa[i - 1] = newton0;
101 fWeight[i - 1] = -std::exp(GammaLogarith
102 GammaLogarithm(
103 (temp * fNumber * temp2);
104 }
105 }
106
107 // -------------------------------------------
108 //
109 // Gauss-Laguerre method for integration of
110 // std::pow(x,alpha)*std::exp(-x)*pFunction(x)
111 // from zero up to infinity. pFunction is eval
112 // for which fAbscissa[i] and fWeight[i] array
113 // G4VGaussianQuadrature(double,int) construct
114
115 G4double G4GaussLaguerreQ::Integral() const
116 {
117 G4double integral = 0.0;
118 for(G4int i = 0; i < fNumber; ++i)
119 {
120 integral += fWeight[i] * fFunction(fAbscis
121 }
122 return integral;
123 }
124