zlnrel.f 1.04 KB
complex(kind(1.d0)) function zlnrel (z)
      implicit none
c april 1977 version.  w. fullerton, c3, los alamos scientific lab.
c
c clnrel(z) = clog(1+z) with relative error accuracy near z = 0.
c let   rho = cabs(z)  and
c       r**2 = cabs(1+z)**2 = (1+x)**2 + y**2 = 1 + 2*x + rho**2 .
c now if rho is small we may evaluate clnrel(z) accurately by
c       clog(1+z) = cmplx (alog(r), carg(1+z))
c                 = cmplx (0.5*alog(r**2), carg(1+z))
c                 = cmplx (0.5*alnrel(2*x+rho**2), carg(1+z))
c
      complex(kind(1.d0)) z
      real(kind(1.d0)) dlnrel,d1mach,zarg,sqeps,x,rho
      external dlnrel, zarg, d1mach
      data sqeps /0.0d0/
c
      if (sqeps.eq.0.) sqeps = sqrt (d1mach(4))
c
      if ( abs(1.+z).lt.sqeps) call seteru (
     1  54hzlnrel  answer lt half precision because z too near -1, 54,
     2  1, 1)
c
      rho = abs(z)
      if (rho.gt.0.375) zlnrel = log (1.0+z)
      if (rho.gt.0.375) return
c
      x = real(z)
      zlnrel = cmplx (0.5*dlnrel(2.*x+rho**2), zarg(1.0+z), kind(1.d0))
c
      return
      end