G4GaussJacobiQ Class Reference

#include <G4GaussJacobiQ.hh>

Inheritance diagram for G4GaussJacobiQ:

Public Member Functions
	G4GaussJacobiQ (const G4GaussJacobiQ &)=delete
	G4GaussJacobiQ (function pFunction, G4double alpha, G4double beta, G4int nJacobi)
G4double	GetAbscissa (G4int index) const
G4int	GetNumber () const
G4double	GetWeight (G4int index) const
G4double	Integral () const
G4GaussJacobiQ &	operator= (const G4GaussJacobiQ &)=delete

Protected Member Functions
G4double	GammaLogarithm (G4double xx)

Protected Attributes
G4double *	fAbscissa = nullptr
function	fFunction
G4int	fNumber = 0
G4double *	fWeight = nullptr

Detailed Description

Definition at line 43 of file G4GaussJacobiQ.hh.

Constructor & Destructor Documentation

G4GaussJacobiQ() [1/2]

G4GaussJacobiQ::G4GaussJacobiQ	(	function	pFunction,
		G4double	alpha,
		G4double	beta,
		G4int	nJacobi
	)

Definition at line 38 of file G4GaussJacobiQ.cc.

  : G4VGaussianQuadrature(pFunction)
 
{
  const G4double tolerance = 1.0e-12;
  const G4double maxNumber = 12;
  G4int i = 1, k = 1;
  G4double root      = 0.;
  G4double alphaBeta = 0.0, alphaReduced = 0.0, betaReduced = 0.0, root1 = 0.0,
           root2 = 0.0, root3 = 0.0;
  G4double a = 0.0, b = 0.0, c = 0.0, newton1 = 0.0, newton2 = 0.0,
           newton3 = 0.0, newton0 = 0.0, temp = 0.0, rootTemp = 0.0;
 
  fNumber   = nJacobi;
  fAbscissa = new G4double[fNumber];
  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;
      newton1 = (alpha - beta + temp * root) / 2.0;
      newton2 = 1.0;
      for(G4int j = 2; j <= nJacobi; ++j)
      {
        newton3 = newton2;
        newton2 = newton1;
        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;
        newton1 = (b * newton2 - c * newton3) / a;
      }
      newton0 = (nJacobi * (alpha - beta - temp * root) * newton1 +
                 2.0 * (nJacobi + alpha) * (nJacobi + beta) * newton2) /
                (temp * (1.0 - root * root));
      rootTemp = root;
      root     = rootTemp - newton1 / newton0;
      if(std::fabs(root - rootTemp) <= tolerance)
      {
        break;
      }
    }
    if(k > maxNumber)
    {
      G4Exception("G4GaussJacobiQ::G4GaussJacobiQ()", "OutOfRange",
                  FatalException, "Too many iterations in constructor.");
    }
    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) / (newton0 * newton2);
  }
}

References alpha, anonymous_namespace{G4PionRadiativeDecayChannel.cc}::beta, G4VGaussianQuadrature::fAbscissa, FatalException, G4VGaussianQuadrature::fNumber, G4VGaussianQuadrature::fWeight, G4Exception(), and G4VGaussianQuadrature::GammaLogarithm().

G4GaussJacobiQ() [2/2]

G4GaussJacobiQ::G4GaussJacobiQ ( const G4GaussJacobiQ &  )

delete

Member Function Documentation

GammaLogarithm()

G4double G4VGaussianQuadrature::GammaLogarithm ( G4double  xx )

protectedinherited

Definition at line 70 of file G4VGaussianQuadrature.cc.

{
  // 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)
 
  static const G4double cof[6] = { 76.18009172947146,     -86.50532032941677,
                                   24.01409824083091,     -1.231739572450155,
                                   0.1208650973866179e-2, -0.5395239384953e-5 };
  G4double x                   = xx - 1.0;
  G4double tmp                 = x + 5.5;
  tmp -= (x + 0.5) * std::log(tmp);
  G4double ser = 1.000000000190015;
 
  for(size_t j = 0; j <= 5; ++j)
  {
    x += 1.0;
    ser += cof[j] / x;
  }
  return -tmp + std::log(2.5066282746310005 * ser);
}

Referenced by G4GaussJacobiQ(), and G4GaussLaguerreQ::G4GaussLaguerreQ().

GetAbscissa()

