| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4GaussJacobiQ class implementation 26 // G4GaussJacobiQ class implementation
27 // 27 //
28 // Author: V.Grichine, 13.05.1997 28 // Author: V.Grichine, 13.05.1997
29 // ------------------------------------------- 29 // --------------------------------------------------------------------
30 30
31 #include "G4GaussJacobiQ.hh" 31 #include "G4GaussJacobiQ.hh"
32 32
33 // ------------------------------------------- 33 // -------------------------------------------------------------
34 // 34 //
35 // Constructor for Gauss-Jacobi integration me 35 // Constructor for Gauss-Jacobi integration method.
36 // 36 //
37 37
38 G4GaussJacobiQ::G4GaussJacobiQ(function pFunct 38 G4GaussJacobiQ::G4GaussJacobiQ(function pFunction, G4double alpha,
39 G4double beta, 39 G4double beta, G4int nJacobi)
40 : G4VGaussianQuadrature(pFunction) 40 : G4VGaussianQuadrature(pFunction)
41 41
42 { 42 {
43 const G4double tolerance = 1.0e-12; 43 const G4double tolerance = 1.0e-12;
44 const G4double maxNumber = 12; 44 const G4double maxNumber = 12;
45 G4int i = 1, k = 1; 45 G4int i = 1, k = 1;
46 G4double root = 0.; 46 G4double root = 0.;
47 G4double alphaBeta = 0.0, alphaReduced = 0.0 47 G4double alphaBeta = 0.0, alphaReduced = 0.0, betaReduced = 0.0, root1 = 0.0,
48 root2 = 0.0, root3 = 0.0; 48 root2 = 0.0, root3 = 0.0;
49 G4double a = 0.0, b = 0.0, c = 0.0, newton1 49 G4double a = 0.0, b = 0.0, c = 0.0, newton1 = 0.0, newton2 = 0.0,
50 newton3 = 0.0, newton0 = 0.0, temp 50 newton3 = 0.0, newton0 = 0.0, temp = 0.0, rootTemp = 0.0;
51 51
52 fNumber = nJacobi; 52 fNumber = nJacobi;
53 fAbscissa = new G4double[fNumber]; 53 fAbscissa = new G4double[fNumber];
54 fWeight = new G4double[fNumber]; 54 fWeight = new G4double[fNumber];
55 55
56 for(i = 1; i <= nJacobi; ++i) 56 for(i = 1; i <= nJacobi; ++i)
57 { 57 {
58 if(i == 1) 58 if(i == 1)
59 { 59 {
60 alphaReduced = alpha / nJacobi; 60 alphaReduced = alpha / nJacobi;
61 betaReduced = beta / nJacobi; 61 betaReduced = beta / nJacobi;
62 root1 = (1.0 + alpha) * (2.78002 62 root1 = (1.0 + alpha) * (2.78002 / (4.0 + nJacobi * nJacobi) +
63 0.767999 * alph 63 0.767999 * alphaReduced / nJacobi);
64 root2 = 1.0 + 1.48 * alphaReduced 64 root2 = 1.0 + 1.48 * alphaReduced + 0.96002 * betaReduced +
65 0.451998 * alphaReduced * alphaR 65 0.451998 * alphaReduced * alphaReduced +
66 0.83001 * alphaReduced * betaRed 66 0.83001 * alphaReduced * betaReduced;
67 root = 1.0 - root1 / root2; 67 root = 1.0 - root1 / root2;
68 } 68 }
69 else if(i == 2) 69 else if(i == 2)
70 { 70 {
71 root1 = (4.1002 + alpha) / ((1.0 + alpha 71 root1 = (4.1002 + alpha) / ((1.0 + alpha) * (1.0 + 0.155998 * alpha));
72 root2 = 1.0 + 0.06 * (nJacobi - 8.0) * ( 72 root2 = 1.0 + 0.06 * (nJacobi - 8.0) * (1.0 + 0.12 * alpha) / nJacobi;
73 root3 = 73 root3 =
74 1.0 + 0.012002 * beta * (1.0 + 0.24997 74 1.0 + 0.012002 * beta * (1.0 + 0.24997 * std::fabs(alpha)) / nJacobi;
75 root -= (1.0 - root) * root1 * root2 * r 75 root -= (1.0 - root) * root1 * root2 * root3;
76 } 76 }
77 else if(i == 3) 77 else if(i == 3)
78 { 78 {
79 root1 = (1.67001 + 0.27998 * alpha) / (1 79 root1 = (1.67001 + 0.27998 * alpha) / (1.0 + 0.37002 * alpha);
80 root2 = 1.0 + 0.22 * (nJacobi - 8.0) / n 80 root2 = 1.0 + 0.22 * (nJacobi - 8.0) / nJacobi;
81 root3 = 1.0 + 8.0 * beta / ((6.28001 + b 81 root3 = 1.0 + 8.0 * beta / ((6.28001 + beta) * nJacobi * nJacobi);
82 root -= (fAbscissa[0] - root) * root1 * 82 root -= (fAbscissa[0] - root) * root1 * root2 * root3;
83 } 83 }
84 else if(i == nJacobi - 1) 84 else if(i == nJacobi - 1)
85 { 85 {
86 root1 = (1.0 + 0.235002 * beta) / (0.766 86 root1 = (1.0 + 0.235002 * beta) / (0.766001 + 0.118998 * beta);
87 root2 = 1.0 / (1.0 + 0.639002 * (nJacobi 87 root2 = 1.0 / (1.0 + 0.639002 * (nJacobi - 4.0) /
88 (1.0 + 0.71001 * 88 (1.0 + 0.71001 * (nJacobi - 4.0)));
89 root3 = 1.0 / (1.0 + 20.0 * alpha / ((7. 89 root3 = 1.0 / (1.0 + 20.0 * alpha / ((7.5 + alpha) * nJacobi * nJacobi));
90 root += (root - fAbscissa[nJacobi - 4]) 90 root += (root - fAbscissa[nJacobi - 4]) * root1 * root2 * root3;
91 } 91 }
92 else if(i == nJacobi) 92 else if(i == nJacobi)
93 { 93 {
94 root1 = (1.0 + 0.37002 * beta) / (1.6700 94 root1 = (1.0 + 0.37002 * beta) / (1.67001 + 0.27998 * beta);
95 root2 = 1.0 / (1.0 + 0.22 * (nJacobi - 8 95 root2 = 1.0 / (1.0 + 0.22 * (nJacobi - 8.0) / nJacobi);
96 root3 = 96 root3 =
97 1.0 / (1.0 + 8.0 * alpha / ((6.28002 + 97 1.0 / (1.0 + 8.0 * alpha / ((6.28002 + alpha) * nJacobi * nJacobi));
98 root += (root - fAbscissa[nJacobi - 3]) 98 root += (root - fAbscissa[nJacobi - 3]) * root1 * root2 * root3;
99 } 99 }
100 else 100 else
101 { 101 {
102 root = 3.0 * fAbscissa[i - 2] - 3.0 * fA 102 root = 3.0 * fAbscissa[i - 2] - 3.0 * fAbscissa[i - 3] + fAbscissa[i - 4];
103 } 103 }
104 alphaBeta = alpha + beta; 104 alphaBeta = alpha + beta;
105 for(k = 1; k <= maxNumber; ++k) 105 for(k = 1; k <= maxNumber; ++k)
106 { 106 {
107 temp = 2.0 + alphaBeta; 107 temp = 2.0 + alphaBeta;
108 newton1 = (alpha - beta + temp * root) / 108 newton1 = (alpha - beta + temp * root) / 2.0;
109 newton2 = 1.0; 109 newton2 = 1.0;
110 for(G4int j = 2; j <= nJacobi; ++j) 110 for(G4int j = 2; j <= nJacobi; ++j)
111 { 111 {
112 newton3 = newton2; 112 newton3 = newton2;
113 newton2 = newton1; 113 newton2 = newton1;
114 temp = 2 * j + alphaBeta; 114 temp = 2 * j + alphaBeta;
115 a = 2 * j * (j + alphaBeta) * (t 115 a = 2 * j * (j + alphaBeta) * (temp - 2.0);
116 b = (temp - 1.0) * 116 b = (temp - 1.0) *
117 (alpha * alpha - beta * beta + tem 117 (alpha * alpha - beta * beta + temp * (temp - 2.0) * root);
118 c = 2.0 * (j - 1 + alpha) * (j - 118 c = 2.0 * (j - 1 + alpha) * (j - 1 + beta) * temp;
119 newton1 = (b * newton2 - c * newton3) 119 newton1 = (b * newton2 - c * newton3) / a;
120 } 120 }
121 newton0 = (nJacobi * (alpha - beta - tem 121 newton0 = (nJacobi * (alpha - beta - temp * root) * newton1 +
122 2.0 * (nJacobi + alpha) * (nJ 122 2.0 * (nJacobi + alpha) * (nJacobi + beta) * newton2) /
123 (temp * (1.0 - root * root)); 123 (temp * (1.0 - root * root));
124 rootTemp = root; 124 rootTemp = root;
125 root = rootTemp - newton1 / newton0; 125 root = rootTemp - newton1 / newton0;
126 if(std::fabs(root - rootTemp) <= toleran 126 if(std::fabs(root - rootTemp) <= tolerance)
127 { 127 {
128 break; 128 break;
129 } 129 }
130 } 130 }
131 if(k > maxNumber) 131 if(k > maxNumber)
132 { 132 {
133 G4Exception("G4GaussJacobiQ::G4GaussJaco 133 G4Exception("G4GaussJacobiQ::G4GaussJacobiQ()", "OutOfRange",
134 FatalException, "Too many it 134 FatalException, "Too many iterations in constructor.");
135 } 135 }
136 fAbscissa[i - 1] = root; 136 fAbscissa[i - 1] = root;
137 fWeight[i - 1] = 137 fWeight[i - 1] =
138 std::exp(GammaLogarithm((G4double)(alpha 138 std::exp(GammaLogarithm((G4double)(alpha + nJacobi)) +
139 GammaLogarithm((G4double)(beta 139 GammaLogarithm((G4double)(beta + nJacobi)) -
140 GammaLogarithm((G4double)(nJaco 140 GammaLogarithm((G4double)(nJacobi + 1.0)) -
141 GammaLogarithm((G4double)(nJaco 141 GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0))) *
142 temp * std::pow(2.0, alphaBeta) / (newto 142 temp * std::pow(2.0, alphaBeta) / (newton0 * newton2);
143 } 143 }
144 } 144 }
145 145
146 // ------------------------------------------- 146 // ----------------------------------------------------------
147 // 147 //
148 // Gauss-Jacobi method for integration of 148 // Gauss-Jacobi method for integration of
149 // ((1-x)^alpha)*((1+x)^beta)*pFunction(x) 149 // ((1-x)^alpha)*((1+x)^beta)*pFunction(x)
150 // from minus unit to plus unit . 150 // from minus unit to plus unit .
151 151
152 G4double G4GaussJacobiQ::Integral() const 152 G4double G4GaussJacobiQ::Integral() const
153 { 153 {
154 G4double integral = 0.0; 154 G4double integral = 0.0;
155 for(G4int i = 0; i < fNumber; ++i) 155 for(G4int i = 0; i < fNumber; ++i)
156 { 156 {
157 integral += fWeight[i] * fFunction(fAbscis 157 integral += fWeight[i] * fFunction(fAbscissa[i]);
158 } 158 }
159 return integral; 159 return integral;
160 } 160 }
161 161