Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/global/HEPNumerics/src/G4GaussLegendreQ.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /global/HEPNumerics/src/G4GaussLegendreQ.cc (Version 11.3.0) and /global/HEPNumerics/src/G4GaussLegendreQ.cc (Version 11.2.1)


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