Blame view
fvn_fnlib/zpsi.f
2.79 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 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 |
complex(8) function zpsi (zin)
implicit none
c may 1978 edition. w. fullerton, c3, los alamos scientific lab.
complex(8) zin, z, z2inv, corr, zcot
dimension bern(13)
real(8) bern,d1mach,pi,bound, dxrel, rmin, rbig
real(8) x,y,cabsz
integer nterm,ndx,n,i
external zcot, d1mach
c
data bern( 1) / .8333333333 3333333 e-1 /
data bern( 2) / -.8333333333 3333333 e-2 /
data bern( 3) / .3968253968 2539683 e-2 /
data bern( 4) / -.4166666666 6666667 e-2 /
data bern( 5) / .7575757575 7575758 e-2 /
data bern( 6) / -.2109279609 2796093 e-1 /
data bern( 7) / .8333333333 3333333 e-1 /
data bern( 8) / -.4432598039 2156863 e0 /
data bern( 9) / .3053954330 2701197 e1 /
data bern(10) / -.2645621212 1212121 e2 /
data bern(11) / .2814601449 2753623 e3 /
data bern(12) / -.3454885393 7728938 e4 /
data bern(13) / .5482758333 3333333 e5 /
c
data pi / 3.141592653 589793 e0 /
data nterm, bound, dxrel, rmin, rbig / 0, 4*0.0 /
c
if (nterm.ne.0) go to 10
nterm = -0.30*log(d1mach(3))
c maybe bound = n*(0.1*eps)**(-1/(2*n-1)) / (pi*exp(1))
bound = 0.1171*dble(nterm) *
1 (0.1*d1mach(3))**(-1.0/(2.0*dble(nterm)-1.0))
dxrel = sqrt(d1mach(4))
rmin = exp (max(log(d1mach(1)), -log(d1mach(2))) + 0.011 )
rbig = 1.0/d1mach(3)
c
10 z = zin
x = real(z)
y = aimag(z)
if (y.lt.0.0) z = conjg(z)
c
corr = (0.0, 0.0)
cabsz = abs(z)
if (x.ge.0.0 .and. cabsz.gt.bound) go to 50
if (x.lt.0.0 .and. abs(y).gt.bound) go to 50
c
if (cabsz.lt.bound) go to 20
c
c use the reflection formula for real(z) negative, cabs(z) large, and
c abs(aimag(y)) small.
c
corr = -pi*zcot(pi*z)
z = 1.0 - z
go to 50
c
c use the recursion relation for cabs(z) small.
c
20 if (cabsz.lt.rmin) call seteru (
1 56hzpsi zpsi called with z so near 0 that zpsi overflows,
1 56, 2, 2)
c
if (x.ge.(-0.5) .or. abs(y).gt.dxrel) go to 30
if ( abs((z-aint(x-0.5))/x).lt.dxrel) call seteru (
1 68hzpsi answer lt half precision because z too near negative
2integer, 68, 1, 1)
if (y.eq.0.0 .and. x.eq.aint(x)) call seteru (
1 31hzpsi z is a negative integer, 31, 3, 2)
c
30 n = sqrt(bound**2-y**2) - x + 1.0
do 40 i=1,n
corr = corr - 1.0/z
z = z + 1.0
40 continue
c
c now evaluate the asymptotic series for suitably large z.
c
50 if (cabsz.gt.rbig) zpsi = log(z) + corr
if (cabsz.gt.rbig) go to 70
c
zpsi = (0.0, 0.0)
z2inv = 1.0/z**2
do 60 i=1,nterm
ndx = nterm + 1 - i
zpsi = bern(ndx) + z2inv*zpsi
60 continue
zpsi = log(z) - 0.5/z - zpsi*z2inv + corr
c
70 if (y.lt.0.0) zpsi = conjg(zpsi)
c
return
end
|