-
git-svn-id: https://lxsd.femto-st.fr/svn/fvn@69 b657c933-2333-4658-acf2-d3c7c2708721
zexprl.f
1.33 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
complex(kind(1.d0)) function zexprl (z)
implicit none
c august 1980 edition. w. fullerton, c3, los alamos scientific lab.
c
c evaluate (cexp(z)-1)/z . for small cabs(z), we use the taylor
c series. we could instead use the expression
c cexprl(z) = (exp(x)*exp(i*y)-1)/z
c = (x*exprel(x) * (1 - 2*sin(y/2)**2) - 2*sin(y/2)**2
c + i*sin(y)*(1+x*exprel(x))) / z
c
complex(kind(1.d0)) z
real(kind(1.d0)) sqeps,r,xn,xln,rbnd,d1mach,alneps
external d1mach
integer nterms,irold,irold2,i
data nterms, rbnd, sqeps / 0, 2*0.0 /
c
if (nterms.ne.0) go to 10
alneps = log(d1mach(3))
xn = 3.72 - 0.3*alneps
xln = log((xn+1.0)/1.36)
nterms = xn - (xn*xln+alneps)/(xln+1.36) + 1.5
rbnd = d1mach(3)
sqeps = 0.0
c
10 r = abs(z)
if (r.le.0.5) go to 20
c
call entsrc (irold, 1)
zexprl = (exp(z) - 1.0) / z
call erroff
call entsrc (irold2, irold)
if (abs(real(z)).lt.sqeps .and. abs(zexprl).lt.sqeps) call
1 seteru (41hzexprl answer lt half prec, z near ni2pi,
2 41, 1, 1)
return
c
20 zexprl = (1.0, 0.0)
if (r.lt.rbnd) return
c
zexprl = (0.0, 0.0)
do 30 i=1,nterms
zexprl = 1.0 + zexprl*z/dble(nterms+2-i)
30 continue
c
return
end