zasin.f 1.31 KB
complex(kind(1.d0)) 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(kind(1.d0)) zinp, z, z2, sqzp1, ci
      real(kind(1.d0)) d1mach,pi2,pi,rmin,r
      integer nterms,i,twoi
      external d1mach

      data pi2 /1.5707963267 9489661923d0/
      data pi /3.1415926535 8979324d0/
      data ci /(0.d0,1.d0)/
      data nterms, rmin / 0, 0.0d0 /
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