casin.f
1.21 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 function casin (zinp)
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 zinp, z, z2, sqzp1, ci, clog, csqrt
external r1mach
data pi2 /1.5707963267 9489661923e0/
data pi /3.1415926535 8979324e0/
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*alog(r1mach(3)) + 1.0
rmin = sqrt (6.0*r1mach(3))
c
10 z = zinp
r = cabs (z)
if (r.gt.0.1) go to 30
c
casin = z
if (r.lt.rmin) return
c
casin = (0.0, 0.0)
z2 = z*z
do 20 i=1,nterms
twoi = 2*(nterms-i) + 1
casin = 1.0/twoi + twoi*casin*z2/(twoi+1.0)
20 continue
casin = z*casin
return
c
30 if (real(zinp).lt.0.0) z = -zinp
c
sqzp1 = csqrt (z+1.0)
if (aimag(sqzp1).lt.0.) sqzp1 = -sqzp1
casin = pi2 - ci * clog (z + sqzp1*csqrt(z-1.0))
c
if (real(casin).gt.pi2) casin = pi - casin
if (real(casin).le.(-pi2)) casin = -pi - casin
if (real(zinp).lt.0.) casin = -casin
c
return
end