Geant4 Cross Reference |
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 // G4GaussJacobiQ class implementation 27 // 28 // Author: V.Grichine, 13.05.1997 29 // -------------------------------------------------------------------- 30 31 #include "G4GaussJacobiQ.hh" 32 33 // ------------------------------------------------------------- 34 // 35 // Constructor for Gauss-Jacobi integration method. 36 // 37 38 G4GaussJacobiQ::G4GaussJacobiQ(function pFunction, G4double alpha, 39 G4double beta, G4int nJacobi) 40 : G4VGaussianQuadrature(pFunction) 41 42 { 43 const G4double tolerance = 1.0e-12; 44 const G4double maxNumber = 12; 45 G4int i = 1, k = 1; 46 G4double root = 0.; 47 G4double alphaBeta = 0.0, alphaReduced = 0.0, betaReduced = 0.0, root1 = 0.0, 48 root2 = 0.0, root3 = 0.0; 49 G4double a = 0.0, b = 0.0, c = 0.0, newton1 = 0.0, newton2 = 0.0, 50 newton3 = 0.0, newton0 = 0.0, temp = 0.0, rootTemp = 0.0; 51 52 fNumber = nJacobi; 53 fAbscissa = new G4double[fNumber]; 54 fWeight = new G4double[fNumber]; 55 56 for(i = 1; i <= nJacobi; ++i) 57 { 58 if(i == 1) 59 { 60 alphaReduced = alpha / nJacobi; 61 betaReduced = beta / nJacobi; 62 root1 = (1.0 + alpha) * (2.78002 / (4.0 + nJacobi * nJacobi) + 63 0.767999 * alphaReduced / nJacobi); 64 root2 = 1.0 + 1.48 * alphaReduced + 0.96002 * betaReduced + 65 0.451998 * alphaReduced * alphaReduced + 66 0.83001 * alphaReduced * betaReduced; 67 root = 1.0 - root1 / root2; 68 } 69 else if(i == 2) 70 { 71 root1 = (4.1002 + alpha) / ((1.0 + alpha) * (1.0 + 0.155998 * alpha)); 72 root2 = 1.0 + 0.06 * (nJacobi - 8.0) * (1.0 + 0.12 * alpha) / nJacobi; 73 root3 = 74 1.0 + 0.012002 * beta * (1.0 + 0.24997 * std::fabs(alpha)) / nJacobi; 75 root -= (1.0 - root) * root1 * root2 * root3; 76 } 77 else if(i == 3) 78 { 79 root1 = (1.67001 + 0.27998 * alpha) / (1.0 + 0.37002 * alpha); 80 root2 = 1.0 + 0.22 * (nJacobi - 8.0) / nJacobi; 81 root3 = 1.0 + 8.0 * beta / ((6.28001 + beta) * nJacobi * nJacobi); 82 root -= (fAbscissa[0] - root) * root1 * root2 * root3; 83 } 84 else if(i == nJacobi - 1) 85 { 86 root1 = (1.0 + 0.235002 * beta) / (0.766001 + 0.118998 * beta); 87 root2 = 1.0 / (1.0 + 0.639002 * (nJacobi - 4.0) / 88 (1.0 + 0.71001 * (nJacobi - 4.0))); 89 root3 = 1.0 / (1.0 + 20.0 * alpha / ((7.5 + alpha) * nJacobi * nJacobi)); 90 root += (root - fAbscissa[nJacobi - 4]) * root1 * root2 * root3; 91 } 92 else if(i == nJacobi) 93 { 94 root1 = (1.0 + 0.37002 * beta) / (1.67001 + 0.27998 * beta); 95 root2 = 1.0 / (1.0 + 0.22 * (nJacobi - 8.0) / nJacobi); 96 root3 = 97 1.0 / (1.0 + 8.0 * alpha / ((6.28002 + alpha) * nJacobi * nJacobi)); 98 root += (root - fAbscissa[nJacobi - 3]) * root1 * root2 * root3; 99 } 100 else 101 { 102 root = 3.0 * fAbscissa[i - 2] - 3.0 * fAbscissa[i - 3] + fAbscissa[i - 4]; 103 } 104 alphaBeta = alpha + beta; 105 for(k = 1; k <= maxNumber; ++k) 106 { 107 temp = 2.0 + alphaBeta; 108 newton1 = (alpha - beta + temp * root) / 2.0; 109 newton2 = 1.0; 110 for(G4int j = 2; j <= nJacobi; ++j) 111 { 112 newton3 = newton2; 113 newton2 = newton1; 114 temp = 2 * j + alphaBeta; 115 a = 2 * j * (j + alphaBeta) * (temp - 2.0); 116 b = (temp - 1.0) * 117 (alpha * alpha - beta * beta + temp * (temp - 2.0) * root); 118 c = 2.0 * (j - 1 + alpha) * (j - 1 + beta) * temp; 119 newton1 = (b * newton2 - c * newton3) / a; 120 } 121 newton0 = (nJacobi * (alpha - beta - temp * root) * newton1 + 122 2.0 * (nJacobi + alpha) * (nJacobi + beta) * newton2) / 123 (temp * (1.0 - root * root)); 124 rootTemp = root; 125 root = rootTemp - newton1 / newton0; 126 if(std::fabs(root - rootTemp) <= tolerance) 127 { 128 break; 129 } 130 } 131 if(k > maxNumber) 132 { 133 G4Exception("G4GaussJacobiQ::G4GaussJacobiQ()", "OutOfRange", 134 FatalException, "Too many iterations in constructor."); 135 } 136 fAbscissa[i - 1] = root; 137 fWeight[i - 1] = 138 std::exp(GammaLogarithm((G4double)(alpha + nJacobi)) + 139 GammaLogarithm((G4double)(beta + nJacobi)) - 140 GammaLogarithm((G4double)(nJacobi + 1.0)) - 141 GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0))) * 142 temp * std::pow(2.0, alphaBeta) / (newton0 * newton2); 143 } 144 } 145 146 // ---------------------------------------------------------- 147 // 148 // Gauss-Jacobi method for integration of 149 // ((1-x)^alpha)*((1+x)^beta)*pFunction(x) 150 // from minus unit to plus unit . 151 152 G4double G4GaussJacobiQ::Integral() const 153 { 154 G4double integral = 0.0; 155 for(G4int i = 0; i < fNumber; ++i) 156 { 157 integral += fWeight[i] * fFunction(fAbscissa[i]); 158 } 159 return integral; 160 } 161