Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lend/src/nf_exponentialIntegral.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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