| Geant4 Cross Reference |
1 /********************************************* 1 /*********************************************************************
2 Returns the exponential integral function 2 Returns the exponential integral function
3 3
4 E_n(x) = int_1^infinity e^( -x * t ) / t^n 4 E_n(x) = int_1^infinity e^( -x * t ) / t^n dt, for x > 0.
5 5
6 C.A. Bertulani May/15/2000 6 C.A. Bertulani May/15/2000
7 ********************************************** 7 *********************************************************************/
8 8
9 #include "nf_specialFunctions.h" 9 #include "nf_specialFunctions.h"
10 10
11 #if defined __cplusplus 11 #if defined __cplusplus
12 #include <cmath> 12 #include <cmath>
13 #include "G4Exp.hh" 13 #include "G4Exp.hh"
14 #include "G4Log.hh" 14 #include "G4Log.hh"
15 namespace GIDI { 15 namespace GIDI {
16 using namespace GIDI; 16 using namespace GIDI;
17 using namespace std; 17 using namespace std;
18 #endif 18 #endif
19 19
20 #define EULER 0.57721566490153286 /* Euler's 20 #define EULER 0.57721566490153286 /* Euler's constant gamma */
21 #define MAXIT 100 /* Maximum 21 #define MAXIT 100 /* Maximum allowed number of iterations. */
22 #define FPMIN 1.0e-300 /* close t 22 #define FPMIN 1.0e-300 /* close to the smallest representable floting-point number. */
23 #define EPS 1.0e-15 /* Desired 23 #define EPS 1.0e-15 /* Desired relative error, not smaller than the machine precision. */
24 24
25 /* 25 /*
26 ********************************************** 26 ************************************************************
27 */ 27 */
28 double nf_exponentialIntegral( int n, double x 28 double nf_exponentialIntegral( int n, double x, nfu_status *status ) {
29 29
30 int i, ii, nm1; 30 int i, ii, nm1;
31 double a, b, c, d, del, fact, h, psi; 31 double a, b, c, d, del, fact, h, psi;
32 double ans = 0.0; 32 double ans = 0.0;
33 33
34 *status = nfu_badInput; 34 *status = nfu_badInput;
35 if( !isfinite( x ) ) return( x ); 35 if( !isfinite( x ) ) return( x );
36 *status = nfu_Okay; 36 *status = nfu_Okay;
37 37
38 nm1 = n - 1; 38 nm1 = n - 1;
39 if( ( n < 0 ) || ( x < 0.0 ) || ( ( x == 0 39 if( ( n < 0 ) || ( x < 0.0 ) || ( ( x == 0.0 ) && ( ( n == 0 ) || ( n == 1 ) ) ) ) {
40 *status = nfu_badInput; } 40 *status = nfu_badInput; }
41 else { 41 else {
42 if( n == 0 ) { 42 if( n == 0 ) {
43 ans = G4Exp( -x ) / x; } 43 ans = G4Exp( -x ) / x; } /* Special case */
44 else if( x == 0.0 ) { 44 else if( x == 0.0 ) {
45 ans = 1.0 / nm1; } 45 ans = 1.0 / nm1; } /* Another special case */
46 else if( x > 1.0 ) { 46 else if( x > 1.0 ) { /* Lentz's algorithm */
47 b = x + n; 47 b = x + n;
48 c = 1.0 / FPMIN; 48 c = 1.0 / FPMIN;
49 d = 1.0 / b; 49 d = 1.0 / b;
50 h = d; 50 h = d;
51 for( i = 1; i <= MAXIT; i++ ) { 51 for( i = 1; i <= MAXIT; i++ ) {
52 a = -i * ( nm1 + i ); 52 a = -i * ( nm1 + i );
53 b += 2.0; 53 b += 2.0;
54 d = 1.0 / ( a * d + b ); 54 d = 1.0 / ( a * d + b ); /* Denominators cannot be zero */
55 c = b + a / c; 55 c = b + a / c;
56 del = c * d; 56 del = c * d;
57 h *= del; 57 h *= del;
58 if( fabs( del - 1.0 ) < EPS ) 58 if( fabs( del - 1.0 ) < EPS ) return( h * G4Exp( -x ) );
59 } 59 }
60 *status = nfu_failedToConverge; } 60 *status = nfu_failedToConverge; }
61 else { 61 else {
62 ans = ( nm1 != 0 ) ? 1.0 / nm1 : - 62 ans = ( nm1 != 0 ) ? 1.0 / nm1 : -G4Log(x) - EULER; /* Set first term */
63 fact = 1.0; 63 fact = 1.0;
64 for( i = 1; i <= MAXIT; i++ ) { 64 for( i = 1; i <= MAXIT; i++ ) {
65 fact *= -x / i; 65 fact *= -x / i;
66 if( i != nm1 ) { 66 if( i != nm1 ) {
67 del = -fact / ( i - nm1 ); 67 del = -fact / ( i - nm1 ); }
68 else { 68 else {
69 psi = -EULER; 69 psi = -EULER; /* Compute psi(n) */
70 for( ii = 1; ii <= nm1; ii 70 for( ii = 1; ii <= nm1; ii++ ) psi += 1.0 / ii;
71 del = fact * ( -G4Log( x ) 71 del = fact * ( -G4Log( x ) + psi );
72 } 72 }
73 ans += del; 73 ans += del;
74 if( fabs( del ) < fabs( ans ) 74 if( fabs( del ) < fabs( ans ) * EPS ) return( ans );
75 } 75 }
76 *status = nfu_failedToConverge; 76 *status = nfu_failedToConverge;
77 } 77 }
78 } 78 }
79 return( ans ); 79 return( ans );
80 } 80 }
81 81
82 #if defined __cplusplus 82 #if defined __cplusplus
83 } 83 }
84 #endif 84 #endif
85 85