1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
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
|
!
! 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
! -------------------------------------------------------------------
|
