complex(kind(1.d0)) function zatan (z)

        implicit none
  c jan 1984 edition.  w. fullerton, los alamos scientific lab.
  c

        complex(kind(1.d0)) z, z2
        real(kind(1.d0)) pi2,sqeps,rmin,rmax,r2sml,r,x,y,r2,xans,yans
        real(kind(1.d0)) d1mach

        integer nterms,i,twoi
        data pi2 /1.5707963267 9489661923d0 /

        data nterms, sqeps, rmin, rmax, r2sml / 0, 4*0.0d0 /
        external d1mach

  c
        if (nterms.ne.0) go to 10
  c nterms = alog(eps)/alog(rbnd) where rbnd = 0.1
        nterms = -0.4343*log(d1mach(3)) + 1.0
        sqeps = sqrt (d1mach(4))
        rmin = sqrt (3.0*d1mach(3))
        rmax = 1.0/d1mach(3)
        r2sml = 10.0*d1mach (4)
  c
   10   r = abs(z)
        if (r.gt.0.1) go to 30
  c
        zatan = z
        if (r.lt.rmin) return
  c
        zatan = (0.0, 0.0)
        z2 = z*z
        do 20 i=1,nterms
          twoi = 2*(nterms-i) + 1
          zatan = 1.0/twoi - z2*zatan
   20   continue
        zatan = z*zatan
        return
  c
   30   if (r.gt.rmax) go to 50
        x = real(z)
        y = aimag(z)
        r2 = r*r
  c
        if (abs(r2-1.0).gt.sqeps) go to 40
        if (r2+2.0*y+1.0.lt.r2sml .or. r2-2.0*y+1.0.lt.r2sml) call
       1  seteru (
       2  55hzatan   no precision because z is too close to +i or -i,
       3  55, 2, 2)

        if (abs(cmplx(1.0d0,0.0d0,kind(1.d0))+z*z).lt.sqeps) call seteru (

       1 50hzatan   answer lt half precision, z**2 close to -1, 50, 1, 1)
  c
   40   xans = 0.5*atan2 (2.0*x, 1.0-r2)
        yans = 0.25*log((r2+2.0*y+1.0)/(r2-2.0*y+1.0))

        zatan = cmplx (xans, yans, kind(1.d0))

        return
  c

   50   zatan = cmplx (pi2, 0.0d0, kind(1.d0))
        if (real(z).lt.0.0) zatan = cmplx (-pi2, 0.0d0, kind(1.d0))

        return
  c
        end