#include <G4Bessel.hh>

Public Member Functions
	G4Bessel ()

	~G4Bessel ()

G4double	I0 (G4double)

G4double	I1 (G4double)

G4double	K0 (G4double)

G4double	K1 (G4double)

G4double	pI0 (G4double)

G4double	pI1 (G4double)

G4double	pK0 (G4double)

G4double	pK1 (G4double)

Detailed Description

Definition at line 66 of file G4Bessel.hh.

Constructor & Destructor Documentation

◆ G4Bessel()

G4Bessel::G4Bessel ( )

Definition at line 72 of file G4Bessel.cc.

73{;}

◆ ~G4Bessel()

G4Bessel::~G4Bessel ( )

Definition at line 76 of file G4Bessel.cc.

77{;}

Member Function Documentation

◆ I0()

G4double G4Bessel::I0 ( G4double x )

Definition at line 80 of file G4Bessel.cc.

{
  const G4double P1 = 1.0,
                 P2 = 3.5156229,
                 P3 = 3.0899424,
                 P4 = 1.2067492,
                 P5 = 0.2659732,
                 P6 = 0.0360768,
                 P7 = 0.0045813;
  const G4double Q1 = 0.39894228,
                 Q2 = 0.01328592,
                 Q3 = 0.00225319,
                 Q4 =-0.00157565,
                 Q5 = 0.00916281,
                 Q6 =-0.02057706,
                 Q7 = 0.02635537,
                 Q8 =-0.01647633,
                 Q9 = 0.00392377;
  
  G4double I = 0.0;
  if (std::abs(x) < 3.75)
  {
    G4double y = G4Pow::GetInstance()->powN(x/3.75, 2);
    I = P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7)))));
  }
  else
  {
    G4double ax = std::abs(x);
    G4double y  = 3.75/ax;
    I  = G4Exp(ax) / std::sqrt(ax) *
      (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*(Q7+y*(Q8+y*Q9))))))));
  }
  return I;
}

Referenced by K0().

◆ I1()

G4double G4Bessel::I1 ( G4double x )

Definition at line 150 of file G4Bessel.cc.

{
  const G4double P1 = 0.5,
                 P2 = 0.87890594,
                 P3 = 0.51498869,
                 P4 = 0.15084934,
                 P5 = 0.02658733,
                 P6 = 0.00301532,
                 P7 = 0.00032411;
  const G4double Q1 = 0.39894228,
                 Q2 =-0.03988024,
                 Q3 =-0.00362018,
                 Q4 = 0.00163801,
                 Q5 =-0.01031555,
                 Q6 = 0.02282967,
                 Q7 =-0.02895312,
                 Q8 = 0.01787654,
                 Q9 =-0.00420059;
 
  G4double I = 0.0;
  if (std::abs(x) < 3.75)
  {
    G4double ax = std::abs(x);
    G4double y = G4Pow::GetInstance()->powN(x/3.75, 2);
    I = ax*(P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7))))));
  }
  else
  {
    G4double ax = std::abs(x);
    G4double y  = 3.75/ax;
    I  = G4Exp(ax) / std::sqrt(ax) *
      (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*(Q7+y*(Q8+y*Q9))))))));
  }
  if (x < 0.0) I = -I;
  return I;
}

Referenced by K1().

◆ K0()

G4double G4Bessel::K0 ( G4double x )

Definition at line 116 of file G4Bessel.cc.

{
  const G4double P1 =-0.57721566,
                 P2 = 0.42278420,
                 P3 = 0.23069756,
                 P4 = 0.03488590,
                 P5 = 0.00262698,
                 P6 = 0.00010750,
                 P7 = 0.00000740;
  const G4double Q1 = 1.25331414,
                 Q2 =-0.07832358,
                 Q3 = 0.02189568,
                 Q4 =-0.01062446,
                 Q5 = 0.00587872,
                 Q6 =-0.00251540,
                 Q7 = 0.00053208;
 
  G4double K = 0.0;
  if (x <= 2.0)
  {
    G4double y = x * x / 4.0;
    K = (-G4Log(x/2.0)) * I0(x) +
      P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7)))));
  }
  else
  {
    G4double y = 2.0 / x;
    K = G4Exp(-x)  / std::sqrt(x) *
      (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*Q7))))));
  }
  return K;
}

Referenced by G4EMDissociationSpectrum::GetGeneralE1Spectrum(), and G4EMDissociationSpectrum::GetGeneralE2Spectrum().

◆ K1()

G4double G4Bessel::K1 ( G4double x )

Definition at line 188 of file G4Bessel.cc.

