Gabriele Cosmo
1478 lines
46 KiB
Plaintext
1478 lines
46 KiB
Plaintext
//
|
|
// ********************************************************************
|
|
// * License and Disclaimer *
|
|
// * *
|
|
// * The Geant4 software is copyright of the Copyright Holders of *
|
|
// * the Geant4 Collaboration. It is provided under the terms and *
|
|
// * conditions of the Geant4 Software License, included in the file *
|
|
// * LICENSE and available at http://cern.ch/geant4/license . These *
|
|
// * include a list of copyright holders. *
|
|
// * *
|
|
// * Neither the authors of this software system, nor their employing *
|
|
// * institutes,nor the agencies providing financial support for this *
|
|
// * work make any representation or warranty, express or implied, *
|
|
// * regarding this software system or assume any liability for its *
|
|
// * use. Please see the license in the file LICENSE and URL above *
|
|
// * for the full disclaimer and the limitation of liability. *
|
|
// * *
|
|
// * This code implementation is the result of the scientific and *
|
|
// * technical work of the GEANT4 collaboration. *
|
|
// * By using, copying, modifying or distributing the software (or *
|
|
// * any work based on the software) you agree to acknowledge its *
|
|
// * use in resulting scientific publications, and indicate your *
|
|
// * acceptance of all terms of the Geant4 Software license. *
|
|
// ********************************************************************
|
|
//
|
|
// G4Integrator inline methods implementation
|
|
//
|
|
// Author: V.Grichine, 04.09.1999 - First implementation based on
|
|
// G4SimpleIntegration class with H.P.Wellisch, G.Cosmo, and
|
|
// E.TCherniaev advises
|
|
// --------------------------------------------------------------------
|
|
|
|
/////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Sympson integration method
|
|
//
|
|
/////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Integration of class member functions T::f by Simpson method.
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Simpson(T& typeT, F f, G4double xInitial,
|
|
G4double xFinal, G4int iterationNumber)
|
|
{
|
|
G4int i;
|
|
G4double step = (xFinal - xInitial) / iterationNumber;
|
|
G4double x = xInitial;
|
|
G4double xPlus = xInitial + 0.5 * step;
|
|
G4double mean = ((typeT.*f)(xInitial) + (typeT.*f)(xFinal)) * 0.5;
|
|
G4double sum = (typeT.*f)(xPlus);
|
|
|
|
for(i = 1; i < iterationNumber; ++i)
|
|
{
|
|
x += step;
|
|
xPlus += step;
|
|
mean += (typeT.*f)(x);
|
|
sum += (typeT.*f)(xPlus);
|
|
}
|
|
mean += 2.0 * sum;
|
|
|
|
return mean * step / 3.0;
|
|
}
|
|
|
|
/////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Integration of class member functions T::f by Simpson method.
|
|
// Convenient to use with 'this' pointer
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Simpson(T* ptrT, F f, G4double xInitial,
|
|
G4double xFinal, G4int iterationNumber)
|
|
{
|
|
G4int i;
|
|
G4double step = (xFinal - xInitial) / iterationNumber;
|
|
G4double x = xInitial;
|
|
G4double xPlus = xInitial + 0.5 * step;
|
|
G4double mean = ((ptrT->*f)(xInitial) + (ptrT->*f)(xFinal)) * 0.5;
|
|
G4double sum = (ptrT->*f)(xPlus);
|
|
|
|
for(i = 1; i < iterationNumber; ++i)
|
|
{
|
|
x += step;
|
|
xPlus += step;
|
|
mean += (ptrT->*f)(x);
|
|
sum += (ptrT->*f)(xPlus);
|
|
}
|
|
mean += 2.0 * sum;
|
|
|
|
return mean * step / 3.0;
|
|
}
|
|
|
|
/////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Integration of class member functions T::f by Simpson method.
|
|
// Convenient to use, when function f is defined in global scope, i.e. in main()
|
|
// program
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Simpson(G4double (*f)(G4double), G4double xInitial,
|
|
G4double xFinal, G4int iterationNumber)
|
|
{
|
|
G4int i;
|
|
G4double step = (xFinal - xInitial) / iterationNumber;
|
|
G4double x = xInitial;
|
|
G4double xPlus = xInitial + 0.5 * step;
|
|
G4double mean = ((*f)(xInitial) + (*f)(xFinal)) * 0.5;
|
|
G4double sum = (*f)(xPlus);
|
|
|
|
for(i = 1; i < iterationNumber; ++i)
|
|
{
|
|
x += step;
|
|
xPlus += step;
|
|
mean += (*f)(x);
|
|
sum += (*f)(xPlus);
|
|
}
|
|
mean += 2.0 * sum;
|
|
|
|
return mean * step / 3.0;
|
|
}
|
|
|
|
//////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Adaptive Gauss method
|
|
//
|
|
//////////////////////////////////////////////////////////////////////////
|
|
//
|
|
//
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Gauss(T& typeT, F f, G4double xInitial,
|
|
G4double xFinal)
|
|
{
|
|
static const G4double root = 1.0 / std::sqrt(3.0);
|
|
|
|
G4double xMean = (xInitial + xFinal) / 2.0;
|
|
G4double Step = (xFinal - xInitial) / 2.0;
|
|
G4double delta = Step * root;
|
|
G4double sum = ((typeT.*f)(xMean + delta) + (typeT.*f)(xMean - delta));
|
|
|
|
return sum * Step;
|
|
}
|
|
|
|
//////////////////////////////////////////////////////////////////////
|
|
//
|
|
//
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Gauss(T* ptrT, F f, G4double a, G4double b)
|
|
{
|
|
return Gauss(*ptrT, f, a, b);
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
//
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Gauss(G4double (*f)(G4double), G4double xInitial,
|
|
G4double xFinal)
|
|
{
|
|
static const G4double root = 1.0 / std::sqrt(3.0);
|
|
|
|
G4double xMean = (xInitial + xFinal) / 2.0;
|
|
G4double Step = (xFinal - xInitial) / 2.0;
|
|
G4double delta = Step * root;
|
|
G4double sum = ((*f)(xMean + delta) + (*f)(xMean - delta));
|
|
|
|
return sum * Step;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////////
|
|
//
|
|
//
|
|
|
|
template <class T, class F>
|
|
void G4Integrator<T, F>::AdaptGauss(T& typeT, F f, G4double xInitial,
|
|
G4double xFinal, G4double fTolerance,
|
|
G4double& sum, G4int& depth)
|
|
{
|
|
if(depth > 100)
|
|
{
|
|
G4cout << "G4Integrator<T,F>::AdaptGauss: WARNING !!!" << G4endl;
|
|
G4cout << "Function varies too rapidly to get stated accuracy in 100 steps "
|
|
<< G4endl;
|
|
|
|
return;
|
|
}
|
|
G4double xMean = (xInitial + xFinal) / 2.0;
|
|
G4double leftHalf = Gauss(typeT, f, xInitial, xMean);
|
|
G4double rightHalf = Gauss(typeT, f, xMean, xFinal);
|
|
G4double full = Gauss(typeT, f, xInitial, xFinal);
|
|
if(std::fabs(leftHalf + rightHalf - full) < fTolerance)
|
|
{
|
|
sum += full;
|
|
}
|
|
else
|
|
{
|
|
++depth;
|
|
AdaptGauss(typeT, f, xInitial, xMean, fTolerance, sum, depth);
|
|
AdaptGauss(typeT, f, xMean, xFinal, fTolerance, sum, depth);
|
|
}
|
|
}
|
|
|
|
template <class T, class F>
|
|
void G4Integrator<T, F>::AdaptGauss(T* ptrT, F f, G4double xInitial,
|
|
G4double xFinal, G4double fTolerance,
|
|
G4double& sum, G4int& depth)
|
|
{
|
|
AdaptGauss(*ptrT, f, xInitial, xFinal, fTolerance, sum, depth);
|
|
}
|
|
|
|
/////////////////////////////////////////////////////////////////////////
|
|
//
|
|
//
|
|
template <class T, class F>
|
|
void G4Integrator<T, F>::AdaptGauss(G4double (*f)(G4double), G4double xInitial,
|
|
G4double xFinal, G4double fTolerance,
|
|
G4double& sum, G4int& depth)
|
|
{
|
|
if(depth > 100)
|
|
{
|
|
G4cout << "G4SimpleIntegration::AdaptGauss: WARNING !!!" << G4endl;
|
|
G4cout << "Function varies too rapidly to get stated accuracy in 100 steps "
|
|
<< G4endl;
|
|
|
|
return;
|
|
}
|
|
G4double xMean = (xInitial + xFinal) / 2.0;
|
|
G4double leftHalf = Gauss(f, xInitial, xMean);
|
|
G4double rightHalf = Gauss(f, xMean, xFinal);
|
|
G4double full = Gauss(f, xInitial, xFinal);
|
|
if(std::fabs(leftHalf + rightHalf - full) < fTolerance)
|
|
{
|
|
sum += full;
|
|
}
|
|
else
|
|
{
|
|
++depth;
|
|
AdaptGauss(f, xInitial, xMean, fTolerance, sum, depth);
|
|
AdaptGauss(f, xMean, xFinal, fTolerance, sum, depth);
|
|
}
|
|
}
|
|
|
|
////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Adaptive Gauss integration with accuracy 'e'
|
|
// Convenient for using with class object typeT
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::AdaptiveGauss(T& typeT, F f, G4double xInitial,
|
|
G4double xFinal, G4double e)
|
|
{
|
|
G4int depth = 0;
|
|
G4double sum = 0.0;
|
|
AdaptGauss(typeT, f, xInitial, xFinal, e, sum, depth);
|
|
return sum;
|
|
}
|
|
|
|
////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Adaptive Gauss integration with accuracy 'e'
|
|
// Convenient for using with 'this' pointer
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::AdaptiveGauss(T* ptrT, F f, G4double xInitial,
|
|
G4double xFinal, G4double e)
|
|
{
|
|
return AdaptiveGauss(*ptrT, f, xInitial, xFinal, e);
|
|
}
|
|
|
|
////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Adaptive Gauss integration with accuracy 'e'
|
|
// Convenient for using with global scope function f
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::AdaptiveGauss(G4double (*f)(G4double),
|
|
G4double xInitial, G4double xFinal,
|
|
G4double e)
|
|
{
|
|
G4int depth = 0;
|
|
G4double sum = 0.0;
|
|
AdaptGauss(f, xInitial, xFinal, e, sum, depth);
|
|
return sum;
|
|
}
|
|
|
|
////////////////////////////////////////////////////////////////////////////
|
|
// Gauss integration methods involving ortogonal polynomials
|
|
////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Methods involving Legendre polynomials
|
|
//
|
|
/////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// The value nLegendre set the accuracy required, i.e the number of points
|
|
// where the function pFunction will be evaluated during integration.
|
|
// The function creates the arrays for abscissas and weights that used
|
|
// in Gauss-Legendre quadrature method.
|
|
// The values a and b are the limits of integration of the function f .
|
|
// nLegendre MUST BE EVEN !!!
|
|
// Returns the integral of the function f between a and b, by 2*fNumber point
|
|
// Gauss-Legendre integration: the function is evaluated exactly
|
|
// 2*fNumber times at interior points in the range of integration.
|
|
// Since the weights and abscissas are, in this case, symmetric around
|
|
// the midpoint of the range of integration, there are actually only
|
|
// fNumber distinct values of each.
|
|
// Convenient for using with some class object dataT
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Legendre(T& typeT, F f, G4double a, G4double b,
|
|
G4int nLegendre)
|
|
{
|
|
G4double nwt, nwt1, temp1, temp2, temp3, temp;
|
|
G4double xDiff, xMean, dx, integral;
|
|
|
|
const G4double tolerance = 1.6e-10;
|
|
G4int i, j, k = nLegendre;
|
|
G4int fNumber = (nLegendre + 1) / 2;
|
|
|
|
if(2 * fNumber != k)
|
|
{
|
|
G4Exception("G4Integrator<T,F>::Legendre(T&,F, ...)", "InvalidCall",
|
|
FatalException, "Invalid (odd) nLegendre in constructor.");
|
|
}
|
|
|
|
G4double* fAbscissa = new G4double[fNumber];
|
|
G4double* fWeight = new G4double[fNumber];
|
|
|
|
for(i = 1; i <= fNumber; ++i) // Loop over the desired roots
|
|
{
|
|
nwt = std::cos(CLHEP::pi * (i - 0.25) /
|
|
(k + 0.5)); // Initial root approximation
|
|
|
|
do // loop of Newton's method
|
|
{
|
|
temp1 = 1.0;
|
|
temp2 = 0.0;
|
|
for(j = 1; j <= k; ++j)
|
|
{
|
|
temp3 = temp2;
|
|
temp2 = temp1;
|
|
temp1 = ((2.0 * j - 1.0) * nwt * temp2 - (j - 1.0) * temp3) / j;
|
|
}
|
|
temp = k * (nwt * temp1 - temp2) / (nwt * nwt - 1.0);
|
|
nwt1 = nwt;
|
|
nwt = nwt1 - temp1 / temp; // Newton's method
|
|
} while(std::fabs(nwt - nwt1) > tolerance);
|
|
|
|
fAbscissa[fNumber - i] = nwt;
|
|
fWeight[fNumber - i] = 2.0 / ((1.0 - nwt * nwt) * temp * temp);
|
|
}
|
|
|
|
//
|
|
// Now we ready to get integral
|
|
//
|
|
|
|
xMean = 0.5 * (a + b);
|
|
xDiff = 0.5 * (b - a);
|
|
integral = 0.0;
|
|
for(i = 0; i < fNumber; ++i)
|
|
{
|
|
dx = xDiff * fAbscissa[i];
|
|
integral += fWeight[i] * ((typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx));
|
|
}
|
|
delete[] fAbscissa;
|
|
delete[] fWeight;
|
|
return integral *= xDiff;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Convenient for using with the pointer 'this'
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Legendre(T* ptrT, F f, G4double a, G4double b,
|
|
G4int nLegendre)
|
|
{
|
|
return Legendre(*ptrT, f, a, b, nLegendre);
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Convenient for using with global scope function f
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Legendre(G4double (*f)(G4double), G4double a,
|
|
G4double b, G4int nLegendre)
|
|
{
|
|
G4double nwt, nwt1, temp1, temp2, temp3, temp;
|
|
G4double xDiff, xMean, dx, integral;
|
|
|
|
const G4double tolerance = 1.6e-10;
|
|
G4int i, j, k = nLegendre;
|
|
G4int fNumber = (nLegendre + 1) / 2;
|
|
|
|
if(2 * fNumber != k)
|
|
{
|
|
G4Exception("G4Integrator<T,F>::Legendre(...)", "InvalidCall",
|
|
FatalException, "Invalid (odd) nLegendre in constructor.");
|
|
}
|
|
|
|
G4double* fAbscissa = new G4double[fNumber];
|
|
G4double* fWeight = new G4double[fNumber];
|
|
|
|
for(i = 1; i <= fNumber; i++) // Loop over the desired roots
|
|
{
|
|
nwt = std::cos(CLHEP::pi * (i - 0.25) /
|
|
(k + 0.5)); // Initial root approximation
|
|
|
|
do // loop of Newton's method
|
|
{
|
|
temp1 = 1.0;
|
|
temp2 = 0.0;
|
|
for(j = 1; j <= k; ++j)
|
|
{
|
|
temp3 = temp2;
|
|
temp2 = temp1;
|
|
temp1 = ((2.0 * j - 1.0) * nwt * temp2 - (j - 1.0) * temp3) / j;
|
|
}
|
|
temp = k * (nwt * temp1 - temp2) / (nwt * nwt - 1.0);
|
|
nwt1 = nwt;
|
|
nwt = nwt1 - temp1 / temp; // Newton's method
|
|
} while(std::fabs(nwt - nwt1) > tolerance);
|
|
|
|
fAbscissa[fNumber - i] = nwt;
|
|
fWeight[fNumber - i] = 2.0 / ((1.0 - nwt * nwt) * temp * temp);
|
|
}
|
|
|
|
//
|
|
// Now we ready to get integral
|
|
//
|
|
|
|
xMean = 0.5 * (a + b);
|
|
xDiff = 0.5 * (b - a);
|
|
integral = 0.0;
|
|
for(i = 0; i < fNumber; ++i)
|
|
{
|
|
dx = xDiff * fAbscissa[i];
|
|
integral += fWeight[i] * ((*f)(xMean + dx) + (*f)(xMean - dx));
|
|
}
|
|
delete[] fAbscissa;
|
|
delete[] fWeight;
|
|
|
|
return integral *= xDiff;
|
|
}
|
|
|
|
////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Returns the integral of the function to be pointed by T::f between a and b,
|
|
// by ten point Gauss-Legendre integration: the function is evaluated exactly
|
|
// ten times at interior points in the range of integration. Since the weights
|
|
// and abscissas are, in this case, symmetric around the midpoint of the
|
|
// range of integration, there are actually only five distinct values of each
|
|
// Convenient for using with class object typeT
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Legendre10(T& typeT, F f, G4double a, G4double b)
|
|
{
|
|
G4int i;
|
|
G4double xDiff, xMean, dx, integral;
|
|
|
|
// From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 916
|
|
|
|
static const G4double abscissa[] = { 0.148874338981631, 0.433395394129247,
|
|
0.679409568299024, 0.865063366688985,
|
|
0.973906528517172 };
|
|
|
|
static const G4double weight[] = { 0.295524224714753, 0.269266719309996,
|
|
0.219086362515982, 0.149451349150581,
|
|
0.066671344308688 };
|
|
xMean = 0.5 * (a + b);
|
|
xDiff = 0.5 * (b - a);
|
|
integral = 0.0;
|
|
for(i = 0; i < 5; ++i)
|
|
{
|
|
dx = xDiff * abscissa[i];
|
|
integral += weight[i] * ((typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx));
|
|
}
|
|
return integral *= xDiff;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Convenient for using with the pointer 'this'
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Legendre10(T* ptrT, F f, G4double a, G4double b)
|
|
{
|
|
return Legendre10(*ptrT, f, a, b);
|
|
}
|
|
|
|
//////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Convenient for using with global scope functions
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Legendre10(G4double (*f)(G4double), G4double a,
|
|
G4double b)
|
|
{
|
|
G4int i;
|
|
G4double xDiff, xMean, dx, integral;
|
|
|
|
// From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 916
|
|
|
|
static const G4double abscissa[] = { 0.148874338981631, 0.433395394129247,
|
|
0.679409568299024, 0.865063366688985,
|
|
0.973906528517172 };
|
|
|
|
static const G4double weight[] = { 0.295524224714753, 0.269266719309996,
|
|
0.219086362515982, 0.149451349150581,
|
|
0.066671344308688 };
|
|
xMean = 0.5 * (a + b);
|
|
xDiff = 0.5 * (b - a);
|
|
integral = 0.0;
|
|
for(i = 0; i < 5; ++i)
|
|
{
|
|
dx = xDiff * abscissa[i];
|
|
integral += weight[i] * ((*f)(xMean + dx) + (*f)(xMean - dx));
|
|
}
|
|
return integral *= xDiff;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Returns the integral of the function to be pointed by T::f between a and b,
|
|
// by 96 point Gauss-Legendre integration: the function is evaluated exactly
|
|
// ten Times at interior points in the range of integration. Since the weights
|
|
// and abscissas are, in this case, symmetric around the midpoint of the
|
|
// range of integration, there are actually only five distinct values of each
|
|
// Convenient for using with some class object typeT
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Legendre96(T& typeT, F f, G4double a, G4double b)
|
|
{
|
|
G4int i;
|
|
G4double xDiff, xMean, dx, integral;
|
|
|
|
// From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919
|
|
|
|
static const G4double abscissa[] = {
|
|
0.016276744849602969579, 0.048812985136049731112,
|
|
0.081297495464425558994, 0.113695850110665920911,
|
|
0.145973714654896941989, 0.178096882367618602759, // 6
|
|
|
|
0.210031310460567203603, 0.241743156163840012328,
|
|
0.273198812591049141487, 0.304364944354496353024,
|
|
0.335208522892625422616, 0.365696861472313635031, // 12
|
|
|
|
0.395797649828908603285, 0.425478988407300545365,
|
|
0.454709422167743008636, 0.483457973920596359768,
|
|
0.511694177154667673586, 0.539388108324357436227, // 18
|
|
|
|
0.566510418561397168404, 0.593032364777572080684,
|
|
0.618925840125468570386, 0.644163403784967106798,
|
|
0.668718310043916153953, 0.692564536642171561344, // 24
|
|
|
|
0.715676812348967626225, 0.738030643744400132851,
|
|
0.759602341176647498703, 0.780369043867433217604,
|
|
0.800308744139140817229, 0.819400310737931675539, // 30
|
|
|
|
0.837623511228187121494, 0.854959033434601455463,
|
|
0.871388505909296502874, 0.886894517402420416057,
|
|
0.901460635315852341319, 0.915071423120898074206, // 36
|
|
|
|
0.927712456722308690965, 0.939370339752755216932,
|
|
0.950032717784437635756, 0.959688291448742539300,
|
|
0.968326828463264212174, 0.975939174585136466453, // 42
|
|
|
|
0.982517263563014677447, 0.988054126329623799481,
|
|
0.992543900323762624572, 0.995981842987209290650,
|
|
0.998364375863181677724, 0.999689503883230766828 // 48
|
|
};
|
|
|
|
static const G4double weight[] = {
|
|
0.032550614492363166242, 0.032516118713868835987,
|
|
0.032447163714064269364, 0.032343822568575928429,
|
|
0.032206204794030250669, 0.032034456231992663218, // 6
|
|
|
|
0.031828758894411006535, 0.031589330770727168558,
|
|
0.031316425596862355813, 0.031010332586313837423,
|
|
0.030671376123669149014, 0.030299915420827593794, // 12
|
|
|
|
0.029896344136328385984, 0.029461089958167905970,
|
|
0.028994614150555236543, 0.028497411065085385646,
|
|
0.027970007616848334440, 0.027412962726029242823, // 18
|
|
|
|
0.026826866725591762198, 0.026212340735672413913,
|
|
0.025570036005349361499, 0.024900633222483610288,
|
|
0.024204841792364691282, 0.023483399085926219842, // 24
|
|
|
|
0.022737069658329374001, 0.021966644438744349195,
|
|
0.021172939892191298988, 0.020356797154333324595,
|
|
0.019519081140145022410, 0.018660679627411467385, // 30
|
|
|
|
0.017782502316045260838, 0.016885479864245172450,
|
|
0.015970562902562291381, 0.015038721026994938006,
|
|
0.014090941772314860916, 0.013128229566961572637, // 36
|
|
|
|
0.012151604671088319635, 0.011162102099838498591,
|
|
0.010160770535008415758, 0.009148671230783386633,
|
|
0.008126876925698759217, 0.007096470791153865269, // 42
|
|
|
|
0.006058545504235961683, 0.005014202742927517693,
|
|
0.003964554338444686674, 0.002910731817934946408,
|
|
0.001853960788946921732, 0.000796792065552012429 // 48
|
|
};
|
|
xMean = 0.5 * (a + b);
|
|
xDiff = 0.5 * (b - a);
|
|
integral = 0.0;
|
|
for(i = 0; i < 48; ++i)
|
|
{
|
|
dx = xDiff * abscissa[i];
|
|
integral += weight[i] * ((typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx));
|
|
}
|
|
return integral *= xDiff;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Convenient for using with the pointer 'this'
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Legendre96(T* ptrT, F f, G4double a, G4double b)
|
|
{
|
|
return Legendre96(*ptrT, f, a, b);
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Convenient for using with global scope function f
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Legendre96(G4double (*f)(G4double), G4double a,
|
|
G4double b)
|
|
{
|
|
G4int i;
|
|
G4double xDiff, xMean, dx, integral;
|
|
|
|
// From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919
|
|
|
|
static const G4double abscissa[] = {
|
|
0.016276744849602969579, 0.048812985136049731112,
|
|
0.081297495464425558994, 0.113695850110665920911,
|
|
0.145973714654896941989, 0.178096882367618602759, // 6
|
|
|
|
0.210031310460567203603, 0.241743156163840012328,
|
|
0.273198812591049141487, 0.304364944354496353024,
|
|
0.335208522892625422616, 0.365696861472313635031, // 12
|
|
|
|
0.395797649828908603285, 0.425478988407300545365,
|
|
0.454709422167743008636, 0.483457973920596359768,
|
|
0.511694177154667673586, 0.539388108324357436227, // 18
|
|
|
|
0.566510418561397168404, 0.593032364777572080684,
|
|
0.618925840125468570386, 0.644163403784967106798,
|
|
0.668718310043916153953, 0.692564536642171561344, // 24
|
|
|
|
0.715676812348967626225, 0.738030643744400132851,
|
|
0.759602341176647498703, 0.780369043867433217604,
|
|
0.800308744139140817229, 0.819400310737931675539, // 30
|
|
|
|
0.837623511228187121494, 0.854959033434601455463,
|
|
0.871388505909296502874, 0.886894517402420416057,
|
|
0.901460635315852341319, 0.915071423120898074206, // 36
|
|
|
|
0.927712456722308690965, 0.939370339752755216932,
|
|
0.950032717784437635756, 0.959688291448742539300,
|
|
0.968326828463264212174, 0.975939174585136466453, // 42
|
|
|
|
0.982517263563014677447, 0.988054126329623799481,
|
|
0.992543900323762624572, 0.995981842987209290650,
|
|
0.998364375863181677724, 0.999689503883230766828 // 48
|
|
};
|
|
|
|
static const G4double weight[] = {
|
|
0.032550614492363166242, 0.032516118713868835987,
|
|
0.032447163714064269364, 0.032343822568575928429,
|
|
0.032206204794030250669, 0.032034456231992663218, // 6
|
|
|
|
0.031828758894411006535, 0.031589330770727168558,
|
|
0.031316425596862355813, 0.031010332586313837423,
|
|
0.030671376123669149014, 0.030299915420827593794, // 12
|
|
|
|
0.029896344136328385984, 0.029461089958167905970,
|
|
0.028994614150555236543, 0.028497411065085385646,
|
|
0.027970007616848334440, 0.027412962726029242823, // 18
|
|
|
|
0.026826866725591762198, 0.026212340735672413913,
|
|
0.025570036005349361499, 0.024900633222483610288,
|
|
0.024204841792364691282, 0.023483399085926219842, // 24
|
|
|
|
0.022737069658329374001, 0.021966644438744349195,
|
|
0.021172939892191298988, 0.020356797154333324595,
|
|
0.019519081140145022410, 0.018660679627411467385, // 30
|
|
|
|
0.017782502316045260838, 0.016885479864245172450,
|
|
0.015970562902562291381, 0.015038721026994938006,
|
|
0.014090941772314860916, 0.013128229566961572637, // 36
|
|
|
|
0.012151604671088319635, 0.011162102099838498591,
|
|
0.010160770535008415758, 0.009148671230783386633,
|
|
0.008126876925698759217, 0.007096470791153865269, // 42
|
|
|
|
0.006058545504235961683, 0.005014202742927517693,
|
|
0.003964554338444686674, 0.002910731817934946408,
|
|
0.001853960788946921732, 0.000796792065552012429 // 48
|
|
};
|
|
xMean = 0.5 * (a + b);
|
|
xDiff = 0.5 * (b - a);
|
|
integral = 0.0;
|
|
for(i = 0; i < 48; ++i)
|
|
{
|
|
dx = xDiff * abscissa[i];
|
|
integral += weight[i] * ((*f)(xMean + dx) + (*f)(xMean - dx));
|
|
}
|
|
return integral *= xDiff;
|
|
}
|
|
|
|
//////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Methods involving Chebyshev polynomials
|
|
//
|
|
///////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Integrates function pointed by T::f from a to b by Gauss-Chebyshev
|
|
// quadrature method.
|
|
// Convenient for using with class object typeT
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Chebyshev(T& typeT, F f, G4double a, G4double b,
|
|
G4int nChebyshev)
|
|
{
|
|
G4int i;
|
|
G4double xDiff, xMean, dx, integral = 0.0;
|
|
|
|
G4int fNumber = nChebyshev; // Try to reduce fNumber twice ??
|
|
G4double cof = CLHEP::pi / fNumber;
|
|
G4double* fAbscissa = new G4double[fNumber];
|
|
G4double* fWeight = new G4double[fNumber];
|
|
for(i = 0; i < fNumber; ++i)
|
|
{
|
|
fAbscissa[i] = std::cos(cof * (i + 0.5));
|
|
fWeight[i] = cof * std::sqrt(1 - fAbscissa[i] * fAbscissa[i]);
|
|
}
|
|
|
|
//
|
|
// Now we ready to estimate the integral
|
|
//
|
|
|
|
xMean = 0.5 * (a + b);
|
|
xDiff = 0.5 * (b - a);
|
|
for(i = 0; i < fNumber; ++i)
|
|
{
|
|
dx = xDiff * fAbscissa[i];
|
|
integral += fWeight[i] * (typeT.*f)(xMean + dx);
|
|
}
|
|
delete[] fAbscissa;
|
|
delete[] fWeight;
|
|
return integral *= xDiff;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Convenient for using with 'this' pointer
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Chebyshev(T* ptrT, F f, G4double a, G4double b,
|
|
G4int n)
|
|
{
|
|
return Chebyshev(*ptrT, f, a, b, n);
|
|
}
|
|
|
|
////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// For use with global scope functions f
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Chebyshev(G4double (*f)(G4double), G4double a,
|
|
G4double b, G4int nChebyshev)
|
|
{
|
|
G4int i;
|
|
G4double xDiff, xMean, dx, integral = 0.0;
|
|
|
|
G4int fNumber = nChebyshev; // Try to reduce fNumber twice ??
|
|
G4double cof = CLHEP::pi / fNumber;
|
|
G4double* fAbscissa = new G4double[fNumber];
|
|
G4double* fWeight = new G4double[fNumber];
|
|
for(i = 0; i < fNumber; ++i)
|
|
{
|
|
fAbscissa[i] = std::cos(cof * (i + 0.5));
|
|
fWeight[i] = cof * std::sqrt(1 - fAbscissa[i] * fAbscissa[i]);
|
|
}
|
|
|
|
//
|
|
// Now we ready to estimate the integral
|
|
//
|
|
|
|
xMean = 0.5 * (a + b);
|
|
xDiff = 0.5 * (b - a);
|
|
for(i = 0; i < fNumber; ++i)
|
|
{
|
|
dx = xDiff * fAbscissa[i];
|
|
integral += fWeight[i] * (*f)(xMean + dx);
|
|
}
|
|
delete[] fAbscissa;
|
|
delete[] fWeight;
|
|
return integral *= xDiff;
|
|
}
|
|
|
|
//////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Method involving Laguerre polynomials
|
|
//
|
|
//////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Integral from zero to infinity of std::pow(x,alpha)*std::exp(-x)*f(x).
|
|
// The value of nLaguerre sets the accuracy.
|
|
// The function creates arrays fAbscissa[0,..,nLaguerre-1] and
|
|
// fWeight[0,..,nLaguerre-1] .
|
|
// Convenient for using with class object 'typeT' and (typeT.*f) function
|
|
// (T::f)
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Laguerre(T& typeT, F f, G4double alpha,
|
|
G4int nLaguerre)
|
|
{
|
|
const G4double tolerance = 1.0e-10;
|
|
const G4int maxNumber = 12;
|
|
G4int i, j, k;
|
|
G4double nwt = 0., nwt1, temp1, temp2, temp3, temp, cofi;
|
|
G4double integral = 0.0;
|
|
|
|
G4int fNumber = nLaguerre;
|
|
G4double* fAbscissa = new G4double[fNumber];
|
|
G4double* fWeight = new G4double[fNumber];
|
|
|
|
for(i = 1; i <= fNumber; ++i) // Loop over the desired roots
|
|
{
|
|
if(i == 1)
|
|
{
|
|
nwt = (1.0 + alpha) * (3.0 + 0.92 * alpha) /
|
|
(1.0 + 2.4 * fNumber + 1.8 * alpha);
|
|
}
|
|
else if(i == 2)
|
|
{
|
|
nwt += (15.0 + 6.25 * alpha) / (1.0 + 0.9 * alpha + 2.5 * fNumber);
|
|
}
|
|
else
|
|
{
|
|
cofi = i - 2;
|
|
nwt += ((1.0 + 2.55 * cofi) / (1.9 * cofi) +
|
|
1.26 * cofi * alpha / (1.0 + 3.5 * cofi)) *
|
|
(nwt - fAbscissa[i - 3]) / (1.0 + 0.3 * alpha);
|
|
}
|
|
for(k = 1; k <= maxNumber; ++k)
|
|
{
|
|
temp1 = 1.0;
|
|
temp2 = 0.0;
|
|
|
|
for(j = 1; j <= fNumber; ++j)
|
|
{
|
|
temp3 = temp2;
|
|
temp2 = temp1;
|
|
temp1 =
|
|
((2 * j - 1 + alpha - nwt) * temp2 - (j - 1 + alpha) * temp3) / j;
|
|
}
|
|
temp = (fNumber * temp1 - (fNumber + alpha) * temp2) / nwt;
|
|
nwt1 = nwt;
|
|
nwt = nwt1 - temp1 / temp;
|
|
|
|
if(std::fabs(nwt - nwt1) <= tolerance)
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
if(k > maxNumber)
|
|
{
|
|
G4Exception("G4Integrator<T,F>::Laguerre(T,F, ...)", "Error",
|
|
FatalException, "Too many (>12) iterations.");
|
|
}
|
|
|
|
fAbscissa[i - 1] = nwt;
|
|
fWeight[i - 1] = -std::exp(GammaLogarithm(alpha + fNumber) -
|
|
GammaLogarithm((G4double) fNumber)) /
|
|
(temp * fNumber * temp2);
|
|
}
|
|
|
|
//
|
|
// Integral evaluation
|
|
//
|
|
|
|
for(i = 0; i < fNumber; ++i)
|
|
{
|
|
integral += fWeight[i] * (typeT.*f)(fAbscissa[i]);
|
|
}
|
|
delete[] fAbscissa;
|
|
delete[] fWeight;
|
|
return integral;
|
|
}
|
|
|
|
//////////////////////////////////////////////////////////////////////
|
|
//
|
|
//
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Laguerre(T* ptrT, F f, G4double alpha,
|
|
G4int nLaguerre)
|
|
{
|
|
return Laguerre(*ptrT, f, alpha, nLaguerre);
|
|
}
|
|
|
|
////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// For use with global scope functions f
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Laguerre(G4double (*f)(G4double), G4double alpha,
|
|
G4int nLaguerre)
|
|
{
|
|
const G4double tolerance = 1.0e-10;
|
|
const G4int maxNumber = 12;
|
|
G4int i, j, k;
|
|
G4double nwt = 0., nwt1, temp1, temp2, temp3, temp, cofi;
|
|
G4double integral = 0.0;
|
|
|
|
G4int fNumber = nLaguerre;
|
|
G4double* fAbscissa = new G4double[fNumber];
|
|
G4double* fWeight = new G4double[fNumber];
|
|
|
|
for(i = 1; i <= fNumber; ++i) // Loop over the desired roots
|
|
{
|
|
if(i == 1)
|
|
{
|
|
nwt = (1.0 + alpha) * (3.0 + 0.92 * alpha) /
|
|
(1.0 + 2.4 * fNumber + 1.8 * alpha);
|
|
}
|
|
else if(i == 2)
|
|
{
|
|
nwt += (15.0 + 6.25 * alpha) / (1.0 + 0.9 * alpha + 2.5 * fNumber);
|
|
}
|
|
else
|
|
{
|
|
cofi = i - 2;
|
|
nwt += ((1.0 + 2.55 * cofi) / (1.9 * cofi) +
|
|
1.26 * cofi * alpha / (1.0 + 3.5 * cofi)) *
|
|
(nwt - fAbscissa[i - 3]) / (1.0 + 0.3 * alpha);
|
|
}
|
|
for(k = 1; k <= maxNumber; ++k)
|
|
{
|
|
temp1 = 1.0;
|
|
temp2 = 0.0;
|
|
|
|
for(j = 1; j <= fNumber; ++j)
|
|
{
|
|
temp3 = temp2;
|
|
temp2 = temp1;
|
|
temp1 =
|
|
((2 * j - 1 + alpha - nwt) * temp2 - (j - 1 + alpha) * temp3) / j;
|
|
}
|
|
temp = (fNumber * temp1 - (fNumber + alpha) * temp2) / nwt;
|
|
nwt1 = nwt;
|
|
nwt = nwt1 - temp1 / temp;
|
|
|
|
if(std::fabs(nwt - nwt1) <= tolerance)
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
if(k > maxNumber)
|
|
{
|
|
G4Exception("G4Integrator<T,F>::Laguerre( ...)", "Error", FatalException,
|
|
"Too many (>12) iterations.");
|
|
}
|
|
|
|
fAbscissa[i - 1] = nwt;
|
|
fWeight[i - 1] = -std::exp(GammaLogarithm(alpha + fNumber) -
|
|
GammaLogarithm((G4double) fNumber)) /
|
|
(temp * fNumber * temp2);
|
|
}
|
|
|
|
//
|
|
// Integral evaluation
|
|
//
|
|
|
|
for(i = 0; i < fNumber; i++)
|
|
{
|
|
integral += fWeight[i] * (*f)(fAbscissa[i]);
|
|
}
|
|
delete[] fAbscissa;
|
|
delete[] fWeight;
|
|
return integral;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Auxiliary function which returns the value of std::log(gamma-function(x))
|
|
// Returns the value ln(Gamma(xx) for xx > 0. Full accuracy is obtained for
|
|
// xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first.
|
|
// (Adapted from Numerical Recipes in C)
|
|
//
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::GammaLogarithm(G4double xx)
|
|
{
|
|
static const G4double cof[6] = { 76.18009172947146, -86.50532032941677,
|
|
24.01409824083091, -1.231739572450155,
|
|
0.1208650973866179e-2, -0.5395239384953e-5 };
|
|
G4int j;
|
|
G4double x = xx - 1.0;
|
|
G4double tmp = x + 5.5;
|
|
tmp -= (x + 0.5) * std::log(tmp);
|
|
G4double ser = 1.000000000190015;
|
|
|
|
for(j = 0; j <= 5; ++j)
|
|
{
|
|
x += 1.0;
|
|
ser += cof[j] / x;
|
|
}
|
|
return -tmp + std::log(2.5066282746310005 * ser);
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Method involving Hermite polynomials
|
|
//
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
//
|
|
// Gauss-Hermite method for integration of std::exp(-x*x)*f(x)
|
|
// from minus infinity to plus infinity .
|
|
//
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Hermite(T& typeT, F f, G4int nHermite)
|
|
{
|
|
const G4double tolerance = 1.0e-12;
|
|
const G4int maxNumber = 12;
|
|
|
|
G4int i, j, k;
|
|
G4double integral = 0.0;
|
|
G4double nwt = 0., nwt1, temp1, temp2, temp3, temp;
|
|
|
|
G4double piInMinusQ =
|
|
std::pow(CLHEP::pi, -0.25); // 1.0/std::sqrt(std::sqrt(pi)) ??
|
|
|
|
G4int fNumber = (nHermite + 1) / 2;
|
|
G4double* fAbscissa = new G4double[fNumber];
|
|
G4double* fWeight = new G4double[fNumber];
|
|
|
|
for(i = 1; i <= fNumber; ++i)
|
|
{
|
|
if(i == 1)
|
|
{
|
|
nwt = std::sqrt((G4double)(2 * nHermite + 1)) -
|
|
1.85575001 * std::pow((G4double)(2 * nHermite + 1), -0.16666999);
|
|
}
|
|
else if(i == 2)
|
|
{
|
|
nwt -= 1.14001 * std::pow((G4double) nHermite, 0.425999) / nwt;
|
|
}
|
|
else if(i == 3)
|
|
{
|
|
nwt = 1.86002 * nwt - 0.86002 * fAbscissa[0];
|
|
}
|
|
else if(i == 4)
|
|
{
|
|
nwt = 1.91001 * nwt - 0.91001 * fAbscissa[1];
|
|
}
|
|
else
|
|
{
|
|
nwt = 2.0 * nwt - fAbscissa[i - 3];
|
|
}
|
|
for(k = 1; k <= maxNumber; ++k)
|
|
{
|
|
temp1 = piInMinusQ;
|
|
temp2 = 0.0;
|
|
|
|
for(j = 1; j <= nHermite; ++j)
|
|
{
|
|
temp3 = temp2;
|
|
temp2 = temp1;
|
|
temp1 = nwt * std::sqrt(2.0 / j) * temp2 -
|
|
std::sqrt(((G4double)(j - 1)) / j) * temp3;
|
|
}
|
|
temp = std::sqrt((G4double) 2 * nHermite) * temp2;
|
|
nwt1 = nwt;
|
|
nwt = nwt1 - temp1 / temp;
|
|
|
|
if(std::fabs(nwt - nwt1) <= tolerance)
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
if(k > maxNumber)
|
|
{
|
|
G4Exception("G4Integrator<T,F>::Hermite(T,F, ...)", "Error",
|
|
FatalException, "Too many (>12) iterations.");
|
|
}
|
|
fAbscissa[i - 1] = nwt;
|
|
fWeight[i - 1] = 2.0 / (temp * temp);
|
|
}
|
|
|
|
//
|
|
// Integral calculation
|
|
//
|
|
|
|
for(i = 0; i < fNumber; ++i)
|
|
{
|
|
integral +=
|
|
fWeight[i] * ((typeT.*f)(fAbscissa[i]) + (typeT.*f)(-fAbscissa[i]));
|
|
}
|
|
delete[] fAbscissa;
|
|
delete[] fWeight;
|
|
return integral;
|
|
}
|
|
|
|
////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// For use with 'this' pointer
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Hermite(T* ptrT, F f, G4int n)
|
|
{
|
|
return Hermite(*ptrT, f, n);
|
|
}
|
|
|
|
////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// For use with global scope f
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Hermite(G4double (*f)(G4double), G4int nHermite)
|
|
{
|
|
const G4double tolerance = 1.0e-12;
|
|
const G4int maxNumber = 12;
|
|
|
|
G4int i, j, k;
|
|
G4double integral = 0.0;
|
|
G4double nwt = 0., nwt1, temp1, temp2, temp3, temp;
|
|
|
|
G4double piInMinusQ =
|
|
std::pow(CLHEP::pi, -0.25); // 1.0/std::sqrt(std::sqrt(pi)) ??
|
|
|
|
G4int fNumber = (nHermite + 1) / 2;
|
|
G4double* fAbscissa = new G4double[fNumber];
|
|
G4double* fWeight = new G4double[fNumber];
|
|
|
|
for(i = 1; i <= fNumber; ++i)
|
|
{
|
|
if(i == 1)
|
|
{
|
|
nwt = std::sqrt((G4double)(2 * nHermite + 1)) -
|
|
1.85575001 * std::pow((G4double)(2 * nHermite + 1), -0.16666999);
|
|
}
|
|
else if(i == 2)
|
|
{
|
|
nwt -= 1.14001 * std::pow((G4double) nHermite, 0.425999) / nwt;
|
|
}
|
|
else if(i == 3)
|
|
{
|
|
nwt = 1.86002 * nwt - 0.86002 * fAbscissa[0];
|
|
}
|
|
else if(i == 4)
|
|
{
|
|
nwt = 1.91001 * nwt - 0.91001 * fAbscissa[1];
|
|
}
|
|
else
|
|
{
|
|
nwt = 2.0 * nwt - fAbscissa[i - 3];
|
|
}
|
|
for(k = 1; k <= maxNumber; ++k)
|
|
{
|
|
temp1 = piInMinusQ;
|
|
temp2 = 0.0;
|
|
|
|
for(j = 1; j <= nHermite; ++j)
|
|
{
|
|
temp3 = temp2;
|
|
temp2 = temp1;
|
|
temp1 = nwt * std::sqrt(2.0 / j) * temp2 -
|
|
std::sqrt(((G4double)(j - 1)) / j) * temp3;
|
|
}
|
|
temp = std::sqrt((G4double) 2 * nHermite) * temp2;
|
|
nwt1 = nwt;
|
|
nwt = nwt1 - temp1 / temp;
|
|
|
|
if(std::fabs(nwt - nwt1) <= tolerance)
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
if(k > maxNumber)
|
|
{
|
|
G4Exception("G4Integrator<T,F>::Hermite(...)", "Error", FatalException,
|
|
"Too many (>12) iterations.");
|
|
}
|
|
fAbscissa[i - 1] = nwt;
|
|
fWeight[i - 1] = 2.0 / (temp * temp);
|
|
}
|
|
|
|
//
|
|
// Integral calculation
|
|
//
|
|
|
|
for(i = 0; i < fNumber; ++i)
|
|
{
|
|
integral += fWeight[i] * ((*f)(fAbscissa[i]) + (*f)(-fAbscissa[i]));
|
|
}
|
|
delete[] fAbscissa;
|
|
delete[] fWeight;
|
|
return integral;
|
|
}
|
|
|
|
////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Method involving Jacobi polynomials
|
|
//
|
|
////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Gauss-Jacobi method for integration of ((1-x)^alpha)*((1+x)^beta)*f(x)
|
|
// from minus unit to plus unit .
|
|
//
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Jacobi(T& typeT, F f, G4double alpha,
|
|
G4double beta, G4int nJacobi)
|
|
{
|
|
const G4double tolerance = 1.0e-12;
|
|
const G4double maxNumber = 12;
|
|
G4int i, k, j;
|
|
G4double alphaBeta, alphaReduced, betaReduced, root1 = 0., root2 = 0.,
|
|
root3 = 0.;
|
|
G4double a, b, c, nwt1, nwt2, nwt3, nwt, temp, root = 0., rootTemp;
|
|
|
|
G4int fNumber = nJacobi;
|
|
G4double* fAbscissa = new G4double[fNumber];
|
|
G4double* fWeight = new G4double[fNumber];
|
|
|
|
for(i = 1; i <= nJacobi; ++i)
|
|
{
|
|
if(i == 1)
|
|
{
|
|
alphaReduced = alpha / nJacobi;
|
|
betaReduced = beta / nJacobi;
|
|
root1 = (1.0 + alpha) * (2.78002 / (4.0 + nJacobi * nJacobi) +
|
|
0.767999 * alphaReduced / nJacobi);
|
|
root2 = 1.0 + 1.48 * alphaReduced + 0.96002 * betaReduced +
|
|
0.451998 * alphaReduced * alphaReduced +
|
|
0.83001 * alphaReduced * betaReduced;
|
|
root = 1.0 - root1 / root2;
|
|
}
|
|
else if(i == 2)
|
|
{
|
|
root1 = (4.1002 + alpha) / ((1.0 + alpha) * (1.0 + 0.155998 * alpha));
|
|
root2 = 1.0 + 0.06 * (nJacobi - 8.0) * (1.0 + 0.12 * alpha) / nJacobi;
|
|
root3 =
|
|
1.0 + 0.012002 * beta * (1.0 + 0.24997 * std::fabs(alpha)) / nJacobi;
|
|
root -= (1.0 - root) * root1 * root2 * root3;
|
|
}
|
|
else if(i == 3)
|
|
{
|
|
root1 = (1.67001 + 0.27998 * alpha) / (1.0 + 0.37002 * alpha);
|
|
root2 = 1.0 + 0.22 * (nJacobi - 8.0) / nJacobi;
|
|
root3 = 1.0 + 8.0 * beta / ((6.28001 + beta) * nJacobi * nJacobi);
|
|
root -= (fAbscissa[0] - root) * root1 * root2 * root3;
|
|
}
|
|
else if(i == nJacobi - 1)
|
|
{
|
|
root1 = (1.0 + 0.235002 * beta) / (0.766001 + 0.118998 * beta);
|
|
root2 = 1.0 / (1.0 + 0.639002 * (nJacobi - 4.0) /
|
|
(1.0 + 0.71001 * (nJacobi - 4.0)));
|
|
root3 = 1.0 / (1.0 + 20.0 * alpha / ((7.5 + alpha) * nJacobi * nJacobi));
|
|
root += (root - fAbscissa[nJacobi - 4]) * root1 * root2 * root3;
|
|
}
|
|
else if(i == nJacobi)
|
|
{
|
|
root1 = (1.0 + 0.37002 * beta) / (1.67001 + 0.27998 * beta);
|
|
root2 = 1.0 / (1.0 + 0.22 * (nJacobi - 8.0) / nJacobi);
|
|
root3 =
|
|
1.0 / (1.0 + 8.0 * alpha / ((6.28002 + alpha) * nJacobi * nJacobi));
|
|
root += (root - fAbscissa[nJacobi - 3]) * root1 * root2 * root3;
|
|
}
|
|
else
|
|
{
|
|
root = 3.0 * fAbscissa[i - 2] - 3.0 * fAbscissa[i - 3] + fAbscissa[i - 4];
|
|
}
|
|
alphaBeta = alpha + beta;
|
|
for(k = 1; k <= maxNumber; ++k)
|
|
{
|
|
temp = 2.0 + alphaBeta;
|
|
nwt1 = (alpha - beta + temp * root) / 2.0;
|
|
nwt2 = 1.0;
|
|
for(j = 2; j <= nJacobi; ++j)
|
|
{
|
|
nwt3 = nwt2;
|
|
nwt2 = nwt1;
|
|
temp = 2 * j + alphaBeta;
|
|
a = 2 * j * (j + alphaBeta) * (temp - 2.0);
|
|
b = (temp - 1.0) *
|
|
(alpha * alpha - beta * beta + temp * (temp - 2.0) * root);
|
|
c = 2.0 * (j - 1 + alpha) * (j - 1 + beta) * temp;
|
|
nwt1 = (b * nwt2 - c * nwt3) / a;
|
|
}
|
|
nwt = (nJacobi * (alpha - beta - temp * root) * nwt1 +
|
|
2.0 * (nJacobi + alpha) * (nJacobi + beta) * nwt2) /
|
|
(temp * (1.0 - root * root));
|
|
rootTemp = root;
|
|
root = rootTemp - nwt1 / nwt;
|
|
if(std::fabs(root - rootTemp) <= tolerance)
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
if(k > maxNumber)
|
|
{
|
|
G4Exception("G4Integrator<T,F>::Jacobi(T,F, ...)", "Error",
|
|
FatalException, "Too many (>12) iterations.");
|
|
}
|
|
fAbscissa[i - 1] = root;
|
|
fWeight[i - 1] =
|
|
std::exp(GammaLogarithm((G4double)(alpha + nJacobi)) +
|
|
GammaLogarithm((G4double)(beta + nJacobi)) -
|
|
GammaLogarithm((G4double)(nJacobi + 1.0)) -
|
|
GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0))) *
|
|
temp * std::pow(2.0, alphaBeta) / (nwt * nwt2);
|
|
}
|
|
|
|
//
|
|
// Calculation of the integral
|
|
//
|
|
|
|
G4double integral = 0.0;
|
|
for(i = 0; i < fNumber; ++i)
|
|
{
|
|
integral += fWeight[i] * (typeT.*f)(fAbscissa[i]);
|
|
}
|
|
delete[] fAbscissa;
|
|
delete[] fWeight;
|
|
return integral;
|
|
}
|
|
|
|
/////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// For use with 'this' pointer
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Jacobi(T* ptrT, F f, G4double alpha, G4double beta,
|
|
G4int n)
|
|
{
|
|
return Jacobi(*ptrT, f, alpha, beta, n);
|
|
}
|
|
|
|
/////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// For use with global scope f
|
|
|
|
template <class T, class F>
|
|
G4double G4Integrator<T, F>::Jacobi(G4double (*f)(G4double), G4double alpha,
|
|
G4double beta, G4int nJacobi)
|
|
{
|
|
const G4double tolerance = 1.0e-12;
|
|
const G4double maxNumber = 12;
|
|
G4int i, k, j;
|
|
G4double alphaBeta, alphaReduced, betaReduced, root1 = 0., root2 = 0.,
|
|
root3 = 0.;
|
|
G4double a, b, c, nwt1, nwt2, nwt3, nwt, temp, root = 0., rootTemp;
|
|
|
|
G4int fNumber = nJacobi;
|
|
G4double* fAbscissa = new G4double[fNumber];
|
|
G4double* fWeight = new G4double[fNumber];
|
|
|
|
for(i = 1; i <= nJacobi; ++i)
|
|
{
|
|
if(i == 1)
|
|
{
|
|
alphaReduced = alpha / nJacobi;
|
|
betaReduced = beta / nJacobi;
|
|
root1 = (1.0 + alpha) * (2.78002 / (4.0 + nJacobi * nJacobi) +
|
|
0.767999 * alphaReduced / nJacobi);
|
|
root2 = 1.0 + 1.48 * alphaReduced + 0.96002 * betaReduced +
|
|
0.451998 * alphaReduced * alphaReduced +
|
|
0.83001 * alphaReduced * betaReduced;
|
|
root = 1.0 - root1 / root2;
|
|
}
|
|
else if(i == 2)
|
|
{
|
|
root1 = (4.1002 + alpha) / ((1.0 + alpha) * (1.0 + 0.155998 * alpha));
|
|
root2 = 1.0 + 0.06 * (nJacobi - 8.0) * (1.0 + 0.12 * alpha) / nJacobi;
|
|
root3 =
|
|
1.0 + 0.012002 * beta * (1.0 + 0.24997 * std::fabs(alpha)) / nJacobi;
|
|
root -= (1.0 - root) * root1 * root2 * root3;
|
|
}
|
|
else if(i == 3)
|
|
{
|
|
root1 = (1.67001 + 0.27998 * alpha) / (1.0 + 0.37002 * alpha);
|
|
root2 = 1.0 + 0.22 * (nJacobi - 8.0) / nJacobi;
|
|
root3 = 1.0 + 8.0 * beta / ((6.28001 + beta) * nJacobi * nJacobi);
|
|
root -= (fAbscissa[0] - root) * root1 * root2 * root3;
|
|
}
|
|
else if(i == nJacobi - 1)
|
|
{
|
|
root1 = (1.0 + 0.235002 * beta) / (0.766001 + 0.118998 * beta);
|
|
root2 = 1.0 / (1.0 + 0.639002 * (nJacobi - 4.0) /
|
|
(1.0 + 0.71001 * (nJacobi - 4.0)));
|
|
root3 = 1.0 / (1.0 + 20.0 * alpha / ((7.5 + alpha) * nJacobi * nJacobi));
|
|
root += (root - fAbscissa[nJacobi - 4]) * root1 * root2 * root3;
|
|
}
|
|
else if(i == nJacobi)
|
|
{
|
|
root1 = (1.0 + 0.37002 * beta) / (1.67001 + 0.27998 * beta);
|
|
root2 = 1.0 / (1.0 + 0.22 * (nJacobi - 8.0) / nJacobi);
|
|
root3 =
|
|
1.0 / (1.0 + 8.0 * alpha / ((6.28002 + alpha) * nJacobi * nJacobi));
|
|
root += (root - fAbscissa[nJacobi - 3]) * root1 * root2 * root3;
|
|
}
|
|
else
|
|
{
|
|
root = 3.0 * fAbscissa[i - 2] - 3.0 * fAbscissa[i - 3] + fAbscissa[i - 4];
|
|
}
|
|
alphaBeta = alpha + beta;
|
|
for(k = 1; k <= maxNumber; ++k)
|
|
{
|
|
temp = 2.0 + alphaBeta;
|
|
nwt1 = (alpha - beta + temp * root) / 2.0;
|
|
nwt2 = 1.0;
|
|
for(j = 2; j <= nJacobi; ++j)
|
|
{
|
|
nwt3 = nwt2;
|
|
nwt2 = nwt1;
|
|
temp = 2 * j + alphaBeta;
|
|
a = 2 * j * (j + alphaBeta) * (temp - 2.0);
|
|
b = (temp - 1.0) *
|
|
(alpha * alpha - beta * beta + temp * (temp - 2.0) * root);
|
|
c = 2.0 * (j - 1 + alpha) * (j - 1 + beta) * temp;
|
|
nwt1 = (b * nwt2 - c * nwt3) / a;
|
|
}
|
|
nwt = (nJacobi * (alpha - beta - temp * root) * nwt1 +
|
|
2.0 * (nJacobi + alpha) * (nJacobi + beta) * nwt2) /
|
|
(temp * (1.0 - root * root));
|
|
rootTemp = root;
|
|
root = rootTemp - nwt1 / nwt;
|
|
if(std::fabs(root - rootTemp) <= tolerance)
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
if(k > maxNumber)
|
|
{
|
|
G4Exception("G4Integrator<T,F>::Jacobi(...)", "Error", FatalException,
|
|
"Too many (>12) iterations.");
|
|
}
|
|
fAbscissa[i - 1] = root;
|
|
fWeight[i - 1] =
|
|
std::exp(GammaLogarithm((G4double)(alpha + nJacobi)) +
|
|
GammaLogarithm((G4double)(beta + nJacobi)) -
|
|
GammaLogarithm((G4double)(nJacobi + 1.0)) -
|
|
GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0))) *
|
|
temp * std::pow(2.0, alphaBeta) / (nwt * nwt2);
|
|
}
|
|
|
|
//
|
|
// Calculation of the integral
|
|
//
|
|
|
|
G4double integral = 0.0;
|
|
for(i = 0; i < fNumber; ++i)
|
|
{
|
|
integral += fWeight[i] * (*f)(fAbscissa[i]);
|
|
}
|
|
delete[] fAbscissa;
|
|
delete[] fWeight;
|
|
return integral;
|
|
}
|
|
|
|
//
|
|
//
|
|
///////////////////////////////////////////////////////////////////
|