-
git-svn-id: https://lxsd.femto-st.fr/svn/fvn@56 b657c933-2333-4658-acf2-d3c7c2708721
ze1.f90
1.31 KB
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
function ze1(z)
!
! ====================================================
! Purpose: Compute complex exponential integral E1(z)
! Input : z --- Argument of E1(z)
! Output: CE1 --- E1(z)
! ====================================================
!
! Déclaration des variables en passage de paramètre
!
implicit none
complex(8), intent(in) :: z
complex(8) :: ze1
!
! Déclaration des variables locales
!
integer(4) :: k
real(8) :: pi,el,x,a0
complex(8) :: cr,ct0,ct
parameter(pi=3.141592653589793D0,el=0.5772156649015328D0)
!
! traitement en fonction des différents cas
! - Z nul entraîne E1 infini
! - module de Z inférieur à 10 ou 20 : formule log+gam+somme
! - module de Z supérieur à 10 ou 20 : formule asymptotique
!
x=dreal(z)
a0=cdabs(z)
if (a0==0.0D0) then
ze1 = dcmplx(1.0D+300,0.0D0)
else if ((a0 <= 10.D0).or.(x <= 0.D0.and.a0 <= 20.D0)) then
ze1 = dcmplx(1.0D0,0.0D0)
cr = dcmplx(1.0D0,0.0D0)
k=0
do while (cdabs(cr)>=cdabs(ze1)*1.0D-15)
k = k+1
cr = -cr*k*z/(k+1.0D0)**2
ze1 = ze1+cr
end do
ze1 = -el-cdlog(z)+z*ze1
else
ct0 = dcmplx(0.0D0,0.0D0)
do k=120,1,-1
ct0 = k/(1.0D0+k/(z+ct0))
end do
ct = 1.0D0/(z+ct0)
ze1 = cdexp(-z)*ct
if (x <= 0.D0.AND.dimag(z) == 0.0) ze1 = ze1-pi*dcmplx(0.D0,1.D0)
end if
!
return
end function