zasin.f
1.28 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
complex(8) function zasin (zinp)
implicit none
c august 1980 edition. w. fullerton, c3, los alamos scientific lab.
c
c ref -- l. l. pennisi, elements of complex variables, holt, rinehart
c and winston, 1963. page 126.
c
complex(8) zinp, z, z2, sqzp1, ci
real(8) d1mach,pi2,pi,rmin,r
integer nterms,i,twoi
external d1mach
data pi2 /1.5707963267 9489661923d0/
data pi /3.1415926535 8979324d0/
data ci /(0.,1.)/
data nterms, rmin / 0, 0.0 /
c
if (nterms.ne.0) go to 10
c nterms = alog(eps)/alog(rmax) where rmax = 0.1
nterms = -0.4343*log(d1mach(3)) + 1.0
rmin = sqrt (6.0*d1mach(3))
c
10 z = zinp
r = abs (z)
if (r.gt.0.1) go to 30
c
zasin = z
if (r.lt.rmin) return
c
zasin = (0.0, 0.0)
z2 = z*z
do 20 i=1,nterms
twoi = 2*(nterms-i) + 1
zasin = 1.0/twoi + twoi*zasin*z2/(twoi+1.0)
20 continue
zasin = z*zasin
return
c
30 if (real(zinp).lt.0.0) z = -zinp
c
sqzp1 = sqrt (z+1.0)
if (aimag(sqzp1).lt.0.) sqzp1 = -sqzp1
zasin = pi2 - ci * log (z + sqzp1*sqrt(z-1.0))
c
if (real(zasin).gt.pi2) zasin = pi - zasin
if (real(zasin).le.(-pi2)) zasin = -pi - zasin
if (real(zinp).lt.0.) zasin = -zasin
c
return
end