expint.f90 2.62 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
!
! 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
! 
dalcorso's avatar
dalcorso committed
18
      USE kinds, ONLY : DP
19 20 21 22
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: n
      REAL(DP), INTENT(IN) :: x
      REAL(DP) :: expint
23 24
      INTEGER, parameter :: maxit=200
      REAL(DP), parameter :: eps=1E-12_DP, big=huge(x)*eps
25 26 27 28 29 30 31
      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
giannozz's avatar
giannozz committed
32
         CALL errore('expint','bad arguments', 1)
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
      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

! -------------------------------------------------------------------