lxsd / fvn

Blame view

fvn_fnlib/zasin.f 1.28 KB
edit raw normal view history
38581db0c   daniau   git-svn-id: https... 
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