complex function catan (z)
c jan 1984 edition.  w. fullerton, los alamos scientific lab.
c
      complex z, z2
      external  r1mach
      data pi2 /1.5707963267 9489661923e0 /
      data nterms, sqeps, rmin, rmax, r2sml / 0, 4*0.0 /
c
      if (nterms.ne.0) go to 10
c nterms = alog(eps)/alog(rbnd) where rbnd = 0.1
      nterms = -0.4343*alog(r1mach(3)) + 1.0
      sqeps = sqrt (r1mach(4))
      rmin = sqrt (3.0*r1mach(3))
      rmax = 1.0/r1mach(3)
      r2sml = 10.0*r1mach (4)
c
 10   r = cabs(z)
      if (r.gt.0.1) go to 30
c
      catan = z
      if (r.lt.rmin) return
c
      catan = (0.0, 0.0)
      z2 = z*z
      do 20 i=1,nterms
        twoi = 2*(nterms-i) + 1
        catan = 1.0/twoi - z2*catan
 20   continue
      catan = z*catan
      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  55hcatan   no precision because z is too close to +i or -i,
     3  55, 2, 2)
      if (cabs(cmplx(1.0,0.0)+z*z).lt.sqeps) call seteru (
     1 50hcatan   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*alog((r2+2.0*y+1.0)/(r2-2.0*y+1.0))
      catan = cmplx (xans, yans)
      return
c
 50   catan = cmplx (pi2, 0.0)
      if (real(z).lt.0.0) catan = cmplx (-pi2, 0.0)
      return
c
      end