| 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 // 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 method: integral from zero to
36 // infinity of std::pow(x,alpha)*std::exp(-x)*f(x).
37 // The value of nLaguerre sets the accuracy.
38 // The constructor creates arrays fAbscissa[0,..,nLaguerre-1] and
39 // fWeight[0,..,nLaguerre-1] .
40 //
41
42 G4GaussLaguerreQ::G4GaussLaguerreQ(function pFunction, G4double alpha,
43 G4int nLaguerre)
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 = 0.0, temp2 = 0.0, temp3 = 0.0,
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 the desired roots
57 {
58 if(i == 1)
59 {
60 newton0 = (1.0 + alpha) * (3.0 + 0.92 * alpha) /
61 (1.0 + 2.4 * fNumber + 1.8 * alpha);
62 }
63 else if(i == 2)
64 {
65 newton0 += (15.0 + 6.25 * alpha) / (1.0 + 0.9 * alpha + 2.5 * fNumber);
66 }
67 else
68 {
69 cofi = i - 2;
70 newton0 += ((1.0 + 2.55 * cofi) / (1.9 * cofi) +
71 1.26 * cofi * alpha / (1.0 + 3.5 * cofi)) *
72 (newton0 - fAbscissa[i - 3]) / (1.0 + 0.3 * alpha);
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) * temp2 - (j - 1 + alpha) * temp3) / j;
84 }
85 temp = (fNumber * temp1 - (fNumber + alpha) * temp2) / newton0;
86 newton1 = newton0;
87 newton0 = newton1 - temp1 / temp;
88 if(std::fabs(newton0 - newton1) <= tolerance)
89 {
90 break;
91 }
92 }
93 if(k > maxNumber)
94 {
95 G4Exception("G4GaussLaguerreQ::G4GaussLaguerreQ()", "OutOfRange",
96 FatalException,
97 "Too many iterations in Gauss-Laguerre constructor");
98 }
99
100 fAbscissa[i - 1] = newton0;
101 fWeight[i - 1] = -std::exp(GammaLogarithm(alpha + fNumber) -
102 GammaLogarithm((G4double) fNumber)) /
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 evaluated in fNumber points
112 // for which fAbscissa[i] and fWeight[i] arrays were created in
113 // G4VGaussianQuadrature(double,int) constructor
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(fAbscissa[i]);
121 }
122 return integral;
123 }
124