ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4GaussLegendreQ.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4GaussLegendreQ.cc
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 //
27 //
28 
29 #include "G4GaussLegendreQ.hh"
30 #include "G4PhysicalConstants.hh"
31 
32 G4GaussLegendreQ::G4GaussLegendreQ( function pFunction )
33  : G4VGaussianQuadrature(pFunction)
34 {
35 }
36 
37 // --------------------------------------------------------------------------
38 //
39 // Constructor for GaussLegendre quadrature method. The value nLegendre sets
40 // the accuracy required, i.e the number of points where the function pFunction
41 // will be evaluated during integration. The constructor creates the arrays for
42 // abscissas and weights that are used in Gauss-Legendre quadrature method.
43 // The values a and b are the limits of integration of the pFunction.
44 // nLegendre MUST BE EVEN !!!
45 
46 G4GaussLegendreQ::G4GaussLegendreQ( function pFunction,
47  G4int nLegendre )
48  : G4VGaussianQuadrature(pFunction)
49 {
50  const G4double tolerance = 1.6e-10 ;
51  G4int k = nLegendre ;
52  fNumber = (nLegendre + 1)/2 ;
53  if(2*fNumber != k)
54  {
55  G4Exception("G4GaussLegendreQ::G4GaussLegendreQ()", "InvalidCall",
56  FatalException, "Invalid nLegendre argument !") ;
57  }
58  G4double newton0=0.0, newton1=0.0,
59  temp1=0.0, temp2=0.0, temp3=0.0, temp=0.0 ;
60 
61  fAbscissa = new G4double[fNumber] ;
62  fWeight = new G4double[fNumber] ;
63 
64  for(G4int i=1;i<=fNumber;i++) // Loop over the desired roots
65  {
66  newton0 = std::cos(pi*(i - 0.25)/(k + 0.5)) ; // Initial root
67  do // approximation
68  { // loop of Newton's method
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*temp2 - (j - 1.0)*temp3)/j ;
76  }
77  temp = k*(newton0*temp1 - temp2)/(newton0*newton0 - 1.0) ;
78  newton1 = newton0 ;
79  newton0 = newton1 - temp1/temp ; // Newton's method
80  }
81  while(std::fabs(newton0 - newton1) > tolerance) ;
82 
83  fAbscissa[fNumber-i] = newton0 ;
84  fWeight[fNumber-i] = 2.0/((1.0 - newton0*newton0)*temp*temp) ;
85  }
86 }
87 
88 // --------------------------------------------------------------------------
89 //
90 // Returns the integral of the function to be pointed by fFunction between a
91 // and b, by 2*fNumber point Gauss-Legendre integration: the function is
92 // evaluated exactly 2*fNumber times at interior points in the range of
93 // integration. Since the weights and abscissas are, in this case, symmetric
94 // around the midpoint of the range of integration, there are actually only
95 // fNumber distinct values of each.
96 
97 G4double
98 G4GaussLegendreQ::Integral(G4double a, G4double b) const
99 {
100  G4double xMean = 0.5*(a + b),
101  xDiff = 0.5*(b - a),
102  integral = 0.0, dx = 0.0 ;
103  for(G4int i=0;i<fNumber;i++)
104  {
105  dx = xDiff*fAbscissa[i] ;
106  integral += fWeight[i]*(fFunction(xMean + dx) + fFunction(xMean - dx)) ;
107  }
108  return integral *= xDiff ;
109 }
110 
111 // --------------------------------------------------------------------------
112 //
113 // Returns the integral of the function to be pointed by fFunction between a
114 // and b, by ten point Gauss-Legendre integration: the function is evaluated
115 // exactly ten times at interior points in the range of integration. Since the
116 // weights and abscissas are, in this case, symmetric around the midpoint of
117 // the range of integration, there are actually only five distinct values of
118 // each.
119 
120 G4double
121  G4GaussLegendreQ::QuickIntegral(G4double a, G4double b) const
122 {
123  // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 916
124 
125  static const G4double abscissa[] = { 0.148874338981631, 0.433395394129247,
126  0.679409568299024, 0.865063366688985,
127  0.973906528517172 } ;
128 
129  static const G4double weight[] = { 0.295524224714753, 0.269266719309996,
130  0.219086362515982, 0.149451349150581,
131  0.066671344308688 } ;
132  G4double xMean = 0.5*(a + b),
133  xDiff = 0.5*(b - a),
134  integral = 0.0, dx = 0.0 ;
135  for(G4int i=0;i<5;i++)
136  {
137  dx = xDiff*abscissa[i] ;
138  integral += weight[i]*(fFunction(xMean + dx) + fFunction(xMean - dx)) ;
139  }
140  return integral *= xDiff ;
141 }
142 
143 // -------------------------------------------------------------------------
144 //
145 // Returns the integral of the function to be pointed by fFunction between a
146 // and b, by 96 point Gauss-Legendre integration: the function is evaluated
147 // exactly ten times at interior points in the range of integration. Since the
148 // weights and abscissas are, in this case, symmetric around the midpoint of
149 // the range of integration, there are actually only five distinct values of
150 // each.
151 
152 G4double
153  G4GaussLegendreQ::AccurateIntegral(G4double a, G4double b) const
154 {
155  // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919
156 
157  static const
158  G4double abscissa[] = {
159  0.016276744849602969579, 0.048812985136049731112,
160  0.081297495464425558994, 0.113695850110665920911,
161  0.145973714654896941989, 0.178096882367618602759, // 6
162 
163  0.210031310460567203603, 0.241743156163840012328,
164  0.273198812591049141487, 0.304364944354496353024,
165  0.335208522892625422616, 0.365696861472313635031, // 12
166 
167  0.395797649828908603285, 0.425478988407300545365,
168  0.454709422167743008636, 0.483457973920596359768,
169  0.511694177154667673586, 0.539388108324357436227, // 18
170 
171  0.566510418561397168404, 0.593032364777572080684,
172  0.618925840125468570386, 0.644163403784967106798,
173  0.668718310043916153953, 0.692564536642171561344, // 24
174 
175  0.715676812348967626225, 0.738030643744400132851,
176  0.759602341176647498703, 0.780369043867433217604,
177  0.800308744139140817229, 0.819400310737931675539, // 30
178 
179  0.837623511228187121494, 0.854959033434601455463,
180  0.871388505909296502874, 0.886894517402420416057,
181  0.901460635315852341319, 0.915071423120898074206, // 36
182 
183  0.927712456722308690965, 0.939370339752755216932,
184  0.950032717784437635756, 0.959688291448742539300,
185  0.968326828463264212174, 0.975939174585136466453, // 42
186 
187  0.982517263563014677447, 0.988054126329623799481,
188  0.992543900323762624572, 0.995981842987209290650,
189  0.998364375863181677724, 0.999689503883230766828 // 48
190  } ;
191 
192  static const
193  G4double weight[] = {
194  0.032550614492363166242, 0.032516118713868835987,
195  0.032447163714064269364, 0.032343822568575928429,
196  0.032206204794030250669, 0.032034456231992663218, // 6
197 
198  0.031828758894411006535, 0.031589330770727168558,
199  0.031316425596862355813, 0.031010332586313837423,
200  0.030671376123669149014, 0.030299915420827593794, // 12
201 
202  0.029896344136328385984, 0.029461089958167905970,
203  0.028994614150555236543, 0.028497411065085385646,
204  0.027970007616848334440, 0.027412962726029242823, // 18
205 
206  0.026826866725591762198, 0.026212340735672413913,
207  0.025570036005349361499, 0.024900633222483610288,
208  0.024204841792364691282, 0.023483399085926219842, // 24
209 
210  0.022737069658329374001, 0.021966644438744349195,
211  0.021172939892191298988, 0.020356797154333324595,
212  0.019519081140145022410, 0.018660679627411467385, // 30
213 
214  0.017782502316045260838, 0.016885479864245172450,
215  0.015970562902562291381, 0.015038721026994938006,
216  0.014090941772314860916, 0.013128229566961572637, // 36
217 
218  0.012151604671088319635, 0.011162102099838498591,
219  0.010160770535008415758, 0.009148671230783386633,
220  0.008126876925698759217, 0.007096470791153865269, // 42
221 
222  0.006058545504235961683, 0.005014202742927517693,
223  0.003964554338444686674, 0.002910731817934946408,
224  0.001853960788946921732, 0.000796792065552012429 // 48
225  } ;
226  G4double xMean = 0.5*(a + b),
227  xDiff = 0.5*(b - a),
228  integral = 0.0, dx = 0.0 ;
229  for(G4int i=0;i<48;i++)
230  {
231  dx = xDiff*abscissa[i] ;
232  integral += weight[i]*(fFunction(xMean + dx) + fFunction(xMean - dx)) ;
233  }
234  return integral *= xDiff ;
235 }