! ! Copyright (C) 2001-2008 Quantum ESPRESSO group ! This file is distributed under the terms of the ! GNU General Public License. See the file `License' ! in the root directory of the present distribution, ! or http://www.gnu.org/copyleft/gpl.txt . ! !----------------------------------------------------------------------- FUNCTION EXPINT(n, x) !----------------------------------------------------------------------- ! ! Evaluates the exponential integral E_n(x) ! Parameters: maxit is the maximum allowed number of iterations, ! eps is the desired relative error, not smaller than the machine precision, ! big is a number near the largest representable floating-point number, ! Inspired from Numerical Recipes ! USE kinds, ONLY : DP IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL(DP), INTENT(IN) :: x REAL(DP) :: expint INTEGER, parameter :: maxit=200 REAL(DP), parameter :: eps=1E-12_DP, big=huge(x)*eps REAL(DP), parameter :: euler = 0.577215664901532860606512_DP ! EPS=1E-9, FPMIN=1E-30 INTEGER :: i, nm1, k REAL(DP) :: a,b,c,d,del,fact,h,iarsum IF (.NOT. ((n >= 0).AND.(x >= 0.0).AND.((x > 0.0).OR.(n > 1)))) THEN CALL errore('expint','bad arguments', 1) END IF IF (n == 0) THEN expint = exp(-x)/x RETURN END IF nm1 = n-1 IF (x == 0.0_DP) THEN expint = 1.0_DP/nm1 ELSE IF (x > 1.0_DP) THEN b = x+n c = big d = 1.0_DP/b h = d DO i=1,maxit a = -i*(nm1+i) b = b+2.0_DP d = 1.0_DP/(a*d+b) c = b+a/c del = c*d h = h*del IF (ABS(del-1.0_DP) <= EPS) EXIT END DO IF (i > maxit) CALL errore('expint','continued fraction failed',1) expint = h*EXP(-x) ELSE IF (nm1 /= 0) THEN expint = 1.0_DP/nm1 ELSE expint = -LOG(x)-euler END IF fact = 1.0_DP do i=1,maxit fact = -fact*x/i IF (i /= nm1) THEN del = -fact/(i-nm1) ELSE iarsum = 0.0_DP do k=1,nm1 iarsum = iarsum + 1.0_DP/k end do del = fact*(-LOG(x)-euler+iarsum) ! del = fact*(-LOG(x)-euler+sum(1.0_DP/arth(1,1,nm1))) END IF expint = expint+del IF (ABS(del) < ABS(expint)*eps) EXIT END DO IF (i > maxit) CALL errore('expint','series failed',1) END IF END FUNCTION EXPINT ! -------------------------------------------------------------------