Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
nf_exponentialIntegral.cc File Reference
Include dependency graph for nf_exponentialIntegral.cc:

Go to the source code of this file.

Macros

#define EULER   0.57721566490153286 /* Euler's constant gamma */
 
#define MAXIT   100 /* Maximum allowed number of iterations. */
 
#define FPMIN   1.0e-300 /* close to the smallest representable floting-point number. */
 
#define EPS   1.0e-15 /* Desired relative error, not smaller than the machine precision. */
 

Functions

double nf_exponentialIntegral (int n, double x, nfu_status *status)
 

Macro Definition Documentation

#define EPS   1.0e-15 /* Desired relative error, not smaller than the machine precision. */

Definition at line 23 of file nf_exponentialIntegral.cc.

#define EULER   0.57721566490153286 /* Euler's constant gamma */

Definition at line 20 of file nf_exponentialIntegral.cc.

#define FPMIN   1.0e-300 /* close to the smallest representable floting-point number. */

Definition at line 22 of file nf_exponentialIntegral.cc.

#define MAXIT   100 /* Maximum allowed number of iterations. */

Definition at line 21 of file nf_exponentialIntegral.cc.

Function Documentation

double nf_exponentialIntegral ( int  n,
double  x,
nfu_status status 
)

Definition at line 28 of file nf_exponentialIntegral.cc.

28  {
29 
30  int i, ii, nm1;
31  double a, b, c, d, del, fact, h, psi;
32  double ans = 0.0;
33 
34  *status = nfu_badInput;
35  if( !isfinite( x ) ) return( x );
36  *status = nfu_Okay;
37 
38  nm1 = n - 1;
39  if( ( n < 0 ) || ( x < 0.0 ) || ( ( x == 0.0 ) && ( ( n == 0 ) || ( n == 1 ) ) ) ) {
40  *status = nfu_badInput; }
41  else {
42  if( n == 0 ) {
43  ans = G4Exp( -x ) / x; } /* Special case */
44  else if( x == 0.0 ) {
45  ans = 1.0 / nm1; } /* Another special case */
46  else if( x > 1.0 ) { /* Lentz's algorithm */
47  b = x + n;
48  c = 1.0 / FPMIN;
49  d = 1.0 / b;
50  h = d;
51  for( i = 1; i <= MAXIT; i++ ) {
52  a = -i * ( nm1 + i );
53  b += 2.0;
54  d = 1.0 / ( a * d + b ); /* Denominators cannot be zero */
55  c = b + a / c;
56  del = c * d;
57  h *= del;
58  if( fabs( del - 1.0 ) < EPS ) return( h * G4Exp( -x ) );
59  }
60  *status = nfu_failedToConverge; }
61  else {
62  ans = ( nm1 != 0 ) ? 1.0 / nm1 : -G4Log(x) - EULER; /* Set first term */
63  fact = 1.0;
64  for( i = 1; i <= MAXIT; i++ ) {
65  fact *= -x / i;
66  if( i != nm1 ) {
67  del = -fact / ( i - nm1 ); }
68  else {
69  psi = -EULER; /* Compute psi(n) */
70  for( ii = 1; ii <= nm1; ii++ ) psi += 1.0 / ii;
71  del = fact * ( -G4Log( x ) + psi );
72  }
73  ans += del;
74  if( fabs( del ) < fabs( ans ) * EPS ) return( ans );
75  }
76  *status = nfu_failedToConverge;
77  }
78  }
79  return( ans );
80 }
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
#define FPMIN
tuple x
Definition: test.py:50
tuple b
Definition: test.py:12
#define MAXIT
#define isfinite
const G4int n
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
#define EULER
#define EPS
tuple c
Definition: test.py:13

Here is the call graph for this function:

Here is the caller graph for this function: