#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 65 of file G4Bessel.cc.

66{;}

◆ ~G4Bessel()

G4Bessel::~G4Bessel ( )

Definition at line 69 of file G4Bessel.cc.

70{;}

Member Function Documentation

◆ I0()

G4double G4Bessel::I0 ( G4double x )

Definition at line 73 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::fabs(x) < 3.75)
  {
    G4double y = std::pow(x/3.75, 2.0);
    I = P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7)))));
  }
  else
  {
    G4double ax = std::fabs(x);
    G4double y  = 3.75/ax;
    I  = std::exp(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 143 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::fabs(x) < 3.75)
  {
    G4double ax = std::fabs(x);
    G4double y = std::pow(x/3.75, 2.0);
    I = ax*(P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7))))));
  }
  else
  {
    G4double ax = std::fabs(x);
    G4double y  = 3.75/ax;
    I  = std::exp(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 109 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 = (-std::log(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 = std::exp(-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 181 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 = std::log(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 = std::exp(-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 215 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 = std::exp(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 258 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 = std::exp(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 303 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 = -std::log(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;
}