Blame view
fvn_fnlib/zexprl.f
1.3 KB
38581db0c git-svn-id: https... |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 |
complex(8) function zexprl (z) implicit none c august 1980 edition. w. fullerton, c3, los alamos scientific lab. c c evaluate (cexp(z)-1)/z . for small cabs(z), we use the taylor c series. we could instead use the expression c cexprl(z) = (exp(x)*exp(i*y)-1)/z c = (x*exprel(x) * (1 - 2*sin(y/2)**2) - 2*sin(y/2)**2 c + i*sin(y)*(1+x*exprel(x))) / z c complex(8) z real(8) sqeps,r,xn,xln,rbnd,d1mach,alneps external d1mach integer nterms,irold,irold2,i data nterms, rbnd, sqeps / 0, 2*0.0 / c if (nterms.ne.0) go to 10 alneps = log(d1mach(3)) xn = 3.72 - 0.3*alneps xln = log((xn+1.0)/1.36) nterms = xn - (xn*xln+alneps)/(xln+1.36) + 1.5 rbnd = d1mach(3) sqeps = 0.0 c 10 r = abs(z) if (r.le.0.5) go to 20 c call entsrc (irold, 1) zexprl = (exp(z) - 1.0) / z call erroff call entsrc (irold2, irold) if (abs(real(z)).lt.sqeps .and. abs(zexprl).lt.sqeps) call 1 seteru (41hzexprl answer lt half prec, z near ni2pi, 2 41, 1, 1) return c 20 zexprl = (1.0, 0.0) if (r.lt.rbnd) return c zexprl = (0.0, 0.0) do 30 i=1,nterms zexprl = 1.0 + zexprl*z/dble(nterms+2-i) 30 continue c return end |