-
git-svn-id: https://lxsd.femto-st.fr/svn/fvn@69 b657c933-2333-4658-acf2-d3c7c2708721
zatan.f
1.63 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
53
54
55
56
57
58
59
complex(kind(1.d0)) function zatan (z)
implicit none
c jan 1984 edition. w. fullerton, los alamos scientific lab.
c
complex(kind(1.d0)) z, z2
real(kind(1.d0)) pi2,sqeps,rmin,rmax,r2sml,r,x,y,r2,xans,yans
real(kind(1.d0)) d1mach
integer nterms,i,twoi
data pi2 /1.5707963267 9489661923d0 /
data nterms, sqeps, rmin, rmax, r2sml / 0, 4*0.0d0 /
external d1mach
c
if (nterms.ne.0) go to 10
c nterms = alog(eps)/alog(rbnd) where rbnd = 0.1
nterms = -0.4343*log(d1mach(3)) + 1.0
sqeps = sqrt (d1mach(4))
rmin = sqrt (3.0*d1mach(3))
rmax = 1.0/d1mach(3)
r2sml = 10.0*d1mach (4)
c
10 r = abs(z)
if (r.gt.0.1) go to 30
c
zatan = z
if (r.lt.rmin) return
c
zatan = (0.0, 0.0)
z2 = z*z
do 20 i=1,nterms
twoi = 2*(nterms-i) + 1
zatan = 1.0/twoi - z2*zatan
20 continue
zatan = z*zatan
return
c
30 if (r.gt.rmax) go to 50
x = real(z)
y = aimag(z)
r2 = r*r
c
if (abs(r2-1.0).gt.sqeps) go to 40
if (r2+2.0*y+1.0.lt.r2sml .or. r2-2.0*y+1.0.lt.r2sml) call
1 seteru (
2 55hzatan no precision because z is too close to +i or -i,
3 55, 2, 2)
if (abs(cmplx(1.0d0,0.0d0,kind(1.d0))+z*z).lt.sqeps) call seteru (
1 50hzatan answer lt half precision, z**2 close to -1, 50, 1, 1)
c
40 xans = 0.5*atan2 (2.0*x, 1.0-r2)
yans = 0.25*log((r2+2.0*y+1.0)/(r2-2.0*y+1.0))
zatan = cmplx (xans, yans, kind(1.d0))
return
c
50 zatan = cmplx (pi2, 0.0d0, kind(1.d0))
if (real(z).lt.0.0) zatan = cmplx (-pi2, 0.0d0, kind(1.d0))
return
c
end