Geant4 Cross Reference |
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 85