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