Geant4  10.02
nf_exponentialIntegral.cc
Go to the documentation of this file.
1 /*********************************************************************
2  Returns the exponential integral function
3 
4  E_n(x) = int_1^infinity e^( -x * t ) / t^n dt, for x > 0.
5 
6  C.A. Bertulani May/15/2000
7 *********************************************************************/
8 
9 #include "nf_specialFunctions.h"
10 
11 #if defined __cplusplus
12 #include <cmath>
13 #include "G4Exp.hh"
14 #include "G4Log.hh"
15 namespace GIDI {
16 using namespace GIDI;
17 using namespace std;
18 #endif
19 
20 #define EULER 0.57721566490153286 /* Euler's constant gamma */
21 #define MAXIT 100 /* Maximum allowed number of iterations. */
22 #define FPMIN 1.0e-300 /* close to the smallest representable floting-point number. */
23 #define EPS 1.0e-15 /* Desired relative error, not smaller than the machine precision. */
24 
25 /*
26 ************************************************************
27 */
28 double nf_exponentialIntegral( int n, double x, nfu_status *status ) {
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 }
81 
82 #if defined __cplusplus
83 }
84 #endif
const G4double fact
#define FPMIN
G4double a
Definition: TRTMaterials.hh:39
#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
const G4double x[NPOINTSGL]
#define EPS
double nf_exponentialIntegral(int n, double x, nfu_status *status)