| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // G4GaussJacobiQ class implementation <<
27 // 23 //
28 // Author: V.Grichine, 13.05.1997 << 24 // $Id: G4GaussJacobiQ.cc,v 1.6 2005/03/15 19:11:35 gcosmo Exp $
29 // ------------------------------------------- << 25 // GEANT4 tag $Name: geant4-08-00 $
30 << 26 //
31 #include "G4GaussJacobiQ.hh" 27 #include "G4GaussJacobiQ.hh"
32 28
>> 29
33 // ------------------------------------------- 30 // -------------------------------------------------------------
34 // 31 //
35 // Constructor for Gauss-Jacobi integration me << 32 // Constructor for Gauss-Jacobi integration method.
36 // 33 //
37 34
38 G4GaussJacobiQ::G4GaussJacobiQ(function pFunct << 35 G4GaussJacobiQ::G4GaussJacobiQ( function pFunction,
39 G4double beta, << 36 G4double alpha,
40 : G4VGaussianQuadrature(pFunction) << 37 G4double beta,
>> 38 G4int nJacobi )
>> 39 : G4VGaussianQuadrature(pFunction)
41 40
42 { 41 {
43 const G4double tolerance = 1.0e-12; << 42 const G4double tolerance = 1.0e-12 ;
44 const G4double maxNumber = 12; << 43 const G4double maxNumber = 12 ;
45 G4int i = 1, k = 1; << 44 G4int i=1, k=1 ;
46 G4double root = 0.; << 45 G4double root=0.;
47 G4double alphaBeta = 0.0, alphaReduced = 0.0 << 46 G4double alphaBeta=0.0, alphaReduced=0.0, betaReduced=0.0,
48 root2 = 0.0, root3 = 0.0; << 47 root1=0.0, root2=0.0, root3=0.0 ;
49 G4double a = 0.0, b = 0.0, c = 0.0, newton1 << 48 G4double a=0.0, b=0.0, c=0.0,
50 newton3 = 0.0, newton0 = 0.0, temp << 49 newton1=0.0, newton2=0.0, newton3=0.0, newton=0.0,
51 << 50 temp=0.0, rootTemp=0.0 ;
52 fNumber = nJacobi; << 51
53 fAbscissa = new G4double[fNumber]; << 52 fNumber = nJacobi ;
54 fWeight = new G4double[fNumber]; << 53 fAbscissa = new G4double[fNumber] ;
>> 54 fWeight = new G4double[fNumber] ;
55 55
56 for(i = 1; i <= nJacobi; ++i) << 56 for (i=1;i<=nJacobi;i++)
57 { 57 {
58 if(i == 1) << 58 if (i == 1)
59 { << 59 {
60 alphaReduced = alpha / nJacobi; << 60 alphaReduced = alpha/nJacobi ;
61 betaReduced = beta / nJacobi; << 61 betaReduced = beta/nJacobi ;
62 root1 = (1.0 + alpha) * (2.78002 << 62 root1 = (1.0+alpha)*(2.78002/(4.0+nJacobi*nJacobi)+
63 0.767999 * alph << 63 0.767999*alphaReduced/nJacobi) ;
64 root2 = 1.0 + 1.48 * alphaReduced << 64 root2 = 1.0+1.48*alphaReduced+0.96002*betaReduced
65 0.451998 * alphaReduced * alphaR << 65 + 0.451998*alphaReduced*alphaReduced
66 0.83001 * alphaReduced * betaRed << 66 + 0.83001*alphaReduced*betaReduced ;
67 root = 1.0 - root1 / root2; << 67 root = 1.0-root1/root2 ;
68 } << 68 }
69 else if(i == 2) << 69 else if (i == 2)
70 { << 70 {
71 root1 = (4.1002 + alpha) / ((1.0 + alpha << 71 root1=(4.1002+alpha)/((1.0+alpha)*(1.0+0.155998*alpha)) ;
72 root2 = 1.0 + 0.06 * (nJacobi - 8.0) * ( << 72 root2=1.0+0.06*(nJacobi-8.0)*(1.0+0.12*alpha)/nJacobi ;
73 root3 = << 73 root3=1.0+0.012002*beta*(1.0+0.24997*std::fabs(alpha))/nJacobi ;
74 1.0 + 0.012002 * beta * (1.0 + 0.24997 << 74 root -= (1.0-root)*root1*root2*root3 ;
75 root -= (1.0 - root) * root1 * root2 * r << 75 }
76 } << 76 else if (i == 3)
77 else if(i == 3) << 77 {
78 { << 78 root1=(1.67001+0.27998*alpha)/(1.0+0.37002*alpha) ;
79 root1 = (1.67001 + 0.27998 * alpha) / (1 << 79 root2=1.0+0.22*(nJacobi-8.0)/nJacobi ;
80 root2 = 1.0 + 0.22 * (nJacobi - 8.0) / n << 80 root3=1.0+8.0*beta/((6.28001+beta)*nJacobi*nJacobi) ;
81 root3 = 1.0 + 8.0 * beta / ((6.28001 + b << 81 root -= (fAbscissa[0]-root)*root1*root2*root3 ;
82 root -= (fAbscissa[0] - root) * root1 * << 82 }
83 } << 83 else if (i == nJacobi-1)
84 else if(i == nJacobi - 1) << 84 {
85 { << 85 root1=(1.0+0.235002*beta)/(0.766001+0.118998*beta) ;
86 root1 = (1.0 + 0.235002 * beta) / (0.766 << 86 root2=1.0/(1.0+0.639002*(nJacobi-4.0)/(1.0+0.71001*(nJacobi-4.0))) ;
87 root2 = 1.0 / (1.0 + 0.639002 * (nJacobi << 87 root3=1.0/(1.0+20.0*alpha/((7.5+alpha)*nJacobi*nJacobi)) ;
88 (1.0 + 0.71001 * << 88 root += (root-fAbscissa[nJacobi-4])*root1*root2*root3 ;
89 root3 = 1.0 / (1.0 + 20.0 * alpha / ((7. << 89 }
90 root += (root - fAbscissa[nJacobi - 4]) << 90 else if (i == nJacobi)
91 } << 91 {
92 else if(i == nJacobi) << 92 root1 = (1.0+0.37002*beta)/(1.67001+0.27998*beta) ;
93 { << 93 root2 = 1.0/(1.0+0.22*(nJacobi-8.0)/nJacobi) ;
94 root1 = (1.0 + 0.37002 * beta) / (1.6700 << 94 root3 = 1.0/(1.0+8.0*alpha/((6.28002+alpha)*nJacobi*nJacobi)) ;
95 root2 = 1.0 / (1.0 + 0.22 * (nJacobi - 8 << 95 root += (root-fAbscissa[nJacobi-3])*root1*root2*root3 ;
96 root3 = << 96 }
97 1.0 / (1.0 + 8.0 * alpha / ((6.28002 + << 97 else
98 root += (root - fAbscissa[nJacobi - 3]) << 98 {
99 } << 99 root = 3.0*fAbscissa[i-2]-3.0*fAbscissa[i-3]+fAbscissa[i-4] ;
100 else << 100 }
101 { << 101 alphaBeta = alpha + beta ;
102 root = 3.0 * fAbscissa[i - 2] - 3.0 * fA << 102 for (k=1;k<=maxNumber;k++)
103 } << 103 {
104 alphaBeta = alpha + beta; << 104 temp = 2.0 + alphaBeta ;
105 for(k = 1; k <= maxNumber; ++k) << 105 newton1 = (alpha-beta+temp*root)/2.0 ;
106 { << 106 newton2 = 1.0 ;
107 temp = 2.0 + alphaBeta; << 107 for (G4int j=2;j<=nJacobi;j++)
108 newton1 = (alpha - beta + temp * root) / << 108 {
109 newton2 = 1.0; << 109 newton3 = newton2 ;
110 for(G4int j = 2; j <= nJacobi; ++j) << 110 newton2 = newton1 ;
111 { << 111 temp = 2*j+alphaBeta ;
112 newton3 = newton2; << 112 a = 2*j*(j+alphaBeta)*(temp-2.0) ;
113 newton2 = newton1; << 113 b = (temp-1.0)*(alpha*alpha-beta*beta+temp*(temp-2.0)*root) ;
114 temp = 2 * j + alphaBeta; << 114 c = 2.0*(j-1+alpha)*(j-1+beta)*temp ;
115 a = 2 * j * (j + alphaBeta) * (t << 115 newton1 = (b*newton2-c*newton3)/a ;
116 b = (temp - 1.0) * << 116 }
117 (alpha * alpha - beta * beta + tem << 117 newton = (nJacobi*(alpha - beta - temp*root)*newton1 +
118 c = 2.0 * (j - 1 + alpha) * (j - << 118 2.0*(nJacobi + alpha)*(nJacobi + beta)*newton2)/
119 newton1 = (b * newton2 - c * newton3) << 119 (temp*(1.0 - root*root)) ;
120 } << 120 rootTemp = root ;
121 newton0 = (nJacobi * (alpha - beta - tem << 121 root = rootTemp - newton1/newton ;
122 2.0 * (nJacobi + alpha) * (nJ << 122 if (std::fabs(root-rootTemp) <= tolerance)
123 (temp * (1.0 - root * root)); << 123 {
124 rootTemp = root; << 124 break ;
125 root = rootTemp - newton1 / newton0; << 125 }
126 if(std::fabs(root - rootTemp) <= toleran << 126 }
127 { << 127 if (k > maxNumber)
128 break; << 128 {
129 } << 129 G4Exception("G4GaussJacobiQ::G4GaussJacobiQ()", "OutOfRange",
130 } << 130 FatalException, "Too many iterations in constructor.") ;
131 if(k > maxNumber) << 131 }
132 { << 132 fAbscissa[i-1] = root ;
133 G4Exception("G4GaussJacobiQ::G4GaussJaco << 133 fWeight[i-1] = std::exp(GammaLogarithm((G4double)(alpha+nJacobi)) +
134 FatalException, "Too many it << 134 GammaLogarithm((G4double)(beta+nJacobi)) -
135 } << 135 GammaLogarithm((G4double)(nJacobi+1.0)) -
136 fAbscissa[i - 1] = root; << 136 GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0)))
137 fWeight[i - 1] = << 137 *temp*std::pow(2.0,alphaBeta)/(newton*newton2) ;
138 std::exp(GammaLogarithm((G4double)(alpha <<
139 GammaLogarithm((G4double)(beta <<
140 GammaLogarithm((G4double)(nJaco <<
141 GammaLogarithm((G4double)(nJaco <<
142 temp * std::pow(2.0, alphaBeta) / (newto <<
143 } 138 }
144 } 139 }
145 140
>> 141
146 // ------------------------------------------- 142 // ----------------------------------------------------------
147 // 143 //
148 // Gauss-Jacobi method for integration of 144 // Gauss-Jacobi method for integration of
149 // ((1-x)^alpha)*((1+x)^beta)*pFunction(x) 145 // ((1-x)^alpha)*((1+x)^beta)*pFunction(x)
150 // from minus unit to plus unit . 146 // from minus unit to plus unit .
151 147
152 G4double G4GaussJacobiQ::Integral() const << 148
>> 149 G4double
>> 150 G4GaussJacobiQ::Integral() const
153 { 151 {
154 G4double integral = 0.0; << 152 G4double integral = 0.0 ;
155 for(G4int i = 0; i < fNumber; ++i) << 153 for(G4int i=0;i<fNumber;i++)
156 { << 154 {
157 integral += fWeight[i] * fFunction(fAbscis << 155 integral += fWeight[i]*fFunction(fAbscissa[i]) ;
158 } << 156 }
159 return integral; << 157 return integral ;
160 } 158 }
>> 159
161 160