lxsd / fvn

Download as
Browse Code ยป

Commit 00055ac08315061e3fd7a85487896e6297dd7726

Authored by kwagner
1 parent 0c3098aed9
Exists in master and in 2 other branches geevx_2020, multi

Define some constants as double precision objects, in order to remove compiler warnings. 

git-svn-id: https://lxsd.femto-st.fr/svn/fvn@70 b657c933-2333-4658-acf2-d3c7c2708721

Showing 16 changed files with 44 additions and 44 deletions Inline Diff

subroutine beskes (xnu, x, nin, bke) 1 1 subroutine beskes (xnu, x, nin, bke)
c july 1980 edition. w. fullerton, c3, los alamos scientific lab. 2 2 c july 1980 edition. w. fullerton, c3, los alamos scientific lab.
dimension bke(1) 3 3 dimension bke(*)
external r1mach 4 4 external r1mach
data alnbig / 0. / 5 5 data alnbig / 0. /
c 6 6 c
if (alnbig.eq.0.) alnbig = alog (r1mach(2)) 7 7 if (alnbig.eq.0.) alnbig = alog (r1mach(2))
c 8 8 c
v = abs(xnu) 9 9 v = abs(xnu)
n = iabs(nin) 10 10 n = iabs(nin)
c 11 11 c
if (v.ge.1.) call seteru (29hbeskes abs(xnu) must be lt 1, 29, 12 12 if (v.ge.1.) call seteru (29hbeskes abs(xnu) must be lt 1, 29,
1 1, 2) 13 13 1 1, 2)
if (x.le.0.) call seteru (17hbeskes x is le 0, 17, 2, 2) 14 14 if (x.le.0.) call seteru (17hbeskes x is le 0, 17, 2, 2)
if (n.eq.0) call seteru ( 15 15 if (n.eq.0) call seteru (
1 41hbeskes n the number in the sequence is 0, 41, 3, 2) 16 16 1 41hbeskes n the number in the sequence is 0, 41, 3, 2)
c 17 17 c
call r9knus (v, x, bke(1), bknu1, iswtch) 18 18 call r9knus (v, x, bke(1), bknu1, iswtch)
if (n.eq.1) return 19 19 if (n.eq.1) return
c 20 20 c
vincr = sign (1.0, float(nin)) 21 21 vincr = sign (1.0, float(nin))
direct = vincr 22 22 direct = vincr
if (xnu.ne.0.) direct = vincr*sign(1.0,xnu) 23 23 if (xnu.ne.0.) direct = vincr*sign(1.0,xnu)
if (iswtch.eq.1 .and. direct.gt.0.) call seteru ( 24 24 if (iswtch.eq.1 .and. direct.gt.0.) call seteru (
1 47hbeskes x so small bessel k-sub-xnu+1 overflows, 47, 4, 2) 25 25 1 47hbeskes x so small bessel k-sub-xnu+1 overflows, 47, 4, 2)
bke(2) = bknu1 26 26 bke(2) = bknu1
c 27 27 c
if (direct.lt.0.) call r9knus (abs(xnu+vincr), x, bke(2), bknu1, 28 28 if (direct.lt.0.) call r9knus (abs(xnu+vincr), x, bke(2), bknu1,
1 iswtch) 29 29 1 iswtch)
if (n.eq.2) return 30 30 if (n.eq.2) return
c 31 31 c
vend = abs(xnu+float(nin)) - 1.0 32 32 vend = abs(xnu+float(nin)) - 1.0
if (x+(vend-0.5)*alog(vend)+.27 - vend*(alog(x)+.306).gt.alnbig) 33 33 if (x+(vend-0.5)*alog(vend)+.27 - vend*(alog(x)+.306).gt.alnbig)
1 call seteru (67hbeskes x so small or abs(nu) so big that bessel 34 34 1 call seteru (67hbeskes x so small or abs(nu) so big that bessel
2 k-sub-nu overflows, 67, 5, 2) 35 35 2 k-sub-nu overflows, 67, 5, 2)
c 36 36 c
v = xnu 37 37 v = xnu
do 10 i=3,n 38 38 do 10 i=3,n
v = v + vincr 39 39 v = v + vincr
bke(i) = 2.0*v*bke(i-1)/x + bke(i-2) 40 40 bke(i) = 2.0*v*bke(i-1)/x + bke(i-2)
10 continue 41 41 10 continue
c 42 42 c
return 43 43 return
end 44 44 end
45 45
subroutine dbskes (xnu, x, nin, bke) 1 1 subroutine dbskes (xnu, x, nin, bke)
c july 1980 edition. w. fullerton, c3, los alamos scientific lab. 2 2 c july 1980 edition. w. fullerton, c3, los alamos scientific lab.
double precision xnu, x, bke(1), bknu1, v, vincr, vend, alnbig, 3 3 double precision xnu, x, bke(*), bknu1, v, vincr, vend, alnbig,
1 d1mach, dlog 4 4 1 d1mach, dlog
external d1mach 5 5 external d1mach
data alnbig / 0.d0 / 6 6 data alnbig / 0.d0 /
c 7 7 c
if (alnbig.eq.0.d0) alnbig = dlog (d1mach(2)) 8 8 if (alnbig.eq.0.d0) alnbig = dlog (d1mach(2))
c 9 9 c
v = dabs(xnu) 10 10 v = dabs(xnu)
n = iabs(nin) 11 11 n = iabs(nin)
c 12 12 c
if (v.ge.1.d0) call seteru (29hdbskes abs(xnu) must be lt 1, 29, 13 13 if (v.ge.1.d0) call seteru (29hdbskes abs(xnu) must be lt 1, 29,
1 1, 2) 14 14 1 1, 2)
if (x.le.0.d0) call seteru (17hdbskes x is le 0, 17, 2, 2) 15 15 if (x.le.0.d0) call seteru (17hdbskes x is le 0, 17, 2, 2)
if (n.eq.0) call seteru ( 16 16 if (n.eq.0) call seteru (
1 41hdbskes n the number in the sequence is 0, 41, 3, 2) 17 17 1 41hdbskes n the number in the sequence is 0, 41, 3, 2)
c 18 18 c
call d9knus (v, x, bke(1), bknu1, iswtch) 19 19 call d9knus (v, x, bke(1), bknu1, iswtch)
if (n.eq.1) return 20 20 if (n.eq.1) return
c 21 21 c
vincr = sign (1.0, float(nin)) 22 22 vincr = sign (1.0, float(nin))
direct = vincr 23 23 direct = vincr
if (xnu.ne.0.d0) direct = vincr*dsign(1.d0, xnu) 24 24 if (xnu.ne.0.d0) direct = vincr*dsign(1.d0, xnu)
if (iswtch.eq.1 .and. direct.gt.0.) call seteru ( 25 25 if (iswtch.eq.1 .and. direct.gt.0.) call seteru (
1 47hdbskes x so small bessel k-sub-xnu+1 overflows, 47, 4, 2) 26 26 1 47hdbskes x so small bessel k-sub-xnu+1 overflows, 47, 4, 2)
bke(2) = bknu1 27 27 bke(2) = bknu1
c 28 28 c
if (direct.lt.0.) call d9knus (dabs(xnu+vincr), x, bke(2), bknu1, 29 29 if (direct.lt.0.) call d9knus (dabs(xnu+vincr), x, bke(2), bknu1,
1 iswtch) 30 30 1 iswtch)
if (n.eq.2) return 31 31 if (n.eq.2) return
c 32 32 c
vend = dabs (xnu+dble(float(nin))) - 1.0d0 33 33 vend = dabs (xnu+dble(float(nin))) - 1.0d0
if (x + (vend-.5d0)*dlog(vend)+0.27d0-vend*(dlog(x)+.306d0) 34 34 if (x + (vend-.5d0)*dlog(vend)+0.27d0-vend*(dlog(x)+.306d0)
1 .gt.alnbig) call seteru ( 35 35 1 .gt.alnbig) call seteru (
2 67hdbskes x so small or abs(nu) so big that bessel 36 36 2 67hdbskes x so small or abs(nu) so big that bessel
3 k-sub-nu overflows, 67, 5, 2) 37 37 3 k-sub-nu overflows, 67, 5, 2)
c 38 38 c
v = xnu 39 39 v = xnu
do 10 i=3,n 40 40 do 10 i=3,n
v = v + vincr 41 41 v = v + vincr
bke(i) = 2.0d0*v*bke(i-1)/x + bke(i-2) 42 42 bke(i) = 2.0d0*v*bke(i-1)/x + bke(i-2)
10 continue 43 43 10 continue
c 44 44 c
return 45 45 return
end 46 46 end
47 47
complex(kind(1.d0)) function z8lgmc (zin) 1 1 complex(kind(1.d0)) function z8lgmc (zin)
implicit none 2 2 implicit none
c may 1978 edition. w. fullerton, c3, los alamos scientific lab. 3 3 c may 1978 edition. w. fullerton, c3, los alamos scientific lab.
c 4 4 c
c compute the log gamma correction term for almost all z so that 5 5 c compute the log gamma correction term for almost all z so that
c clog(cgamma(z)) = 0.5*alog(2.*pi) + (z-0.5)*clog(z) - z + c8lgmc(z) 6 6 c clog(cgamma(z)) = 0.5*alog(2.*pi) + (z-0.5)*clog(z) - z + c8lgmc(z)
c this routine is valid even when real(z) is negative or when cabs(z) 7 7 c this routine is valid even when real(z) is negative or when cabs(z)
c is small. 8 8 c is small.
c when real(z) is negative, c8lgmc merely returns a correction which 9 9 c when real(z) is negative, c8lgmc merely returns a correction which
c may be wrong by a multiple of 2*pi*i. 10 10 c may be wrong by a multiple of 2*pi*i.
c 11 11 c
complex(kind(1.d0)) zin, z, corr, z9lgmc, z0lgmc, zexp, zlnrel 12 12 complex(kind(1.d0)) zin, z, corr, z9lgmc, z0lgmc, zlnrel
real(kind(1.d0)) d1mach,pi,bound,sqeps,eps,x,y,absz,test 13 13 real(kind(1.d0)) d1mach,pi,bound,sqeps,eps,x,y,absz,test
integer nterm,n,i,irold,ir 14 14 integer nterm,n,i,irold,ir
c complex(kind(1.d0)) tmp_arg 15 15 c complex(kind(1.d0)) tmp_arg
external z0lgmc, z9lgmc, zlnrel, 16 16 external z0lgmc, z9lgmc, zlnrel,
1 d1mach 17 17 1 d1mach
data pi / 3.1415926535 8979324d0 / 18 18 data pi / 3.1415926535 8979324d0 /
data bound, sqeps, eps / 3*0.0 / 19 19 data bound, sqeps, eps / 3*0.0d0 /
c 20 20 c
if (bound.ne.0.0) go to 10 21 21 if (bound.ne.0.0) go to 10
nterm = -0.30*log(d1mach(3)) 22 22 nterm = -0.30*log(d1mach(3))
bound = 0.1170*dble(nterm)* 23 23 bound = 0.1170*dble(nterm)*
1 (0.1*d1mach(3))**(-1.0/(2.0*dble(nterm)-1.0)) 24 24 1 (0.1*d1mach(3))**(-1.0/(2.0*dble(nterm)-1.0))
eps = 10.0*d1mach(4) 25 25 eps = 10.0*d1mach(4)
sqeps = sqrt (d1mach(4)) 26 26 sqeps = sqrt (d1mach(4))
c 27 27 c
10 z = zin 28 28 10 z = zin
x = real(z) 29 29 x = real(z)
y = aimag(z) 30 30 y = aimag(z)
if (x.lt.0.0) z = -z 31 31 if (x.lt.0.0) z = -z
c 32 32 c
corr = (0.0, 0.0) 33 33 corr = (0.0, 0.0)
if (x.gt.bound .or. abs(y).ge.bound) go to 30 34 34 if (x.gt.bound .or. abs(y).ge.bound) go to 30
absz = abs(z) 35 35 absz = abs(z)
if (absz.ge.bound) go to 30 36 36 if (absz.ge.bound) go to 30
c 37 37 c
if (x.eq.0.0 .and. y.eq.0.0) call seteru ( 38 38 if (x.eq.0.0 .and. y.eq.0.0) call seteru (
1 57hz8lgmc log gamma correction term is infinite for z = 0.0, 39 39 1 57hz8lgmc log gamma correction term is infinite for z = 0.0,
2 57, 3, 2) 40 40 2 57, 3, 2)
c 41 41 c
n = sqrt (bound**2-y**2) - abs(x) + 1.0 42 42 n = sqrt (bound**2-y**2) - abs(x) + 1.0
do 20 i=1,n 43 43 do 20 i=1,n
corr = corr + z0lgmc(z) 44 44 corr = corr + z0lgmc(z)
z = z + 1.0 45 45 z = z + 1.0
20 continue 46 46 20 continue
c 47 47 c
30 z8lgmc = z9lgmc(z) + corr 48 48 30 z8lgmc = z9lgmc(z) + corr
if (x.gt.0.0) return 49 49 if (x.gt.0.0) return
c 50 50 c
call entsrc (irold, 1) 51 51 call entsrc (irold, 1)
test = 1.0 52 52 test = 1.0
if (x.lt.(-0.5)) test = abs ((z-aint(x-0.5))/z) 53 53 if (x.lt.(-0.5)) test = abs ((z-aint(x-0.5))/z)
if (test.lt.eps) call seteru ( 54 54 if (test.lt.eps) call seteru (
1 31hz8lgmc z is a negative integer, 31, 2, 2) 55 55 1 31hz8lgmc z is a negative integer, 31, 2, 2)
c 56 56 c
z = zin 57 57 z = zin
if (y.lt.0.0) z = conjg(z) 58 58 if (y.lt.0.0) z = conjg(z)
c tmp_arg = -exp(cmplx(0.,2.*pi,kind(1.d0))*z) 59 59 c tmp_arg = -exp(cmplx(0.,2.*pi,kind(1.d0))*z)
corr = -cmplx(0.0,pi,kind(1.d0)) + 60 60 corr = -cmplx(0.0,pi,kind(1.d0)) +
1 zlnrel(-exp(cmplx(0.,2.*pi,kind(1.d0))*z)) 61 61 1 zlnrel(-exp(cmplx(0.,2.*pi,kind(1.d0))*z))
if (y.lt.0.0) corr = conjg(corr) 62 62 if (y.lt.0.0) corr = conjg(corr)
c 63 63 c
z8lgmc = corr - z8lgmc 64 64 z8lgmc = corr - z8lgmc
c 65 65 c
call erroff 66 66 call erroff
call entsrc (ir, irold) 67 67 call entsrc (ir, irold)
if (test.lt.sqeps) call seteru ( 63hc8lgmc answer lt half precisi 68 68 if (test.lt.sqeps) call seteru ( 63hc8lgmc answer lt half precisi
1on, z too near a negative integer, 63, 1, 1) 69 69 1on, z too near a negative integer, 63, 1, 1)
c 70 70 c
return 71 71 return
end 72 72 end
73 73
complex(kind(1.d0)) function z9lgmc (zin) 1 1 complex(kind(1.d0)) function z9lgmc (zin)
implicit none 2 2 implicit none
c april 1978 edition. w. fullerton c3, los alamos scientific lab. 3 3 c april 1978 edition. w. fullerton c3, los alamos scientific lab.
c 4 4 c
c compute the log gamma correction term for large cabs(z) when real(z) 5 5 c compute the log gamma correction term for large cabs(z) when real(z)
c .ge. 0.0 and for large abs(aimag(y)) when real(z) .lt. 0.0. we find 6 6 c .ge. 0.0 and for large abs(aimag(y)) when real(z) .lt. 0.0. we find
c c9lgmc so that 7 7 c c9lgmc so that
c clog(cgamma(z)) = 0.5*alog(2.*pi) + (z-0.5)*clog(z) - z + c9lgmc(z). 8 8 c clog(cgamma(z)) = 0.5*alog(2.*pi) + (z-0.5)*clog(z) - z + c9lgmc(z).
c 9 9 c
complex(kind(1.d0)) zin, z, z2inv 10 10 complex(kind(1.d0)) zin, z, z2inv
real(kind(1.d0)) d1mach,bern,xbig,xmax,bound,cabsz,x,y 11 11 real(kind(1.d0)) d1mach,bern,xbig,xmax,bound,cabsz,x,y
integer nterm,i,ndx 12 12 integer nterm,i,ndx
external d1mach 13 13 external d1mach
c 14 14 c
dimension bern(11) 15 15 dimension bern(11)
data bern( 1) / .08333333333 3333333e0 / 16 16 data bern( 1) / .08333333333 3333333d0 /
data bern( 2) / -.002777777777 7777778e0 / 17 17 data bern( 2) / -.002777777777 7777778d0 /
data bern( 3) / .0007936507936 5079365e0 / 18 18 data bern( 3) / .0007936507936 5079365d0 /
data bern( 4) / -.0005952380952 3809524e0 / 19 19 data bern( 4) / -.0005952380952 3809524d0 /
data bern( 5) / .0008417508417 5084175e0 / 20 20 data bern( 5) / .0008417508417 5084175d0 /
data bern( 6) / -.001917526917 5269175e0 / 21 21 data bern( 6) / -.001917526917 5269175d0 /
data bern( 7) / .006410256410 2564103e0 / 22 22 data bern( 7) / .006410256410 2564103d0 /
data bern( 8) / -.02955065359 4771242e0 / 23 23 data bern( 8) / -.02955065359 4771242d0 /
data bern( 9) / .1796443723 6883057e0 / 24 24 data bern( 9) / .1796443723 6883057d0 /
data bern(10) / -1.392432216 9059011e0 / 25 25 data bern(10) / -1.392432216 9059011d0 /
data bern(11) / 13.40286404 4168392e0 / 26 26 data bern(11) / 13.40286404 4168392d0 /
c 27 27 c
data nterm, bound, xbig, xmax / 0, 3*0.0 / 28 28 data nterm, bound, xbig, xmax / 0d0, 3*0.0d0 /
c 29 29 c
if (nterm.ne.0) go to 10 30 30 if (nterm.ne.0) go to 10
c 31 31 c
nterm = -0.30*log(d1mach(3)) 32 32 nterm = -0.30*log(d1mach(3))
bound = 0.1170*dble(nterm)* 33 33 bound = 0.1170*dble(nterm)*
1 (0.1*d1mach(3))**(-1./(2.*dble(nterm)-1.)) 34 34 1 (0.1*d1mach(3))**(-1./(2.*dble(nterm)-1.))
xbig = 1.0/sqrt(d1mach(3)) 35 35 xbig = 1.0/sqrt(d1mach(3))
xmax = exp (min(log(d1mach(2)/12.0), -log(12.*d1mach(1))) ) 36 36 xmax = exp (min(log(d1mach(2)/12.0), -log(12.*d1mach(1))) )
c 37 37 c
10 z = zin 38 38 10 z = zin
x = real (z) 39 39 x = real (z)
y = aimag(z) 40 40 y = aimag(z)
cabsz = abs(z) 41 41 cabsz = abs(z)
c 42 42 c
if (x.lt.0.0 .and. abs(y).lt.bound) call seteru ( 69hc9lgmc c9lgm 43 43 if (x.lt.0.0 .and. abs(y).lt.bound) call seteru ( 69hc9lgmc c9lgm
1c not valid for negative real(z) and small abs(aimag(z)), 69, 2,2) 44 44 1c not valid for negative real(z) and small abs(aimag(z)), 69, 2,2)
if (cabsz.lt.bound) call seteru ( 45 45 if (cabsz.lt.bound) call seteru (
1 42hz9lgmc z9lgmc not valid for small cabs(z), 42, 3, 2) 46 46 1 42hz9lgmc z9lgmc not valid for small cabs(z), 42, 3, 2)
c 47 47 c
if (cabsz.ge.xmax) go to 50 48 48 if (cabsz.ge.xmax) go to 50
c 49 49 c
if (cabsz.ge.xbig) z9lgmc = 1.0/(12.0*z) 50 50 if (cabsz.ge.xbig) z9lgmc = 1.0/(12.0*z)
if (cabsz.ge.xbig) return 51 51 if (cabsz.ge.xbig) return
c 52 52 c
z2inv = 1.0/z**2 53 53 z2inv = 1.0/z**2
z9lgmc = (0.0, 0.0) 54 54 z9lgmc = (0.0, 0.0)
do 40 i=1,nterm 55 55 do 40 i=1,nterm
ndx = nterm + 1 - i 56 56 ndx = nterm + 1 - i
z9lgmc = bern(ndx) + z9lgmc*z2inv 57 57 z9lgmc = bern(ndx) + z9lgmc*z2inv
40 continue 58 58 40 continue
c 59 59 c
z9lgmc = z9lgmc/z 60 60 z9lgmc = z9lgmc/z
return 61 61 return
c 62 62 c
50 z9lgmc = (0.0, 0.0) 63 63 50 z9lgmc = (0.0, 0.0)
call seteru (34hz9lgmc z so big z9lgmc underflows, 34, 1, 0) 64 64 call seteru (34hz9lgmc z so big z9lgmc underflows, 34, 1, 0)
return 65 65 return
c 66 66 c
end 67 67 end
68 68
complex(kind(1.d0)) function zasin (zinp) 1 1 complex(kind(1.d0)) function zasin (zinp)
implicit none 2 2 implicit none
c august 1980 edition. w. fullerton, c3, los alamos scientific lab. 3 3 c august 1980 edition. w. fullerton, c3, los alamos scientific lab.
c 4 4 c
c ref -- l. l. pennisi, elements of complex variables, holt, rinehart 5 5 c ref -- l. l. pennisi, elements of complex variables, holt, rinehart
c and winston, 1963. page 126. 6 6 c and winston, 1963. page 126.
c 7 7 c
complex(kind(1.d0)) zinp, z, z2, sqzp1, ci 8 8 complex(kind(1.d0)) zinp, z, z2, sqzp1, ci
real(kind(1.d0)) d1mach,pi2,pi,rmin,r 9 9 real(kind(1.d0)) d1mach,pi2,pi,rmin,r
integer nterms,i,twoi 10 10 integer nterms,i,twoi
external d1mach 11 11 external d1mach
12 12
data pi2 /1.5707963267 9489661923d0/ 13 13 data pi2 /1.5707963267 9489661923d0/
data pi /3.1415926535 8979324d0/ 14 14 data pi /3.1415926535 8979324d0/
data ci /(0.,1.)/ 15 15 data ci /(0.d0,1.d0)/
data nterms, rmin / 0, 0.0 / 16 16 data nterms, rmin / 0, 0.0d0 /
c 17 17 c
if (nterms.ne.0) go to 10 18 18 if (nterms.ne.0) go to 10
c nterms = alog(eps)/alog(rmax) where rmax = 0.1 19 19 c nterms = alog(eps)/alog(rmax) where rmax = 0.1
nterms = -0.4343*log(d1mach(3)) + 1.0 20 20 nterms = -0.4343*log(d1mach(3)) + 1.0
rmin = sqrt (6.0*d1mach(3)) 21 21 rmin = sqrt (6.0*d1mach(3))
c 22 22 c
10 z = zinp 23 23 10 z = zinp
r = abs (z) 24 24 r = abs (z)
if (r.gt.0.1) go to 30 25 25 if (r.gt.0.1) go to 30
c 26 26 c
zasin = z 27 27 zasin = z
if (r.lt.rmin) return 28 28 if (r.lt.rmin) return
c 29 29 c
zasin = (0.0, 0.0) 30 30 zasin = (0.0, 0.0)
z2 = z*z 31 31 z2 = z*z
do 20 i=1,nterms 32 32 do 20 i=1,nterms
twoi = 2*(nterms-i) + 1 33 33 twoi = 2*(nterms-i) + 1
zasin = 1.0/twoi + twoi*zasin*z2/(twoi+1.0) 34 34 zasin = 1.0/twoi + twoi*zasin*z2/(twoi+1.0)
20 continue 35 35 20 continue
zasin = z*zasin 36 36 zasin = z*zasin
return 37 37 return
c 38 38 c
30 if (real(zinp).lt.0.0) z = -zinp 39 39 30 if (real(zinp).lt.0.0) z = -zinp
c 40 40 c
sqzp1 = sqrt (z+1.0) 41 41 sqzp1 = sqrt (z+1.0)
if (aimag(sqzp1).lt.0.) sqzp1 = -sqzp1 42 42 if (aimag(sqzp1).lt.0.) sqzp1 = -sqzp1
zasin = pi2 - ci * log (z + sqzp1*sqrt(z-1.0)) 43 43 zasin = pi2 - ci * log (z + sqzp1*sqrt(z-1.0))
c 44 44 c
if (real(zasin).gt.pi2) zasin = pi - zasin 45 45 if (real(zasin).gt.pi2) zasin = pi - zasin
if (real(zasin).le.(-pi2)) zasin = -pi - zasin 46 46 if (real(zasin).le.(-pi2)) zasin = -pi - zasin
if (real(zinp).lt.0.) zasin = -zasin 47 47 if (real(zinp).lt.0.) zasin = -zasin
c 48 48 c
return 49 49 return
end 50 50 end
51 51
complex(kind(1.d0)) function zatanh (z) 1 1 complex(kind(1.d0)) function zatanh (z)
implicit none 2 2 implicit none
c april 1977 version. w. fullerton, c3, los alamos scientific lab. 3 3 c april 1977 version. w. fullerton, c3, los alamos scientific lab.
complex(kind(1.d0)) z, ci, zatan 4 4 complex(kind(1.d0)) z, ci, zatan
external zatan 5 5 external zatan
data ci /(0.,1.)/ 6 6 data ci /(0.d0,1.d0)/
c 7 7 c
zatanh = -ci*zatan(ci*z) 8 8 zatanh = -ci*zatan(ci*z)
c 9 9 c
return 10 10 return
end 11 11 end
12 12
complex(kind(1.d0)) function zbeta (a, b) 1 1 complex(kind(1.d0)) function zbeta (a, b)
implicit none 2 2 implicit none
c july 1977 edition. w. fullerton, c3, los alamos scientific lab. 3 3 c july 1977 edition. w. fullerton, c3, los alamos scientific lab.
complex(kind(1.d0)) a, b, zgamma, zlbeta 4 4 complex(kind(1.d0)) a, b, zgamma, zlbeta
real(kind(1.d0)) xmax,xmin 5 5 real(kind(1.d0)) xmax,xmin
external zgamma, zlbeta 6 6 external zgamma, zlbeta
data xmax / 0.0 / 7 7 data xmax / 0.0d0 /
c 8 8 c
if (xmax.eq.0.0) call d9gaml (xmin, xmax) 9 9 if (xmax.eq.0.0) call d9gaml (xmin, xmax)
c 10 10 c
if (real(a).le.0.0 .or. real(b).le.0.0) call seteru ( 11 11 if (real(a).le.0.0 .or. real(b).le.0.0) call seteru (
1 48hzbeta real part of both arguments must be gt 0, 48, 1, 2) 12 12 1 48hzbeta real part of both arguments must be gt 0, 48, 1, 2)
c 13 13 c
if (real(a)+real(b).lt.xmax) zbeta = zgamma(a) * (zgamma(b)/ 14 14 if (real(a)+real(b).lt.xmax) zbeta = zgamma(a) * (zgamma(b)/
1 zgamma(a+b) ) 15 15 1 zgamma(a+b) )
if (real(a)+real(b).lt.xmax) return 16 16 if (real(a)+real(b).lt.xmax) return
c 17 17 c
zbeta = exp (zlbeta(a, b)) 18 18 zbeta = exp (zlbeta(a, b))
c 19 19 c
return 20 20 return
end 21 21 end
22 22
complex(kind(1.d0)) function zcosh (z) 1 1 complex(kind(1.d0)) function zcosh (z)
implicit none 2 2 implicit none
c april 1977 version. w. fullerton, c3, los alamos scientific lab. 3 3 c april 1977 version. w. fullerton, c3, los alamos scientific lab.
complex(kind(1.d0)) z, ci 4 4 complex(kind(1.d0)) z, ci
data ci /(0.,1.)/ 5 5 data ci /(0.d0,1.d0)/
c 6 6 c
zcosh = cos (ci*z) 7 7 zcosh = cos (ci*z)
c 8 8 c
return 9 9 return
end 10 10 end
11 11
complex(kind(1.d0)) function zcot (z) 1 1 complex(kind(1.d0)) function zcot (z)
implicit none 2 2 implicit none
c march 1979 edition. w. fullerton, c3, los alamos scientific lab. 3 3 c march 1979 edition. w. fullerton, c3, los alamos scientific lab.
complex(kind(1.d0)) z 4 4 complex(kind(1.d0)) z
real(kind(1.d0)) d1mach 5 5 real(kind(1.d0)) d1mach
real(kind(1.d0)) eps, xmax, ylarge, ybig, rmin, ymin 6 6 real(kind(1.d0)) eps, xmax, ylarge, ybig, rmin, ymin
real(kind(1.d0)) x,y,x2,y2,sn2x,den 7 7 real(kind(1.d0)) x,y,x2,y2,sn2x,den
integer irold,irold2 8 8 integer irold,irold2
external d1mach 9 9 external d1mach
data eps, xmax, ylarge, ybig, rmin, ymin / 5*0.0, 1.5 / 10 10 data eps, xmax, ylarge, ybig, rmin, ymin / 5*0.0d0, 1.5d0 /
c 11 11 c
if (eps.ne.0.0) go to 10 12 12 if (eps.ne.0.0) go to 10
eps = d1mach(4) 13 13 eps = d1mach(4)
xmax = 1.0/eps 14 14 xmax = 1.0/eps
ylarge = -0.5*log(0.5*d1mach(3)) 15 15 ylarge = -0.5*log(0.5*d1mach(3))
ybig = -0.5*log(0.5*sqrt(d1mach(3))) 16 16 ybig = -0.5*log(0.5*sqrt(d1mach(3)))
rmin = exp (max(log(d1mach(1)), -log(d1mach(2))) + 0.01) 17 17 rmin = exp (max(log(d1mach(1)), -log(d1mach(2))) + 0.01)
c 18 18 c
10 x = real (z) 19 19 10 x = real (z)
y = aimag (z) 20 20 y = aimag (z)
if (abs(y).gt.ylarge) go to 30 21 21 if (abs(y).gt.ylarge) go to 30
c 22 22 c
x2 = 2.0*x 23 23 x2 = 2.0*x
y2 = 2.0*y 24 24 y2 = 2.0*y
if (abs(x2).gt.xmax) go to 20 25 25 if (abs(x2).gt.xmax) go to 20
if ( abs(z).lt.rmin) call seteru (55hzcot abs(z) is zero or so 26 26 if ( abs(z).lt.rmin) call seteru (55hzcot abs(z) is zero or so
1 small that zcot overflows, 55, 5, 2) 27 27 1 small that zcot overflows, 55, 5, 2)
c 28 28 c
call entsrc (irold, 1) 29 29 call entsrc (irold, 1)
sn2x = sin(x2) 30 30 sn2x = sin(x2)
call erroff 31 31 call erroff
den = cosh(y2) - cos(x2) 32 32 den = cosh(y2) - cos(x2)
call erroff 33 33 call erroff
call entsrc (irold2, irold) 34 34 call entsrc (irold2, irold)
c 35 35 c
if (den.lt.x*x2*eps*eps) call seteru (74hzcot cot is nearly sin 36 36 if (den.lt.x*x2*eps*eps) call seteru (74hzcot cot is nearly sin
1gular, x is near a multiple of pi and y is near 0, 74, 3, 2) 37 37 1gular, x is near a multiple of pi and y is near 0, 74, 3, 2)
if (den.lt.x*x2*eps .and. abs(x).gt.0.5) call seteru (65hzcot a 38 38 if (den.lt.x*x2*eps .and. abs(x).gt.0.5) call seteru (65hzcot a
1nswer is lt half precision, x is near pi and y is near 0, 65, 1,1) 39 39 1nswer is lt half precision, x is near pi and y is near 0, 65, 1,1)
c 40 40 c
zcot = cmplx (sn2x/den, -sinh(y2)/den, kind(1.d0)) 41 41 zcot = cmplx (sn2x/den, -sinh(y2)/den, kind(1.d0))
return 42 42 return
c 43 43 c
20 if (abs(y).lt.ymin) call seteru (75hzcot answer would have no p 44 44 20 if (abs(y).lt.ymin) call seteru (75hzcot answer would have no p
1recision, abs(x) is very big and abs(y) small, 75, 4, 2) 45 45 1recision, abs(x) is very big and abs(y) small, 75, 4, 2)
if (abs(y).lt.ybig) call seteru (69hzcot answer lt half precisi 46 46 if (abs(y).lt.ybig) call seteru (69hzcot answer lt half precisi
1on, abs(x) is very big and abs(y) small, 69, 2, 1) 47 47 1on, abs(x) is very big and abs(y) small, 69, 2, 1)
c 48 48 c
zcot = cmplx (0.0, 1.0/tanh(y2), kind(1.d0)) 49 49 zcot = cmplx (0.0, 1.0/tanh(y2), kind(1.d0))
return 50 50 return
c 51 51 c
30 zcot = cmplx (0.0, sign (1.0d0, y), kind(1.d0)) 52 52 30 zcot = cmplx (0.0, sign (1.0d0, y), kind(1.d0))
return 53 53 return
c 54 54 c
end 55 55 end
56 56
complex(kind(1.d0)) function zexprl (z) 1 1 complex(kind(1.d0)) function zexprl (z)
implicit none 2 2 implicit none
c august 1980 edition. w. fullerton, c3, los alamos scientific lab. 3 3 c august 1980 edition. w. fullerton, c3, los alamos scientific lab.
c 4 4 c
c evaluate (cexp(z)-1)/z . for small cabs(z), we use the taylor 5 5 c evaluate (cexp(z)-1)/z . for small cabs(z), we use the taylor
c series. we could instead use the expression 6 6 c series. we could instead use the expression
c cexprl(z) = (exp(x)*exp(i*y)-1)/z 7 7 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 8 8 c = (x*exprel(x) * (1 - 2*sin(y/2)**2) - 2*sin(y/2)**2
c + i*sin(y)*(1+x*exprel(x))) / z 9 9 c + i*sin(y)*(1+x*exprel(x))) / z
c 10 10 c
complex(kind(1.d0)) z 11 11 complex(kind(1.d0)) z
real(kind(1.d0)) sqeps,r,xn,xln,rbnd,d1mach,alneps 12 12 real(kind(1.d0)) sqeps,r,xn,xln,rbnd,d1mach,alneps
external d1mach 13 13 external d1mach
integer nterms,irold,irold2,i 14 14 integer nterms,irold,irold2,i
data nterms, rbnd, sqeps / 0, 2*0.0 / 15 15 data nterms, rbnd, sqeps / 0d0, 2*0.0d0 /
c 16 16 c
if (nterms.ne.0) go to 10 17 17 if (nterms.ne.0) go to 10
alneps = log(d1mach(3)) 18 18 alneps = log(d1mach(3))
xn = 3.72 - 0.3*alneps 19 19 xn = 3.72 - 0.3*alneps
xln = log((xn+1.0)/1.36) 20 20 xln = log((xn+1.0)/1.36)
nterms = xn - (xn*xln+alneps)/(xln+1.36) + 1.5 21 21 nterms = xn - (xn*xln+alneps)/(xln+1.36) + 1.5
rbnd = d1mach(3) 22 22 rbnd = d1mach(3)
sqeps = 0.0 23 23 sqeps = 0.0
c 24 24 c
10 r = abs(z) 25 25 10 r = abs(z)
if (r.le.0.5) go to 20 26 26 if (r.le.0.5) go to 20
c 27 27 c
call entsrc (irold, 1) 28 28 call entsrc (irold, 1)
zexprl = (exp(z) - 1.0) / z 29 29 zexprl = (exp(z) - 1.0) / z
call erroff 30 30 call erroff
call entsrc (irold2, irold) 31 31 call entsrc (irold2, irold)
if (abs(real(z)).lt.sqeps .and. abs(zexprl).lt.sqeps) call 32 32 if (abs(real(z)).lt.sqeps .and. abs(zexprl).lt.sqeps) call
1 seteru (41hzexprl answer lt half prec, z near ni2pi, 33 33 1 seteru (41hzexprl answer lt half prec, z near ni2pi,
2 41, 1, 1) 34 34 2 41, 1, 1)
return 35 35 return
c 36 36 c
20 zexprl = (1.0, 0.0) 37 37 20 zexprl = (1.0, 0.0)
if (r.lt.rbnd) return 38 38 if (r.lt.rbnd) return
c 39 39 c
zexprl = (0.0, 0.0) 40 40 zexprl = (0.0, 0.0)
do 30 i=1,nterms 41 41 do 30 i=1,nterms
zexprl = 1.0 + zexprl*z/dble(nterms+2-i) 42 42 zexprl = 1.0 + zexprl*z/dble(nterms+2-i)
30 continue 43 43 30 continue
c 44 44 c
return 45 45 return
end 46 46 end
47 47
complex(kind(1.d0)) function zlngam (zin) 1 1 complex(kind(1.d0)) function zlngam (zin)
implicit none 2 2 implicit none
c august 1980 edition. w. fullerton c3, los alamos scientific lab. 3 3 c august 1980 edition. w. fullerton c3, los alamos scientific lab.
c eventually clngam should make use of c8lgmc for all z except for 4 4 c eventually clngam should make use of c8lgmc for all z except for
c z in the vicinity of 1 and 2. 5 5 c z in the vicinity of 1 and 2.
complex(kind(1.d0)) zin, z, corr, zlnrel, z9lgmc 6 6 complex(kind(1.d0)) zin, z, corr, zlnrel, z9lgmc
real(kind(1.d0)) d1mach,pi,sq2pil,bound,dxrel,rmax,cabsz 7 7 real(kind(1.d0)) d1mach,pi,sq2pil,bound,dxrel,rmax,cabsz
real(kind(1.d0)) x,y,argsum,zarg 8 8 real(kind(1.d0)) x,y,argsum,zarg
integer irold,ir,i,n 9 9 integer irold,ir,i,n
external z9lgmc, zarg, zlnrel,d1mach 10 10 external z9lgmc, zarg, zlnrel,d1mach
data pi / 3.1415926535 8979324d0 / 11 11 data pi / 3.1415926535 8979324d0 /
data sq2pil / 0.9189385332 0467274d0 / 12 12 data sq2pil / 0.9189385332 0467274d0 /
c 13 13 c
data bound, dxrel, rmax / 3*0.0 / 14 14 data bound, dxrel, rmax / 3*0.0d0 /
c 15 15 c
if (bound.ne.0.) go to 10 16 16 if (bound.ne.0.) go to 10
n = -0.30*log(d1mach(3)) 17 17 n = -0.30*log(d1mach(3))
c bound = n*(0.1*eps)**(-1/(2*n-1))/(pi*exp(1)) 18 18 c bound = n*(0.1*eps)**(-1/(2*n-1))/(pi*exp(1))
bound = 0.1171*dble(n)*(0.1*d1mach(3))**(-1./(2.*dble(n)-1.)) 19 19 bound = 0.1171*dble(n)*(0.1*d1mach(3))**(-1./(2.*dble(n)-1.))
dxrel = sqrt (d1mach(4)) 20 20 dxrel = sqrt (d1mach(4))
rmax = d1mach(2)/log(d1mach(2)) 21 21 rmax = d1mach(2)/log(d1mach(2))
c 22 22 c
10 z = zin 23 23 10 z = zin
x = real(zin) 24 24 x = real(zin)
y = aimag(zin) 25 25 y = aimag(zin)
c 26 26 c
corr = (0.0, 0.0) 27 27 corr = (0.0, 0.0)
cabsz = abs(z) 28 28 cabsz = abs(z)
if (cabsz.gt.rmax) call seteru ( 29 29 if (cabsz.gt.rmax) call seteru (
1 44hzlngam abs(z) so large result may overflow, 44, 2, 2) 30 30 1 44hzlngam abs(z) so large result may overflow, 44, 2, 2)
if (x.ge.0.0 .and. cabsz.gt.bound) go to 50 31 31 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 32 32 if (x.lt.0.0 .and. abs(y).gt.bound) go to 50
c 33 33 c
if (cabsz.lt.bound) go to 20 34 34 if (cabsz.lt.bound) go to 20
c 35 35 c
c use the reflection formula for real(z) negative, cabs(z) large, and 36 36 c use the reflection formula for real(z) negative, cabs(z) large, and
c abs(aimag(y)) small. 37 37 c abs(aimag(y)) small.
c 38 38 c
call entsrc (irold, 1) 39 39 call entsrc (irold, 1)
if (y.gt.0.0) z = conjg (z) 40 40 if (y.gt.0.0) z = conjg (z)
corr = exp (-cmplx(0.0,2.0*pi,kind(1.d0))*z) 41 41 corr = exp (-cmplx(0.0,2.0*pi,kind(1.d0))*z)
if (real(corr).eq.1.0 .and. aimag(corr).eq.0.0) call seteru ( 42 42 if (real(corr).eq.1.0 .and. aimag(corr).eq.0.0) call seteru (
1 31hzlngam z is a negative integer, 31, 3, 2) 43 43 1 31hzlngam z is a negative integer, 31, 3, 2)
c 44 44 c
zlngam = sq2pil + 1.0 - cmplx(0.0,pi,kind(1.d0))*(z-0.5) - 45 45 zlngam = sq2pil + 1.0 - cmplx(0.0,pi,kind(1.d0))*(z-0.5) -
1 zlnrel(-corr) + (z-0.5)*log(1.0-z) - z - z9lgmc(1.0-z) 46 46 1 zlnrel(-corr) + (z-0.5)*log(1.0-z) - z - z9lgmc(1.0-z)
if (y.gt.0.0) zlngam = conjg (zlngam) 47 47 if (y.gt.0.0) zlngam = conjg (zlngam)
c 48 48 c
call erroff 49 49 call erroff
call entsrc (ir, irold) 50 50 call entsrc (ir, irold)
if (abs(y).gt.dxrel) return 51 51 if (abs(y).gt.dxrel) return
c X in sequence field on next line is to preserve trailing blank 52 52 c X in sequence field on next line is to preserve trailing blank
if (0.5* abs((1.-corr)*zlngam/z).lt.dxrel) call seteru (68hzlngam X 53 53 if (0.5* abs((1.-corr)*zlngam/z).lt.dxrel) call seteru (68hzlngam X
1 answer lt half precision because z too near negative integer, 68, 54 54 1 answer lt half precision because z too near negative integer, 68,
2 1, 1) 55 55 2 1, 1)
return 56 56 return
c 57 57 c
c use the recursion relation for cabs(z) small. 58 58 c use the recursion relation for cabs(z) small.
c 59 59 c
20 if (x.ge.(-0.5) .or. abs(y).gt.dxrel) go to 30 60 60 20 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 ( 68hzlngam ans 61 61 if ( abs((z-aint(x-0.5))/x).lt.dxrel) call seteru ( 68hzlngam ans
1wer lt half precision because z too near negative integer, 68,1,1) 62 62 1wer lt half precision because z too near negative integer, 68,1,1)
c 63 63 c
30 n = sqrt (bound**2 - y**2) - x + 1.0 64 64 30 n = sqrt (bound**2 - y**2) - x + 1.0
argsum = 0.0 65 65 argsum = 0.0
corr = (1.0, 0.0) 66 66 corr = (1.0, 0.0)
do 40 i=1,n 67 67 do 40 i=1,n
argsum = argsum + zarg(z) 68 68 argsum = argsum + zarg(z)
corr = z*corr 69 69 corr = z*corr
z = 1.0 + z 70 70 z = 1.0 + z
40 continue 71 71 40 continue
c 72 72 c
if (real(corr).eq.0.0 .and. aimag(corr).eq.0.0) call seteru ( 73 73 if (real(corr).eq.0.0 .and. aimag(corr).eq.0.0) call seteru (
1 31hzlngam z is a negative integer, 31, 3, 2) 74 74 1 31hzlngam z is a negative integer, 31, 3, 2)
corr = -cmplx (log(abs(corr)), argsum, kind(1.d0)) 75 75 corr = -cmplx (log(abs(corr)), argsum, kind(1.d0))
c 76 76 c
c use stirling-s approximation for large z. 77 77 c use stirling-s approximation for large z.
c 78 78 c
50 zlngam = sq2pil + (z-0.5)*log(z) - z + corr + z9lgmc(z) 79 79 50 zlngam = sq2pil + (z-0.5)*log(z) - z + corr + z9lgmc(z)
return 80 80 return
c 81 81 c
end 82 82 end
83 83
complex(kind(1.d0)) function zlnrel (z) 1 1 complex(kind(1.d0)) function zlnrel (z)
implicit none 2 2 implicit none
c april 1977 version. w. fullerton, c3, los alamos scientific lab. 3 3 c april 1977 version. w. fullerton, c3, los alamos scientific lab.
c 4 4 c
c clnrel(z) = clog(1+z) with relative error accuracy near z = 0. 5 5 c clnrel(z) = clog(1+z) with relative error accuracy near z = 0.
c let rho = cabs(z) and 6 6 c let rho = cabs(z) and
c r**2 = cabs(1+z)**2 = (1+x)**2 + y**2 = 1 + 2*x + rho**2 . 7 7 c r**2 = cabs(1+z)**2 = (1+x)**2 + y**2 = 1 + 2*x + rho**2 .
c now if rho is small we may evaluate clnrel(z) accurately by 8 8 c now if rho is small we may evaluate clnrel(z) accurately by
c clog(1+z) = cmplx (alog(r), carg(1+z)) 9 9 c clog(1+z) = cmplx (alog(r), carg(1+z))
c = cmplx (0.5*alog(r**2), carg(1+z)) 10 10 c = cmplx (0.5*alog(r**2), carg(1+z))
c = cmplx (0.5*alnrel(2*x+rho**2), carg(1+z)) 11 11 c = cmplx (0.5*alnrel(2*x+rho**2), carg(1+z))
c 12 12 c
complex(kind(1.d0)) z 13 13 complex(kind(1.d0)) z
real(kind(1.d0)) dlnrel,d1mach,zarg,sqeps,x,rho 14 14 real(kind(1.d0)) dlnrel,d1mach,zarg,sqeps,x,rho
external dlnrel, zarg, d1mach 15 15 external dlnrel, zarg, d1mach
data sqeps /0.0/ 16 16 data sqeps /0.0d0/
c 17 17 c
if (sqeps.eq.0.) sqeps = sqrt (d1mach(4)) 18 18 if (sqeps.eq.0.) sqeps = sqrt (d1mach(4))
c 19 19 c
if ( abs(1.+z).lt.sqeps) call seteru ( 20 20 if ( abs(1.+z).lt.sqeps) call seteru (
1 54hzlnrel answer lt half precision because z too near -1, 54, 21 21 1 54hzlnrel answer lt half precision because z too near -1, 54,
2 1, 1) 22 22 2 1, 1)
c 23 23 c
rho = abs(z) 24 24 rho = abs(z)
if (rho.gt.0.375) zlnrel = log (1.0+z) 25 25 if (rho.gt.0.375) zlnrel = log (1.0+z)
if (rho.gt.0.375) return 26 26 if (rho.gt.0.375) return
c 27 27 c
x = real(z) 28 28 x = real(z)
zlnrel = cmplx (0.5*dlnrel(2.*x+rho**2), zarg(1.0+z), kind(1.d0)) 29 29 zlnrel = cmplx (0.5*dlnrel(2.*x+rho**2), zarg(1.0+z), kind(1.d0))
c 30 30 c
return 31 31 return
end 32 32 end
33 33
complex(kind(1.d0)) function zpsi (zin) 1 1 complex(kind(1.d0)) function zpsi (zin)
implicit none 2 2 implicit none
c may 1978 edition. w. fullerton, c3, los alamos scientific lab. 3 3 c may 1978 edition. w. fullerton, c3, los alamos scientific lab.
complex(kind(1.d0)) zin, z, z2inv, corr, zcot 4 4 complex(kind(1.d0)) zin, z, z2inv, corr, zcot
dimension bern(13) 5 5 dimension bern(13)
real(kind(1.d0)) bern,d1mach,pi,bound, dxrel, rmin, rbig 6 6 real(kind(1.d0)) bern,d1mach,pi,bound, dxrel, rmin, rbig
real(kind(1.d0)) x,y,cabsz 7 7 real(kind(1.d0)) x,y,cabsz
integer nterm,ndx,n,i 8 8 integer nterm,ndx,n,i
external zcot, d1mach 9 9 external zcot, d1mach
c 10 10 c
data bern( 1) / .8333333333 3333333 e-1 / 11 11 data bern( 1) / .8333333333 3333333 d-1 /
data bern( 2) / -.8333333333 3333333 e-2 / 12 12 data bern( 2) / -.8333333333 3333333 d-2 /
data bern( 3) / .3968253968 2539683 e-2 / 13 13 data bern( 3) / .3968253968 2539683 d-2 /
data bern( 4) / -.4166666666 6666667 e-2 / 14 14 data bern( 4) / -.4166666666 6666667 d-2 /
data bern( 5) / .7575757575 7575758 e-2 / 15 15 data bern( 5) / .7575757575 7575758 d-2 /
data bern( 6) / -.2109279609 2796093 e-1 / 16 16 data bern( 6) / -.2109279609 2796093 d-1 /
data bern( 7) / .8333333333 3333333 e-1 / 17 17 data bern( 7) / .8333333333 3333333 d-1 /
data bern( 8) / -.4432598039 2156863 e0 / 18 18 data bern( 8) / -.4432598039 2156863 d0 /
data bern( 9) / .3053954330 2701197 e1 / 19 19 data bern( 9) / .3053954330 2701197 d1 /
data bern(10) / -.2645621212 1212121 e2 / 20 20 data bern(10) / -.2645621212 1212121 d2 /
data bern(11) / .2814601449 2753623 e3 / 21 21 data bern(11) / .2814601449 2753623 d3 /
data bern(12) / -.3454885393 7728938 e4 / 22 22 data bern(12) / -.3454885393 7728938 d4 /
data bern(13) / .5482758333 3333333 e5 / 23 23 data bern(13) / .5482758333 3333333 d5 /
c 24 24 c
data pi / 3.141592653 589793 e0 / 25 25 data pi / 3.141592653 589793 d0 /
data nterm, bound, dxrel, rmin, rbig / 0, 4*0.0 / 26 26 data nterm, bound, dxrel, rmin, rbig / 0d0, 4*0.0d0 /
c 27 27 c
if (nterm.ne.0) go to 10 28 28 if (nterm.ne.0) go to 10
nterm = -0.30*log(d1mach(3)) 29 29 nterm = -0.30*log(d1mach(3))
c maybe bound = n*(0.1*eps)**(-1/(2*n-1)) / (pi*exp(1)) 30 30 c maybe bound = n*(0.1*eps)**(-1/(2*n-1)) / (pi*exp(1))
bound = 0.1171*dble(nterm) * 31 31 bound = 0.1171*dble(nterm) *
1 (0.1*d1mach(3))**(-1.0/(2.0*dble(nterm)-1.0)) 32 32 1 (0.1*d1mach(3))**(-1.0/(2.0*dble(nterm)-1.0))
dxrel = sqrt(d1mach(4)) 33 33 dxrel = sqrt(d1mach(4))
rmin = exp (max(log(d1mach(1)), -log(d1mach(2))) + 0.011 ) 34 34 rmin = exp (max(log(d1mach(1)), -log(d1mach(2))) + 0.011 )
rbig = 1.0/d1mach(3) 35 35 rbig = 1.0/d1mach(3)
c 36 36 c
10 z = zin 37 37 10 z = zin
x = real(z) 38 38 x = real(z)
y = aimag(z) 39 39 y = aimag(z)
if (y.lt.0.0) z = conjg(z) 40 40 if (y.lt.0.0) z = conjg(z)
c 41 41 c
corr = (0.0, 0.0) 42 42 corr = (0.0, 0.0)
cabsz = abs(z) 43 43 cabsz = abs(z)
if (x.ge.0.0 .and. cabsz.gt.bound) go to 50 44 44 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 45 45 if (x.lt.0.0 .and. abs(y).gt.bound) go to 50
c 46 46 c
if (cabsz.lt.bound) go to 20 47 47 if (cabsz.lt.bound) go to 20
c 48 48 c
c use the reflection formula for real(z) negative, cabs(z) large, and 49 49 c use the reflection formula for real(z) negative, cabs(z) large, and
c abs(aimag(y)) small. 50 50 c abs(aimag(y)) small.
c 51 51 c
corr = -pi*zcot(pi*z) 52 52 corr = -pi*zcot(pi*z)
z = 1.0 - z 53 53 z = 1.0 - z
go to 50 54 54 go to 50
c 55 55 c
c use the recursion relation for cabs(z) small. 56 56 c use the recursion relation for cabs(z) small.
c 57 57 c
20 if (cabsz.lt.rmin) call seteru ( 58 58 20 if (cabsz.lt.rmin) call seteru (
1 56hzpsi zpsi called with z so near 0 that zpsi overflows, 59 59 1 56hzpsi zpsi called with z so near 0 that zpsi overflows,
1 56, 2, 2) 60 60 1 56, 2, 2)
c 61 61 c
if (x.ge.(-0.5) .or. abs(y).gt.dxrel) go to 30 62 62 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 ( 63 63 if ( abs((z-aint(x-0.5))/x).lt.dxrel) call seteru (
1 68hzpsi answer lt half precision because z too near negative 64 64 1 68hzpsi answer lt half precision because z too near negative
2integer, 68, 1, 1) 65 65 2integer, 68, 1, 1)
if (y.eq.0.0 .and. x.eq.aint(x)) call seteru ( 66 66 if (y.eq.0.0 .and. x.eq.aint(x)) call seteru (
1 31hzpsi z is a negative integer, 31, 3, 2) 67 67 1 31hzpsi z is a negative integer, 31, 3, 2)
c 68 68 c
30 n = sqrt(bound**2-y**2) - x + 1.0 69 69 30 n = sqrt(bound**2-y**2) - x + 1.0
do 40 i=1,n 70 70 do 40 i=1,n
corr = corr - 1.0/z 71 71 corr = corr - 1.0/z
z = z + 1.0 72 72 z = z + 1.0
40 continue 73 73 40 continue
c 74 74 c
c now evaluate the asymptotic series for suitably large z. 75 75 c now evaluate the asymptotic series for suitably large z.
c 76 76 c
50 if (cabsz.gt.rbig) zpsi = log(z) + corr 77 77 50 if (cabsz.gt.rbig) zpsi = log(z) + corr
if (cabsz.gt.rbig) go to 70 78 78 if (cabsz.gt.rbig) go to 70
c 79 79 c
zpsi = (0.0, 0.0) 80 80 zpsi = (0.0, 0.0)
z2inv = 1.0/z**2 81 81 z2inv = 1.0/z**2
do 60 i=1,nterm 82 82 do 60 i=1,nterm
ndx = nterm + 1 - i 83 83 ndx = nterm + 1 - i
zpsi = bern(ndx) + z2inv*zpsi 84 84 zpsi = bern(ndx) + z2inv*zpsi
60 continue 85 85 60 continue
zpsi = log(z) - 0.5/z - zpsi*z2inv + corr 86 86 zpsi = log(z) - 0.5/z - zpsi*z2inv + corr
c 87 87 c
70 if (y.lt.0.0) zpsi = conjg(zpsi) 88 88 70 if (y.lt.0.0) zpsi = conjg(zpsi)
c 89 89 c
return 90 90 return
end 91 91 end
92 92
complex(kind(1.d0)) function zsinh (z) 1 1 complex(kind(1.d0)) function zsinh (z)
implicit none 2 2 implicit none
c april 1977 version. w. fullerton, c3, los alamos scientific lab. 3 3 c april 1977 version. w. fullerton, c3, los alamos scientific lab.
complex(kind(1.d0)) z, ci 4 4 complex(kind(1.d0)) z, ci
data ci /(0.,1.)/ 5 5 data ci /(0.d0,1.d0)/
c 6 6 c
zsinh = -ci*sin(ci*z) 7 7 zsinh = -ci*sin(ci*z)
c 8 8 c
return 9 9 return
end 10 10 end
11 11
complex(kind(1.d0)) function ztan (z) 1 1 complex(kind(1.d0)) function ztan (z)
implicit none 2 2 implicit none
c march 1979 edition. w. fullerton, c3, los alamos scientific lab. 3 3 c march 1979 edition. w. fullerton, c3, los alamos scientific lab.
complex(kind(1.d0)) z 4 4 complex(kind(1.d0)) z
real(kind(1.d0)) d1mach,eps, xmax, ylarge, ybig, ymin,sn2x,den 5 5 real(kind(1.d0)) d1mach,eps, xmax, ylarge, ybig, ymin,sn2x,den
real(kind(1.d0)) x,y,x2,y2 6 6 real(kind(1.d0)) x,y,x2,y2
integer irold,irold2 7 7 integer irold,irold2
external d1mach 8 8 external d1mach
data eps, xmax, ylarge, ybig, ymin / 4*0.0, 1.50 / 9 9 data eps, xmax, ylarge, ybig, ymin / 4*0.0d0, 1.50d0 /
c 10 10 c
if (eps.ne.0.0) go to 10 11 11 if (eps.ne.0.0) go to 10
eps = d1mach(4) 12 12 eps = d1mach(4)
xmax = 1.0/eps 13 13 xmax = 1.0/eps
ylarge = -0.5*log(0.5*d1mach(3)) 14 14 ylarge = -0.5*log(0.5*d1mach(3))
ybig = -0.5*log(0.5*sqrt(d1mach(3))) 15 15 ybig = -0.5*log(0.5*sqrt(d1mach(3)))
c 16 16 c
10 x = real(z) 17 17 10 x = real(z)
y = aimag(z) 18 18 y = aimag(z)
if (abs(y).gt.ylarge) go to 30 19 19 if (abs(y).gt.ylarge) go to 30
c 20 20 c
x2 = 2.0*x 21 21 x2 = 2.0*x
y2 = 2.0*y 22 22 y2 = 2.0*y
if (abs(x2).gt.xmax) go to 20 23 23 if (abs(x2).gt.xmax) go to 20
c 24 24 c
call entsrc (irold, 1) 25 25 call entsrc (irold, 1)
sn2x = sin(x2) 26 26 sn2x = sin(x2)
call erroff 27 27 call erroff
den = cos(x2) + cosh(y2) 28 28 den = cos(x2) + cosh(y2)
call erroff 29 29 call erroff
call entsrc (irold2, irold) 30 30 call entsrc (irold2, irold)
c 31 31 c
if (den.lt.x*x2*eps*eps) call seteru ( 72hztan tan is nearly si 32 32 if (den.lt.x*x2*eps*eps) call seteru ( 72hztan tan is nearly si
1ngular, x is near pi/2 or 3*pi/2 and y is near 0, 72, 3, 2) 33 33 1ngular, x is near pi/2 or 3*pi/2 and y is near 0, 72, 3, 2)
if (den.lt.x*x2*eps) call seteru (74hztan answer lt half precis 34 34 if (den.lt.x*x2*eps) call seteru (74hztan answer lt half precis
1ion, x is near pi/2 or 3*pi/2 and y is near 0, 74, 1, 1) 35 35 1ion, x is near pi/2 or 3*pi/2 and y is near 0, 74, 1, 1)
c 36 36 c
ztan = cmplx (sn2x/den, sinh(y2)/den, kind(1.d0)) 37 37 ztan = cmplx (sn2x/den, sinh(y2)/den, kind(1.d0))
return 38 38 return
c 39 39 c
20 if (abs(y).lt.ymin) call seteru (75hztan answer would have no p 40 40 20 if (abs(y).lt.ymin) call seteru (75hztan answer would have no p
1recision, abs(x) is very big and abs(y) small, 75, 4, 2) 41 41 1recision, abs(x) is very big and abs(y) small, 75, 4, 2)
if (abs(y).lt.ybig) call seteru (69hztan answer lt half precisi 42 42 if (abs(y).lt.ybig) call seteru (69hztan answer lt half precisi
1on, abs(x) is very big and abs(y) small, 69, 2, 1) 43 43 1on, abs(x) is very big and abs(y) small, 69, 2, 1)
c 44 44 c
ztan = cmplx (0.0, tanh(y2), kind(1.d0)) 45 45 ztan = cmplx (0.0, tanh(y2), kind(1.d0))
return 46 46 return
c 47 47 c
30 ztan = cmplx (0.0, -sign(1.0d0, y), kind(1.d0)) 48 48 30 ztan = cmplx (0.0, -sign(1.0d0, y), kind(1.d0))
return 49 49 return
c 50 50 c
end 51 51 end
52 52
complex(kind(1.d0)) function ztanh (z) 1 1 complex(kind(1.d0)) function ztanh (z)
c april 1977 version. w. fullerton, c3, los alamos scientific lab. 2 2 c april 1977 version. w. fullerton, c3, los alamos scientific lab.
complex(kind(1.d0)) z, ci, ztan 3 3 complex(kind(1.d0)) z, ci, ztan
external ztan 4 4 external ztan
data ci /(0.,1.)/ 5 5 data ci /(0.d0,1.d0)/
c 6 6 c
ztanh = -ci*ztan(ci*z) 7 7 ztanh = -ci*ztan(ci*z)
c 8 8 c
return 9 9 return
end 10 10 end
11 11