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 9.5.p1)


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