| Geant4 Cross Reference |
1 /********************************************* 1
2 Returns the exponential integral function
3
4 E_n(x) = int_1^infinity e^( -x * t ) / t^n
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
21 #define MAXIT 100 /* Maximum
22 #define FPMIN 1.0e-300 /* close t
23 #define EPS 1.0e-15 /* Desired
24
25 /*
26 **********************************************
27 */
28 double nf_exponentialIntegral( int n, double x
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
40 *status = nfu_badInput; }
41 else {
42 if( n == 0 ) {
43 ans = G4Exp( -x ) / x; }
44 else if( x == 0.0 ) {
45 ans = 1.0 / nm1; }
46 else if( x > 1.0 ) {
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 );
55 c = b + a / c;
56 del = c * d;
57 h *= del;
58 if( fabs( del - 1.0 ) < EPS )
59 }
60 *status = nfu_failedToConverge; }
61 else {
62 ans = ( nm1 != 0 ) ? 1.0 / nm1 : -
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;
70 for( ii = 1; ii <= nm1; ii
71 del = fact * ( -G4Log( x )
72 }
73 ans += del;
74 if( fabs( del ) < fabs( ans )
75 }
76 *status = nfu_failedToConverge;
77 }
78 }
79 return( ans );
80 }
81
82 #if defined __cplusplus
83 }
84 #endif
85