G4double G4VGaussianQuadrature::GetAbscissa ( G4int  index ) const

inherited

Definition at line 53 of file G4VGaussianQuadrature.cc.

{
  return fAbscissa[index];
}

References G4VGaussianQuadrature::fAbscissa.

GetNumber()

G4int G4VGaussianQuadrature::GetNumber ( ) const

inherited

Definition at line 63 of file G4VGaussianQuadrature.cc.

{ return fNumber; }

References G4VGaussianQuadrature::fNumber.

GetWeight()

G4double G4VGaussianQuadrature::GetWeight ( G4int  index ) const

inherited

Definition at line 58 of file G4VGaussianQuadrature.cc.

{
  return fWeight[index];
}

References G4VGaussianQuadrature::fWeight.

Integral()

G4double G4GaussJacobiQ::Integral

(

)

const

Definition at line 152 of file G4GaussJacobiQ.cc.

{
  G4double integral = 0.0;
  for(G4int i = 0; i < fNumber; ++i)
  {
    integral += fWeight[i] * fFunction(fAbscissa[i]);
  }
  return integral;
}

References G4VGaussianQuadrature::fAbscissa, G4VGaussianQuadrature::fFunction, G4VGaussianQuadrature::fNumber, and G4VGaussianQuadrature::fWeight.

operator=()

G4GaussJacobiQ & G4GaussJacobiQ::operator= ( const G4GaussJacobiQ &  )

delete

Field Documentation

fAbscissa

G4double* G4VGaussianQuadrature::fAbscissa = nullptr

protectedinherited

Definition at line 71 of file G4VGaussianQuadrature.hh.

Referenced by G4GaussChebyshevQ::G4GaussChebyshevQ(), G4GaussHermiteQ::G4GaussHermiteQ(), G4GaussJacobiQ(), G4GaussLaguerreQ::G4GaussLaguerreQ(), G4GaussLegendreQ::G4GaussLegendreQ(), G4VGaussianQuadrature::GetAbscissa(), G4GaussHermiteQ::Integral(), Integral(), G4GaussLaguerreQ::Integral(), G4GaussChebyshevQ::Integral(), G4GaussLegendreQ::Integral(), and G4VGaussianQuadrature::~G4VGaussianQuadrature().

fFunction

function G4VGaussianQuadrature::fFunction

protectedinherited

Definition at line 70 of file G4VGaussianQuadrature.hh.

Referenced by G4GaussLegendreQ::AccurateIntegral(), G4GaussHermiteQ::Integral(), Integral(), G4GaussLaguerreQ::Integral(), G4GaussChebyshevQ::Integral(), G4GaussLegendreQ::Integral(), and G4GaussLegendreQ::QuickIntegral().

fNumber

G4int G4VGaussianQuadrature::fNumber = 0

protectedinherited

Definition at line 73 of file G4VGaussianQuadrature.hh.

Referenced by G4GaussChebyshevQ::G4GaussChebyshevQ(), G4GaussHermiteQ::G4GaussHermiteQ(), G4GaussJacobiQ(), G4GaussLaguerreQ::G4GaussLaguerreQ(), G4GaussLegendreQ::G4GaussLegendreQ(), G4VGaussianQuadrature::GetNumber(), G4GaussHermiteQ::Integral(), Integral(), G4GaussLaguerreQ::Integral(), G4GaussChebyshevQ::Integral(), and G4GaussLegendreQ::Integral().

fWeight

G4double* G4VGaussianQuadrature::fWeight = nullptr

protectedinherited

Definition at line 72 of file G4VGaussianQuadrature.hh.

Referenced by G4GaussChebyshevQ::G4GaussChebyshevQ(), G4GaussHermiteQ::G4GaussHermiteQ(), G4GaussJacobiQ(), G4GaussLaguerreQ::G4GaussLaguerreQ(), G4GaussLegendreQ::G4GaussLegendreQ(), G4VGaussianQuadrature::GetWeight(), G4GaussHermiteQ::Integral(), Integral(), G4GaussLaguerreQ::Integral(), G4GaussChebyshevQ::Integral(), G4GaussLegendreQ::Integral(), and G4VGaussianQuadrature::~G4VGaussianQuadrature().

---

The documentation for this class was generated from the following files:

• source/global/HEPNumerics/include/G4GaussJacobiQ.hh
• source/global/HEPNumerics/src/G4GaussJacobiQ.cc