nf_exponentialIntegral.cc File Reference

#include "nf_specialFunctions.h"

Go to the source code of this file.

Macros
#define	EULER   0.57721566490153286 /* Euler's constant gamma */
#define	MAXIT   100 /* Maximum allowed number of iterations. */
#define	FPMIN   1.0e-300 /* close to the smallest representable floting-point number. */
#define	EPS   1.0e-15 /* Desired relative error, not smaller than the machine precision. */

Functions
double	nf_exponentialIntegral (int n, double x, nfu_status *status)

Macro Definition Documentation

EPS

#define EPS   1.0e-15 /* Desired relative error, not smaller than the machine precision. */

Definition at line 23 of file nf_exponentialIntegral.cc.

EULER

#define EULER   0.57721566490153286 /* Euler's constant gamma */

Definition at line 20 of file nf_exponentialIntegral.cc.

FPMIN

#define FPMIN   1.0e-300 /* close to the smallest representable floting-point number. */

Definition at line 22 of file nf_exponentialIntegral.cc.

MAXIT

#define MAXIT   100 /* Maximum allowed number of iterations. */

Definition at line 21 of file nf_exponentialIntegral.cc.

Function Documentation

nf_exponentialIntegral()

double nf_exponentialIntegral	(	int	n,
		double	x,
		nfu_status *	status
	)

Definition at line 28 of file nf_exponentialIntegral.cc.

                                                                     {
 
    int i, ii, nm1;
    double a, b, c, d, del, fact, h, psi;
    double ans = 0.0;
 
    *status = nfu_badInput;
    if( !isfinite( x ) ) return( x );
    *status = nfu_Okay;
 
    nm1 = n - 1;
    if( ( n < 0 ) || ( x < 0.0 ) || ( ( x == 0.0 ) && ( ( n == 0 ) || ( n == 1 ) ) ) ) {
        *status = nfu_badInput; }
    else {
        if( n == 0 ) {
            ans = G4Exp( -x ) / x; }                  /* Special case */
        else if( x == 0.0 ) {
            ans = 1.0 / nm1; }                      /* Another special case */
        else if( x > 1.0 ) {                        /* Lentz's algorithm */
            b = x + n;
            c = 1.0 / FPMIN;
            d = 1.0 / b;
            h = d;
            for( i = 1; i <= MAXIT; i++ ) {
                a = -i * ( nm1 + i );
                b += 2.0;
                d = 1.0 / ( a * d + b );            /* Denominators cannot be zero */
                c = b + a / c;
                del = c * d;
                h *= del;
                if( fabs( del - 1.0 ) < EPS ) return( h * G4Exp( -x ) );
            }
            *status = nfu_failedToConverge; }
        else {
            ans = ( nm1 != 0 ) ? 1.0 / nm1 : -G4Log(x) - EULER;   /* Set first term */
            fact = 1.0;
            for( i = 1; i <= MAXIT; i++ ) {
                fact *= -x / i;
                if( i != nm1 ) {
                    del = -fact / ( i - nm1 ); }
                else {
                    psi = -EULER;                   /* Compute psi(n) */
                    for( ii = 1; ii <= nm1; ii++ ) psi += 1.0 / ii;
                    del = fact * ( -G4Log( x ) + psi );
                }
                ans += del;
                if( fabs( del ) < fabs( ans ) * EPS ) return( ans );
            }
            *status = nfu_failedToConverge;
        }
    }
    return( ans );
}