Geant4  10.02.p01
G4GaussLaguerreQ.cc
Go to the documentation of this file.
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 // $Id: G4GaussLaguerreQ.cc 67970 2013-03-13 10:10:06Z gcosmo $
28 //
29 #include "G4GaussLaguerreQ.hh"
30 
31 
32 
33 // ------------------------------------------------------------
34 //
35 // Constructor for Gauss-Laguerre quadrature method: integral from zero to
36 // infinity of std::pow(x,alpha)*std::exp(-x)*f(x).
37 // The value of nLaguerre sets the accuracy.
38 // The constructor creates arrays fAbscissa[0,..,nLaguerre-1] and
39 // fWeight[0,..,nLaguerre-1] .
40 //
41 
44  G4int nLaguerre )
45  : G4VGaussianQuadrature(pFunction)
46 {
47  const G4double tolerance = 1.0e-10 ;
48  const G4int maxNumber = 12 ;
49  G4int i=1, k=1 ;
50  G4double newton0=0.0, newton1=0.0,
51  temp1=0.0, temp2=0.0, temp3=0.0, temp=0.0, cofi=0.0 ;
52 
53  fNumber = nLaguerre ;
54  fAbscissa = new G4double[fNumber] ;
55  fWeight = new G4double[fNumber] ;
56 
57  for(i=1;i<=fNumber;i++) // Loop over the desired roots
58  {
59  if(i == 1)
60  {
61  newton0 = (1.0 + alpha)*(3.0 + 0.92*alpha)
62  / (1.0 + 2.4*fNumber + 1.8*alpha) ;
63  }
64  else if(i == 2)
65  {
66  newton0 += (15.0 + 6.25*alpha)/(1.0 + 0.9*alpha + 2.5*fNumber) ;
67  }
68  else
69  {
70  cofi = i - 2 ;
71  newton0 += ((1.0+2.55*cofi)/(1.9*cofi)
72  + 1.26*cofi*alpha/(1.0+3.5*cofi))
73  * (newton0 - fAbscissa[i-3])/(1.0 + 0.3*alpha) ;
74  }
75  for(k=1;k<=maxNumber;k++)
76  {
77  temp1 = 1.0 ;
78  temp2 = 0.0 ;
79  for(G4int j=1;j<=fNumber;j++)
80  {
81  temp3 = temp2 ;
82  temp2 = temp1 ;
83  temp1 = ((2*j - 1 + alpha - newton0)*temp2
84  - (j - 1 + alpha)*temp3)/j ;
85  }
86  temp = (fNumber*temp1 - (fNumber +alpha)*temp2)/newton0 ;
87  newton1 = newton0 ;
88  newton0 = newton1 - temp1/temp ;
89  if(std::fabs(newton0 - newton1) <= tolerance)
90  {
91  break ;
92  }
93  }
94  if(k > maxNumber)
95  {
96  G4Exception("G4GaussLaguerreQ::G4GaussLaguerreQ()",
97  "OutOfRange", FatalException,
98  "Too many iterations in Gauss-Laguerre constructor") ;
99  }
100 
101  fAbscissa[i-1] = newton0 ;
102  fWeight[i-1] = -std::exp(GammaLogarithm(alpha + fNumber)
103  - GammaLogarithm((G4double)fNumber))/(temp*fNumber*temp2) ;
104  }
105 }
106 
107 // -----------------------------------------------------------------
108 //
109 // Gauss-Laguerre method for integration of
110 // std::pow(x,alpha)*std::exp(-x)*pFunction(x)
111 // from zero up to infinity. pFunction is evaluated in fNumber points
112 // for which fAbscissa[i] and fWeight[i] arrays were created in
113 // G4VGaussianQuadrature(double,int) constructor
114 
115 G4double
117 {
118  G4double integral = 0.0 ;
119  for(G4int i=0;i<fNumber;i++)
120  {
121  integral += fWeight[i]*fFunction(fAbscissa[i]) ;
122  }
123  return integral ;
124 }
G4GaussLaguerreQ(function pFunction, G4double alpha, G4int nLaguerre)
static const G4float tolerance
int G4int
Definition: G4Types.hh:78
G4double GammaLogarithm(G4double xx)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double G4double
Definition: G4Types.hh:76
static const G4double alpha
G4double Integral() const