Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // G4Integrator inline methods implementation << 27 // 26 // 28 // Author: V.Grichine, 04.09.1999 - First impl << 27 // 29 // G4SimpleIntegration class with H.P. << 28 // Implementation of G4Integrator methods. 30 // E.TCherniaev advises << 29 // 31 // ------------------------------------------- << 30 // 32 31 33 ////////////////////////////////////////////// 32 ///////////////////////////////////////////////////////////////////// 34 // 33 // 35 // Sympson integration method 34 // Sympson integration method 36 // 35 // 37 ////////////////////////////////////////////// 36 ///////////////////////////////////////////////////////////////////// 38 // 37 // 39 // Integration of class member functions T::f << 38 // Integration of class member functions T::f by Simpson method. 40 39 41 template <class T, class F> << 40 template <class T, class F> 42 G4double G4Integrator<T, F>::Simpson(T& typeT, << 41 G4double G4Integrator<T,F>::Simpson( T& typeT, 43 G4double << 42 F f, 44 { << 43 G4double xInitial, 45 G4int i; << 44 G4double xFinal, 46 G4double step = (xFinal - xInitial) / itera << 45 G4int iterationNumber ) 47 G4double x = xInitial; << 46 { 48 G4double xPlus = xInitial + 0.5 * step; << 47 G4int i ; 49 G4double mean = ((typeT.*f)(xInitial) + (ty << 48 G4double step = (xFinal - xInitial)/iterationNumber ; 50 G4double sum = (typeT.*f)(xPlus); << 49 G4double x = xInitial ; 51 << 50 G4double xPlus = xInitial + 0.5*step ; 52 for(i = 1; i < iterationNumber; ++i) << 51 G4double mean = ( (typeT.*f)(xInitial) + (typeT.*f)(xFinal) )*0.5 ; 53 { << 52 G4double sum = (typeT.*f)(xPlus) ; 54 x += step; << 53 55 xPlus += step; << 54 for(i=1;i<iterationNumber;i++) 56 mean += (typeT.*f)(x); << 55 { 57 sum += (typeT.*f)(xPlus); << 56 x += step ; 58 } << 57 xPlus += step ; 59 mean += 2.0 * sum; << 58 mean += (typeT.*f)(x) ; >> 59 sum += (typeT.*f)(xPlus) ; >> 60 } >> 61 mean += 2.0*sum ; 60 62 61 return mean * step / 3.0; << 63 return mean*step/3.0 ; 62 } 64 } 63 65 64 ////////////////////////////////////////////// 66 ///////////////////////////////////////////////////////////////////// 65 // 67 // 66 // Integration of class member functions T::f 68 // Integration of class member functions T::f by Simpson method. 67 // Convenient to use with 'this' pointer 69 // Convenient to use with 'this' pointer 68 70 69 template <class T, class F> << 71 template <class T, class F> 70 G4double G4Integrator<T, F>::Simpson(T* ptrT, << 72 G4double G4Integrator<T,F>::Simpson( T* ptrT, 71 G4double << 73 F f, 72 { << 74 G4double xInitial, 73 G4int i; << 75 G4double xFinal, 74 G4double step = (xFinal - xInitial) / itera << 76 G4int iterationNumber ) 75 G4double x = xInitial; << 77 { 76 G4double xPlus = xInitial + 0.5 * step; << 78 G4int i ; 77 G4double mean = ((ptrT->*f)(xInitial) + (pt << 79 G4double step = (xFinal - xInitial)/iterationNumber ; 78 G4double sum = (ptrT->*f)(xPlus); << 80 G4double x = xInitial ; >> 81 G4double xPlus = xInitial + 0.5*step ; >> 82 G4double mean = ( (ptrT->*f)(xInitial) + (ptrT->*f)(xFinal) )*0.5 ; >> 83 G4double sum = (ptrT->*f)(xPlus) ; >> 84 >> 85 for(i=1;i<iterationNumber;i++) >> 86 { >> 87 x += step ; >> 88 xPlus += step ; >> 89 mean += (ptrT->*f)(x) ; >> 90 sum += (ptrT->*f)(xPlus) ; >> 91 } >> 92 mean += 2.0*sum ; 79 93 80 for(i = 1; i < iterationNumber; ++i) << 94 return mean*step/3.0 ; 81 { << 82 x += step; << 83 xPlus += step; << 84 mean += (ptrT->*f)(x); << 85 sum += (ptrT->*f)(xPlus); << 86 } << 87 mean += 2.0 * sum; << 88 << 89 return mean * step / 3.0; << 90 } 95 } 91 96 92 ////////////////////////////////////////////// 97 ///////////////////////////////////////////////////////////////////// 93 // 98 // 94 // Integration of class member functions T::f 99 // Integration of class member functions T::f by Simpson method. 95 // Convenient to use, when function f is defin 100 // Convenient to use, when function f is defined in global scope, i.e. in main() 96 // program 101 // program 97 102 98 template <class T, class F> << 103 template <class T, class F> 99 G4double G4Integrator<T, F>::Simpson(G4double << 104 G4double G4Integrator<T,F>::Simpson( G4double (*f)(G4double), 100 G4double << 105 G4double xInitial, 101 { << 106 G4double xFinal, 102 G4int i; << 107 G4int iterationNumber ) 103 G4double step = (xFinal - xInitial) / itera << 108 { 104 G4double x = xInitial; << 109 G4int i ; 105 G4double xPlus = xInitial + 0.5 * step; << 110 G4double step = (xFinal - xInitial)/iterationNumber ; 106 G4double mean = ((*f)(xInitial) + (*f)(xFin << 111 G4double x = xInitial ; 107 G4double sum = (*f)(xPlus); << 112 G4double xPlus = xInitial + 0.5*step ; 108 << 113 G4double mean = ( (*f)(xInitial) + (*f)(xFinal) )*0.5 ; 109 for(i = 1; i < iterationNumber; ++i) << 114 G4double sum = (*f)(xPlus) ; 110 { << 115 111 x += step; << 116 for(i=1;i<iterationNumber;i++) 112 xPlus += step; << 117 { 113 mean += (*f)(x); << 118 x += step ; 114 sum += (*f)(xPlus); << 119 xPlus += step ; 115 } << 120 mean += (*f)(x) ; 116 mean += 2.0 * sum; << 121 sum += (*f)(xPlus) ; >> 122 } >> 123 mean += 2.0*sum ; 117 124 118 return mean * step / 3.0; << 125 return mean*step/3.0 ; 119 } 126 } 120 127 121 ////////////////////////////////////////////// 128 ////////////////////////////////////////////////////////////////////////// 122 // 129 // 123 // Adaptive Gauss method 130 // Adaptive Gauss method 124 // 131 // 125 ////////////////////////////////////////////// 132 ////////////////////////////////////////////////////////////////////////// 126 // 133 // 127 // 134 // 128 135 129 template <class T, class F> << 136 template <class T, class F> 130 G4double G4Integrator<T, F>::Gauss(T& typeT, F << 137 G4double G4Integrator<T,F>::Gauss( T& typeT, F f, 131 G4double xF << 138 G4double xInitial, G4double xFinal ) 132 { << 139 { 133 static const G4double root = 1.0 / std::sqrt << 140 static const G4double root = 1.0/std::sqrt(3.0) ; 134 << 141 135 G4double xMean = (xInitial + xFinal) / 2.0; << 142 G4double xMean = (xInitial + xFinal)/2.0 ; 136 G4double Step = (xFinal - xInitial) / 2.0; << 143 G4double Step = (xFinal - xInitial)/2.0 ; 137 G4double delta = Step * root; << 144 G4double delta = Step*root ; 138 G4double sum = ((typeT.*f)(xMean + delta) << 145 G4double sum = ((typeT.*f)(xMean + delta) + 139 << 146 (typeT.*f)(xMean - delta)) ; 140 return sum * Step; << 147 >> 148 return sum*Step ; 141 } 149 } 142 150 143 ////////////////////////////////////////////// 151 ////////////////////////////////////////////////////////////////////// 144 // 152 // 145 // 153 // 146 154 147 template <class T, class F> << 155 template <class T, class F> G4double 148 G4double G4Integrator<T, F>::Gauss(T* ptrT, F << 156 G4Integrator<T,F>::Gauss( T* ptrT, F f, G4double a, G4double b ) 149 { 157 { 150 return Gauss(*ptrT, f, a, b); << 158 return Gauss(*ptrT,f,a,b) ; 151 } 159 } 152 160 153 ////////////////////////////////////////////// 161 /////////////////////////////////////////////////////////////////////// 154 // 162 // 155 // 163 // 156 164 157 template <class T, class F> 165 template <class T, class F> 158 G4double G4Integrator<T, F>::Gauss(G4double (* << 166 G4double G4Integrator<T,F>::Gauss( G4double (*f)(G4double), 159 G4double xF << 167 G4double xInitial, G4double xFinal) 160 { 168 { 161 static const G4double root = 1.0 / std::sqrt << 169 static const G4double root = 1.0/std::sqrt(3.0) ; 162 170 163 G4double xMean = (xInitial + xFinal) / 2.0; << 171 G4double xMean = (xInitial + xFinal)/2.0 ; 164 G4double Step = (xFinal - xInitial) / 2.0; << 172 G4double Step = (xFinal - xInitial)/2.0 ; 165 G4double delta = Step * root; << 173 G4double delta = Step*root ; 166 G4double sum = ((*f)(xMean + delta) + (*f) << 174 G4double sum = ( (*f)(xMean + delta) + (*f)(xMean - delta) ) ; 167 << 175 168 return sum * Step; << 176 return sum*Step ; 169 } 177 } 170 178 171 ////////////////////////////////////////////// 179 /////////////////////////////////////////////////////////////////////////// 172 // 180 // 173 // 181 // 174 182 175 template <class T, class F> << 183 template <class T, class F> 176 void G4Integrator<T, F>::AdaptGauss(T& typeT, << 184 void G4Integrator<T,F>::AdaptGauss( T& typeT, F f, G4double xInitial, 177 G4double x << 185 G4double xFinal, G4double fTolerance, 178 G4double& << 186 G4double& sum, >> 187 G4int& depth ) >> 188 { >> 189 if(depth > 100) >> 190 { >> 191 G4cout<<"G4Integrator<T,F>::AdaptGauss: WARNING !!!"<<G4endl ; >> 192 G4cout<<"Function varies too rapidly to get stated accuracy in 100 steps " >> 193 <<G4endl ; >> 194 >> 195 return ; >> 196 } >> 197 G4double xMean = (xInitial + xFinal)/2.0 ; >> 198 G4double leftHalf = Gauss(typeT,f,xInitial,xMean) ; >> 199 G4double rightHalf = Gauss(typeT,f,xMean,xFinal) ; >> 200 G4double full = Gauss(typeT,f,xInitial,xFinal) ; >> 201 if(std::fabs(leftHalf+rightHalf-full) < fTolerance) >> 202 { >> 203 sum += full ; >> 204 } >> 205 else >> 206 { >> 207 depth++ ; >> 208 AdaptGauss(typeT,f,xInitial,xMean,fTolerance,sum,depth) ; >> 209 AdaptGauss(typeT,f,xMean,xFinal,fTolerance,sum,depth) ; >> 210 } >> 211 } >> 212 >> 213 template <class T, class F> >> 214 void G4Integrator<T,F>::AdaptGauss( T* ptrT, F f, G4double xInitial, >> 215 G4double xFinal, G4double fTolerance, >> 216 G4double& sum, >> 217 G4int& depth ) 179 { 218 { 180 if(depth > 100) << 219 AdaptGauss(*ptrT,f,xInitial,xFinal,fTolerance,sum,depth) ; 181 { << 182 G4cout << "G4Integrator<T,F>::AdaptGauss: << 183 G4cout << "Function varies too rapidly to << 184 << G4endl; << 185 << 186 return; << 187 } << 188 G4double xMean = (xInitial + xFinal) / 2 << 189 G4double leftHalf = Gauss(typeT, f, xInitia << 190 G4double rightHalf = Gauss(typeT, f, xMean, << 191 G4double full = Gauss(typeT, f, xInitia << 192 if(std::fabs(leftHalf + rightHalf - full) < << 193 { << 194 sum += full; << 195 } << 196 else << 197 { << 198 ++depth; << 199 AdaptGauss(typeT, f, xInitial, xMean, fTol << 200 AdaptGauss(typeT, f, xMean, xFinal, fToler << 201 } << 202 } << 203 << 204 template <class T, class F> << 205 void G4Integrator<T, F>::AdaptGauss(T* ptrT, F << 206 G4double x << 207 G4double& << 208 { << 209 AdaptGauss(*ptrT, f, xInitial, xFinal, fTole << 210 } 220 } 211 221 212 ////////////////////////////////////////////// 222 ///////////////////////////////////////////////////////////////////////// 213 // 223 // 214 // 224 // 215 template <class T, class F> 225 template <class T, class F> 216 void G4Integrator<T, F>::AdaptGauss(G4double ( << 226 void G4Integrator<T,F>::AdaptGauss( G4double (*f)(G4double), 217 G4double x << 227 G4double xInitial, G4double xFinal, 218 G4double& << 228 G4double fTolerance, G4double& sum, 219 { << 229 G4int& depth ) 220 if(depth > 100) << 230 { 221 { << 231 if(depth > 100) 222 G4cout << "G4SimpleIntegration::AdaptGauss << 232 { 223 G4cout << "Function varies too rapidly to << 233 G4cout<<"G4SimpleIntegration::AdaptGauss: WARNING !!!"<<G4endl ; 224 << G4endl; << 234 G4cout<<"Function varies too rapidly to get stated accuracy in 100 steps " 225 << 235 <<G4endl ; 226 return; << 236 227 } << 237 return ; 228 G4double xMean = (xInitial + xFinal) / 2 << 238 } 229 G4double leftHalf = Gauss(f, xInitial, xMea << 239 G4double xMean = (xInitial + xFinal)/2.0 ; 230 G4double rightHalf = Gauss(f, xMean, xFinal) << 240 G4double leftHalf = Gauss(f,xInitial,xMean) ; 231 G4double full = Gauss(f, xInitial, xFin << 241 G4double rightHalf = Gauss(f,xMean,xFinal) ; 232 if(std::fabs(leftHalf + rightHalf - full) < << 242 G4double full = Gauss(f,xInitial,xFinal) ; 233 { << 243 if(std::fabs(leftHalf+rightHalf-full) < fTolerance) 234 sum += full; << 244 { 235 } << 245 sum += full ; 236 else << 246 } 237 { << 247 else 238 ++depth; << 248 { 239 AdaptGauss(f, xInitial, xMean, fTolerance, << 249 depth++ ; 240 AdaptGauss(f, xMean, xFinal, fTolerance, s << 250 AdaptGauss(f,xInitial,xMean,fTolerance,sum,depth) ; 241 } << 251 AdaptGauss(f,xMean,xFinal,fTolerance,sum,depth) ; >> 252 } 242 } 253 } 243 254 244 ////////////////////////////////////////////// 255 //////////////////////////////////////////////////////////////////////// 245 // 256 // 246 // Adaptive Gauss integration with accuracy 'e 257 // Adaptive Gauss integration with accuracy 'e' 247 // Convenient for using with class object type 258 // Convenient for using with class object typeT 248 << 259 249 template <class T, class F> << 260 template<class T, class F> 250 G4double G4Integrator<T, F>::AdaptiveGauss(T& << 261 G4double G4Integrator<T,F>::AdaptiveGauss( T& typeT, F f, G4double xInitial, 251 G4d << 262 G4double xFinal, G4double e ) 252 { << 263 { 253 G4int depth = 0; << 264 G4int depth = 0 ; 254 G4double sum = 0.0; << 265 G4double sum = 0.0 ; 255 AdaptGauss(typeT, f, xInitial, xFinal, e, su << 266 AdaptGauss(typeT,f,xInitial,xFinal,e,sum,depth) ; 256 return sum; << 267 return sum ; 257 } 268 } 258 269 259 ////////////////////////////////////////////// 270 //////////////////////////////////////////////////////////////////////// 260 // 271 // 261 // Adaptive Gauss integration with accuracy 'e 272 // Adaptive Gauss integration with accuracy 'e' 262 // Convenient for using with 'this' pointer 273 // Convenient for using with 'this' pointer 263 << 274 264 template <class T, class F> << 275 template<class T, class F> 265 G4double G4Integrator<T, F>::AdaptiveGauss(T* << 276 G4double G4Integrator<T,F>::AdaptiveGauss( T* ptrT, F f, G4double xInitial, 266 G4d << 277 G4double xFinal, G4double e ) 267 { 278 { 268 return AdaptiveGauss(*ptrT, f, xInitial, xFi << 279 return AdaptiveGauss(*ptrT,f,xInitial,xFinal,e) ; 269 } 280 } 270 281 271 ////////////////////////////////////////////// 282 //////////////////////////////////////////////////////////////////////// 272 // 283 // 273 // Adaptive Gauss integration with accuracy 'e 284 // Adaptive Gauss integration with accuracy 'e' 274 // Convenient for using with global scope func 285 // Convenient for using with global scope function f 275 << 286 276 template <class T, class F> 287 template <class T, class F> 277 G4double G4Integrator<T, F>::AdaptiveGauss(G4d << 288 G4double G4Integrator<T,F>::AdaptiveGauss( G4double (*f)(G4double), 278 G4d << 289 G4double xInitial, G4double xFinal, G4double e ) 279 G4d << 290 { 280 { << 291 G4int depth = 0 ; 281 G4int depth = 0; << 292 G4double sum = 0.0 ; 282 G4double sum = 0.0; << 293 AdaptGauss(f,xInitial,xFinal,e,sum,depth) ; 283 AdaptGauss(f, xInitial, xFinal, e, sum, dept << 294 return sum ; 284 return sum; << 285 } 295 } 286 296 287 ////////////////////////////////////////////// 297 //////////////////////////////////////////////////////////////////////////// 288 // Gauss integration methods involving ortogon 298 // Gauss integration methods involving ortogonal polynomials 289 ////////////////////////////////////////////// 299 //////////////////////////////////////////////////////////////////////////// 290 // 300 // 291 // Methods involving Legendre polynomials << 301 // Methods involving Legendre polynomials 292 // 302 // 293 ////////////////////////////////////////////// 303 ///////////////////////////////////////////////////////////////////////// 294 // 304 // 295 // The value nLegendre set the accuracy requir 305 // The value nLegendre set the accuracy required, i.e the number of points 296 // where the function pFunction will be evalua 306 // where the function pFunction will be evaluated during integration. 297 // The function creates the arrays for absciss << 307 // The function creates the arrays for abscissas and weights that used 298 // in Gauss-Legendre quadrature method. << 308 // in Gauss-Legendre quadrature method. 299 // The values a and b are the limits of integr 309 // The values a and b are the limits of integration of the function f . 300 // nLegendre MUST BE EVEN !!! 310 // nLegendre MUST BE EVEN !!! 301 // Returns the integral of the function f betw << 311 // Returns the integral of the function f between a and b, by 2*fNumber point 302 // Gauss-Legendre integration: the function is 312 // Gauss-Legendre integration: the function is evaluated exactly 303 // 2*fNumber times at interior points in the r << 313 // 2*fNumber times at interior points in the range of integration. 304 // Since the weights and abscissas are, in thi << 314 // Since the weights and abscissas are, in this case, symmetric around 305 // the midpoint of the range of integration, t << 315 // the midpoint of the range of integration, there are actually only 306 // fNumber distinct values of each. 316 // fNumber distinct values of each. 307 // Convenient for using with some class object 317 // Convenient for using with some class object dataT 308 318 309 template <class T, class F> 319 template <class T, class F> 310 G4double G4Integrator<T, F>::Legendre(T& typeT << 320 G4double G4Integrator<T,F>::Legendre( T& typeT, F f, G4double a, G4double b, 311 G4int nL << 321 G4int nLegendre ) 312 { 322 { 313 G4double nwt, nwt1, temp1, temp2, temp3, tem << 323 G4double nwt, nwt1, temp1, temp2, temp3, temp ; 314 G4double xDiff, xMean, dx, integral; << 324 G4double xDiff, xMean, dx, integral ; 315 << 316 const G4double tolerance = 1.6e-10; << 317 G4int i, j, k = nLegendre; << 318 G4int fNumber = (nLegendre + 1) / 2; << 319 << 320 if(2 * fNumber != k) << 321 { << 322 G4Exception("G4Integrator<T,F>::Legendre(T << 323 FatalException, "Invalid (odd) << 324 } << 325 325 326 G4double* fAbscissa = new G4double[fNumber]; << 326 const G4double tolerance = 1.6e-10 ; 327 G4double* fWeight = new G4double[fNumber]; << 327 G4int i, j, k = nLegendre ; 328 << 328 G4int fNumber = (nLegendre + 1)/2 ; 329 for(i = 1; i <= fNumber; ++i) // Loop over << 329 330 { << 330 if(2*fNumber != k) 331 nwt = std::cos(CLHEP::pi * (i - 0.25) / << 331 { 332 (k + 0.5)); // Initial roo << 332 G4Exception("G4Integrator<T,F>::Legendre(T&,F, ...)", "InvalidCall", 333 << 333 FatalException, "Invalid (odd) nLegendre in constructor."); 334 do // loop of Newton's method << 334 } 335 { << 335 336 temp1 = 1.0; << 336 G4double* fAbscissa = new G4double[fNumber] ; 337 temp2 = 0.0; << 337 G4double* fWeight = new G4double[fNumber] ; 338 for(j = 1; j <= k; ++j) << 338 339 { << 339 for(i=1;i<=fNumber;i++) // Loop over the desired roots 340 temp3 = temp2; << 340 { 341 temp2 = temp1; << 341 nwt = std::cos(CLHEP::pi*(i - 0.25)/(k + 0.5)) ; // Initial root approximation 342 temp1 = ((2.0 * j - 1.0) * nwt * temp2 << 342 343 } << 343 do // loop of Newton's method 344 temp = k * (nwt * temp1 - temp2) / (nwt << 344 { 345 nwt1 = nwt; << 345 temp1 = 1.0 ; 346 nwt = nwt1 - temp1 / temp; // Newton's << 346 temp2 = 0.0 ; 347 } while(std::fabs(nwt - nwt1) > tolerance) << 347 for(j=1;j<=k;j++) 348 << 348 { 349 fAbscissa[fNumber - i] = nwt; << 349 temp3 = temp2 ; 350 fWeight[fNumber - i] = 2.0 / ((1.0 - nwt << 350 temp2 = temp1 ; 351 } << 351 temp1 = ((2.0*j - 1.0)*nwt*temp2 - (j - 1.0)*temp3)/j ; 352 << 352 } 353 // << 353 temp = k*(nwt*temp1 - temp2)/(nwt*nwt - 1.0) ; 354 // Now we ready to get integral << 354 nwt1 = nwt ; 355 // << 355 nwt = nwt1 - temp1/temp ; // Newton's method 356 << 356 } 357 xMean = 0.5 * (a + b); << 357 while(std::fabs(nwt - nwt1) > tolerance) ; 358 xDiff = 0.5 * (b - a); << 358 359 integral = 0.0; << 359 fAbscissa[fNumber-i] = nwt ; 360 for(i = 0; i < fNumber; ++i) << 360 fWeight[fNumber-i] = 2.0/((1.0 - nwt*nwt)*temp*temp) ; 361 { << 361 } 362 dx = xDiff * fAbscissa[i]; << 362 363 integral += fWeight[i] * ((typeT.*f)(xMean << 363 // 364 } << 364 // Now we ready to get integral 365 delete[] fAbscissa; << 365 // 366 delete[] fWeight; << 366 367 return integral *= xDiff; << 367 xMean = 0.5*(a + b) ; 368 } << 368 xDiff = 0.5*(b - a) ; >> 369 integral = 0.0 ; >> 370 for(i=0;i<fNumber;i++) >> 371 { >> 372 dx = xDiff*fAbscissa[i] ; >> 373 integral += fWeight[i]*( (typeT.*f)(xMean + dx) + >> 374 (typeT.*f)(xMean - dx) ) ; >> 375 } >> 376 delete[] fAbscissa; >> 377 delete[] fWeight; >> 378 return integral *= xDiff ; >> 379 } 369 380 370 ////////////////////////////////////////////// 381 /////////////////////////////////////////////////////////////////////// 371 // 382 // 372 // Convenient for using with the pointer 'this 383 // Convenient for using with the pointer 'this' 373 384 374 template <class T, class F> 385 template <class T, class F> 375 G4double G4Integrator<T, F>::Legendre(T* ptrT, << 386 G4double G4Integrator<T,F>::Legendre( T* ptrT, F f, G4double a, 376 G4int nL << 387 G4double b, G4int nLegendre ) 377 { 388 { 378 return Legendre(*ptrT, f, a, b, nLegendre); << 389 return Legendre(*ptrT,f,a,b,nLegendre) ; 379 } 390 } 380 391 381 ////////////////////////////////////////////// 392 /////////////////////////////////////////////////////////////////////// 382 // 393 // 383 // Convenient for using with global scope func 394 // Convenient for using with global scope function f 384 395 385 template <class T, class F> 396 template <class T, class F> 386 G4double G4Integrator<T, F>::Legendre(G4double << 397 G4double G4Integrator<T,F>::Legendre( G4double (*f)(G4double), 387 G4double << 398 G4double a, G4double b, G4int nLegendre) 388 { 399 { 389 G4double nwt, nwt1, temp1, temp2, temp3, tem << 400 G4double nwt, nwt1, temp1, temp2, temp3, temp ; 390 G4double xDiff, xMean, dx, integral; << 401 G4double xDiff, xMean, dx, integral ; 391 << 392 const G4double tolerance = 1.6e-10; << 393 G4int i, j, k = nLegendre; << 394 G4int fNumber = (nLegendre + 1) / 2; << 395 << 396 if(2 * fNumber != k) << 397 { << 398 G4Exception("G4Integrator<T,F>::Legendre(. << 399 FatalException, "Invalid (odd) << 400 } << 401 402 402 G4double* fAbscissa = new G4double[fNumber]; << 403 const G4double tolerance = 1.6e-10 ; 403 G4double* fWeight = new G4double[fNumber]; << 404 G4int i, j, k = nLegendre ; 404 << 405 G4int fNumber = (nLegendre + 1)/2 ; 405 for(i = 1; i <= fNumber; i++) // Loop over << 406 406 { << 407 if(2*fNumber != k) 407 nwt = std::cos(CLHEP::pi * (i - 0.25) / << 408 { 408 (k + 0.5)); // Initial roo << 409 G4Exception("G4Integrator<T,F>::Legendre(...)", "InvalidCall", 409 << 410 FatalException, "Invalid (odd) nLegendre in constructor."); 410 do // loop of Newton's method << 411 } 411 { << 412 412 temp1 = 1.0; << 413 G4double* fAbscissa = new G4double[fNumber] ; 413 temp2 = 0.0; << 414 G4double* fWeight = new G4double[fNumber] ; 414 for(j = 1; j <= k; ++j) << 415 415 { << 416 for(i=1;i<=fNumber;i++) // Loop over the desired roots 416 temp3 = temp2; << 417 { 417 temp2 = temp1; << 418 nwt = std::cos(CLHEP::pi*(i - 0.25)/(k + 0.5)) ; // Initial root approximation 418 temp1 = ((2.0 * j - 1.0) * nwt * temp2 << 419 419 } << 420 do // loop of Newton's method 420 temp = k * (nwt * temp1 - temp2) / (nwt << 421 { 421 nwt1 = nwt; << 422 temp1 = 1.0 ; 422 nwt = nwt1 - temp1 / temp; // Newton's << 423 temp2 = 0.0 ; 423 } while(std::fabs(nwt - nwt1) > tolerance) << 424 for(j=1;j<=k;j++) 424 << 425 { 425 fAbscissa[fNumber - i] = nwt; << 426 temp3 = temp2 ; 426 fWeight[fNumber - i] = 2.0 / ((1.0 - nwt << 427 temp2 = temp1 ; 427 } << 428 temp1 = ((2.0*j - 1.0)*nwt*temp2 - (j - 1.0)*temp3)/j ; 428 << 429 } 429 // << 430 temp = k*(nwt*temp1 - temp2)/(nwt*nwt - 1.0) ; 430 // Now we ready to get integral << 431 nwt1 = nwt ; 431 // << 432 nwt = nwt1 - temp1/temp ; // Newton's method 432 << 433 } 433 xMean = 0.5 * (a + b); << 434 while(std::fabs(nwt - nwt1) > tolerance) ; 434 xDiff = 0.5 * (b - a); << 435 435 integral = 0.0; << 436 fAbscissa[fNumber-i] = nwt ; 436 for(i = 0; i < fNumber; ++i) << 437 fWeight[fNumber-i] = 2.0/((1.0 - nwt*nwt)*temp*temp) ; 437 { << 438 } 438 dx = xDiff * fAbscissa[i]; << 439 439 integral += fWeight[i] * ((*f)(xMean + dx) << 440 // 440 } << 441 // Now we ready to get integral 441 delete[] fAbscissa; << 442 // 442 delete[] fWeight; << 443 >> 444 xMean = 0.5*(a + b) ; >> 445 xDiff = 0.5*(b - a) ; >> 446 integral = 0.0 ; >> 447 for(i=0;i<fNumber;i++) >> 448 { >> 449 dx = xDiff*fAbscissa[i] ; >> 450 integral += fWeight[i]*( (*f)(xMean + dx) + (*f)(xMean - dx) ) ; >> 451 } >> 452 delete[] fAbscissa; >> 453 delete[] fWeight; 443 454 444 return integral *= xDiff; << 455 return integral *= xDiff ; 445 } << 456 } 446 457 447 ////////////////////////////////////////////// 458 //////////////////////////////////////////////////////////////////////////// 448 // 459 // 449 // Returns the integral of the function to be 460 // Returns the integral of the function to be pointed by T::f between a and b, 450 // by ten point Gauss-Legendre integration: th 461 // by ten point Gauss-Legendre integration: the function is evaluated exactly 451 // ten times at interior points in the range o 462 // ten times at interior points in the range of integration. Since the weights 452 // and abscissas are, in this case, symmetric << 463 // and abscissas are, in this case, symmetric around the midpoint of the 453 // range of integration, there are actually on 464 // range of integration, there are actually only five distinct values of each 454 // Convenient for using with class object type 465 // Convenient for using with class object typeT 455 466 456 template <class T, class F> << 467 template <class T, class F> 457 G4double G4Integrator<T, F>::Legendre10(T& typ << 468 G4double G4Integrator<T,F>::Legendre10( T& typeT, F f,G4double a, G4double b) 458 { 469 { 459 G4int i; << 470 G4int i ; 460 G4double xDiff, xMean, dx, integral; << 471 G4double xDiff, xMean, dx, integral ; 461 << 472 462 // From Abramowitz M., Stegan I.A. 1964 , Ha << 473 // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 916 463 << 474 464 static const G4double abscissa[] = { 0.14887 << 475 static const G4double abscissa[] = { 0.148874338981631, 0.433395394129247, 465 0.67940 << 476 0.679409568299024, 0.865063366688985, 466 0.97390 << 477 0.973906528517172 } ; 467 << 478 468 static const G4double weight[] = { 0.2955242 << 479 static const G4double weight[] = { 0.295524224714753, 0.269266719309996, 469 0.2190863 << 480 0.219086362515982, 0.149451349150581, 470 0.0666713 << 481 0.066671344308688 } ; 471 xMean = 0.5 * (a + << 482 xMean = 0.5*(a + b) ; 472 xDiff = 0.5 * (b - << 483 xDiff = 0.5*(b - a) ; 473 integral = 0.0; << 484 integral = 0.0 ; 474 for(i = 0; i < 5; ++i) << 485 for(i=0;i<5;i++) 475 { << 486 { 476 dx = xDiff * abscissa[i]; << 487 dx = xDiff*abscissa[i] ; 477 integral += weight[i] * ((typeT.*f)(xMean << 488 integral += weight[i]*( (typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx)) ; 478 } << 489 } 479 return integral *= xDiff; << 490 return integral *= xDiff ; 480 } 491 } 481 492 482 ////////////////////////////////////////////// 493 /////////////////////////////////////////////////////////////////////////// 483 // 494 // 484 // Convenient for using with the pointer 'this 495 // Convenient for using with the pointer 'this' 485 496 486 template <class T, class F> << 497 template <class T, class F> 487 G4double G4Integrator<T, F>::Legendre10(T* ptr << 498 G4double G4Integrator<T,F>::Legendre10( T* ptrT, F f,G4double a, G4double b) 488 { 499 { 489 return Legendre10(*ptrT, f, a, b); << 500 return Legendre10(*ptrT,f,a,b) ; 490 } << 501 } 491 502 492 ////////////////////////////////////////////// 503 ////////////////////////////////////////////////////////////////////////// 493 // 504 // 494 // Convenient for using with global scope func 505 // Convenient for using with global scope functions 495 506 496 template <class T, class F> 507 template <class T, class F> 497 G4double G4Integrator<T, F>::Legendre10(G4doub << 508 G4double G4Integrator<T,F>::Legendre10( G4double (*f)(G4double), 498 G4doub << 509 G4double a, G4double b ) 499 { 510 { 500 G4int i; << 511 G4int i ; 501 G4double xDiff, xMean, dx, integral; << 512 G4double xDiff, xMean, dx, integral ; 502 << 513 503 // From Abramowitz M., Stegan I.A. 1964 , Ha << 514 // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 916 504 << 515 505 static const G4double abscissa[] = { 0.14887 << 516 static const G4double abscissa[] = { 0.148874338981631, 0.433395394129247, 506 0.67940 << 517 0.679409568299024, 0.865063366688985, 507 0.97390 << 518 0.973906528517172 } ; 508 << 519 509 static const G4double weight[] = { 0.2955242 << 520 static const G4double weight[] = { 0.295524224714753, 0.269266719309996, 510 0.2190863 << 521 0.219086362515982, 0.149451349150581, 511 0.0666713 << 522 0.066671344308688 } ; 512 xMean = 0.5 * (a + << 523 xMean = 0.5*(a + b) ; 513 xDiff = 0.5 * (b - << 524 xDiff = 0.5*(b - a) ; 514 integral = 0.0; << 525 integral = 0.0 ; 515 for(i = 0; i < 5; ++i) << 526 for(i=0;i<5;i++) 516 { << 527 { 517 dx = xDiff * abscissa[i]; << 528 dx = xDiff*abscissa[i] ; 518 integral += weight[i] * ((*f)(xMean + dx) << 529 integral += weight[i]*( (*f)(xMean + dx) + (*f)(xMean - dx)) ; 519 } << 530 } 520 return integral *= xDiff; << 531 return integral *= xDiff ; 521 } 532 } 522 533 523 ////////////////////////////////////////////// 534 /////////////////////////////////////////////////////////////////////// 524 // 535 // 525 // Returns the integral of the function to be 536 // Returns the integral of the function to be pointed by T::f between a and b, 526 // by 96 point Gauss-Legendre integration: the 537 // by 96 point Gauss-Legendre integration: the function is evaluated exactly 527 // ten Times at interior points in the range o 538 // ten Times at interior points in the range of integration. Since the weights 528 // and abscissas are, in this case, symmetric << 539 // and abscissas are, in this case, symmetric around the midpoint of the 529 // range of integration, there are actually on 540 // range of integration, there are actually only five distinct values of each 530 // Convenient for using with some class object 541 // Convenient for using with some class object typeT 531 542 532 template <class T, class F> << 543 template <class T, class F> 533 G4double G4Integrator<T, F>::Legendre96(T& typ << 544 G4double G4Integrator<T,F>::Legendre96( T& typeT, F f,G4double a, G4double b) 534 { 545 { 535 G4int i; << 546 G4int i ; 536 G4double xDiff, xMean, dx, integral; << 547 G4double xDiff, xMean, dx, integral ; 537 << 548 538 // From Abramowitz M., Stegan I.A. 1964 , Ha << 549 // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919 539 << 550 540 static const G4double abscissa[] = { << 551 static const G4double 541 0.016276744849602969579, 0.048812985136049 << 552 abscissa[] = { 542 0.081297495464425558994, 0.113695850110665 << 553 0.016276744849602969579, 0.048812985136049731112, 543 0.145973714654896941989, 0.178096882367618 << 554 0.081297495464425558994, 0.113695850110665920911, 544 << 555 0.145973714654896941989, 0.178096882367618602759, // 6 545 0.210031310460567203603, 0.241743156163840 << 556 546 0.273198812591049141487, 0.304364944354496 << 557 0.210031310460567203603, 0.241743156163840012328, 547 0.335208522892625422616, 0.365696861472313 << 558 0.273198812591049141487, 0.304364944354496353024, 548 << 559 0.335208522892625422616, 0.365696861472313635031, // 12 549 0.395797649828908603285, 0.425478988407300 << 560 550 0.454709422167743008636, 0.483457973920596 << 561 0.395797649828908603285, 0.425478988407300545365, 551 0.511694177154667673586, 0.539388108324357 << 562 0.454709422167743008636, 0.483457973920596359768, 552 << 563 0.511694177154667673586, 0.539388108324357436227, // 18 553 0.566510418561397168404, 0.593032364777572 << 564 554 0.618925840125468570386, 0.644163403784967 << 565 0.566510418561397168404, 0.593032364777572080684, 555 0.668718310043916153953, 0.692564536642171 << 566 0.618925840125468570386, 0.644163403784967106798, 556 << 567 0.668718310043916153953, 0.692564536642171561344, // 24 557 0.715676812348967626225, 0.738030643744400 << 568 558 0.759602341176647498703, 0.780369043867433 << 569 0.715676812348967626225, 0.738030643744400132851, 559 0.800308744139140817229, 0.819400310737931 << 570 0.759602341176647498703, 0.780369043867433217604, 560 << 571 0.800308744139140817229, 0.819400310737931675539, // 30 561 0.837623511228187121494, 0.854959033434601 << 572 562 0.871388505909296502874, 0.886894517402420 << 573 0.837623511228187121494, 0.854959033434601455463, 563 0.901460635315852341319, 0.915071423120898 << 574 0.871388505909296502874, 0.886894517402420416057, 564 << 575 0.901460635315852341319, 0.915071423120898074206, // 36 565 0.927712456722308690965, 0.939370339752755 << 576 566 0.950032717784437635756, 0.959688291448742 << 577 0.927712456722308690965, 0.939370339752755216932, 567 0.968326828463264212174, 0.975939174585136 << 578 0.950032717784437635756, 0.959688291448742539300, 568 << 579 0.968326828463264212174, 0.975939174585136466453, // 42 569 0.982517263563014677447, 0.988054126329623 << 580 570 0.992543900323762624572, 0.995981842987209 << 581 0.982517263563014677447, 0.988054126329623799481, 571 0.998364375863181677724, 0.999689503883230 << 582 0.992543900323762624572, 0.995981842987209290650, 572 }; << 583 0.998364375863181677724, 0.999689503883230766828 // 48 573 << 584 } ; 574 static const G4double weight[] = { << 585 575 0.032550614492363166242, 0.032516118713868 << 586 static const G4double 576 0.032447163714064269364, 0.032343822568575 << 587 weight[] = { 577 0.032206204794030250669, 0.032034456231992 << 588 0.032550614492363166242, 0.032516118713868835987, 578 << 589 0.032447163714064269364, 0.032343822568575928429, 579 0.031828758894411006535, 0.031589330770727 << 590 0.032206204794030250669, 0.032034456231992663218, // 6 580 0.031316425596862355813, 0.031010332586313 << 591 581 0.030671376123669149014, 0.030299915420827 << 592 0.031828758894411006535, 0.031589330770727168558, 582 << 593 0.031316425596862355813, 0.031010332586313837423, 583 0.029896344136328385984, 0.029461089958167 << 594 0.030671376123669149014, 0.030299915420827593794, // 12 584 0.028994614150555236543, 0.028497411065085 << 595 585 0.027970007616848334440, 0.027412962726029 << 596 0.029896344136328385984, 0.029461089958167905970, 586 << 597 0.028994614150555236543, 0.028497411065085385646, 587 0.026826866725591762198, 0.026212340735672 << 598 0.027970007616848334440, 0.027412962726029242823, // 18 588 0.025570036005349361499, 0.024900633222483 << 599 589 0.024204841792364691282, 0.023483399085926 << 600 0.026826866725591762198, 0.026212340735672413913, 590 << 601 0.025570036005349361499, 0.024900633222483610288, 591 0.022737069658329374001, 0.021966644438744 << 602 0.024204841792364691282, 0.023483399085926219842, // 24 592 0.021172939892191298988, 0.020356797154333 << 603 593 0.019519081140145022410, 0.018660679627411 << 604 0.022737069658329374001, 0.021966644438744349195, 594 << 605 0.021172939892191298988, 0.020356797154333324595, 595 0.017782502316045260838, 0.016885479864245 << 606 0.019519081140145022410, 0.018660679627411467385, // 30 596 0.015970562902562291381, 0.015038721026994 << 607 597 0.014090941772314860916, 0.013128229566961 << 608 0.017782502316045260838, 0.016885479864245172450, 598 << 609 0.015970562902562291381, 0.015038721026994938006, 599 0.012151604671088319635, 0.011162102099838 << 610 0.014090941772314860916, 0.013128229566961572637, // 36 600 0.010160770535008415758, 0.009148671230783 << 611 601 0.008126876925698759217, 0.007096470791153 << 612 0.012151604671088319635, 0.011162102099838498591, 602 << 613 0.010160770535008415758, 0.009148671230783386633, 603 0.006058545504235961683, 0.005014202742927 << 614 0.008126876925698759217, 0.007096470791153865269, // 42 604 0.003964554338444686674, 0.002910731817934 << 615 605 0.001853960788946921732, 0.000796792065552 << 616 0.006058545504235961683, 0.005014202742927517693, 606 }; << 617 0.003964554338444686674, 0.002910731817934946408, 607 xMean = 0.5 * (a + b); << 618 0.001853960788946921732, 0.000796792065552012429 // 48 608 xDiff = 0.5 * (b - a); << 619 } ; 609 integral = 0.0; << 620 xMean = 0.5*(a + b) ; 610 for(i = 0; i < 48; ++i) << 621 xDiff = 0.5*(b - a) ; 611 { << 622 integral = 0.0 ; 612 dx = xDiff * abscissa[i]; << 623 for(i=0;i<48;i++) 613 integral += weight[i] * ((typeT.*f)(xMean << 624 { 614 } << 625 dx = xDiff*abscissa[i] ; 615 return integral *= xDiff; << 626 integral += weight[i]*((typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx)) ; >> 627 } >> 628 return integral *= xDiff ; 616 } 629 } 617 630 618 ////////////////////////////////////////////// 631 /////////////////////////////////////////////////////////////////////// 619 // 632 // 620 // Convenient for using with the pointer 'this 633 // Convenient for using with the pointer 'this' 621 634 622 template <class T, class F> << 635 template <class T, class F> 623 G4double G4Integrator<T, F>::Legendre96(T* ptr << 636 G4double G4Integrator<T,F>::Legendre96( T* ptrT, F f,G4double a, G4double b) 624 { 637 { 625 return Legendre96(*ptrT, f, a, b); << 638 return Legendre96(*ptrT,f,a,b) ; 626 } << 639 } 627 640 628 ////////////////////////////////////////////// 641 /////////////////////////////////////////////////////////////////////// 629 // 642 // 630 // Convenient for using with global scope func << 643 // Convenient for using with global scope function f 631 644 632 template <class T, class F> 645 template <class T, class F> 633 G4double G4Integrator<T, F>::Legendre96(G4doub << 646 G4double G4Integrator<T,F>::Legendre96( G4double (*f)(G4double), 634 G4doub << 647 G4double a, G4double b ) 635 { 648 { 636 G4int i; << 649 G4int i ; 637 G4double xDiff, xMean, dx, integral; << 650 G4double xDiff, xMean, dx, integral ; 638 << 651 639 // From Abramowitz M., Stegan I.A. 1964 , Ha << 652 // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919 640 << 653 641 static const G4double abscissa[] = { << 654 static const G4double 642 0.016276744849602969579, 0.048812985136049 << 655 abscissa[] = { 643 0.081297495464425558994, 0.113695850110665 << 656 0.016276744849602969579, 0.048812985136049731112, 644 0.145973714654896941989, 0.178096882367618 << 657 0.081297495464425558994, 0.113695850110665920911, 645 << 658 0.145973714654896941989, 0.178096882367618602759, // 6 646 0.210031310460567203603, 0.241743156163840 << 659 647 0.273198812591049141487, 0.304364944354496 << 660 0.210031310460567203603, 0.241743156163840012328, 648 0.335208522892625422616, 0.365696861472313 << 661 0.273198812591049141487, 0.304364944354496353024, 649 << 662 0.335208522892625422616, 0.365696861472313635031, // 12 650 0.395797649828908603285, 0.425478988407300 << 663 651 0.454709422167743008636, 0.483457973920596 << 664 0.395797649828908603285, 0.425478988407300545365, 652 0.511694177154667673586, 0.539388108324357 << 665 0.454709422167743008636, 0.483457973920596359768, 653 << 666 0.511694177154667673586, 0.539388108324357436227, // 18 654 0.566510418561397168404, 0.593032364777572 << 667 655 0.618925840125468570386, 0.644163403784967 << 668 0.566510418561397168404, 0.593032364777572080684, 656 0.668718310043916153953, 0.692564536642171 << 669 0.618925840125468570386, 0.644163403784967106798, 657 << 670 0.668718310043916153953, 0.692564536642171561344, // 24 658 0.715676812348967626225, 0.738030643744400 << 671 659 0.759602341176647498703, 0.780369043867433 << 672 0.715676812348967626225, 0.738030643744400132851, 660 0.800308744139140817229, 0.819400310737931 << 673 0.759602341176647498703, 0.780369043867433217604, 661 << 674 0.800308744139140817229, 0.819400310737931675539, // 30 662 0.837623511228187121494, 0.854959033434601 << 675 663 0.871388505909296502874, 0.886894517402420 << 676 0.837623511228187121494, 0.854959033434601455463, 664 0.901460635315852341319, 0.915071423120898 << 677 0.871388505909296502874, 0.886894517402420416057, 665 << 678 0.901460635315852341319, 0.915071423120898074206, // 36 666 0.927712456722308690965, 0.939370339752755 << 679 667 0.950032717784437635756, 0.959688291448742 << 680 0.927712456722308690965, 0.939370339752755216932, 668 0.968326828463264212174, 0.975939174585136 << 681 0.950032717784437635756, 0.959688291448742539300, 669 << 682 0.968326828463264212174, 0.975939174585136466453, // 42 670 0.982517263563014677447, 0.988054126329623 << 683 671 0.992543900323762624572, 0.995981842987209 << 684 0.982517263563014677447, 0.988054126329623799481, 672 0.998364375863181677724, 0.999689503883230 << 685 0.992543900323762624572, 0.995981842987209290650, 673 }; << 686 0.998364375863181677724, 0.999689503883230766828 // 48 674 << 687 } ; 675 static const G4double weight[] = { << 688 676 0.032550614492363166242, 0.032516118713868 << 689 static const G4double 677 0.032447163714064269364, 0.032343822568575 << 690 weight[] = { 678 0.032206204794030250669, 0.032034456231992 << 691 0.032550614492363166242, 0.032516118713868835987, 679 << 692 0.032447163714064269364, 0.032343822568575928429, 680 0.031828758894411006535, 0.031589330770727 << 693 0.032206204794030250669, 0.032034456231992663218, // 6 681 0.031316425596862355813, 0.031010332586313 << 694 682 0.030671376123669149014, 0.030299915420827 << 695 0.031828758894411006535, 0.031589330770727168558, 683 << 696 0.031316425596862355813, 0.031010332586313837423, 684 0.029896344136328385984, 0.029461089958167 << 697 0.030671376123669149014, 0.030299915420827593794, // 12 685 0.028994614150555236543, 0.028497411065085 << 698 686 0.027970007616848334440, 0.027412962726029 << 699 0.029896344136328385984, 0.029461089958167905970, 687 << 700 0.028994614150555236543, 0.028497411065085385646, 688 0.026826866725591762198, 0.026212340735672 << 701 0.027970007616848334440, 0.027412962726029242823, // 18 689 0.025570036005349361499, 0.024900633222483 << 702 690 0.024204841792364691282, 0.023483399085926 << 703 0.026826866725591762198, 0.026212340735672413913, 691 << 704 0.025570036005349361499, 0.024900633222483610288, 692 0.022737069658329374001, 0.021966644438744 << 705 0.024204841792364691282, 0.023483399085926219842, // 24 693 0.021172939892191298988, 0.020356797154333 << 706 694 0.019519081140145022410, 0.018660679627411 << 707 0.022737069658329374001, 0.021966644438744349195, 695 << 708 0.021172939892191298988, 0.020356797154333324595, 696 0.017782502316045260838, 0.016885479864245 << 709 0.019519081140145022410, 0.018660679627411467385, // 30 697 0.015970562902562291381, 0.015038721026994 << 710 698 0.014090941772314860916, 0.013128229566961 << 711 0.017782502316045260838, 0.016885479864245172450, 699 << 712 0.015970562902562291381, 0.015038721026994938006, 700 0.012151604671088319635, 0.011162102099838 << 713 0.014090941772314860916, 0.013128229566961572637, // 36 701 0.010160770535008415758, 0.009148671230783 << 714 702 0.008126876925698759217, 0.007096470791153 << 715 0.012151604671088319635, 0.011162102099838498591, 703 << 716 0.010160770535008415758, 0.009148671230783386633, 704 0.006058545504235961683, 0.005014202742927 << 717 0.008126876925698759217, 0.007096470791153865269, // 42 705 0.003964554338444686674, 0.002910731817934 << 718 706 0.001853960788946921732, 0.000796792065552 << 719 0.006058545504235961683, 0.005014202742927517693, 707 }; << 720 0.003964554338444686674, 0.002910731817934946408, 708 xMean = 0.5 * (a + b); << 721 0.001853960788946921732, 0.000796792065552012429 // 48 709 xDiff = 0.5 * (b - a); << 722 } ; 710 integral = 0.0; << 723 xMean = 0.5*(a + b) ; 711 for(i = 0; i < 48; ++i) << 724 xDiff = 0.5*(b - a) ; 712 { << 725 integral = 0.0 ; 713 dx = xDiff * abscissa[i]; << 726 for(i=0;i<48;i++) 714 integral += weight[i] * ((*f)(xMean + dx) << 727 { 715 } << 728 dx = xDiff*abscissa[i] ; 716 return integral *= xDiff; << 729 integral += weight[i]*((*f)(xMean + dx) + (*f)(xMean - dx)) ; >> 730 } >> 731 return integral *= xDiff ; 717 } 732 } 718 733 719 ////////////////////////////////////////////// 734 ////////////////////////////////////////////////////////////////////////////// 720 // 735 // 721 // Methods involving Chebyshev polynomials << 736 // Methods involving Chebyshev polynomials 722 // 737 // 723 ////////////////////////////////////////////// 738 /////////////////////////////////////////////////////////////////////////// 724 // 739 // 725 // Integrates function pointed by T::f from a << 740 // Integrates function pointed by T::f from a to b by Gauss-Chebyshev 726 // quadrature method. 741 // quadrature method. 727 // Convenient for using with class object type 742 // Convenient for using with class object typeT 728 743 729 template <class T, class F> 744 template <class T, class F> 730 G4double G4Integrator<T, F>::Chebyshev(T& type << 745 G4double G4Integrator<T,F>::Chebyshev( T& typeT, F f, G4double a, 731 G4int n << 746 G4double b, G4int nChebyshev ) 732 { 747 { 733 G4int i; << 748 G4int i ; 734 G4double xDiff, xMean, dx, integral = 0.0; << 749 G4double xDiff, xMean, dx, integral = 0.0 ; 735 << 750 736 G4int fNumber = nChebyshev; // Try to << 751 G4int fNumber = nChebyshev ; // Try to reduce fNumber twice ?? 737 G4double cof = CLHEP::pi / fNumber; << 752 G4double cof = CLHEP::pi/fNumber ; 738 G4double* fAbscissa = new G4double[fNumber]; << 753 G4double* fAbscissa = new G4double[fNumber] ; 739 G4double* fWeight = new G4double[fNumber]; << 754 G4double* fWeight = new G4double[fNumber] ; 740 for(i = 0; i < fNumber; ++i) << 755 for(i=0;i<fNumber;i++) 741 { << 756 { 742 fAbscissa[i] = std::cos(cof * (i + 0.5)); << 757 fAbscissa[i] = std::cos(cof*(i + 0.5)) ; 743 fWeight[i] = cof * std::sqrt(1 - fAbscis << 758 fWeight[i] = cof*std::sqrt(1 - fAbscissa[i]*fAbscissa[i]) ; 744 } << 759 } 745 << 760 746 // << 761 // 747 // Now we ready to estimate the integral << 762 // Now we ready to estimate the integral 748 // << 763 // 749 << 764 750 xMean = 0.5 * (a + b); << 765 xMean = 0.5*(a + b) ; 751 xDiff = 0.5 * (b - a); << 766 xDiff = 0.5*(b - a) ; 752 for(i = 0; i < fNumber; ++i) << 767 for(i=0;i<fNumber;i++) 753 { << 768 { 754 dx = xDiff * fAbscissa[i]; << 769 dx = xDiff*fAbscissa[i] ; 755 integral += fWeight[i] * (typeT.*f)(xMean << 770 integral += fWeight[i]*(typeT.*f)(xMean + dx) ; 756 } << 771 } 757 delete[] fAbscissa; << 772 delete[] fAbscissa; 758 delete[] fWeight; << 773 delete[] fWeight; 759 return integral *= xDiff; << 774 return integral *= xDiff ; 760 } 775 } 761 776 762 ////////////////////////////////////////////// 777 /////////////////////////////////////////////////////////////////////// 763 // 778 // 764 // Convenient for using with 'this' pointer 779 // Convenient for using with 'this' pointer 765 780 766 template <class T, class F> 781 template <class T, class F> 767 G4double G4Integrator<T, F>::Chebyshev(T* ptrT << 782 G4double G4Integrator<T,F>::Chebyshev( T* ptrT, F f, G4double a, 768 G4int n << 783 G4double b, G4int n ) 769 { 784 { 770 return Chebyshev(*ptrT, f, a, b, n); << 785 return Chebyshev(*ptrT,f,a,b,n) ; 771 } << 786 } 772 787 773 ////////////////////////////////////////////// 788 //////////////////////////////////////////////////////////////////////// 774 // 789 // 775 // For use with global scope functions f << 790 // For use with global scope functions f 776 791 777 template <class T, class F> 792 template <class T, class F> 778 G4double G4Integrator<T, F>::Chebyshev(G4doubl << 793 G4double G4Integrator<T,F>::Chebyshev( G4double (*f)(G4double), 779 G4doubl << 794 G4double a, G4double b, G4int nChebyshev ) 780 { 795 { 781 G4int i; << 796 G4int i ; 782 G4double xDiff, xMean, dx, integral = 0.0; << 797 G4double xDiff, xMean, dx, integral = 0.0 ; 783 << 798 784 G4int fNumber = nChebyshev; // Try to << 799 G4int fNumber = nChebyshev ; // Try to reduce fNumber twice ?? 785 G4double cof = CLHEP::pi / fNumber; << 800 G4double cof = CLHEP::pi/fNumber ; 786 G4double* fAbscissa = new G4double[fNumber]; << 801 G4double* fAbscissa = new G4double[fNumber] ; 787 G4double* fWeight = new G4double[fNumber]; << 802 G4double* fWeight = new G4double[fNumber] ; 788 for(i = 0; i < fNumber; ++i) << 803 for(i=0;i<fNumber;i++) 789 { << 804 { 790 fAbscissa[i] = std::cos(cof * (i + 0.5)); << 805 fAbscissa[i] = std::cos(cof*(i + 0.5)) ; 791 fWeight[i] = cof * std::sqrt(1 - fAbscis << 806 fWeight[i] = cof*std::sqrt(1 - fAbscissa[i]*fAbscissa[i]) ; 792 } << 807 } 793 << 808 794 // << 809 // 795 // Now we ready to estimate the integral << 810 // Now we ready to estimate the integral 796 // << 811 // 797 << 812 798 xMean = 0.5 * (a + b); << 813 xMean = 0.5*(a + b) ; 799 xDiff = 0.5 * (b - a); << 814 xDiff = 0.5*(b - a) ; 800 for(i = 0; i < fNumber; ++i) << 815 for(i=0;i<fNumber;i++) 801 { << 816 { 802 dx = xDiff * fAbscissa[i]; << 817 dx = xDiff*fAbscissa[i] ; 803 integral += fWeight[i] * (*f)(xMean + dx); << 818 integral += fWeight[i]*(*f)(xMean + dx) ; 804 } << 819 } 805 delete[] fAbscissa; << 820 delete[] fAbscissa; 806 delete[] fWeight; << 821 delete[] fWeight; 807 return integral *= xDiff; << 822 return integral *= xDiff ; 808 } 823 } 809 824 810 ////////////////////////////////////////////// 825 ////////////////////////////////////////////////////////////////////// 811 // 826 // 812 // Method involving Laguerre polynomials 827 // Method involving Laguerre polynomials 813 // 828 // 814 ////////////////////////////////////////////// 829 ////////////////////////////////////////////////////////////////////// 815 // 830 // 816 // Integral from zero to infinity of std::pow( << 831 // Integral from zero to infinity of std::pow(x,alpha)*std::exp(-x)*f(x). 817 // The value of nLaguerre sets the accuracy. 832 // The value of nLaguerre sets the accuracy. 818 // The function creates arrays fAbscissa[0,.., << 833 // The function creates arrays fAbscissa[0,..,nLaguerre-1] and 819 // fWeight[0,..,nLaguerre-1] . << 834 // fWeight[0,..,nLaguerre-1] . 820 // Convenient for using with class object 'typ 835 // Convenient for using with class object 'typeT' and (typeT.*f) function 821 // (T::f) 836 // (T::f) 822 837 823 template <class T, class F> 838 template <class T, class F> 824 G4double G4Integrator<T, F>::Laguerre(T& typeT << 839 G4double G4Integrator<T,F>::Laguerre( T& typeT, F f, G4double alpha, 825 G4int nL << 840 G4int nLaguerre ) 826 { 841 { 827 const G4double tolerance = 1.0e-10; << 842 const G4double tolerance = 1.0e-10 ; 828 const G4int maxNumber = 12; << 843 const G4int maxNumber = 12 ; 829 G4int i, j, k; << 844 G4int i, j, k ; 830 G4double nwt = 0., nwt1, temp1, temp2, << 845 G4double nwt=0., nwt1, temp1, temp2, temp3, temp, cofi ; 831 G4double integral = 0.0; << 846 G4double integral = 0.0 ; 832 << 847 833 G4int fNumber = nLaguerre; << 848 G4int fNumber = nLaguerre ; 834 G4double* fAbscissa = new G4double[fNumber]; << 849 G4double* fAbscissa = new G4double[fNumber] ; 835 G4double* fWeight = new G4double[fNumber]; << 850 G4double* fWeight = new G4double[fNumber] ; >> 851 >> 852 for(i=1;i<=fNumber;i++) // Loop over the desired roots >> 853 { >> 854 if(i == 1) >> 855 { >> 856 nwt = (1.0 + alpha)*(3.0 + 0.92*alpha) >> 857 / (1.0 + 2.4*fNumber + 1.8*alpha) ; >> 858 } >> 859 else if(i == 2) >> 860 { >> 861 nwt += (15.0 + 6.25*alpha)/(1.0 + 0.9*alpha + 2.5*fNumber) ; >> 862 } >> 863 else >> 864 { >> 865 cofi = i - 2 ; >> 866 nwt += ((1.0+2.55*cofi)/(1.9*cofi) >> 867 + 1.26*cofi*alpha/(1.0+3.5*cofi)) >> 868 * (nwt - fAbscissa[i-3])/(1.0 + 0.3*alpha) ; >> 869 } >> 870 for(k=1;k<=maxNumber;k++) >> 871 { >> 872 temp1 = 1.0 ; >> 873 temp2 = 0.0 ; 836 874 837 for(i = 1; i <= fNumber; ++i) // Loop over << 875 for(j=1;j<=fNumber;j++) 838 { << 876 { 839 if(i == 1) << 877 temp3 = temp2 ; 840 { << 878 temp2 = temp1 ; 841 nwt = (1.0 + alpha) * (3.0 + 0.92 * alph << 879 temp1 = ((2*j - 1 + alpha - nwt)*temp2 - (j - 1 + alpha)*temp3)/j ; 842 (1.0 + 2.4 * fNumber + 1.8 * alpha << 880 } 843 } << 881 temp = (fNumber*temp1 - (fNumber +alpha)*temp2)/nwt ; 844 else if(i == 2) << 882 nwt1 = nwt ; 845 { << 883 nwt = nwt1 - temp1/temp ; 846 nwt += (15.0 + 6.25 * alpha) / (1.0 + 0. << 884 847 } << 885 if(std::fabs(nwt - nwt1) <= tolerance) 848 else << 886 { 849 { << 887 break ; 850 cofi = i - 2; << 888 } 851 nwt += ((1.0 + 2.55 * cofi) / (1.9 * cof << 889 } 852 1.26 * cofi * alpha / (1.0 + 3.5 << 890 if(k > maxNumber) 853 (nwt - fAbscissa[i - 3]) / (1.0 + << 891 { 854 } << 892 G4Exception("G4Integrator<T,F>::Laguerre(T,F, ...)", "Error", 855 for(k = 1; k <= maxNumber; ++k) << 893 FatalException, "Too many (>12) iterations."); 856 { << 894 } 857 temp1 = 1.0; << 895 858 temp2 = 0.0; << 896 fAbscissa[i-1] = nwt ; 859 << 897 fWeight[i-1] = -std::exp(GammaLogarithm(alpha + fNumber) - 860 for(j = 1; j <= fNumber; ++j) << 898 GammaLogarithm((G4double)fNumber))/(temp*fNumber*temp2) ; 861 { << 899 } 862 temp3 = temp2; << 900 863 temp2 = temp1; << 901 // 864 temp1 = << 902 // Integral evaluation 865 ((2 * j - 1 + alpha - nwt) * temp2 - << 903 // 866 } << 904 867 temp = (fNumber * temp1 - (fNumber + alp << 905 for(i=0;i<fNumber;i++) 868 nwt1 = nwt; << 906 { 869 nwt = nwt1 - temp1 / temp; << 907 integral += fWeight[i]*(typeT.*f)(fAbscissa[i]) ; 870 << 908 } 871 if(std::fabs(nwt - nwt1) <= tolerance) << 909 delete[] fAbscissa; 872 { << 910 delete[] fWeight; 873 break; << 911 return integral ; 874 } << 912 } 875 } << 876 if(k > maxNumber) << 877 { << 878 G4Exception("G4Integrator<T,F>::Laguerre << 879 FatalException, "Too many (> << 880 } << 881 << 882 fAbscissa[i - 1] = nwt; << 883 fWeight[i - 1] = -std::exp(GammaLogarith << 884 GammaLogarithm( << 885 (temp * fNumber * temp2); << 886 } << 887 913 888 // << 889 // Integral evaluation << 890 // << 891 914 892 for(i = 0; i < fNumber; ++i) << 893 { << 894 integral += fWeight[i] * (typeT.*f)(fAbsci << 895 } << 896 delete[] fAbscissa; << 897 delete[] fWeight; << 898 return integral; << 899 } << 900 915 901 ////////////////////////////////////////////// 916 ////////////////////////////////////////////////////////////////////// 902 // 917 // 903 // 918 // 904 919 905 template <class T, class F> << 920 template <class T, class F> G4double 906 G4double G4Integrator<T, F>::Laguerre(T* ptrT, << 921 G4Integrator<T,F>::Laguerre( T* ptrT, F f, G4double alpha, G4int nLaguerre ) 907 G4int nL << 908 { 922 { 909 return Laguerre(*ptrT, f, alpha, nLaguerre); << 923 return Laguerre(*ptrT,f,alpha,nLaguerre) ; 910 } 924 } 911 925 912 ////////////////////////////////////////////// 926 //////////////////////////////////////////////////////////////////////// 913 // 927 // 914 // For use with global scope functions f << 928 // For use with global scope functions f 915 << 916 template <class T, class F> << 917 G4double G4Integrator<T, F>::Laguerre(G4double << 918 G4int nL << 919 { << 920 const G4double tolerance = 1.0e-10; << 921 const G4int maxNumber = 12; << 922 G4int i, j, k; << 923 G4double nwt = 0., nwt1, temp1, temp2, << 924 G4double integral = 0.0; << 925 << 926 G4int fNumber = nLaguerre; << 927 G4double* fAbscissa = new G4double[fNumber]; << 928 G4double* fWeight = new G4double[fNumber]; << 929 << 930 for(i = 1; i <= fNumber; ++i) // Loop over << 931 { << 932 if(i == 1) << 933 { << 934 nwt = (1.0 + alpha) * (3.0 + 0.92 * alph << 935 (1.0 + 2.4 * fNumber + 1.8 * alpha << 936 } << 937 else if(i == 2) << 938 { << 939 nwt += (15.0 + 6.25 * alpha) / (1.0 + 0. << 940 } << 941 else << 942 { << 943 cofi = i - 2; << 944 nwt += ((1.0 + 2.55 * cofi) / (1.9 * cof << 945 1.26 * cofi * alpha / (1.0 + 3.5 << 946 (nwt - fAbscissa[i - 3]) / (1.0 + << 947 } << 948 for(k = 1; k <= maxNumber; ++k) << 949 { << 950 temp1 = 1.0; << 951 temp2 = 0.0; << 952 << 953 for(j = 1; j <= fNumber; ++j) << 954 { << 955 temp3 = temp2; << 956 temp2 = temp1; << 957 temp1 = << 958 ((2 * j - 1 + alpha - nwt) * temp2 - << 959 } << 960 temp = (fNumber * temp1 - (fNumber + alp << 961 nwt1 = nwt; << 962 nwt = nwt1 - temp1 / temp; << 963 << 964 if(std::fabs(nwt - nwt1) <= tolerance) << 965 { << 966 break; << 967 } << 968 } << 969 if(k > maxNumber) << 970 { << 971 G4Exception("G4Integrator<T,F>::Laguerre << 972 "Too many (>12) iterations." << 973 } << 974 << 975 fAbscissa[i - 1] = nwt; << 976 fWeight[i - 1] = -std::exp(GammaLogarith << 977 GammaLogarithm( << 978 (temp * fNumber * temp2); << 979 } << 980 929 981 // << 930 template <class T, class F> G4double 982 // Integral evaluation << 931 G4Integrator<T,F>::Laguerre( G4double (*f)(G4double), 983 // << 932 G4double alpha, G4int nLaguerre ) >> 933 { >> 934 const G4double tolerance = 1.0e-10 ; >> 935 const G4int maxNumber = 12 ; >> 936 G4int i, j, k ; >> 937 G4double nwt=0., nwt1, temp1, temp2, temp3, temp, cofi ; >> 938 G4double integral = 0.0 ; >> 939 >> 940 G4int fNumber = nLaguerre ; >> 941 G4double* fAbscissa = new G4double[fNumber] ; >> 942 G4double* fWeight = new G4double[fNumber] ; >> 943 >> 944 for(i=1;i<=fNumber;i++) // Loop over the desired roots >> 945 { >> 946 if(i == 1) >> 947 { >> 948 nwt = (1.0 + alpha)*(3.0 + 0.92*alpha) >> 949 / (1.0 + 2.4*fNumber + 1.8*alpha) ; >> 950 } >> 951 else if(i == 2) >> 952 { >> 953 nwt += (15.0 + 6.25*alpha)/(1.0 + 0.9*alpha + 2.5*fNumber) ; >> 954 } >> 955 else >> 956 { >> 957 cofi = i - 2 ; >> 958 nwt += ((1.0+2.55*cofi)/(1.9*cofi) >> 959 + 1.26*cofi*alpha/(1.0+3.5*cofi)) >> 960 * (nwt - fAbscissa[i-3])/(1.0 + 0.3*alpha) ; >> 961 } >> 962 for(k=1;k<=maxNumber;k++) >> 963 { >> 964 temp1 = 1.0 ; >> 965 temp2 = 0.0 ; 984 966 985 for(i = 0; i < fNumber; i++) << 967 for(j=1;j<=fNumber;j++) 986 { << 968 { 987 integral += fWeight[i] * (*f)(fAbscissa[i] << 969 temp3 = temp2 ; 988 } << 970 temp2 = temp1 ; 989 delete[] fAbscissa; << 971 temp1 = ((2*j - 1 + alpha - nwt)*temp2 - (j - 1 + alpha)*temp3)/j ; 990 delete[] fWeight; << 972 } 991 return integral; << 973 temp = (fNumber*temp1 - (fNumber +alpha)*temp2)/nwt ; >> 974 nwt1 = nwt ; >> 975 nwt = nwt1 - temp1/temp ; >> 976 >> 977 if(std::fabs(nwt - nwt1) <= tolerance) >> 978 { >> 979 break ; >> 980 } >> 981 } >> 982 if(k > maxNumber) >> 983 { >> 984 G4Exception("G4Integrator<T,F>::Laguerre( ...)", "Error", >> 985 FatalException, "Too many (>12) iterations."); >> 986 } >> 987 >> 988 fAbscissa[i-1] = nwt ; >> 989 fWeight[i-1] = -std::exp(GammaLogarithm(alpha + fNumber) - >> 990 GammaLogarithm((G4double)fNumber))/(temp*fNumber*temp2) ; >> 991 } >> 992 >> 993 // >> 994 // Integral evaluation >> 995 // >> 996 >> 997 for(i=0;i<fNumber;i++) >> 998 { >> 999 integral += fWeight[i]*(*f)(fAbscissa[i]) ; >> 1000 } >> 1001 delete[] fAbscissa; >> 1002 delete[] fWeight; >> 1003 return integral ; 992 } 1004 } 993 1005 994 ////////////////////////////////////////////// 1006 /////////////////////////////////////////////////////////////////////// 995 // 1007 // 996 // Auxiliary function which returns the value 1008 // Auxiliary function which returns the value of std::log(gamma-function(x)) 997 // Returns the value ln(Gamma(xx) for xx > 0. << 1009 // Returns the value ln(Gamma(xx) for xx > 0. Full accuracy is obtained for 998 // xx > 1. For 0 < xx < 1. the reflection form 1010 // xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first. 999 // (Adapted from Numerical Recipes in C) 1011 // (Adapted from Numerical Recipes in C) 1000 // 1012 // 1001 1013 1002 template <class T, class F> 1014 template <class T, class F> 1003 G4double G4Integrator<T, F>::GammaLogarithm(G << 1015 G4double G4Integrator<T,F>::GammaLogarithm(G4double xx) 1004 { 1016 { 1005 static const G4double cof[6] = { 76.1800917 1017 static const G4double cof[6] = { 76.18009172947146, -86.50532032941677, 1006 24.0140982 << 1018 24.01409824083091, -1.231739572450155, 1007 0.12086509 1019 0.1208650973866179e-2, -0.5395239384953e-5 }; 1008 G4int j; 1020 G4int j; 1009 G4double x = xx - 1.0; << 1021 G4double x = xx - 1.0 ; 1010 G4double tmp = x + 5.5; << 1022 G4double tmp = x + 5.5 ; 1011 tmp -= (x + 0.5) * std::log(tmp); << 1023 tmp -= (x + 0.5) * std::log(tmp) ; 1012 G4double ser = 1.000000000190015; << 1024 G4double ser = 1.000000000190015 ; 1013 1025 1014 for(j = 0; j <= 5; ++j) << 1026 for ( j = 0; j <= 5; j++ ) 1015 { 1027 { 1016 x += 1.0; << 1028 x += 1.0 ; 1017 ser += cof[j] / x; << 1029 ser += cof[j]/x ; 1018 } 1030 } 1019 return -tmp + std::log(2.5066282746310005 * << 1031 return -tmp + std::log(2.5066282746310005*ser) ; 1020 } 1032 } 1021 1033 1022 ///////////////////////////////////////////// 1034 /////////////////////////////////////////////////////////////////////// 1023 // 1035 // 1024 // Method involving Hermite polynomials 1036 // Method involving Hermite polynomials 1025 // 1037 // 1026 ///////////////////////////////////////////// 1038 /////////////////////////////////////////////////////////////////////// 1027 // 1039 // 1028 // 1040 // 1029 // Gauss-Hermite method for integration of st << 1041 // Gauss-Hermite method for integration of std::exp(-x*x)*f(x) 1030 // from minus infinity to plus infinity . << 1042 // from minus infinity to plus infinity . 1031 // 1043 // 1032 1044 1033 template <class T, class F> << 1045 template <class T, class F> 1034 G4double G4Integrator<T, F>::Hermite(T& typeT << 1046 G4double G4Integrator<T,F>::Hermite( T& typeT, F f, G4int nHermite ) 1035 { 1047 { 1036 const G4double tolerance = 1.0e-12; << 1048 const G4double tolerance = 1.0e-12 ; 1037 const G4int maxNumber = 12; << 1049 const G4int maxNumber = 12 ; 1038 << 1050 1039 G4int i, j, k; << 1051 G4int i, j, k ; 1040 G4double integral = 0.0; << 1052 G4double integral = 0.0 ; 1041 G4double nwt = 0., nwt1, temp1, temp2, << 1053 G4double nwt=0., nwt1, temp1, temp2, temp3, temp ; 1042 1054 1043 G4double piInMinusQ = << 1055 G4double piInMinusQ = std::pow(CLHEP::pi,-0.25) ; // 1.0/std::sqrt(std::sqrt(pi)) ?? 1044 std::pow(CLHEP::pi, -0.25); // 1.0/std:: << 1045 1056 1046 G4int fNumber = (nHermite + 1) / 2; << 1057 G4int fNumber = (nHermite +1)/2 ; 1047 G4double* fAbscissa = new G4double[fNumber] << 1058 G4double* fAbscissa = new G4double[fNumber] ; 1048 G4double* fWeight = new G4double[fNumber] << 1059 G4double* fWeight = new G4double[fNumber] ; 1049 1060 1050 for(i = 1; i <= fNumber; ++i) << 1061 for(i=1;i<=fNumber;i++) 1051 { << 1062 { 1052 if(i == 1) << 1063 if(i == 1) 1053 { << 1064 { 1054 nwt = std::sqrt((G4double)(2 * nHermite << 1065 nwt = std::sqrt((G4double)(2*nHermite + 1)) - 1055 1.85575001 * std::pow((G4double)( << 1066 1.85575001*std::pow((G4double)(2*nHermite + 1),-0.16666999) ; 1056 } << 1067 } 1057 else if(i == 2) << 1068 else if(i == 2) 1058 { << 1069 { 1059 nwt -= 1.14001 * std::pow((G4double) nH << 1070 nwt -= 1.14001*std::pow((G4double)nHermite,0.425999)/nwt ; 1060 } << 1071 } 1061 else if(i == 3) << 1072 else if(i == 3) 1062 { << 1073 { 1063 nwt = 1.86002 * nwt - 0.86002 * fAbscis << 1074 nwt = 1.86002*nwt - 0.86002*fAbscissa[0] ; 1064 } << 1075 } 1065 else if(i == 4) << 1076 else if(i == 4) 1066 { << 1077 { 1067 nwt = 1.91001 * nwt - 0.91001 * fAbscis << 1078 nwt = 1.91001*nwt - 0.91001*fAbscissa[1] ; 1068 } << 1079 } 1069 else << 1080 else 1070 { << 1081 { 1071 nwt = 2.0 * nwt - fAbscissa[i - 3]; << 1082 nwt = 2.0*nwt - fAbscissa[i - 3] ; 1072 } << 1083 } 1073 for(k = 1; k <= maxNumber; ++k) << 1084 for(k=1;k<=maxNumber;k++) 1074 { << 1085 { 1075 temp1 = piInMinusQ; << 1086 temp1 = piInMinusQ ; 1076 temp2 = 0.0; << 1087 temp2 = 0.0 ; 1077 << 1078 for(j = 1; j <= nHermite; ++j) << 1079 { << 1080 temp3 = temp2; << 1081 temp2 = temp1; << 1082 temp1 = nwt * std::sqrt(2.0 / j) * te << 1083 std::sqrt(((G4double)(j - 1)) << 1084 } << 1085 temp = std::sqrt((G4double) 2 * nHermit << 1086 nwt1 = nwt; << 1087 nwt = nwt1 - temp1 / temp; << 1088 << 1089 if(std::fabs(nwt - nwt1) <= tolerance) << 1090 { << 1091 break; << 1092 } << 1093 } << 1094 if(k > maxNumber) << 1095 { << 1096 G4Exception("G4Integrator<T,F>::Hermite << 1097 FatalException, "Too many ( << 1098 } << 1099 fAbscissa[i - 1] = nwt; << 1100 fWeight[i - 1] = 2.0 / (temp * temp); << 1101 } << 1102 << 1103 // << 1104 // Integral calculation << 1105 // << 1106 1088 1107 for(i = 0; i < fNumber; ++i) << 1089 for(j=1;j<=nHermite;j++) 1108 { << 1090 { 1109 integral += << 1091 temp3 = temp2 ; 1110 fWeight[i] * ((typeT.*f)(fAbscissa[i]) << 1092 temp2 = temp1 ; 1111 } << 1093 temp1 = nwt*std::sqrt(2.0/j)*temp2 - 1112 delete[] fAbscissa; << 1094 std::sqrt(((G4double)(j - 1))/j)*temp3 ; 1113 delete[] fWeight; << 1095 } 1114 return integral; << 1096 temp = std::sqrt((G4double)2*nHermite)*temp2 ; >> 1097 nwt1 = nwt ; >> 1098 nwt = nwt1 - temp1/temp ; >> 1099 >> 1100 if(std::fabs(nwt - nwt1) <= tolerance) >> 1101 { >> 1102 break ; >> 1103 } >> 1104 } >> 1105 if(k > maxNumber) >> 1106 { >> 1107 G4Exception("G4Integrator<T,F>::Hermite(T,F, ...)", "Error", >> 1108 FatalException, "Too many (>12) iterations."); >> 1109 } >> 1110 fAbscissa[i-1] = nwt ; >> 1111 fWeight[i-1] = 2.0/(temp*temp) ; >> 1112 } >> 1113 >> 1114 // >> 1115 // Integral calculation >> 1116 // >> 1117 >> 1118 for(i=0;i<fNumber;i++) >> 1119 { >> 1120 integral += fWeight[i]*( (typeT.*f)(fAbscissa[i]) + >> 1121 (typeT.*f)(-fAbscissa[i]) ) ; >> 1122 } >> 1123 delete[] fAbscissa; >> 1124 delete[] fWeight; >> 1125 return integral ; 1115 } 1126 } 1116 1127 >> 1128 1117 ///////////////////////////////////////////// 1129 //////////////////////////////////////////////////////////////////////// 1118 // 1130 // 1119 // For use with 'this' pointer 1131 // For use with 'this' pointer 1120 1132 1121 template <class T, class F> << 1133 template <class T, class F> 1122 G4double G4Integrator<T, F>::Hermite(T* ptrT, << 1134 G4double G4Integrator<T,F>::Hermite( T* ptrT, F f, G4int n ) 1123 { 1135 { 1124 return Hermite(*ptrT, f, n); << 1136 return Hermite(*ptrT,f,n) ; 1125 } << 1137 } 1126 1138 1127 ///////////////////////////////////////////// 1139 //////////////////////////////////////////////////////////////////////// 1128 // 1140 // 1129 // For use with global scope f 1141 // For use with global scope f 1130 1142 1131 template <class T, class F> 1143 template <class T, class F> 1132 G4double G4Integrator<T, F>::Hermite(G4double << 1144 G4double G4Integrator<T,F>::Hermite( G4double (*f)(G4double), G4int nHermite) 1133 { 1145 { 1134 const G4double tolerance = 1.0e-12; << 1146 const G4double tolerance = 1.0e-12 ; 1135 const G4int maxNumber = 12; << 1147 const G4int maxNumber = 12 ; 1136 << 1148 1137 G4int i, j, k; << 1149 G4int i, j, k ; 1138 G4double integral = 0.0; << 1150 G4double integral = 0.0 ; 1139 G4double nwt = 0., nwt1, temp1, temp2, << 1151 G4double nwt=0., nwt1, temp1, temp2, temp3, temp ; 1140 << 1152 1141 G4double piInMinusQ = << 1153 G4double piInMinusQ = std::pow(CLHEP::pi,-0.25) ; // 1.0/std::sqrt(std::sqrt(pi)) ?? 1142 std::pow(CLHEP::pi, -0.25); // 1.0/std:: << 1154 1143 << 1155 G4int fNumber = (nHermite +1)/2 ; 1144 G4int fNumber = (nHermite + 1) / 2; << 1156 G4double* fAbscissa = new G4double[fNumber] ; 1145 G4double* fAbscissa = new G4double[fNumber] << 1157 G4double* fWeight = new G4double[fNumber] ; 1146 G4double* fWeight = new G4double[fNumber] << 1158 1147 << 1159 for(i=1;i<=fNumber;i++) 1148 for(i = 1; i <= fNumber; ++i) << 1160 { 1149 { << 1161 if(i == 1) 1150 if(i == 1) << 1162 { 1151 { << 1163 nwt = std::sqrt((G4double)(2*nHermite + 1)) - 1152 nwt = std::sqrt((G4double)(2 * nHermite << 1164 1.85575001*std::pow((G4double)(2*nHermite + 1),-0.16666999) ; 1153 1.85575001 * std::pow((G4double)( << 1165 } 1154 } << 1166 else if(i == 2) 1155 else if(i == 2) << 1167 { 1156 { << 1168 nwt -= 1.14001*std::pow((G4double)nHermite,0.425999)/nwt ; 1157 nwt -= 1.14001 * std::pow((G4double) nH << 1169 } 1158 } << 1170 else if(i == 3) 1159 else if(i == 3) << 1171 { 1160 { << 1172 nwt = 1.86002*nwt - 0.86002*fAbscissa[0] ; 1161 nwt = 1.86002 * nwt - 0.86002 * fAbscis << 1173 } 1162 } << 1174 else if(i == 4) 1163 else if(i == 4) << 1175 { 1164 { << 1176 nwt = 1.91001*nwt - 0.91001*fAbscissa[1] ; 1165 nwt = 1.91001 * nwt - 0.91001 * fAbscis << 1177 } 1166 } << 1178 else 1167 else << 1179 { 1168 { << 1180 nwt = 2.0*nwt - fAbscissa[i - 3] ; 1169 nwt = 2.0 * nwt - fAbscissa[i - 3]; << 1181 } 1170 } << 1182 for(k=1;k<=maxNumber;k++) 1171 for(k = 1; k <= maxNumber; ++k) << 1183 { 1172 { << 1184 temp1 = piInMinusQ ; 1173 temp1 = piInMinusQ; << 1185 temp2 = 0.0 ; 1174 temp2 = 0.0; << 1175 << 1176 for(j = 1; j <= nHermite; ++j) << 1177 { << 1178 temp3 = temp2; << 1179 temp2 = temp1; << 1180 temp1 = nwt * std::sqrt(2.0 / j) * te << 1181 std::sqrt(((G4double)(j - 1)) << 1182 } << 1183 temp = std::sqrt((G4double) 2 * nHermit << 1184 nwt1 = nwt; << 1185 nwt = nwt1 - temp1 / temp; << 1186 << 1187 if(std::fabs(nwt - nwt1) <= tolerance) << 1188 { << 1189 break; << 1190 } << 1191 } << 1192 if(k > maxNumber) << 1193 { << 1194 G4Exception("G4Integrator<T,F>::Hermite << 1195 "Too many (>12) iterations. << 1196 } << 1197 fAbscissa[i - 1] = nwt; << 1198 fWeight[i - 1] = 2.0 / (temp * temp); << 1199 } << 1200 << 1201 // << 1202 // Integral calculation << 1203 // << 1204 1186 1205 for(i = 0; i < fNumber; ++i) << 1187 for(j=1;j<=nHermite;j++) 1206 { << 1188 { 1207 integral += fWeight[i] * ((*f)(fAbscissa[ << 1189 temp3 = temp2 ; 1208 } << 1190 temp2 = temp1 ; 1209 delete[] fAbscissa; << 1191 temp1 = nwt*std::sqrt(2.0/j)*temp2 - 1210 delete[] fWeight; << 1192 std::sqrt(((G4double)(j - 1))/j)*temp3 ; 1211 return integral; << 1193 } >> 1194 temp = std::sqrt((G4double)2*nHermite)*temp2 ; >> 1195 nwt1 = nwt ; >> 1196 nwt = nwt1 - temp1/temp ; >> 1197 >> 1198 if(std::fabs(nwt - nwt1) <= tolerance) >> 1199 { >> 1200 break ; >> 1201 } >> 1202 } >> 1203 if(k > maxNumber) >> 1204 { >> 1205 G4Exception("G4Integrator<T,F>::Hermite(...)", "Error", >> 1206 FatalException, "Too many (>12) iterations."); >> 1207 } >> 1208 fAbscissa[i-1] = nwt ; >> 1209 fWeight[i-1] = 2.0/(temp*temp) ; >> 1210 } >> 1211 >> 1212 // >> 1213 // Integral calculation >> 1214 // >> 1215 >> 1216 for(i=0;i<fNumber;i++) >> 1217 { >> 1218 integral += fWeight[i]*( (*f)(fAbscissa[i]) + (*f)(-fAbscissa[i]) ) ; >> 1219 } >> 1220 delete[] fAbscissa; >> 1221 delete[] fWeight; >> 1222 return integral ; 1212 } 1223 } 1213 1224 1214 ///////////////////////////////////////////// 1225 //////////////////////////////////////////////////////////////////////////// 1215 // 1226 // 1216 // Method involving Jacobi polynomials 1227 // Method involving Jacobi polynomials 1217 // 1228 // 1218 ///////////////////////////////////////////// 1229 //////////////////////////////////////////////////////////////////////////// 1219 // 1230 // 1220 // Gauss-Jacobi method for integration of ((1 1231 // Gauss-Jacobi method for integration of ((1-x)^alpha)*((1+x)^beta)*f(x) 1221 // from minus unit to plus unit . 1232 // from minus unit to plus unit . 1222 // 1233 // 1223 1234 1224 template <class T, class F> << 1235 template <class T, class F> 1225 G4double G4Integrator<T, F>::Jacobi(T& typeT, << 1236 G4double G4Integrator<T,F>::Jacobi( T& typeT, F f, G4double alpha, 1226 G4double << 1237 G4double beta, G4int nJacobi) 1227 { << 1238 { 1228 const G4double tolerance = 1.0e-12; << 1239 const G4double tolerance = 1.0e-12 ; 1229 const G4double maxNumber = 12; << 1240 const G4double maxNumber = 12 ; 1230 G4int i, k, j; << 1241 G4int i, k, j ; 1231 G4double alphaBeta, alphaReduced, betaReduc << 1242 G4double alphaBeta, alphaReduced, betaReduced, root1=0., root2=0., root3=0. ; 1232 << 1243 G4double a, b, c, nwt1, nwt2, nwt3, nwt, temp, root=0., rootTemp ; 1233 G4double a, b, c, nwt1, nwt2, nwt3, nwt, te << 1244 1234 << 1245 G4int fNumber = nJacobi ; 1235 G4int fNumber = nJacobi; << 1246 G4double* fAbscissa = new G4double[fNumber] ; 1236 G4double* fAbscissa = new G4double[fNumber] << 1247 G4double* fWeight = new G4double[fNumber] ; 1237 G4double* fWeight = new G4double[fNumber] << 1248 1238 << 1249 for (i=1;i<=nJacobi;i++) 1239 for(i = 1; i <= nJacobi; ++i) << 1250 { 1240 { << 1251 if (i == 1) 1241 if(i == 1) << 1252 { 1242 { << 1253 alphaReduced = alpha/nJacobi ; 1243 alphaReduced = alpha / nJacobi; << 1254 betaReduced = beta/nJacobi ; 1244 betaReduced = beta / nJacobi; << 1255 root1 = (1.0+alpha)*(2.78002/(4.0+nJacobi*nJacobi)+ 1245 root1 = (1.0 + alpha) * (2.78002 << 1256 0.767999*alphaReduced/nJacobi) ; 1246 0.767999 * alp << 1257 root2 = 1.0+1.48*alphaReduced+0.96002*betaReduced + 1247 root2 = 1.0 + 1.48 * alphaReduce << 1258 0.451998*alphaReduced*alphaReduced + 1248 0.451998 * alphaReduced * alpha << 1259 0.83001*alphaReduced*betaReduced ; 1249 0.83001 * alphaReduced * betaRe << 1260 root = 1.0-root1/root2 ; 1250 root = 1.0 - root1 / root2; << 1261 } 1251 } << 1262 else if (i == 2) 1252 else if(i == 2) << 1263 { 1253 { << 1264 root1=(4.1002+alpha)/((1.0+alpha)*(1.0+0.155998*alpha)) ; 1254 root1 = (4.1002 + alpha) / ((1.0 + alph << 1265 root2=1.0+0.06*(nJacobi-8.0)*(1.0+0.12*alpha)/nJacobi ; 1255 root2 = 1.0 + 0.06 * (nJacobi - 8.0) * << 1266 root3=1.0+0.012002*beta*(1.0+0.24997*std::fabs(alpha))/nJacobi ; 1256 root3 = << 1267 root -= (1.0-root)*root1*root2*root3 ; 1257 1.0 + 0.012002 * beta * (1.0 + 0.2499 << 1268 } 1258 root -= (1.0 - root) * root1 * root2 * << 1269 else if (i == 3) 1259 } << 1270 { 1260 else if(i == 3) << 1271 root1=(1.67001+0.27998*alpha)/(1.0+0.37002*alpha) ; 1261 { << 1272 root2=1.0+0.22*(nJacobi-8.0)/nJacobi ; 1262 root1 = (1.67001 + 0.27998 * alpha) / ( << 1273 root3=1.0+8.0*beta/((6.28001+beta)*nJacobi*nJacobi) ; 1263 root2 = 1.0 + 0.22 * (nJacobi - 8.0) / << 1274 root -= (fAbscissa[0]-root)*root1*root2*root3 ; 1264 root3 = 1.0 + 8.0 * beta / ((6.28001 + << 1275 } 1265 root -= (fAbscissa[0] - root) * root1 * << 1276 else if (i == nJacobi-1) 1266 } << 1277 { 1267 else if(i == nJacobi - 1) << 1278 root1=(1.0+0.235002*beta)/(0.766001+0.118998*beta) ; 1268 { << 1279 root2=1.0/(1.0+0.639002*(nJacobi-4.0)/(1.0+0.71001*(nJacobi-4.0))) ; 1269 root1 = (1.0 + 0.235002 * beta) / (0.76 << 1280 root3=1.0/(1.0+20.0*alpha/((7.5+alpha)*nJacobi*nJacobi)) ; 1270 root2 = 1.0 / (1.0 + 0.639002 * (nJacob << 1281 root += (root-fAbscissa[nJacobi-4])*root1*root2*root3 ; 1271 (1.0 + 0.71001 * << 1282 } 1272 root3 = 1.0 / (1.0 + 20.0 * alpha / ((7 << 1283 else if (i == nJacobi) 1273 root += (root - fAbscissa[nJacobi - 4]) << 1284 { 1274 } << 1285 root1 = (1.0+0.37002*beta)/(1.67001+0.27998*beta) ; 1275 else if(i == nJacobi) << 1286 root2 = 1.0/(1.0+0.22*(nJacobi-8.0)/nJacobi) ; 1276 { << 1287 root3 = 1.0/(1.0+8.0*alpha/((6.28002+alpha)*nJacobi*nJacobi)) ; 1277 root1 = (1.0 + 0.37002 * beta) / (1.670 << 1288 root += (root-fAbscissa[nJacobi-3])*root1*root2*root3 ; 1278 root2 = 1.0 / (1.0 + 0.22 * (nJacobi - << 1289 } 1279 root3 = << 1290 else 1280 1.0 / (1.0 + 8.0 * alpha / ((6.28002 << 1291 { 1281 root += (root - fAbscissa[nJacobi - 3]) << 1292 root = 3.0*fAbscissa[i-2]-3.0*fAbscissa[i-3]+fAbscissa[i-4] ; 1282 } << 1293 } 1283 else << 1294 alphaBeta = alpha + beta ; 1284 { << 1295 for (k=1;k<=maxNumber;k++) 1285 root = 3.0 * fAbscissa[i - 2] - 3.0 * f << 1296 { 1286 } << 1297 temp = 2.0 + alphaBeta ; 1287 alphaBeta = alpha + beta; << 1298 nwt1 = (alpha-beta+temp*root)/2.0 ; 1288 for(k = 1; k <= maxNumber; ++k) << 1299 nwt2 = 1.0 ; 1289 { << 1300 for (j=2;j<=nJacobi;j++) 1290 temp = 2.0 + alphaBeta; << 1301 { 1291 nwt1 = (alpha - beta + temp * root) / 2 << 1302 nwt3 = nwt2 ; 1292 nwt2 = 1.0; << 1303 nwt2 = nwt1 ; 1293 for(j = 2; j <= nJacobi; ++j) << 1304 temp = 2*j+alphaBeta ; 1294 { << 1305 a = 2*j*(j+alphaBeta)*(temp-2.0) ; 1295 nwt3 = nwt2; << 1306 b = (temp-1.0)*(alpha*alpha-beta*beta+temp*(temp-2.0)*root) ; 1296 nwt2 = nwt1; << 1307 c = 2.0*(j-1+alpha)*(j-1+beta)*temp ; 1297 temp = 2 * j + alphaBeta; << 1308 nwt1 = (b*nwt2-c*nwt3)/a ; 1298 a = 2 * j * (j + alphaBeta) * (tem << 1309 } 1299 b = (temp - 1.0) * << 1310 nwt = (nJacobi*(alpha - beta - temp*root)*nwt1 + 1300 (alpha * alpha - beta * beta + te << 1311 2.0*(nJacobi + alpha)*(nJacobi + beta)*nwt2)/ 1301 c = 2.0 * (j - 1 + alpha) * (j - 1 << 1312 (temp*(1.0 - root*root)) ; 1302 nwt1 = (b * nwt2 - c * nwt3) / a; << 1313 rootTemp = root ; 1303 } << 1314 root = rootTemp - nwt1/nwt ; 1304 nwt = (nJacobi * (alpha - beta - temp * << 1315 if (std::fabs(root-rootTemp) <= tolerance) 1305 2.0 * (nJacobi + alpha) * (nJaco << 1316 { 1306 (temp * (1.0 - root * root)); << 1317 break ; 1307 rootTemp = root; << 1318 } 1308 root = rootTemp - nwt1 / nwt; << 1319 } 1309 if(std::fabs(root - rootTemp) <= tolera << 1320 if (k > maxNumber) 1310 { << 1321 { 1311 break; << 1322 G4Exception("G4Integrator<T,F>::Jacobi(T,F, ...)", "Error", 1312 } << 1323 FatalException, "Too many (>12) iterations."); 1313 } << 1324 } 1314 if(k > maxNumber) << 1325 fAbscissa[i-1] = root ; 1315 { << 1326 fWeight[i-1] = std::exp(GammaLogarithm((G4double)(alpha+nJacobi)) + 1316 G4Exception("G4Integrator<T,F>::Jacobi( << 1327 GammaLogarithm((G4double)(beta+nJacobi)) - 1317 FatalException, "Too many ( << 1328 GammaLogarithm((G4double)(nJacobi+1.0)) - 1318 } << 1329 GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0))) 1319 fAbscissa[i - 1] = root; << 1330 *temp*std::pow(2.0,alphaBeta)/(nwt*nwt2) ; 1320 fWeight[i - 1] = << 1331 } 1321 std::exp(GammaLogarithm((G4double)(alph << 1332 1322 GammaLogarithm((G4double)(beta << 1333 // 1323 GammaLogarithm((G4double)(nJac << 1334 // Calculation of the integral 1324 GammaLogarithm((G4double)(nJac << 1335 // 1325 temp * std::pow(2.0, alphaBeta) / (nwt << 1336 1326 } << 1337 G4double integral = 0.0 ; 1327 << 1338 for(i=0;i<fNumber;i++) 1328 // << 1339 { 1329 // Calculation of the integral << 1340 integral += fWeight[i]*(typeT.*f)(fAbscissa[i]) ; 1330 // << 1341 } 1331 << 1342 delete[] fAbscissa; 1332 G4double integral = 0.0; << 1343 delete[] fWeight; 1333 for(i = 0; i < fNumber; ++i) << 1344 return integral ; 1334 { << 1335 integral += fWeight[i] * (typeT.*f)(fAbsc << 1336 } << 1337 delete[] fAbscissa; << 1338 delete[] fWeight; << 1339 return integral; << 1340 } 1345 } 1341 1346 >> 1347 1342 ///////////////////////////////////////////// 1348 ///////////////////////////////////////////////////////////////////////// 1343 // 1349 // 1344 // For use with 'this' pointer 1350 // For use with 'this' pointer 1345 1351 1346 template <class T, class F> << 1352 template <class T, class F> 1347 G4double G4Integrator<T, F>::Jacobi(T* ptrT, << 1353 G4double G4Integrator<T,F>::Jacobi( T* ptrT, F f, G4double alpha, 1348 G4int n) << 1354 G4double beta, G4int n) 1349 { 1355 { 1350 return Jacobi(*ptrT, f, alpha, beta, n); << 1356 return Jacobi(*ptrT,f,alpha,beta,n) ; 1351 } << 1357 } 1352 1358 1353 ///////////////////////////////////////////// 1359 ///////////////////////////////////////////////////////////////////////// 1354 // 1360 // 1355 // For use with global scope f << 1361 // For use with global scope f 1356 1362 1357 template <class T, class F> 1363 template <class T, class F> 1358 G4double G4Integrator<T, F>::Jacobi(G4double << 1364 G4double G4Integrator<T,F>::Jacobi( G4double (*f)(G4double), G4double alpha, 1359 G4double << 1365 G4double beta, G4int nJacobi) 1360 { 1366 { 1361 const G4double tolerance = 1.0e-12; << 1367 const G4double tolerance = 1.0e-12 ; 1362 const G4double maxNumber = 12; << 1368 const G4double maxNumber = 12 ; 1363 G4int i, k, j; << 1369 G4int i, k, j ; 1364 G4double alphaBeta, alphaReduced, betaReduc << 1370 G4double alphaBeta, alphaReduced, betaReduced, root1=0., root2=0., root3=0. ; 1365 << 1371 G4double a, b, c, nwt1, nwt2, nwt3, nwt, temp, root=0., rootTemp ; 1366 G4double a, b, c, nwt1, nwt2, nwt3, nwt, te << 1372 1367 << 1373 G4int fNumber = nJacobi ; 1368 G4int fNumber = nJacobi; << 1374 G4double* fAbscissa = new G4double[fNumber] ; 1369 G4double* fAbscissa = new G4double[fNumber] << 1375 G4double* fWeight = new G4double[fNumber] ; 1370 G4double* fWeight = new G4double[fNumber] << 1376 1371 << 1377 for (i=1;i<=nJacobi;i++) 1372 for(i = 1; i <= nJacobi; ++i) << 1378 { 1373 { << 1379 if (i == 1) 1374 if(i == 1) << 1380 { 1375 { << 1381 alphaReduced = alpha/nJacobi ; 1376 alphaReduced = alpha / nJacobi; << 1382 betaReduced = beta/nJacobi ; 1377 betaReduced = beta / nJacobi; << 1383 root1 = (1.0+alpha)*(2.78002/(4.0+nJacobi*nJacobi)+ 1378 root1 = (1.0 + alpha) * (2.78002 << 1384 0.767999*alphaReduced/nJacobi) ; 1379 0.767999 * alp << 1385 root2 = 1.0+1.48*alphaReduced+0.96002*betaReduced + 1380 root2 = 1.0 + 1.48 * alphaReduce << 1386 0.451998*alphaReduced*alphaReduced + 1381 0.451998 * alphaReduced * alpha << 1387 0.83001*alphaReduced*betaReduced ; 1382 0.83001 * alphaReduced * betaRe << 1388 root = 1.0-root1/root2 ; 1383 root = 1.0 - root1 / root2; << 1389 } 1384 } << 1390 else if (i == 2) 1385 else if(i == 2) << 1391 { 1386 { << 1392 root1=(4.1002+alpha)/((1.0+alpha)*(1.0+0.155998*alpha)) ; 1387 root1 = (4.1002 + alpha) / ((1.0 + alph << 1393 root2=1.0+0.06*(nJacobi-8.0)*(1.0+0.12*alpha)/nJacobi ; 1388 root2 = 1.0 + 0.06 * (nJacobi - 8.0) * << 1394 root3=1.0+0.012002*beta*(1.0+0.24997*std::fabs(alpha))/nJacobi ; 1389 root3 = << 1395 root -= (1.0-root)*root1*root2*root3 ; 1390 1.0 + 0.012002 * beta * (1.0 + 0.2499 << 1396 } 1391 root -= (1.0 - root) * root1 * root2 * << 1397 else if (i == 3) 1392 } << 1398 { 1393 else if(i == 3) << 1399 root1=(1.67001+0.27998*alpha)/(1.0+0.37002*alpha) ; 1394 { << 1400 root2=1.0+0.22*(nJacobi-8.0)/nJacobi ; 1395 root1 = (1.67001 + 0.27998 * alpha) / ( << 1401 root3=1.0+8.0*beta/((6.28001+beta)*nJacobi*nJacobi) ; 1396 root2 = 1.0 + 0.22 * (nJacobi - 8.0) / << 1402 root -= (fAbscissa[0]-root)*root1*root2*root3 ; 1397 root3 = 1.0 + 8.0 * beta / ((6.28001 + << 1403 } 1398 root -= (fAbscissa[0] - root) * root1 * << 1404 else if (i == nJacobi-1) 1399 } << 1405 { 1400 else if(i == nJacobi - 1) << 1406 root1=(1.0+0.235002*beta)/(0.766001+0.118998*beta) ; 1401 { << 1407 root2=1.0/(1.0+0.639002*(nJacobi-4.0)/(1.0+0.71001*(nJacobi-4.0))) ; 1402 root1 = (1.0 + 0.235002 * beta) / (0.76 << 1408 root3=1.0/(1.0+20.0*alpha/((7.5+alpha)*nJacobi*nJacobi)) ; 1403 root2 = 1.0 / (1.0 + 0.639002 * (nJacob << 1409 root += (root-fAbscissa[nJacobi-4])*root1*root2*root3 ; 1404 (1.0 + 0.71001 * << 1410 } 1405 root3 = 1.0 / (1.0 + 20.0 * alpha / ((7 << 1411 else if (i == nJacobi) 1406 root += (root - fAbscissa[nJacobi - 4]) << 1412 { 1407 } << 1413 root1 = (1.0+0.37002*beta)/(1.67001+0.27998*beta) ; 1408 else if(i == nJacobi) << 1414 root2 = 1.0/(1.0+0.22*(nJacobi-8.0)/nJacobi) ; 1409 { << 1415 root3 = 1.0/(1.0+8.0*alpha/((6.28002+alpha)*nJacobi*nJacobi)) ; 1410 root1 = (1.0 + 0.37002 * beta) / (1.670 << 1416 root += (root-fAbscissa[nJacobi-3])*root1*root2*root3 ; 1411 root2 = 1.0 / (1.0 + 0.22 * (nJacobi - << 1417 } 1412 root3 = << 1418 else 1413 1.0 / (1.0 + 8.0 * alpha / ((6.28002 << 1419 { 1414 root += (root - fAbscissa[nJacobi - 3]) << 1420 root = 3.0*fAbscissa[i-2]-3.0*fAbscissa[i-3]+fAbscissa[i-4] ; 1415 } << 1421 } 1416 else << 1422 alphaBeta = alpha + beta ; 1417 { << 1423 for (k=1;k<=maxNumber;k++) 1418 root = 3.0 * fAbscissa[i - 2] - 3.0 * f << 1424 { 1419 } << 1425 temp = 2.0 + alphaBeta ; 1420 alphaBeta = alpha + beta; << 1426 nwt1 = (alpha-beta+temp*root)/2.0 ; 1421 for(k = 1; k <= maxNumber; ++k) << 1427 nwt2 = 1.0 ; 1422 { << 1428 for (j=2;j<=nJacobi;j++) 1423 temp = 2.0 + alphaBeta; << 1429 { 1424 nwt1 = (alpha - beta + temp * root) / 2 << 1430 nwt3 = nwt2 ; 1425 nwt2 = 1.0; << 1431 nwt2 = nwt1 ; 1426 for(j = 2; j <= nJacobi; ++j) << 1432 temp = 2*j+alphaBeta ; 1427 { << 1433 a = 2*j*(j+alphaBeta)*(temp-2.0) ; 1428 nwt3 = nwt2; << 1434 b = (temp-1.0)*(alpha*alpha-beta*beta+temp*(temp-2.0)*root) ; 1429 nwt2 = nwt1; << 1435 c = 2.0*(j-1+alpha)*(j-1+beta)*temp ; 1430 temp = 2 * j + alphaBeta; << 1436 nwt1 = (b*nwt2-c*nwt3)/a ; 1431 a = 2 * j * (j + alphaBeta) * (tem << 1437 } 1432 b = (temp - 1.0) * << 1438 nwt = (nJacobi*(alpha - beta - temp*root)*nwt1 + 1433 (alpha * alpha - beta * beta + te << 1439 2.0*(nJacobi + alpha)*(nJacobi + beta)*nwt2) / 1434 c = 2.0 * (j - 1 + alpha) * (j - 1 << 1440 (temp*(1.0 - root*root)) ; 1435 nwt1 = (b * nwt2 - c * nwt3) / a; << 1441 rootTemp = root ; 1436 } << 1442 root = rootTemp - nwt1/nwt ; 1437 nwt = (nJacobi * (alpha - beta - temp * << 1443 if (std::fabs(root-rootTemp) <= tolerance) 1438 2.0 * (nJacobi + alpha) * (nJaco << 1444 { 1439 (temp * (1.0 - root * root)); << 1445 break ; 1440 rootTemp = root; << 1446 } 1441 root = rootTemp - nwt1 / nwt; << 1447 } 1442 if(std::fabs(root - rootTemp) <= tolera << 1448 if (k > maxNumber) 1443 { << 1449 { 1444 break; << 1450 G4Exception("G4Integrator<T,F>::Jacobi(...)", "Error", 1445 } << 1451 FatalException, "Too many (>12) iterations."); 1446 } << 1452 } 1447 if(k > maxNumber) << 1453 fAbscissa[i-1] = root ; 1448 { << 1454 fWeight[i-1] = 1449 G4Exception("G4Integrator<T,F>::Jacobi( << 1455 std::exp(GammaLogarithm((G4double)(alpha+nJacobi)) + 1450 "Too many (>12) iterations. << 1456 GammaLogarithm((G4double)(beta+nJacobi)) - 1451 } << 1457 GammaLogarithm((G4double)(nJacobi+1.0)) - 1452 fAbscissa[i - 1] = root; << 1458 GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0))) 1453 fWeight[i - 1] = << 1459 *temp*std::pow(2.0,alphaBeta)/(nwt*nwt2); 1454 std::exp(GammaLogarithm((G4double)(alph << 1460 } 1455 GammaLogarithm((G4double)(beta << 1461 1456 GammaLogarithm((G4double)(nJac << 1462 // 1457 GammaLogarithm((G4double)(nJac << 1463 // Calculation of the integral 1458 temp * std::pow(2.0, alphaBeta) / (nwt << 1464 // 1459 } << 1465 1460 << 1466 G4double integral = 0.0 ; 1461 // << 1467 for(i=0;i<fNumber;i++) 1462 // Calculation of the integral << 1468 { 1463 // << 1469 integral += fWeight[i]*(*f)(fAbscissa[i]) ; 1464 << 1470 } 1465 G4double integral = 0.0; << 1471 delete[] fAbscissa; 1466 for(i = 0; i < fNumber; ++i) << 1472 delete[] fWeight; 1467 { << 1473 return integral ; 1468 integral += fWeight[i] * (*f)(fAbscissa[i << 1469 } << 1470 delete[] fAbscissa; << 1471 delete[] fWeight; << 1472 return integral; << 1473 } 1474 } 1474 1475 1475 // 1476 // 1476 // 1477 // 1477 ///////////////////////////////////////////// 1478 /////////////////////////////////////////////////////////////////// 1478 1479