ze1.f90 1.31 KB
edit raw blame history



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


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