#include <G4GaussHermiteQ.hh>

Inheritance diagram for G4GaussHermiteQ:

Public Member Functions
	G4GaussHermiteQ (function pFunction, G4int nHermite)

	G4GaussHermiteQ (const G4GaussHermiteQ &)=delete

G4GaussHermiteQ &	operator= (const G4GaussHermiteQ &)=delete

G4double	Integral () const

Public Member Functions inherited from G4VGaussianQuadrature
	G4VGaussianQuadrature (function pFunction)

virtual	~G4VGaussianQuadrature ()

	G4VGaussianQuadrature (const G4VGaussianQuadrature &)=delete

G4VGaussianQuadrature &	operator= (const G4VGaussianQuadrature &)=delete

G4double	GetAbscissa (G4int index) const

G4double	GetWeight (G4int index) const

G4int	GetNumber () const

Additional Inherited Members
Protected Member Functions inherited from G4VGaussianQuadrature
G4double	GammaLogarithm (G4double xx)

Protected Attributes inherited from G4VGaussianQuadrature
function	fFunction

G4double *	fAbscissa = nullptr

G4double *	fWeight = nullptr

G4int	fNumber = 0

Detailed Description

Definition at line 43 of file G4GaussHermiteQ.hh.

Constructor & Destructor Documentation

◆ G4GaussHermiteQ() [1/2]

G4GaussHermiteQ::G4GaussHermiteQ	(	function	pFunction,
		G4int	nHermite )

Definition at line 40 of file G4GaussHermiteQ.cc.

  : G4VGaussianQuadrature(pFunction)
{
  const G4double tolerance = 1.0e-12;
  const G4int maxNumber    = 12;
 
  G4int i = 1, j = 1, k = 1;
  G4double newton0 = 0.;
  G4double newton1 = 0.0, temp1 = 0.0, temp2 = 0.0, temp3 = 0.0, temp = 0.0;
  G4double piInMinusQ = std::pow(pi, -0.25);  // 1.0/std::sqrt(std::sqrt(pi)) ??
 
  fNumber   = (nHermite + 1) / 2;
  fAbscissa = new G4double[fNumber];
  fWeight   = new G4double[fNumber];
 
  for(i = 1; i <= fNumber; ++i)
  {
    if(i == 1)
    {
      newton0 =
        std::sqrt((G4double)(2 * nHermite + 1)) -
        1.85575001 * std::pow((G4double)(2 * nHermite + 1), -0.16666999);
    }
    else if(i == 2)
    {
      newton0 -= 1.14001 * std::pow((G4double) nHermite, 0.425999) / newton0;
    }
    else if(i == 3)
    {
      newton0 = 1.86002 * newton0 - 0.86002 * fAbscissa[0];
    }
    else if(i == 4)
    {
      newton0 = 1.91001 * newton0 - 0.91001 * fAbscissa[1];
    }
    else
    {
      newton0 = 2.0 * newton0 - 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 = newton0 * std::sqrt(2.0 / j) * temp2 -
                std::sqrt(((G4double)(j - 1)) / j) * temp3;
      }
      temp    = std::sqrt((G4double) 2 * nHermite) * temp2;
      newton1 = newton0;
      G4double ratio = std::numeric_limits<G4double>::max();
      if(temp > 0.0)
      {
        ratio = temp1 / temp;
      }
      newton0 = newton1 - ratio;
      if(std::fabs(newton0 - newton1) <= tolerance)
      {
        break;
      }
    }
    if(k > maxNumber)
    {
      G4Exception("G4GaussHermiteQ::G4GaussHermiteQ()", "OutOfRange",
                  FatalException,
                  "Too many iterations in Gauss-Hermite constructor.");
    }
    fAbscissa[i - 1] = newton0;
    fWeight[i - 1]   = 2.0 / (temp * temp);
  }
}

Referenced by G4GaussHermiteQ(), and operator=().

◆ G4GaussHermiteQ() [2/2]

G4GaussHermiteQ::G4GaussHermiteQ ( const G4GaussHermiteQ & )

delete

Member Function Documentation

◆ Integral()

G4double G4GaussHermiteQ::Integral ( ) const

Definition at line 119 of file G4GaussHermiteQ.cc.

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

◆ operator=()

G4GaussHermiteQ & G4GaussHermiteQ::operator= ( const G4GaussHermiteQ & )