include/Geant4/G4Integrator.icc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // G4Integrator inline methods implementation
0027 //
0028 // Author: V.Grichine, 04.09.1999 - First implementation based on
0029 //         G4SimpleIntegration class with H.P.Wellisch, G.Cosmo, and
0030 //         E.TCherniaev advises
0031 // --------------------------------------------------------------------
0032
0033 /////////////////////////////////////////////////////////////////////
0034 //
0035 // Sympson integration method
0036 //
0037 /////////////////////////////////////////////////////////////////////
0038 //
0039 // Integration of class member functions T::f by Simpson method.
0040
0041 template <class T, class F>
0042 G4double G4Integrator<T, F>::Simpson(T& typeT, F f, G4double xInitial,
0043                                      G4double xFinal, G4int iterationNumber)
0044 {
0045   G4int i;
0046   G4double step  = (xFinal - xInitial) / iterationNumber;
0047   G4double x     = xInitial;
0048   G4double xPlus = xInitial + 0.5 * step;
0049   G4double mean  = ((typeT.*f)(xInitial) + (typeT.*f)(xFinal)) * 0.5;
0050   G4double sum   = (typeT.*f)(xPlus);
0051
0052   for(i = 1; i < iterationNumber; ++i)
0053   {
0054     x += step;
0055     xPlus += step;
0056     mean += (typeT.*f)(x);
0057     sum += (typeT.*f)(xPlus);
0058   }
0059   mean += 2.0 * sum;
0060
0061   return mean * step / 3.0;
0062 }
0063
0064 /////////////////////////////////////////////////////////////////////
0065 //
0066 // Integration of class member functions T::f by Simpson method.
0067 // Convenient to use with 'this' pointer
0068
0069 template <class T, class F>
0070 G4double G4Integrator<T, F>::Simpson(T* ptrT, F f, G4double xInitial,
0071                                      G4double xFinal, G4int iterationNumber)
0072 {
0073   G4int i;
0074   G4double step  = (xFinal - xInitial) / iterationNumber;
0075   G4double x     = xInitial;
0076   G4double xPlus = xInitial + 0.5 * step;
0077   G4double mean  = ((ptrT->*f)(xInitial) + (ptrT->*f)(xFinal)) * 0.5;
0078   G4double sum   = (ptrT->*f)(xPlus);
0079
0080   for(i = 1; i < iterationNumber; ++i)
0081   {
0082     x += step;
0083     xPlus += step;
0084     mean += (ptrT->*f)(x);
0085     sum += (ptrT->*f)(xPlus);
0086   }
0087   mean += 2.0 * sum;
0088
0089   return mean * step / 3.0;
0090 }
0091
0092 /////////////////////////////////////////////////////////////////////
0093 //
0094 // Integration of class member functions T::f by Simpson method.
0095 // Convenient to use, when function f is defined in global scope, i.e. in main()
0096 // program
0097
0098 template <class T, class F>
0099 G4double G4Integrator<T, F>::Simpson(G4double (*f)(G4double), G4double xInitial,
0100                                      G4double xFinal, G4int iterationNumber)
0101 {
0102   G4int i;
0103   G4double step  = (xFinal - xInitial) / iterationNumber;
0104   G4double x     = xInitial;
0105   G4double xPlus = xInitial + 0.5 * step;
0106   G4double mean  = ((*f)(xInitial) + (*f)(xFinal)) * 0.5;
0107   G4double sum   = (*f)(xPlus);
0108
0109   for(i = 1; i < iterationNumber; ++i)
0110   {
0111     x += step;
0112     xPlus += step;
0113     mean += (*f)(x);
0114     sum += (*f)(xPlus);
0115   }
0116   mean += 2.0 * sum;
0117
0118   return mean * step / 3.0;
0119 }
0120
0121 //////////////////////////////////////////////////////////////////////////
0122 //
0123 // Adaptive Gauss method
0124 //
0125 //////////////////////////////////////////////////////////////////////////
0126 //
0127 //
0128
0129 template <class T, class F>
0130 G4double G4Integrator<T, F>::Gauss(T& typeT, F f, G4double xInitial,
0131                                    G4double xFinal)
0132 {
0133   static const G4double root = 1.0 / std::sqrt(3.0);
0134
0135   G4double xMean = (xInitial + xFinal) / 2.0;
0136   G4double Step  = (xFinal - xInitial) / 2.0;
0137   G4double delta = Step * root;
0138   G4double sum   = ((typeT.*f)(xMean + delta) + (typeT.*f)(xMean - delta));
0139
0140   return sum * Step;
0141 }
0142
0143 //////////////////////////////////////////////////////////////////////
0144 //
0145 //
0146
0147 template <class T, class F>
0148 G4double G4Integrator<T, F>::Gauss(T* ptrT, F f, G4double a, G4double b)
0149 {
0150   return Gauss(*ptrT, f, a, b);
0151 }
0152
0153 ///////////////////////////////////////////////////////////////////////
0154 //
0155 //
0156
0157 template <class T, class F>
0158 G4double G4Integrator<T, F>::Gauss(G4double (*f)(G4double), G4double xInitial,
0159                                    G4double xFinal)
0160 {
0161   static const G4double root = 1.0 / std::sqrt(3.0);
0162
0163   G4double xMean = (xInitial + xFinal) / 2.0;
0164   G4double Step  = (xFinal - xInitial) / 2.0;
0165   G4double delta = Step * root;
0166   G4double sum   = ((*f)(xMean + delta) + (*f)(xMean - delta));
0167
0168   return sum * Step;
0169 }
0170
0171 ///////////////////////////////////////////////////////////////////////////
0172 //
0173 //
0174
0175 template <class T, class F>
0176 void G4Integrator<T, F>::AdaptGauss(T& typeT, F f, G4double xInitial,
0177                                     G4double xFinal, G4double fTolerance,
0178                                     G4double& sum, G4int& depth)
0179 {
0180   if(depth > 100)
0181   {
0182     G4cout << "G4Integrator<T,F>::AdaptGauss: WARNING !!!" << G4endl;
0183     G4cout << "Function varies too rapidly to get stated accuracy in 100 steps "
0184            << G4endl;
0185
0186     return;
0187   }
0188   G4double xMean     = (xInitial + xFinal) / 2.0;
0189   G4double leftHalf  = Gauss(typeT, f, xInitial, xMean);
0190   G4double rightHalf = Gauss(typeT, f, xMean, xFinal);
0191   G4double full      = Gauss(typeT, f, xInitial, xFinal);
0192   if(std::fabs(leftHalf + rightHalf - full) < fTolerance)
0193   {
0194     sum += full;
0195   }
0196   else
0197   {
0198     ++depth;
0199     AdaptGauss(typeT, f, xInitial, xMean, fTolerance, sum, depth);
0200     AdaptGauss(typeT, f, xMean, xFinal, fTolerance, sum, depth);
0201   }
0202 }
0203
0204 template <class T, class F>
0205 void G4Integrator<T, F>::AdaptGauss(T* ptrT, F f, G4double xInitial,
0206                                     G4double xFinal, G4double fTolerance,
0207                                     G4double& sum, G4int& depth)
0208 {
0209   AdaptGauss(*ptrT, f, xInitial, xFinal, fTolerance, sum, depth);
0210 }
0211
0212 /////////////////////////////////////////////////////////////////////////
0213 //
0214 //
0215 template <class T, class F>
0216 void G4Integrator<T, F>::AdaptGauss(G4double (*f)(G4double), G4double xInitial,
0217                                     G4double xFinal, G4double fTolerance,
0218                                     G4double& sum, G4int& depth)
0219 {
0220   if(depth > 100)
0221   {
0222     G4cout << "G4SimpleIntegration::AdaptGauss: WARNING !!!" << G4endl;
0223     G4cout << "Function varies too rapidly to get stated accuracy in 100 steps "
0224            << G4endl;
0225
0226     return;
0227   }
0228   G4double xMean     = (xInitial + xFinal) / 2.0;
0229   G4double leftHalf  = Gauss(f, xInitial, xMean);
0230   G4double rightHalf = Gauss(f, xMean, xFinal);
0231   G4double full      = Gauss(f, xInitial, xFinal);
0232   if(std::fabs(leftHalf + rightHalf - full) < fTolerance)
0233   {
0234     sum += full;
0235   }
0236   else
0237   {
0238     ++depth;
0239     AdaptGauss(f, xInitial, xMean, fTolerance, sum, depth);
0240     AdaptGauss(f, xMean, xFinal, fTolerance, sum, depth);
0241   }
0242 }
0243
0244 ////////////////////////////////////////////////////////////////////////
0245 //
0246 // Adaptive Gauss integration with accuracy 'e'
0247 // Convenient for using with class object typeT
0248
0249 template <class T, class F>
0250 G4double G4Integrator<T, F>::AdaptiveGauss(T& typeT, F f, G4double xInitial,
0251                                            G4double xFinal, G4double e)
0252 {
0253   G4int depth  = 0;
0254   G4double sum = 0.0;
0255   AdaptGauss(typeT, f, xInitial, xFinal, e, sum, depth);
0256   return sum;
0257 }
0258
0259 ////////////////////////////////////////////////////////////////////////
0260 //
0261 // Adaptive Gauss integration with accuracy 'e'
0262 // Convenient for using with 'this' pointer
0263
0264 template <class T, class F>
0265 G4double G4Integrator<T, F>::AdaptiveGauss(T* ptrT, F f, G4double xInitial,
0266                                            G4double xFinal, G4double e)
0267 {
0268   return AdaptiveGauss(*ptrT, f, xInitial, xFinal, e);
0269 }
0270
0271 ////////////////////////////////////////////////////////////////////////
0272 //
0273 // Adaptive Gauss integration with accuracy 'e'
0274 // Convenient for using with global scope function f
0275
0276 template <class T, class F>
0277 G4double G4Integrator<T, F>::AdaptiveGauss(G4double (*f)(G4double),
0278                                            G4double xInitial, G4double xFinal,
0279                                            G4double e)
0280 {
0281   G4int depth  = 0;
0282   G4double sum = 0.0;
0283   AdaptGauss(f, xInitial, xFinal, e, sum, depth);
0284   return sum;
0285 }
0286
0287 ////////////////////////////////////////////////////////////////////////////
0288 // Gauss integration methods involving ortogonal polynomials
0289 ////////////////////////////////////////////////////////////////////////////
0290 //
0291 // Methods involving Legendre polynomials
0292 //
0293 /////////////////////////////////////////////////////////////////////////
0294 //
0295 // The value nLegendre set the accuracy required, i.e the number of points
0296 // where the function pFunction will be evaluated during integration.
0297 // The function creates the arrays for abscissas and weights that used
0298 // in Gauss-Legendre quadrature method.
0299 // The values a and b are the limits of integration of the function  f .
0300 // nLegendre MUST BE EVEN !!!
0301 // Returns the integral of the function f between a and b, by 2*fNumber point
0302 // Gauss-Legendre integration: the function is evaluated exactly
0303 // 2*fNumber times at interior points in the range of integration.
0304 // Since the weights and abscissas are, in this case, symmetric around
0305 // the midpoint of the range of integration, there are actually only
0306 // fNumber distinct values of each.
0307 // Convenient for using with some class object dataT
0308
0309 template <class T, class F>
0310 G4double G4Integrator<T, F>::Legendre(T& typeT, F f, G4double a, G4double b,
0311                                       G4int nLegendre)
0312 {
0313   G4double nwt, nwt1, temp1, temp2, temp3, temp;
0314   G4double xDiff, xMean, dx, integral;
0315
0316   const G4double tolerance = 1.6e-10;
0317   G4int i, j, k = nLegendre;
0318   G4int fNumber = (nLegendre + 1) / 2;
0319
0320   if(2 * fNumber != k)
0321   {
0322     G4Exception("G4Integrator<T,F>::Legendre(T&,F, ...)", "InvalidCall",
0323                 FatalException, "Invalid (odd) nLegendre in constructor.");
0324   }
0325
0326   G4double* fAbscissa = new G4double[fNumber];
0327   G4double* fWeight   = new G4double[fNumber];
0328
0329   for(i = 1; i <= fNumber; ++i)  // Loop over the desired roots
0330   {
0331     nwt = std::cos(CLHEP::pi * (i - 0.25) /
0332                    (k + 0.5));  // Initial root approximation
0333
0334     do  // loop of Newton's method
0335     {
0336       temp1 = 1.0;
0337       temp2 = 0.0;
0338       for(j = 1; j <= k; ++j)
0339       {
0340         temp3 = temp2;
0341         temp2 = temp1;
0342         temp1 = ((2.0 * j - 1.0) * nwt * temp2 - (j - 1.0) * temp3) / j;
0343       }
0344       temp = k * (nwt * temp1 - temp2) / (nwt * nwt - 1.0);
0345       nwt1 = nwt;
0346       nwt  = nwt1 - temp1 / temp;  // Newton's method
0347     } while(std::fabs(nwt - nwt1) > tolerance);
0348
0349     fAbscissa[fNumber - i] = nwt;
0350     fWeight[fNumber - i]   = 2.0 / ((1.0 - nwt * nwt) * temp * temp);
0351   }
0352
0353   //
0354   // Now we ready to get integral
0355   //
0356
0357   xMean    = 0.5 * (a + b);
0358   xDiff    = 0.5 * (b - a);
0359   integral = 0.0;
0360   for(i = 0; i < fNumber; ++i)
0361   {
0362     dx = xDiff * fAbscissa[i];
0363     integral += fWeight[i] * ((typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx));
0364   }
0365   delete[] fAbscissa;
0366   delete[] fWeight;
0367   return integral *= xDiff;
0368 }
0369
0370 ///////////////////////////////////////////////////////////////////////
0371 //
0372 // Convenient for using with the pointer 'this'
0373
0374 template <class T, class F>
0375 G4double G4Integrator<T, F>::Legendre(T* ptrT, F f, G4double a, G4double b,
0376                                       G4int nLegendre)
0377 {
0378   return Legendre(*ptrT, f, a, b, nLegendre);
0379 }
0380
0381 ///////////////////////////////////////////////////////////////////////
0382 //
0383 // Convenient for using with global scope function f
0384
0385 template <class T, class F>
0386 G4double G4Integrator<T, F>::Legendre(G4double (*f)(G4double), G4double a,
0387                                       G4double b, G4int nLegendre)
0388 {
0389   G4double nwt, nwt1, temp1, temp2, temp3, temp;
0390   G4double xDiff, xMean, dx, integral;
0391
0392   const G4double tolerance = 1.6e-10;
0393   G4int i, j, k = nLegendre;
0394   G4int fNumber = (nLegendre + 1) / 2;
0395
0396   if(2 * fNumber != k)
0397   {
0398     G4Exception("G4Integrator<T,F>::Legendre(...)", "InvalidCall",
0399                 FatalException, "Invalid (odd) nLegendre in constructor.");
0400   }
0401
0402   G4double* fAbscissa = new G4double[fNumber];
0403   G4double* fWeight   = new G4double[fNumber];
0404
0405   for(i = 1; i <= fNumber; i++)  // Loop over the desired roots
0406   {
0407     nwt = std::cos(CLHEP::pi * (i - 0.25) /
0408                    (k + 0.5));  // Initial root approximation
0409
0410     do  // loop of Newton's method
0411     {
0412       temp1 = 1.0;
0413       temp2 = 0.0;
0414       for(j = 1; j <= k; ++j)
0415       {
0416         temp3 = temp2;
0417         temp2 = temp1;
0418         temp1 = ((2.0 * j - 1.0) * nwt * temp2 - (j - 1.0) * temp3) / j;
0419       }
0420       temp = k * (nwt * temp1 - temp2) / (nwt * nwt - 1.0);
0421       nwt1 = nwt;
0422       nwt  = nwt1 - temp1 / temp;  // Newton's method
0423     } while(std::fabs(nwt - nwt1) > tolerance);
0424
0425     fAbscissa[fNumber - i] = nwt;
0426     fWeight[fNumber - i]   = 2.0 / ((1.0 - nwt * nwt) * temp * temp);
0427   }
0428
0429   //
0430   // Now we ready to get integral
0431   //
0432
0433   xMean    = 0.5 * (a + b);
0434   xDiff    = 0.5 * (b - a);
0435   integral = 0.0;
0436   for(i = 0; i < fNumber; ++i)
0437   {
0438     dx = xDiff * fAbscissa[i];
0439     integral += fWeight[i] * ((*f)(xMean + dx) + (*f)(xMean - dx));
0440   }
0441   delete[] fAbscissa;
0442   delete[] fWeight;
0443
0444   return integral *= xDiff;
0445 }
0446
0447 ////////////////////////////////////////////////////////////////////////////
0448 //
0449 // Returns the integral of the function to be pointed by T::f between a and b,
0450 // by ten point Gauss-Legendre integration: the function is evaluated exactly
0451 // ten times at interior points in the range of integration. Since the weights
0452 // and abscissas are, in this case, symmetric around the midpoint of the
0453 // range of integration, there are actually only five distinct values of each
0454 // Convenient for using with class object typeT
0455
0456 template <class T, class F>
0457 G4double G4Integrator<T, F>::Legendre10(T& typeT, F f, G4double a, G4double b)
0458 {
0459   G4int i;
0460   G4double xDiff, xMean, dx, integral;
0461
0462   // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 916
0463
0464   static const G4double abscissa[] = { 0.148874338981631, 0.433395394129247,
0465                                        0.679409568299024, 0.865063366688985,
0466                                        0.973906528517172 };
0467
0468   static const G4double weight[] = { 0.295524224714753, 0.269266719309996,
0469                                      0.219086362515982, 0.149451349150581,
0470                                      0.066671344308688 };
0471   xMean                          = 0.5 * (a + b);
0472   xDiff                          = 0.5 * (b - a);
0473   integral                       = 0.0;
0474   for(i = 0; i < 5; ++i)
0475   {
0476     dx = xDiff * abscissa[i];
0477     integral += weight[i] * ((typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx));
0478   }
0479   return integral *= xDiff;
0480 }
0481
0482 ///////////////////////////////////////////////////////////////////////////
0483 //
0484 // Convenient for using with the pointer 'this'
0485
0486 template <class T, class F>
0487 G4double G4Integrator<T, F>::Legendre10(T* ptrT, F f, G4double a, G4double b)
0488 {
0489   return Legendre10(*ptrT, f, a, b);
0490 }
0491
0492 //////////////////////////////////////////////////////////////////////////
0493 //
0494 // Convenient for using with global scope functions
0495
0496 template <class T, class F>
0497 G4double G4Integrator<T, F>::Legendre10(G4double (*f)(G4double), G4double a,
0498                                         G4double b)
0499 {
0500   G4int i;
0501   G4double xDiff, xMean, dx, integral;
0502
0503   // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 916
0504
0505   static const G4double abscissa[] = { 0.148874338981631, 0.433395394129247,
0506                                        0.679409568299024, 0.865063366688985,
0507                                        0.973906528517172 };
0508
0509   static const G4double weight[] = { 0.295524224714753, 0.269266719309996,
0510                                      0.219086362515982, 0.149451349150581,
0511                                      0.066671344308688 };
0512   xMean                          = 0.5 * (a + b);
0513   xDiff                          = 0.5 * (b - a);
0514   integral                       = 0.0;
0515   for(i = 0; i < 5; ++i)
0516   {
0517     dx = xDiff * abscissa[i];
0518     integral += weight[i] * ((*f)(xMean + dx) + (*f)(xMean - dx));
0519   }
0520   return integral *= xDiff;
0521 }
0522
0523 ///////////////////////////////////////////////////////////////////////
0524 //
0525 // Returns the integral of the function to be pointed by T::f between a and b,
0526 // by 96 point Gauss-Legendre integration: the function is evaluated exactly
0527 // ten Times at interior points in the range of integration. Since the weights
0528 // and abscissas are, in this case, symmetric around the midpoint of the
0529 // range of integration, there are actually only five distinct values of each
0530 // Convenient for using with some class object typeT
0531
0532 template <class T, class F>
0533 G4double G4Integrator<T, F>::Legendre96(T& typeT, F f, G4double a, G4double b)
0534 {
0535   G4int i;
0536   G4double xDiff, xMean, dx, integral;
0537
0538   // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919
0539
0540   static const G4double abscissa[] = {
0541     0.016276744849602969579, 0.048812985136049731112,
0542     0.081297495464425558994, 0.113695850110665920911,
0543     0.145973714654896941989, 0.178096882367618602759,  // 6
0544
0545     0.210031310460567203603, 0.241743156163840012328,
0546     0.273198812591049141487, 0.304364944354496353024,
0547     0.335208522892625422616, 0.365696861472313635031,  // 12
0548
0549     0.395797649828908603285, 0.425478988407300545365,
0550     0.454709422167743008636, 0.483457973920596359768,
0551     0.511694177154667673586, 0.539388108324357436227,  // 18
0552
0553     0.566510418561397168404, 0.593032364777572080684,
0554     0.618925840125468570386, 0.644163403784967106798,
0555     0.668718310043916153953, 0.692564536642171561344,  // 24
0556
0557     0.715676812348967626225, 0.738030643744400132851,
0558     0.759602341176647498703, 0.780369043867433217604,
0559     0.800308744139140817229, 0.819400310737931675539,  // 30
0560
0561     0.837623511228187121494, 0.854959033434601455463,
0562     0.871388505909296502874, 0.886894517402420416057,
0563     0.901460635315852341319, 0.915071423120898074206,  // 36
0564
0565     0.927712456722308690965, 0.939370339752755216932,
0566     0.950032717784437635756, 0.959688291448742539300,
0567     0.968326828463264212174, 0.975939174585136466453,  // 42
0568
0569     0.982517263563014677447, 0.988054126329623799481,
0570     0.992543900323762624572, 0.995981842987209290650,
0571     0.998364375863181677724, 0.999689503883230766828  // 48
0572   };
0573
0574   static const G4double weight[] = {
0575     0.032550614492363166242, 0.032516118713868835987,
0576     0.032447163714064269364, 0.032343822568575928429,
0577     0.032206204794030250669, 0.032034456231992663218,  // 6
0578
0579     0.031828758894411006535, 0.031589330770727168558,
0580     0.031316425596862355813, 0.031010332586313837423,
0581     0.030671376123669149014, 0.030299915420827593794,  // 12
0582
0583     0.029896344136328385984, 0.029461089958167905970,
0584     0.028994614150555236543, 0.028497411065085385646,
0585     0.027970007616848334440, 0.027412962726029242823,  // 18
0586
0587     0.026826866725591762198, 0.026212340735672413913,
0588     0.025570036005349361499, 0.024900633222483610288,
0589     0.024204841792364691282, 0.023483399085926219842,  // 24
0590
0591     0.022737069658329374001, 0.021966644438744349195,
0592     0.021172939892191298988, 0.020356797154333324595,
0593     0.019519081140145022410, 0.018660679627411467385,  // 30
0594
0595     0.017782502316045260838, 0.016885479864245172450,
0596     0.015970562902562291381, 0.015038721026994938006,
0597     0.014090941772314860916, 0.013128229566961572637,  // 36
0598
0599     0.012151604671088319635, 0.011162102099838498591,
0600     0.010160770535008415758, 0.009148671230783386633,
0601     0.008126876925698759217, 0.007096470791153865269,  // 42
0602
0603     0.006058545504235961683, 0.005014202742927517693,
0604     0.003964554338444686674, 0.002910731817934946408,
0605     0.001853960788946921732, 0.000796792065552012429  // 48
0606   };
0607   xMean    = 0.5 * (a + b);
0608   xDiff    = 0.5 * (b - a);
0609   integral = 0.0;
0610   for(i = 0; i < 48; ++i)
0611   {
0612     dx = xDiff * abscissa[i];
0613     integral += weight[i] * ((typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx));
0614   }
0615   return integral *= xDiff;
0616 }
0617
0618 ///////////////////////////////////////////////////////////////////////
0619 //
0620 // Convenient for using with the pointer 'this'
0621
0622 template <class T, class F>
0623 G4double G4Integrator<T, F>::Legendre96(T* ptrT, F f, G4double a, G4double b)
0624 {
0625   return Legendre96(*ptrT, f, a, b);
0626 }
0627
0628 ///////////////////////////////////////////////////////////////////////
0629 //
0630 // Convenient for using with global scope function f
0631
0632 template <class T, class F>
0633 G4double G4Integrator<T, F>::Legendre96(G4double (*f)(G4double), G4double a,
0634                                         G4double b)
0635 {
0636   G4int i;
0637   G4double xDiff, xMean, dx, integral;
0638
0639   // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919
0640
0641   static const G4double abscissa[] = {
0642     0.016276744849602969579, 0.048812985136049731112,
0643     0.081297495464425558994, 0.113695850110665920911,
0644     0.145973714654896941989, 0.178096882367618602759,  // 6
0645
0646     0.210031310460567203603, 0.241743156163840012328,
0647     0.273198812591049141487, 0.304364944354496353024,
0648     0.335208522892625422616, 0.365696861472313635031,  // 12
0649
0650     0.395797649828908603285, 0.425478988407300545365,
0651     0.454709422167743008636, 0.483457973920596359768,
0652     0.511694177154667673586, 0.539388108324357436227,  // 18
0653
0654     0.566510418561397168404, 0.593032364777572080684,
0655     0.618925840125468570386, 0.644163403784967106798,
0656     0.668718310043916153953, 0.692564536642171561344,  // 24
0657
0658     0.715676812348967626225, 0.738030643744400132851,
0659     0.759602341176647498703, 0.780369043867433217604,
0660     0.800308744139140817229, 0.819400310737931675539,  // 30
0661
0662     0.837623511228187121494, 0.854959033434601455463,
0663     0.871388505909296502874, 0.886894517402420416057,
0664     0.901460635315852341319, 0.915071423120898074206,  // 36
0665
0666     0.927712456722308690965, 0.939370339752755216932,
0667     0.950032717784437635756, 0.959688291448742539300,
0668     0.968326828463264212174, 0.975939174585136466453,  // 42
0669
0670     0.982517263563014677447, 0.988054126329623799481,
0671     0.992543900323762624572, 0.995981842987209290650,
0672     0.998364375863181677724, 0.999689503883230766828  // 48
0673   };
0674
0675   static const G4double weight[] = {
0676     0.032550614492363166242, 0.032516118713868835987,
0677     0.032447163714064269364, 0.032343822568575928429,
0678     0.032206204794030250669, 0.032034456231992663218,  // 6
0679
0680     0.031828758894411006535, 0.031589330770727168558,
0681     0.031316425596862355813, 0.031010332586313837423,
0682     0.030671376123669149014, 0.030299915420827593794,  // 12
0683
0684     0.029896344136328385984, 0.029461089958167905970,
0685     0.028994614150555236543, 0.028497411065085385646,
0686     0.027970007616848334440, 0.027412962726029242823,  // 18
0687
0688     0.026826866725591762198, 0.026212340735672413913,
0689     0.025570036005349361499, 0.024900633222483610288,
0690     0.024204841792364691282, 0.023483399085926219842,  // 24
0691
0692     0.022737069658329374001, 0.021966644438744349195,
0693     0.021172939892191298988, 0.020356797154333324595,
0694     0.019519081140145022410, 0.018660679627411467385,  // 30
0695
0696     0.017782502316045260838, 0.016885479864245172450,
0697     0.015970562902562291381, 0.015038721026994938006,
0698     0.014090941772314860916, 0.013128229566961572637,  // 36
0699
0700     0.012151604671088319635, 0.011162102099838498591,
0701     0.010160770535008415758, 0.009148671230783386633,
0702     0.008126876925698759217, 0.007096470791153865269,  // 42
0703
0704     0.006058545504235961683, 0.005014202742927517693,
0705     0.003964554338444686674, 0.002910731817934946408,
0706     0.001853960788946921732, 0.000796792065552012429  // 48
0707   };
0708   xMean    = 0.5 * (a + b);
0709   xDiff    = 0.5 * (b - a);
0710   integral = 0.0;
0711   for(i = 0; i < 48; ++i)
0712   {
0713     dx = xDiff * abscissa[i];
0714     integral += weight[i] * ((*f)(xMean + dx) + (*f)(xMean - dx));
0715   }
0716   return integral *= xDiff;
0717 }
0718
0719 //////////////////////////////////////////////////////////////////////////////
0720 //
0721 // Methods involving Chebyshev polynomials
0722 //
0723 ///////////////////////////////////////////////////////////////////////////
0724 //
0725 // Integrates function pointed by T::f from a to b by Gauss-Chebyshev
0726 // quadrature method.
0727 // Convenient for using with class object typeT
0728
0729 template <class T, class F>
0730 G4double G4Integrator<T, F>::Chebyshev(T& typeT, F f, G4double a, G4double b,
0731                                        G4int nChebyshev)
0732 {
0733   G4int i;
0734   G4double xDiff, xMean, dx, integral = 0.0;
0735
0736   G4int fNumber       = nChebyshev;  // Try to reduce fNumber twice ??
0737   G4double cof        = CLHEP::pi / fNumber;
0738   G4double* fAbscissa = new G4double[fNumber];
0739   G4double* fWeight   = new G4double[fNumber];
0740   for(i = 0; i < fNumber; ++i)
0741   {
0742     fAbscissa[i] = std::cos(cof * (i + 0.5));
0743     fWeight[i]   = cof * std::sqrt(1 - fAbscissa[i] * fAbscissa[i]);
0744   }
0745
0746   //
0747   // Now we ready to estimate the integral
0748   //
0749
0750   xMean = 0.5 * (a + b);
0751   xDiff = 0.5 * (b - a);
0752   for(i = 0; i < fNumber; ++i)
0753   {
0754     dx = xDiff * fAbscissa[i];
0755     integral += fWeight[i] * (typeT.*f)(xMean + dx);
0756   }
0757   delete[] fAbscissa;
0758   delete[] fWeight;
0759   return integral *= xDiff;
0760 }
0761
0762 ///////////////////////////////////////////////////////////////////////
0763 //
0764 // Convenient for using with 'this' pointer
0765
0766 template <class T, class F>
0767 G4double G4Integrator<T, F>::Chebyshev(T* ptrT, F f, G4double a, G4double b,
0768                                        G4int n)
0769 {
0770   return Chebyshev(*ptrT, f, a, b, n);
0771 }
0772
0773 ////////////////////////////////////////////////////////////////////////
0774 //
0775 // For use with global scope functions f
0776
0777 template <class T, class F>
0778 G4double G4Integrator<T, F>::Chebyshev(G4double (*f)(G4double), G4double a,
0779                                        G4double b, G4int nChebyshev)
0780 {
0781   G4int i;
0782   G4double xDiff, xMean, dx, integral = 0.0;
0783
0784   G4int fNumber       = nChebyshev;  // Try to reduce fNumber twice ??
0785   G4double cof        = CLHEP::pi / fNumber;
0786   G4double* fAbscissa = new G4double[fNumber];
0787   G4double* fWeight   = new G4double[fNumber];
0788   for(i = 0; i < fNumber; ++i)
0789   {
0790     fAbscissa[i] = std::cos(cof * (i + 0.5));
0791     fWeight[i]   = cof * std::sqrt(1 - fAbscissa[i] * fAbscissa[i]);
0792   }
0793
0794   //
0795   // Now we ready to estimate the integral
0796   //
0797
0798   xMean = 0.5 * (a + b);
0799   xDiff = 0.5 * (b - a);
0800   for(i = 0; i < fNumber; ++i)
0801   {
0802     dx = xDiff * fAbscissa[i];
0803     integral += fWeight[i] * (*f)(xMean + dx);
0804   }
0805   delete[] fAbscissa;
0806   delete[] fWeight;
0807   return integral *= xDiff;
0808 }
0809
0810 //////////////////////////////////////////////////////////////////////
0811 //
0812 // Method involving Laguerre polynomials
0813 //
0814 //////////////////////////////////////////////////////////////////////
0815 //
0816 // Integral from zero to infinity of std::pow(x,alpha)*std::exp(-x)*f(x).
0817 // The value of nLaguerre sets the accuracy.
0818 // The function creates arrays fAbscissa[0,..,nLaguerre-1] and
0819 // fWeight[0,..,nLaguerre-1] .
0820 // Convenient for using with class object 'typeT' and (typeT.*f) function
0821 // (T::f)
0822
0823 template <class T, class F>
0824 G4double G4Integrator<T, F>::Laguerre(T& typeT, F f, G4double alpha,
0825                                       G4int nLaguerre)
0826 {
0827   const G4double tolerance = 1.0e-10;
0828   const G4int maxNumber    = 12;
0829   G4int i, j, k;
0830   G4double nwt      = 0., nwt1, temp1, temp2, temp3, temp, cofi;
0831   G4double integral = 0.0;
0832
0833   G4int fNumber       = nLaguerre;
0834   G4double* fAbscissa = new G4double[fNumber];
0835   G4double* fWeight   = new G4double[fNumber];
0836
0837   for(i = 1; i <= fNumber; ++i)  // Loop over the desired roots
0838   {
0839     if(i == 1)
0840     {
0841       nwt = (1.0 + alpha) * (3.0 + 0.92 * alpha) /
0842             (1.0 + 2.4 * fNumber + 1.8 * alpha);
0843     }
0844     else if(i == 2)
0845     {
0846       nwt += (15.0 + 6.25 * alpha) / (1.0 + 0.9 * alpha + 2.5 * fNumber);
0847     }
0848     else
0849     {
0850       cofi = i - 2;
0851       nwt += ((1.0 + 2.55 * cofi) / (1.9 * cofi) +
0852               1.26 * cofi * alpha / (1.0 + 3.5 * cofi)) *
0853              (nwt - fAbscissa[i - 3]) / (1.0 + 0.3 * alpha);
0854     }
0855     for(k = 1; k <= maxNumber; ++k)
0856     {
0857       temp1 = 1.0;
0858       temp2 = 0.0;
0859
0860       for(j = 1; j <= fNumber; ++j)
0861       {
0862         temp3 = temp2;
0863         temp2 = temp1;
0864         temp1 =
0865           ((2 * j - 1 + alpha - nwt) * temp2 - (j - 1 + alpha) * temp3) / j;
0866       }
0867       temp = (fNumber * temp1 - (fNumber + alpha) * temp2) / nwt;
0868       nwt1 = nwt;
0869       nwt  = nwt1 - temp1 / temp;
0870
0871       if(std::fabs(nwt - nwt1) <= tolerance)
0872       {
0873         break;
0874       }
0875     }
0876     if(k > maxNumber)
0877     {
0878       G4Exception("G4Integrator<T,F>::Laguerre(T,F, ...)", "Error",
0879                   FatalException, "Too many (>12) iterations.");
0880     }
0881
0882     fAbscissa[i - 1] = nwt;
0883     fWeight[i - 1]   = -std::exp(GammaLogarithm(alpha + fNumber) -
0884                                GammaLogarithm((G4double) fNumber)) /
0885                      (temp * fNumber * temp2);
0886   }
0887
0888   //
0889   // Integral evaluation
0890   //
0891
0892   for(i = 0; i < fNumber; ++i)
0893   {
0894     integral += fWeight[i] * (typeT.*f)(fAbscissa[i]);
0895   }
0896   delete[] fAbscissa;
0897   delete[] fWeight;
0898   return integral;
0899 }
0900
0901 //////////////////////////////////////////////////////////////////////
0902 //
0903 //
0904
0905 template <class T, class F>
0906 G4double G4Integrator<T, F>::Laguerre(T* ptrT, F f, G4double alpha,
0907                                       G4int nLaguerre)
0908 {
0909   return Laguerre(*ptrT, f, alpha, nLaguerre);
0910 }
0911
0912 ////////////////////////////////////////////////////////////////////////
0913 //
0914 // For use with global scope functions f
0915
0916 template <class T, class F>
0917 G4double G4Integrator<T, F>::Laguerre(G4double (*f)(G4double), G4double alpha,
0918                                       G4int nLaguerre)
0919 {
0920   const G4double tolerance = 1.0e-10;
0921   const G4int maxNumber    = 12;
0922   G4int i, j, k;
0923   G4double nwt      = 0., nwt1, temp1, temp2, temp3, temp, cofi;
0924   G4double integral = 0.0;
0925
0926   G4int fNumber       = nLaguerre;
0927   G4double* fAbscissa = new G4double[fNumber];
0928   G4double* fWeight   = new G4double[fNumber];
0929
0930   for(i = 1; i <= fNumber; ++i)  // Loop over the desired roots
0931   {
0932     if(i == 1)
0933     {
0934       nwt = (1.0 + alpha) * (3.0 + 0.92 * alpha) /
0935             (1.0 + 2.4 * fNumber + 1.8 * alpha);
0936     }
0937     else if(i == 2)
0938     {
0939       nwt += (15.0 + 6.25 * alpha) / (1.0 + 0.9 * alpha + 2.5 * fNumber);
0940     }
0941     else
0942     {
0943       cofi = i - 2;
0944       nwt += ((1.0 + 2.55 * cofi) / (1.9 * cofi) +
0945               1.26 * cofi * alpha / (1.0 + 3.5 * cofi)) *
0946              (nwt - fAbscissa[i - 3]) / (1.0 + 0.3 * alpha);
0947     }
0948     for(k = 1; k <= maxNumber; ++k)
0949     {
0950       temp1 = 1.0;
0951       temp2 = 0.0;
0952
0953       for(j = 1; j <= fNumber; ++j)
0954       {
0955         temp3 = temp2;
0956         temp2 = temp1;
0957         temp1 =
0958           ((2 * j - 1 + alpha - nwt) * temp2 - (j - 1 + alpha) * temp3) / j;
0959       }
0960       temp = (fNumber * temp1 - (fNumber + alpha) * temp2) / nwt;
0961       nwt1 = nwt;
0962       nwt  = nwt1 - temp1 / temp;
0963
0964       if(std::fabs(nwt - nwt1) <= tolerance)
0965       {
0966         break;
0967       }
0968     }
0969     if(k > maxNumber)
0970     {
0971       G4Exception("G4Integrator<T,F>::Laguerre( ...)", "Error", FatalException,
0972                   "Too many (>12) iterations.");
0973     }
0974
0975     fAbscissa[i - 1] = nwt;
0976     fWeight[i - 1]   = -std::exp(GammaLogarithm(alpha + fNumber) -
0977                                GammaLogarithm((G4double) fNumber)) /
0978                      (temp * fNumber * temp2);
0979   }
0980
0981   //
0982   // Integral evaluation
0983   //
0984
0985   for(i = 0; i < fNumber; i++)
0986   {
0987     integral += fWeight[i] * (*f)(fAbscissa[i]);
0988   }
0989   delete[] fAbscissa;
0990   delete[] fWeight;
0991   return integral;
0992 }
0993
0994 ///////////////////////////////////////////////////////////////////////
0995 //
0996 // Auxiliary function which returns the value of std::log(gamma-function(x))
0997 // Returns the value ln(Gamma(xx) for xx > 0.  Full accuracy is obtained for
0998 // xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first.
0999 // (Adapted from Numerical Recipes in C)
1000 //
1001
1002 template <class T, class F>
1003 G4double G4Integrator<T, F>::GammaLogarithm(G4double xx)
1004 {
1005   static const G4double cof[6] = { 76.18009172947146,     -86.50532032941677,
1006                                    24.01409824083091,     -1.231739572450155,
1007                                    0.1208650973866179e-2, -0.5395239384953e-5 };
1008   G4int j;
1009   G4double x   = xx - 1.0;
1010   G4double tmp = x + 5.5;
1011   tmp -= (x + 0.5) * std::log(tmp);
1012   G4double ser = 1.000000000190015;
1013
1014   for(j = 0; j <= 5; ++j)
1015   {
1016     x += 1.0;
1017     ser += cof[j] / x;
1018   }
1019   return -tmp + std::log(2.5066282746310005 * ser);
1020 }
1021
1022 ///////////////////////////////////////////////////////////////////////
1023 //
1024 // Method involving Hermite polynomials
1025 //
1026 ///////////////////////////////////////////////////////////////////////
1027 //
1028 //
1029 // Gauss-Hermite method for integration of std::exp(-x*x)*f(x)
1030 // from minus infinity to plus infinity .
1031 //
1032
1033 template <class T, class F>
1034 G4double G4Integrator<T, F>::Hermite(T& typeT, F f, G4int nHermite)
1035 {
1036   const G4double tolerance = 1.0e-12;
1037   const G4int maxNumber    = 12;
1038
1039   G4int i, j, k;
1040   G4double integral = 0.0;
1041   G4double nwt      = 0., nwt1, temp1, temp2, temp3, temp;
1042
1043   G4double piInMinusQ =
1044     std::pow(CLHEP::pi, -0.25);  // 1.0/std::sqrt(std::sqrt(pi)) ??
1045
1046   G4int fNumber       = (nHermite + 1) / 2;
1047   G4double* fAbscissa = new G4double[fNumber];
1048   G4double* fWeight   = new G4double[fNumber];
1049
1050   for(i = 1; i <= fNumber; ++i)
1051   {
1052     if(i == 1)
1053     {
1054       nwt = std::sqrt((G4double)(2 * nHermite + 1)) -
1055             1.85575001 * std::pow((G4double)(2 * nHermite + 1), -0.16666999);
1056     }
1057     else if(i == 2)
1058     {
1059       nwt -= 1.14001 * std::pow((G4double) nHermite, 0.425999) / nwt;
1060     }
1061     else if(i == 3)
1062     {
1063       nwt = 1.86002 * nwt - 0.86002 * fAbscissa[0];
1064     }
1065     else if(i == 4)
1066     {
1067       nwt = 1.91001 * nwt - 0.91001 * fAbscissa[1];
1068     }
1069     else
1070     {
1071       nwt = 2.0 * nwt - fAbscissa[i - 3];
1072     }
1073     for(k = 1; k <= maxNumber; ++k)
1074     {
1075       temp1 = piInMinusQ;
1076       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) * temp2 -
1083                 std::sqrt(((G4double)(j - 1)) / j) * temp3;
1084       }
1085       temp = std::sqrt((G4double) 2 * nHermite) * temp2;
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(T,F, ...)", "Error",
1097                   FatalException, "Too many (>12) iterations.");
1098     }
1099     fAbscissa[i - 1] = nwt;
1100     fWeight[i - 1]   = 2.0 / (temp * temp);
1101   }
1102
1103   //
1104   // Integral calculation
1105   //
1106
1107   for(i = 0; i < fNumber; ++i)
1108   {
1109     integral +=
1110       fWeight[i] * ((typeT.*f)(fAbscissa[i]) + (typeT.*f)(-fAbscissa[i]));
1111   }
1112   delete[] fAbscissa;
1113   delete[] fWeight;
1114   return integral;
1115 }
1116
1117 ////////////////////////////////////////////////////////////////////////
1118 //
1119 // For use with 'this' pointer
1120
1121 template <class T, class F>
1122 G4double G4Integrator<T, F>::Hermite(T* ptrT, F f, G4int n)
1123 {
1124   return Hermite(*ptrT, f, n);
1125 }
1126
1127 ////////////////////////////////////////////////////////////////////////
1128 //
1129 // For use with global scope f
1130
1131 template <class T, class F>
1132 G4double G4Integrator<T, F>::Hermite(G4double (*f)(G4double), G4int nHermite)
1133 {
1134   const G4double tolerance = 1.0e-12;
1135   const G4int maxNumber    = 12;
1136
1137   G4int i, j, k;
1138   G4double integral = 0.0;
1139   G4double nwt      = 0., nwt1, temp1, temp2, temp3, temp;
1140
1141   G4double piInMinusQ =
1142     std::pow(CLHEP::pi, -0.25);  // 1.0/std::sqrt(std::sqrt(pi)) ??
1143
1144   G4int fNumber       = (nHermite + 1) / 2;
1145   G4double* fAbscissa = new G4double[fNumber];
1146   G4double* fWeight   = new G4double[fNumber];
1147
1148   for(i = 1; i <= fNumber; ++i)
1149   {
1150     if(i == 1)
1151     {
1152       nwt = std::sqrt((G4double)(2 * nHermite + 1)) -
1153             1.85575001 * std::pow((G4double)(2 * nHermite + 1), -0.16666999);
1154     }
1155     else if(i == 2)
1156     {
1157       nwt -= 1.14001 * std::pow((G4double) nHermite, 0.425999) / nwt;
1158     }
1159     else if(i == 3)
1160     {
1161       nwt = 1.86002 * nwt - 0.86002 * fAbscissa[0];
1162     }
1163     else if(i == 4)
1164     {
1165       nwt = 1.91001 * nwt - 0.91001 * fAbscissa[1];
1166     }
1167     else
1168     {
1169       nwt = 2.0 * nwt - fAbscissa[i - 3];
1170     }
1171     for(k = 1; k <= maxNumber; ++k)
1172     {
1173       temp1 = piInMinusQ;
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) * temp2 -
1181                 std::sqrt(((G4double)(j - 1)) / j) * temp3;
1182       }
1183       temp = std::sqrt((G4double) 2 * nHermite) * temp2;
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(...)", "Error", FatalException,
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
1205   for(i = 0; i < fNumber; ++i)
1206   {
1207     integral += fWeight[i] * ((*f)(fAbscissa[i]) + (*f)(-fAbscissa[i]));
1208   }
1209   delete[] fAbscissa;
1210   delete[] fWeight;
1211   return integral;
1212 }
1213
1214 ////////////////////////////////////////////////////////////////////////////
1215 //
1216 // Method involving Jacobi polynomials
1217 //
1218 ////////////////////////////////////////////////////////////////////////////
1219 //
1220 // Gauss-Jacobi method for integration of ((1-x)^alpha)*((1+x)^beta)*f(x)
1221 // from minus unit to plus unit .
1222 //
1223
1224 template <class T, class F>
1225 G4double G4Integrator<T, F>::Jacobi(T& typeT, F f, G4double alpha,
1226                                     G4double beta, G4int nJacobi)
1227 {
1228   const G4double tolerance = 1.0e-12;
1229   const G4double maxNumber = 12;
1230   G4int i, k, j;
1231   G4double alphaBeta, alphaReduced, betaReduced, root1 = 0., root2 = 0.,
1232                                                  root3 = 0.;
1233   G4double a, b, c, nwt1, nwt2, nwt3, nwt, temp, root = 0., rootTemp;
1234
1235   G4int fNumber       = nJacobi;
1236   G4double* fAbscissa = new G4double[fNumber];
1237   G4double* fWeight   = new G4double[fNumber];
1238
1239   for(i = 1; i <= nJacobi; ++i)
1240   {
1241     if(i == 1)
1242     {
1243       alphaReduced = alpha / nJacobi;
1244       betaReduced  = beta / nJacobi;
1245       root1        = (1.0 + alpha) * (2.78002 / (4.0 + nJacobi * nJacobi) +
1246                                0.767999 * alphaReduced / nJacobi);
1247       root2        = 1.0 + 1.48 * alphaReduced + 0.96002 * betaReduced +
1248               0.451998 * alphaReduced * alphaReduced +
1249               0.83001 * alphaReduced * betaReduced;
1250       root = 1.0 - root1 / root2;
1251     }
1252     else if(i == 2)
1253     {
1254       root1 = (4.1002 + alpha) / ((1.0 + alpha) * (1.0 + 0.155998 * alpha));
1255       root2 = 1.0 + 0.06 * (nJacobi - 8.0) * (1.0 + 0.12 * alpha) / nJacobi;
1256       root3 =
1257         1.0 + 0.012002 * beta * (1.0 + 0.24997 * std::fabs(alpha)) / nJacobi;
1258       root -= (1.0 - root) * root1 * root2 * root3;
1259     }
1260     else if(i == 3)
1261     {
1262       root1 = (1.67001 + 0.27998 * alpha) / (1.0 + 0.37002 * alpha);
1263       root2 = 1.0 + 0.22 * (nJacobi - 8.0) / nJacobi;
1264       root3 = 1.0 + 8.0 * beta / ((6.28001 + beta) * nJacobi * nJacobi);
1265       root -= (fAbscissa[0] - root) * root1 * root2 * root3;
1266     }
1267     else if(i == nJacobi - 1)
1268     {
1269       root1 = (1.0 + 0.235002 * beta) / (0.766001 + 0.118998 * beta);
1270       root2 = 1.0 / (1.0 + 0.639002 * (nJacobi - 4.0) /
1271                              (1.0 + 0.71001 * (nJacobi - 4.0)));
1272       root3 = 1.0 / (1.0 + 20.0 * alpha / ((7.5 + alpha) * nJacobi * nJacobi));
1273       root += (root - fAbscissa[nJacobi - 4]) * root1 * root2 * root3;
1274     }
1275     else if(i == nJacobi)
1276     {
1277       root1 = (1.0 + 0.37002 * beta) / (1.67001 + 0.27998 * beta);
1278       root2 = 1.0 / (1.0 + 0.22 * (nJacobi - 8.0) / nJacobi);
1279       root3 =
1280         1.0 / (1.0 + 8.0 * alpha / ((6.28002 + alpha) * nJacobi * nJacobi));
1281       root += (root - fAbscissa[nJacobi - 3]) * root1 * root2 * root3;
1282     }
1283     else
1284     {
1285       root = 3.0 * fAbscissa[i - 2] - 3.0 * fAbscissa[i - 3] + fAbscissa[i - 4];
1286     }
1287     alphaBeta = alpha + beta;
1288     for(k = 1; k <= maxNumber; ++k)
1289     {
1290       temp = 2.0 + alphaBeta;
1291       nwt1 = (alpha - beta + temp * root) / 2.0;
1292       nwt2 = 1.0;
1293       for(j = 2; j <= nJacobi; ++j)
1294       {
1295         nwt3 = nwt2;
1296         nwt2 = nwt1;
1297         temp = 2 * j + alphaBeta;
1298         a    = 2 * j * (j + alphaBeta) * (temp - 2.0);
1299         b    = (temp - 1.0) *
1300             (alpha * alpha - beta * beta + temp * (temp - 2.0) * root);
1301         c    = 2.0 * (j - 1 + alpha) * (j - 1 + beta) * temp;
1302         nwt1 = (b * nwt2 - c * nwt3) / a;
1303       }
1304       nwt = (nJacobi * (alpha - beta - temp * root) * nwt1 +
1305              2.0 * (nJacobi + alpha) * (nJacobi + beta) * nwt2) /
1306             (temp * (1.0 - root * root));
1307       rootTemp = root;
1308       root     = rootTemp - nwt1 / nwt;
1309       if(std::fabs(root - rootTemp) <= tolerance)
1310       {
1311         break;
1312       }
1313     }
1314     if(k > maxNumber)
1315     {
1316       G4Exception("G4Integrator<T,F>::Jacobi(T,F, ...)", "Error",
1317                   FatalException, "Too many (>12) iterations.");
1318     }
1319     fAbscissa[i - 1] = root;
1320     fWeight[i - 1] =
1321       std::exp(GammaLogarithm((G4double)(alpha + nJacobi)) +
1322                GammaLogarithm((G4double)(beta + nJacobi)) -
1323                GammaLogarithm((G4double)(nJacobi + 1.0)) -
1324                GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0))) *
1325       temp * std::pow(2.0, alphaBeta) / (nwt * nwt2);
1326   }
1327
1328   //
1329   // Calculation of the integral
1330   //
1331
1332   G4double integral = 0.0;
1333   for(i = 0; i < fNumber; ++i)
1334   {
1335     integral += fWeight[i] * (typeT.*f)(fAbscissa[i]);
1336   }
1337   delete[] fAbscissa;
1338   delete[] fWeight;
1339   return integral;
1340 }
1341
1342 /////////////////////////////////////////////////////////////////////////
1343 //
1344 // For use with 'this' pointer
1345
1346 template <class T, class F>
1347 G4double G4Integrator<T, F>::Jacobi(T* ptrT, F f, G4double alpha, G4double beta,
1348                                     G4int n)
1349 {
1350   return Jacobi(*ptrT, f, alpha, beta, n);
1351 }
1352
1353 /////////////////////////////////////////////////////////////////////////
1354 //
1355 // For use with global scope f
1356
1357 template <class T, class F>
1358 G4double G4Integrator<T, F>::Jacobi(G4double (*f)(G4double), G4double alpha,
1359                                     G4double beta, G4int nJacobi)
1360 {
1361   const G4double tolerance = 1.0e-12;
1362   const G4double maxNumber = 12;
1363   G4int i, k, j;
1364   G4double alphaBeta, alphaReduced, betaReduced, root1 = 0., root2 = 0.,
1365                                                  root3 = 0.;
1366   G4double a, b, c, nwt1, nwt2, nwt3, nwt, temp, root = 0., rootTemp;
1367
1368   G4int fNumber       = nJacobi;
1369   G4double* fAbscissa = new G4double[fNumber];
1370   G4double* fWeight   = new G4double[fNumber];
1371
1372   for(i = 1; i <= nJacobi; ++i)
1373   {
1374     if(i == 1)
1375     {
1376       alphaReduced = alpha / nJacobi;
1377       betaReduced  = beta / nJacobi;
1378       root1        = (1.0 + alpha) * (2.78002 / (4.0 + nJacobi * nJacobi) +
1379                                0.767999 * alphaReduced / nJacobi);
1380       root2        = 1.0 + 1.48 * alphaReduced + 0.96002 * betaReduced +
1381               0.451998 * alphaReduced * alphaReduced +
1382               0.83001 * alphaReduced * betaReduced;
1383       root = 1.0 - root1 / root2;
1384     }
1385     else if(i == 2)
1386     {
1387       root1 = (4.1002 + alpha) / ((1.0 + alpha) * (1.0 + 0.155998 * alpha));
1388       root2 = 1.0 + 0.06 * (nJacobi - 8.0) * (1.0 + 0.12 * alpha) / nJacobi;
1389       root3 =
1390         1.0 + 0.012002 * beta * (1.0 + 0.24997 * std::fabs(alpha)) / nJacobi;
1391       root -= (1.0 - root) * root1 * root2 * root3;
1392     }
1393     else if(i == 3)
1394     {
1395       root1 = (1.67001 + 0.27998 * alpha) / (1.0 + 0.37002 * alpha);
1396       root2 = 1.0 + 0.22 * (nJacobi - 8.0) / nJacobi;
1397       root3 = 1.0 + 8.0 * beta / ((6.28001 + beta) * nJacobi * nJacobi);
1398       root -= (fAbscissa[0] - root) * root1 * root2 * root3;
1399     }
1400     else if(i == nJacobi - 1)
1401     {
1402       root1 = (1.0 + 0.235002 * beta) / (0.766001 + 0.118998 * beta);
1403       root2 = 1.0 / (1.0 + 0.639002 * (nJacobi - 4.0) /
1404                              (1.0 + 0.71001 * (nJacobi - 4.0)));
1405       root3 = 1.0 / (1.0 + 20.0 * alpha / ((7.5 + alpha) * nJacobi * nJacobi));
1406       root += (root - fAbscissa[nJacobi - 4]) * root1 * root2 * root3;
1407     }
1408     else if(i == nJacobi)
1409     {
1410       root1 = (1.0 + 0.37002 * beta) / (1.67001 + 0.27998 * beta);
1411       root2 = 1.0 / (1.0 + 0.22 * (nJacobi - 8.0) / nJacobi);
1412       root3 =
1413         1.0 / (1.0 + 8.0 * alpha / ((6.28002 + alpha) * nJacobi * nJacobi));
1414       root += (root - fAbscissa[nJacobi - 3]) * root1 * root2 * root3;
1415     }
1416     else
1417     {
1418       root = 3.0 * fAbscissa[i - 2] - 3.0 * fAbscissa[i - 3] + fAbscissa[i - 4];
1419     }
1420     alphaBeta = alpha + beta;
1421     for(k = 1; k <= maxNumber; ++k)
1422     {
1423       temp = 2.0 + alphaBeta;
1424       nwt1 = (alpha - beta + temp * root) / 2.0;
1425       nwt2 = 1.0;
1426       for(j = 2; j <= nJacobi; ++j)
1427       {
1428         nwt3 = nwt2;
1429         nwt2 = nwt1;
1430         temp = 2 * j + alphaBeta;
1431         a    = 2 * j * (j + alphaBeta) * (temp - 2.0);
1432         b    = (temp - 1.0) *
1433             (alpha * alpha - beta * beta + temp * (temp - 2.0) * root);
1434         c    = 2.0 * (j - 1 + alpha) * (j - 1 + beta) * temp;
1435         nwt1 = (b * nwt2 - c * nwt3) / a;
1436       }
1437       nwt = (nJacobi * (alpha - beta - temp * root) * nwt1 +
1438              2.0 * (nJacobi + alpha) * (nJacobi + beta) * nwt2) /
1439             (temp * (1.0 - root * root));
1440       rootTemp = root;
1441       root     = rootTemp - nwt1 / nwt;
1442       if(std::fabs(root - rootTemp) <= tolerance)
1443       {
1444         break;
1445       }
1446     }
1447     if(k > maxNumber)
1448     {
1449       G4Exception("G4Integrator<T,F>::Jacobi(...)", "Error", FatalException,
1450                   "Too many (>12) iterations.");
1451     }
1452     fAbscissa[i - 1] = root;
1453     fWeight[i - 1] =
1454       std::exp(GammaLogarithm((G4double)(alpha + nJacobi)) +
1455                GammaLogarithm((G4double)(beta + nJacobi)) -
1456                GammaLogarithm((G4double)(nJacobi + 1.0)) -
1457                GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0))) *
1458       temp * std::pow(2.0, alphaBeta) / (nwt * nwt2);
1459   }
1460
1461   //
1462   // Calculation of the integral
1463   //
1464
1465   G4double integral = 0.0;
1466   for(i = 0; i < fNumber; ++i)
1467   {
1468     integral += fWeight[i] * (*f)(fAbscissa[i]);
1469   }
1470   delete[] fAbscissa;
1471   delete[] fWeight;
1472   return integral;
1473 }
1474
1475 //
1476 //
1477 ///////////////////////////////////////////////////////////////////