{
  const G4double P1 = 1.0,
                 P2 = 0.15443144,
                 P3 =-0.67278579,
                 P4 =-0.18156897,
                 P5 =-0.01919402,
                 P6 =-0.00110404,
                 P7 =-0.00004686;
  const G4double Q1 = 1.25331414,
                 Q2 = 0.23498619,
                 Q3 =-0.03655620,
                 Q4 = 0.01504268,
                 Q5 =-0.00780353,
                 Q6 = 0.00325614,
                 Q7 =-0.00068245;
 
  G4double K = 0.0;
  if (x <= 2.0)
  {
    G4double y = x * x / 4.0;
    K = G4Log(x/2.0)*I1(x) + 1.0/x *
      (P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7))))));
  }
  else
  {
    G4double y = 2.0 / x;
    K = G4Exp(-x) / std::sqrt(x) *
      (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*Q7))))));
  }
  return K;
}

Referenced by G4EMDissociationSpectrum::GetGeneralE1Spectrum(), and G4EMDissociationSpectrum::GetGeneralE2Spectrum().

◆ pI0()

G4double G4Bessel::pI0 ( G4double x )

Definition at line 222 of file G4Bessel.cc.

{
  const G4double A0  = 0.1250000000000E+00,
                 A1  = 7.0312500000000E-02,
                 A2  = 7.3242187500000E-02,
                 A3  = 1.1215209960938E-01,
                 A4  = 2.2710800170898E-01,
                 A5  = 5.7250142097473E-01,
                 A6  = 1.7277275025845E+00,
                 A7  = 6.0740420012735E+00,
                 A8  = 2.4380529699556E+01,
                 A9  = 1.1001714026925E+02,
                 A10 = 5.5133589612202E+02,
                 A11 = 3.0380905109224E+03;
 
  G4double I = 0.0;
  if (x == 0.0)
  {
    I = 1.0;
  }
  else if (x < 18.0)
  {
    I          = 1.0;
    G4double y = x * x;
    G4double q = 1.0;
    for (G4int i=1; i<101; i++)
    {
      q *= 0.25 * y / i / i;
      I += q;
      if (std::abs(q/I) < 1.0E-15) break;
    }
  }
  else
  {
    G4double y = 1.0 / x;
    I = G4Exp(x) / std::sqrt(twopi*x) *
      (1.0 + y*(A0+y*(A1+y*(A2+y*(A3+y*(A4+y*(A5+y*(A6+y*(A7+y*(A8+y*(A9+y*(A10+y*A11))))))))))));
  }
 
  return I;
}

Referenced by pK0(), and pK1().

◆ pI1()

G4double G4Bessel::pI1 ( G4double x )

Definition at line 265 of file G4Bessel.cc.

{
  const G4double A0  = -0.3750000000000E+00,
                 A1  = -1.1718750000000E-01,
                 A2  = -1.0253906250000E-01,
                 A3  = -1.4419555664063E-01,
                 A4  = -2.775764465332E-01,
                 A5  = -6.7659258842468E-01,
                 A6  = -1.9935317337513E+00,
                 A7  = -6.8839142681099E+00,
                 A8  = -2.7248827311269E+01,
                 A9  = -1.2159789187654E+02,
                 A10 = -6.0384407670507E+02,
                 A11 = -3.3022722944809E+03;
 
  G4double I = 0.0;
  if (x == 0.0)
  {
    I = 0.0;
  }
  else if (x < 18.0)
  {
    I          = 1.0;
    G4double y = x * x;
    G4double q = 1.0;
    for (G4int i=1; i<101; i++)
    {
      q *= 0.25 * y / i / (i+1.0);
      I += q;
      if (std::abs(q/I) < 1.0E-15) break;
    }
    I *= 0.5 * x;
    
  }
  else
  {
    G4double y = 1.0 / x;
    I = G4Exp(x) / std::sqrt(twopi*x) *
      (1.0 + y*(A0+y*(A1+y*(A2+y*(A3+y*(A4+y*(A5+y*(A6+y*(A7+y*(A8+y*(A9+y*(A10+y*A11))))))))))));
  }
 
  return I;
}

Referenced by pK1().

◆ pK0()

G4double G4Bessel::pK0 ( G4double x )

Definition at line 310 of file G4Bessel.cc.

{
  const G4double A0 = 0.1250000000000E+00,
                 A1 = 0.2109375000000E+00,
                 A2 = 1.0986328125000E+00,
                 A3 = 1.1775970458984E+01,
                 A4 = 2.1461706161499E+02,
                 A5 = 5.9511522710323E+03,
                 A6 = 2.3347645606175E+05,
                 A7 = 1.2312234987631E+07;
 
  G4double K = 0.0;
  if (x == 0.0)
  {
    K = 1.0E+307;
  }
  else if (x < 9.0)
  {
    G4double y = x * x;
    G4double C = -G4Log(x/2.0) - 0.5772156649015329;
    G4double q = 1.0;
    G4double t = 0.0;
    for (G4int i=1; i<51; i++)
    {
      q *= 0.25 * y / i / i;
      t += 1.0 / i ;
      K += q * (t+C);
    }
    K += C;
  }
  else
  {
    G4double y = 1.0 / x / x;
    K = 0.5 / x / pI0(x) *
      (1.0 + y*(A0+y*(A1+y*(A2+y*(A3+y*(A4+y*(A5+y*(A6+y*A7))))))));
  }
  
  return K;
}