| 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 // G4GaussLegendreQ class implementation
27 //
28 // Author: V.Grichine, 13.05.1997
29 // -------------------------------------------
30
31 #include "G4GaussLegendreQ.hh"
32 #include "G4PhysicalConstants.hh"
33
34 G4GaussLegendreQ::G4GaussLegendreQ(function pF
35 : G4VGaussianQuadrature(pFunction)
36 {}
37
38 // -------------------------------------------
39 //
40 // Constructor for GaussLegendre quadrature me
41 // the accuracy required, i.e the number of po
42 // will be evaluated during integration. The c
43 // abscissas and weights that are used in Gaus
44 // The values a and b are the limits of integr
45 // nLegendre MUST BE EVEN !!!
46
47 G4GaussLegendreQ::G4GaussLegendreQ(function pF
48 : G4VGaussianQuadrature(pFunction)
49 {
50 const G4double tolerance = 1.6e-10;
51 G4int k = nLegendre;
52 fNumber = (nLegendre + 1) /
53 if(2 * fNumber != k)
54 {
55 G4Exception("G4GaussLegendreQ::G4GaussLege
56 FatalException, "Invalid nLege
57 }
58 G4double newton0 = 0.0, newton1 = 0.0, temp1
59 temp = 0.0;
60
61 fAbscissa = new G4double[fNumber];
62 fWeight = new G4double[fNumber];
63
64 for(G4int i = 1; i <= fNumber; ++i) // Loop
65 {
66 newton0 = std::cos(pi * (i - 0.25) / (k +
67 do
68 {
69 temp1 = 1.0;
70 temp2 = 0.0;
71 for(G4int j = 1; j <= k; ++j)
72 {
73 temp3 = temp2;
74 temp2 = temp1;
75 temp1 = ((2.0 * j - 1.0) * newton0 * t
76 }
77 temp = k * (newton0 * temp1 - temp2)
78 newton1 = newton0;
79 newton0 = newton1 - temp1 / temp; // Ne
80 } while(std::fabs(newton0 - newton1) > tol
81
82 fAbscissa[fNumber - i] = newton0;
83 fWeight[fNumber - i] = 2.0 / ((1.0 - new
84 }
85 }
86
87 // -------------------------------------------
88 //
89 // Returns the integral of the function to be
90 // and b, by 2*fNumber point Gauss-Legendre in
91 // evaluated exactly 2*fNumber times at interi
92 // integration. Since the weights and abscissa
93 // around the midpoint of the range of integra
94 // fNumber distinct values of each.
95
96 G4double G4GaussLegendreQ::Integral(G4double a
97 {
98 G4double xMean = 0.5 * (a + b), xDiff = 0.5
99 dx = 0.0;
100 for(G4int i = 0; i < fNumber; ++i)
101 {
102 dx = xDiff * fAbscissa[i];
103 integral += fWeight[i] * (fFunction(xMean
104 }
105 return integral *= xDiff;
106 }
107
108 // -------------------------------------------
109 //
110 // Returns the integral of the function to be
111 // and b, by ten point Gauss-Legendre integrat
112 // exactly ten times at interior points in the
113 // weights and abscissas are, in this case, sy
114 // the range of integration, there are actuall
115 // each.
116
117 G4double G4GaussLegendreQ::QuickIntegral(G4dou
118 {
119 // From Abramowitz M., Stegan I.A. 1964 , Ha
120
121 static const G4double abscissa[] = { 0.14887
122 0.67940
123 0.97390
124
125 static const G4double weight[] = { 0.2955242
126 0.2190863
127 0.0666713
128 G4double xMean = 0.5 * (a + b), xDiff = 0.5
129 dx = 0.0;
130 for(G4int i = 0; i < 5; ++i)
131 {
132 dx = xDiff * abscissa[i];
133 integral += weight[i] * (fFunction(xMean +
134 }
135 return integral *= xDiff;
136 }
137
138 // -------------------------------------------
139 //
140 // Returns the integral of the function to be
141 // and b, by 96 point Gauss-Legendre integrati
142 // exactly ten times at interior points in the
143 // weights and abscissas are, in this case, sy
144 // the range of integration, there are actuall
145 // each.
146
147 G4double G4GaussLegendreQ::AccurateIntegral(G4
148 {
149 // From Abramowitz M., Stegan I.A. 1964 , Ha
150
151 static const G4double abscissa[] = {
152 0.016276744849602969579, 0.048812985136049
153 0.081297495464425558994, 0.113695850110665
154 0.145973714654896941989, 0.178096882367618
155
156 0.210031310460567203603, 0.241743156163840
157 0.273198812591049141487, 0.304364944354496
158 0.335208522892625422616, 0.365696861472313
159
160 0.395797649828908603285, 0.425478988407300
161 0.454709422167743008636, 0.483457973920596
162 0.511694177154667673586, 0.539388108324357
163
164 0.566510418561397168404, 0.593032364777572
165 0.618925840125468570386, 0.644163403784967
166 0.668718310043916153953, 0.692564536642171
167
168 0.715676812348967626225, 0.738030643744400
169 0.759602341176647498703, 0.780369043867433
170 0.800308744139140817229, 0.819400310737931
171
172 0.837623511228187121494, 0.854959033434601
173 0.871388505909296502874, 0.886894517402420
174 0.901460635315852341319, 0.915071423120898
175
176 0.927712456722308690965, 0.939370339752755
177 0.950032717784437635756, 0.959688291448742
178 0.968326828463264212174, 0.975939174585136
179
180 0.982517263563014677447, 0.988054126329623
181 0.992543900323762624572, 0.995981842987209
182 0.998364375863181677724, 0.999689503883230
183 };
184
185 static const G4double weight[] = {
186 0.032550614492363166242, 0.032516118713868
187 0.032447163714064269364, 0.032343822568575
188 0.032206204794030250669, 0.032034456231992
189
190 0.031828758894411006535, 0.031589330770727
191 0.031316425596862355813, 0.031010332586313
192 0.030671376123669149014, 0.030299915420827
193
194 0.029896344136328385984, 0.029461089958167
195 0.028994614150555236543, 0.028497411065085
196 0.027970007616848334440, 0.027412962726029
197
198 0.026826866725591762198, 0.026212340735672
199 0.025570036005349361499, 0.024900633222483
200 0.024204841792364691282, 0.023483399085926
201
202 0.022737069658329374001, 0.021966644438744
203 0.021172939892191298988, 0.020356797154333
204 0.019519081140145022410, 0.018660679627411
205
206 0.017782502316045260838, 0.016885479864245
207 0.015970562902562291381, 0.015038721026994
208 0.014090941772314860916, 0.013128229566961
209
210 0.012151604671088319635, 0.011162102099838
211 0.010160770535008415758, 0.009148671230783
212 0.008126876925698759217, 0.007096470791153
213
214 0.006058545504235961683, 0.005014202742927
215 0.003964554338444686674, 0.002910731817934
216 0.001853960788946921732, 0.000796792065552
217 };
218 G4double xMean = 0.5 * (a + b), xDiff = 0.5
219 dx = 0.0;
220 for(G4int i = 0; i < 48; ++i)
221 {
222 dx = xDiff * abscissa[i];
223 integral += weight[i] * (fFunction(xMean +
224 }
225 return integral *= xDiff;
226 }
227