complex(8) function zatan (z) implicit none c jan 1984 edition. w. fullerton, los alamos scientific lab. c complex(8) z, z2 real(8) pi2,sqeps,rmin,rmax,r2sml,r,d1mach,x,y,r2,xans,yans external d1mach integer nterms,i,twoi data pi2 /1.5707963267 9489661923d0 / 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*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(dcmplx(1.0,0.0)+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 = dcmplx (xans, yans) return c 50 zatan = dcmplx (pi2, 0.0) if (real(z).lt.0.0) zatan = dcmplx (-pi2, 0.0) return c end