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 // G4SimpleIntegration class implementation << 8 // $Id: G4SimpleIntegration.cc,v 1.2 1999/11/16 17:31:11 gcosmo Exp $ >> 9 // GEANT4 tag $Name: geant4-02-00 $ >> 10 // >> 11 // Implementation file for simple integration methods 27 // 12 // 28 // Author: V.Grichine, 26.03.1997 << 29 // ------------------------------------------- << 30 13 31 #include "G4SimpleIntegration.hh" 14 #include "G4SimpleIntegration.hh" 32 #include "globals.hh" << 33 15 34 G4SimpleIntegration::G4SimpleIntegration(funct << 16 35 : fFunction(pFunction) << 17 G4int G4SimpleIntegration::fMaxDepth = 100 ; 36 {} << 18 37 << 19 38 G4SimpleIntegration::G4SimpleIntegration(funct << 20 G4SimpleIntegration::G4SimpleIntegration( function pFunction ) 39 G4dou << 21 { 40 : fFunction(pFunction) << 22 fFunction = pFunction ; 41 , fTolerance(pTolerance) << 23 } 42 {} << 24 43 << 25 G4SimpleIntegration::G4SimpleIntegration( function pFunction, 44 // Simple integration methods << 26 G4double pTolerance) 45 << 27 { 46 G4double G4SimpleIntegration::Trapezoidal(G4do << 28 fFunction = pFunction ; 47 G4in << 29 fTolerance = pTolerance ; 48 { << 30 } 49 G4double Step = (xFinal - xInitial) / iterat << 31 50 G4double mean = (fFunction(xInitial) + fFunc << 32 51 G4double x = xInitial; << 33 G4SimpleIntegration::~G4SimpleIntegration() 52 for(G4int i = 1; i < iterationNumber; ++i) << 34 { 53 { << 35 ; 54 x += Step; << 36 } 55 mean += fFunction(x); << 37 56 } << 38 // Simple integration methods 57 return mean * Step; << 39 58 } << 40 G4double 59 << 41 G4SimpleIntegration::Trapezoidal(G4double xInitial, 60 G4double G4SimpleIntegration::MidPoint(G4doubl << 42 G4double xFinal, 61 G4int i << 43 G4int iterationNumber ) 62 { << 44 { 63 G4double Step = (xFinal - xInitial) / iterat << 45 G4int i ; 64 G4double x = xInitial + 0.5 * Step; << 46 G4double Step = (xFinal - xInitial)/iterationNumber ; 65 G4double mean = fFunction(x); << 47 G4double mean = (fFunction(xInitial) + fFunction(xFinal))*0.5 ; 66 for(G4int i = 1; i < iterationNumber; ++i) << 48 G4double x = xInitial ; 67 { << 49 for(i=1;i<iterationNumber;i++) 68 x += Step; << 50 { 69 mean += fFunction(x); << 51 x += Step ; 70 } << 52 mean += fFunction(x) ; 71 return mean * Step; << 53 } 72 } << 54 return mean*Step ; 73 << 55 } 74 G4double G4SimpleIntegration::Gauss(G4double x << 56 75 G4int iter << 57 G4double 76 { << 58 G4SimpleIntegration::MidPoint(G4double xInitial, 77 G4double x = 0.; << 59 G4double xFinal, 78 static const G4double root = 1.0 / std::sqrt << 60 G4int iterationNumber ) 79 G4double Step = (xFinal - xInit << 61 { 80 G4double delta = Step * root; << 62 G4int i ; 81 G4double mean = 0.0; << 63 G4double Step = (xFinal - xInitial)/iterationNumber ; 82 for(G4int i = 0; i < iterationNumber; ++i) << 64 G4double x = xInitial + 0.5*Step; 83 { << 65 G4double mean = fFunction(x) ; 84 x = (2 * i + 1) * Step; << 66 for(i=1;i<iterationNumber;i++) 85 mean += (fFunction(x + delta) + fFunction( << 67 { 86 } << 68 x += Step ; 87 return mean * Step; << 69 mean += fFunction(x) ; 88 } << 70 } 89 << 71 return mean*Step ; 90 G4double G4SimpleIntegration::Simpson(G4double << 72 } 91 G4int it << 73 92 { << 74 G4double 93 G4double Step = (xFinal - xInitial) / itera << 75 G4SimpleIntegration::Gauss(G4double xInitial, 94 G4double x = xInitial; << 76 G4double xFinal, 95 G4double xPlus = xInitial + 0.5 * Step; << 77 G4int iterationNumber ) 96 G4double mean = (fFunction(xInitial) + fFun << 78 { 97 G4double sum = fFunction(xPlus); << 79 G4int i ; 98 for(G4int i = 1; i < iterationNumber; ++i) << 80 G4double x ; 99 { << 81 static G4double root = 1.0/sqrt(3.0) ; 100 x += Step; << 82 G4double Step = (xFinal - xInitial)/(2.0*iterationNumber) ; 101 xPlus += Step; << 83 G4double delta = Step*root ; 102 mean += fFunction(x); << 84 G4double mean = 0.0 ; 103 sum += fFunction(xPlus); << 85 for(i=0;i<iterationNumber;i++) 104 } << 86 { 105 mean += 2.0 * sum; << 87 x = (2*i + 1)*Step ; 106 return mean * Step / 3.0; << 88 mean += (fFunction(x+delta) + fFunction(x-delta)) ; 107 } << 89 } 108 << 90 return mean*Step ; 109 // Adaptive Gauss integration << 91 } 110 << 92 111 G4double G4SimpleIntegration::AdaptGaussIntegr << 93 G4double 112 << 94 G4SimpleIntegration::Simpson(G4double xInitial, 113 { << 95 G4double xFinal, 114 G4int depth = 0; << 96 G4int iterationNumber ) 115 G4double sum = 0.0; << 97 { 116 AdaptGauss(xInitial, xFinal, sum, depth); << 98 G4int i ; 117 return sum; << 99 G4double Step = (xFinal - xInitial)/iterationNumber ; 118 } << 100 G4double x = xInitial ; 119 << 101 G4double xPlus = xInitial + 0.5*Step ; 120 G4double G4SimpleIntegration::Gauss(G4double x << 102 G4double mean = (fFunction(xInitial) + fFunction(xFinal))*0.5 ; 121 { << 103 G4double sum = fFunction(xPlus) ; 122 static const G4double root = 1.0 / std::sqrt << 104 for(i=1;i<iterationNumber;i++) 123 << 105 { 124 G4double xMean = (xInitial + xFinal) / 2.0; << 106 x += Step ; 125 G4double Step = (xFinal - xInitial) / 2.0; << 107 xPlus += Step ; 126 G4double delta = Step * root; << 108 mean += fFunction(x) ; 127 G4double sum = (fFunction(xMean + delta) + << 109 sum += fFunction(xPlus) ; 128 << 110 } 129 return sum * Step; << 111 mean += 2.0*sum ; 130 } << 112 return mean*Step/3.0 ; 131 << 113 } 132 void G4SimpleIntegration::AdaptGauss(G4double << 114 133 G4double& << 115 134 { << 116 135 if(depth > fMaxDepth) << 117 // Adaptive Gauss integration 136 { << 118 137 G4Exception("G4SimpleIntegration::AdaptGau << 119 G4double 138 "Function varies too rapidly ! << 120 G4SimpleIntegration::AdaptGaussIntegration( G4double xInitial, 139 } << 121 G4double xFinal ) 140 G4double xMean = (xInitial + xFinal) / 2 << 122 { 141 G4double leftHalf = Gauss(xInitial, xMean); << 123 G4int depth = 0 ; 142 G4double rightHalf = Gauss(xMean, xFinal); << 124 G4double sum = 0.0 ; 143 G4double full = Gauss(xInitial, xFinal) << 125 AdaptGauss(xInitial,xFinal,sum,depth) ; 144 if(std::fabs(leftHalf + rightHalf - full) < << 126 return sum ; 145 { << 127 } 146 sum += full; << 128 147 } << 129 148 else << 130 G4double 149 { << 131 G4SimpleIntegration::Gauss( G4double xInitial, 150 ++depth; << 132 G4double xFinal ) 151 AdaptGauss(xInitial, xMean, sum, depth); << 133 { 152 AdaptGauss(xMean, xFinal, sum, depth); << 134 static G4double root = 1.0/sqrt(3.0) ; 153 } << 135 >> 136 G4double xMean = (xInitial + xFinal)/2.0 ; >> 137 G4double Step = (xFinal - xInitial)/2.0 ; >> 138 G4double delta = Step*root ; >> 139 G4double sum = (fFunction(xMean + delta) + fFunction(xMean - delta)) ; >> 140 >> 141 return sum*Step ; >> 142 } >> 143 >> 144 >> 145 void >> 146 G4SimpleIntegration::AdaptGauss( G4double xInitial, >> 147 G4double xFinal, >> 148 G4double& sum, >> 149 G4int& depth ) >> 150 { >> 151 if(depth >fMaxDepth) >> 152 { >> 153 G4Exception("Function varies too rapidly in G4SimpleIntegration::AdaptGauss") ; >> 154 } >> 155 G4double xMean = (xInitial + xFinal)/2.0 ; >> 156 G4double leftHalf = Gauss(xInitial,xMean) ; >> 157 G4double rightHalf = Gauss(xMean,xFinal) ; >> 158 G4double full = Gauss(xInitial,xFinal) ; >> 159 if(fabs(leftHalf+rightHalf-full) < fTolerance) >> 160 { >> 161 sum += full ; >> 162 } >> 163 else >> 164 { >> 165 depth++ ; >> 166 AdaptGauss(xInitial,xMean,sum,depth) ; >> 167 AdaptGauss(xMean,xFinal,sum,depth) ; >> 168 } 154 } 169 } 155 170