Commit 0c3098aed95e9e568f02887c9ae6145422e8d857
1 parent
e711bb807c
Exists in
master
and in
2 other branches
ChW 02/2010 for type definitions and dcmplx replacements in accordance to ANSI
git-svn-id: https://lxsd.femto-st.fr/svn/fvn@69 b657c933-2333-4658-acf2-d3c7c2708721
Showing 27 changed files with 107 additions and 100 deletions Inline Diff
- fvn_fnlib/z0lgmc.f
- fvn_fnlib/z8lgmc.f
- fvn_fnlib/z9lgmc.f
- fvn_fnlib/z9ln2r.f
- fvn_fnlib/zacos.f
- fvn_fnlib/zacosh.f
- fvn_fnlib/zarg.f
- fvn_fnlib/zasin.f
- fvn_fnlib/zasinh.f
- fvn_fnlib/zatan.f
- fvn_fnlib/zatan2.f
- fvn_fnlib/zatanh.f
- fvn_fnlib/zbeta.f
- fvn_fnlib/zcbrt.f
- fvn_fnlib/zcosh.f
- fvn_fnlib/zcot.f
- fvn_fnlib/zexprl.f
- fvn_fnlib/zgamma.f
- fvn_fnlib/zgamr.f
- fvn_fnlib/zlbeta.f
- fvn_fnlib/zlngam.f
- fvn_fnlib/zlnrel.f
- fvn_fnlib/zlog10.f
- fvn_fnlib/zpsi.f
- fvn_fnlib/zsinh.f
- fvn_fnlib/ztan.f
- fvn_fnlib/ztanh.f
fvn_fnlib/z0lgmc.f
| complex(8) function z0lgmc (z) | 1 | 1 | complex(kind(1.d0)) function z0lgmc (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 (z+0.5)*clog((z+1.0)/z) - 1.0 with relative error accuracy. | 5 | 5 | c evaluate (z+0.5)*clog((z+1.0)/z) - 1.0 with relative error accuracy. | |
| c let q = 1.0/z so that | 6 | 6 | c let q = 1.0/z so that | |
| c (z+0.5)*clog(1+1/z) - 1 = (z+0.5)*(clog(1+q) - q + q*q/2) - q*q/4 | 7 | 7 | c (z+0.5)*clog(1+1/z) - 1 = (z+0.5)*(clog(1+q) - q + q*q/2) - q*q/4 | |
| c = (z+0.5)*q**3*c9ln2r(q) - q**2/4, | 8 | 8 | c = (z+0.5)*q**3*c9ln2r(q) - q**2/4, | |
| c where c9ln2r is (clog(1+q) - q + 0.5*q**2) / q**3. | 9 | 9 | c where c9ln2r is (clog(1+q) - q + 0.5*q**2) / q**3. | |
| c | 10 | 10 | c | |
| complex(8) z, q, z9ln2r | 11 | 11 | complex(kind(1.d0)) z, q, z9ln2r | |
| real(8) zabsz | 12 | 12 | real(kind(1.d0)) zabsz | |
| 13 | 13 | |||
| external z9ln2r | 14 | 14 | external z9ln2r | |
| c | 15 | 15 | c | |
| zabsz = abs(z) | 16 | 16 | zabsz = abs(z) | |
| c | 17 | 17 | c | |
| q = 1.0/z | 18 | 18 | q = 1.0/z | |
| if (zabsz.le.1.23) z0lgmc = (z+0.5)*log(1.0+q) - 1.0 | 19 | 19 | if (zabsz.le.1.23) z0lgmc = (z+0.5)*log(1.0+q) - 1.0 | |
| if (zabsz.gt.1.23) z0lgmc = ((1.+.5*q)*z9ln2r(q) - .25) * q**2 | 20 | 20 | if (zabsz.gt.1.23) z0lgmc = ((1.+.5*q)*z9ln2r(q) - .25) * q**2 | |
| c | 21 | 21 | c | |
| return | 22 | 22 | return | |
| end | 23 | 23 | end | |
| 24 | 24 | |||
fvn_fnlib/z8lgmc.f
| complex(8) 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(8) zin, z, corr, z9lgmc, z0lgmc, zexp, zlnrel | 12 | 12 | complex(kind(1.d0)) zin, z, corr, z9lgmc, z0lgmc, zexp, zlnrel | |
| real(8) 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 | |
| 15 | c complex(kind(1.d0)) tmp_arg | |||
| external z0lgmc, z9lgmc, zlnrel, | 15 | 16 | external z0lgmc, z9lgmc, zlnrel, | |
| 1 d1mach | 16 | 17 | 1 d1mach | |
| data pi / 3.1415926535 8979324d0 / | 17 | 18 | data pi / 3.1415926535 8979324d0 / | |
| data bound, sqeps, eps / 3*0.0 / | 18 | 19 | data bound, sqeps, eps / 3*0.0 / | |
| c | 19 | 20 | c | |
| if (bound.ne.0.0) go to 10 | 20 | 21 | if (bound.ne.0.0) go to 10 | |
| nterm = -0.30*log(d1mach(3)) | 21 | 22 | nterm = -0.30*log(d1mach(3)) | |
| bound = 0.1170*dble(nterm)* | 22 | 23 | bound = 0.1170*dble(nterm)* | |
| 1 (0.1*d1mach(3))**(-1.0/(2.0*dble(nterm)-1.0)) | 23 | 24 | 1 (0.1*d1mach(3))**(-1.0/(2.0*dble(nterm)-1.0)) | |
| eps = 10.0*d1mach(4) | 24 | 25 | eps = 10.0*d1mach(4) | |
| sqeps = sqrt (d1mach(4)) | 25 | 26 | sqeps = sqrt (d1mach(4)) | |
| c | 26 | 27 | c | |
| 10 z = zin | 27 | 28 | 10 z = zin | |
| x = real(z) | 28 | 29 | x = real(z) | |
| y = aimag(z) | 29 | 30 | y = aimag(z) | |
| if (x.lt.0.0) z = -z | 30 | 31 | if (x.lt.0.0) z = -z | |
| c | 31 | 32 | c | |
| corr = (0.0, 0.0) | 32 | 33 | corr = (0.0, 0.0) | |
| if (x.gt.bound .or. abs(y).ge.bound) go to 30 | 33 | 34 | if (x.gt.bound .or. abs(y).ge.bound) go to 30 | |
| absz = abs(z) | 34 | 35 | absz = abs(z) | |
| if (absz.ge.bound) go to 30 | 35 | 36 | if (absz.ge.bound) go to 30 | |
| c | 36 | 37 | c | |
| if (x.eq.0.0 .and. y.eq.0.0) call seteru ( | 37 | 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, | 38 | 39 | 1 57hz8lgmc log gamma correction term is infinite for z = 0.0, | |
| 2 57, 3, 2) | 39 | 40 | 2 57, 3, 2) | |
| c | 40 | 41 | c | |
| n = sqrt (bound**2-y**2) - abs(x) + 1.0 | 41 | 42 | n = sqrt (bound**2-y**2) - abs(x) + 1.0 | |
| do 20 i=1,n | 42 | 43 | do 20 i=1,n | |
| corr = corr + z0lgmc(z) | 43 | 44 | corr = corr + z0lgmc(z) | |
| z = z + 1.0 | 44 | 45 | z = z + 1.0 | |
| 20 continue | 45 | 46 | 20 continue | |
| c | 46 | 47 | c | |
| 30 z8lgmc = z9lgmc(z) + corr | 47 | 48 | 30 z8lgmc = z9lgmc(z) + corr | |
| if (x.gt.0.0) return | 48 | 49 | if (x.gt.0.0) return | |
| c | 49 | 50 | c | |
| call entsrc (irold, 1) | 50 | 51 | call entsrc (irold, 1) | |
| test = 1.0 | 51 | 52 | test = 1.0 | |
| if (x.lt.(-0.5)) test = abs ((z-aint(x-0.5))/z) | 52 | 53 | if (x.lt.(-0.5)) test = abs ((z-aint(x-0.5))/z) | |
| if (test.lt.eps) call seteru ( | 53 | 54 | if (test.lt.eps) call seteru ( | |
| 1 31hz8lgmc z is a negative integer, 31, 2, 2) | 54 | 55 | 1 31hz8lgmc z is a negative integer, 31, 2, 2) | |
| c | 55 | 56 | c | |
| z = zin | 56 | 57 | z = zin | |
| if (y.lt.0.0) z = conjg(z) | 57 | 58 | if (y.lt.0.0) z = conjg(z) | |
| corr = -dcmplx(0.0,pi) + zlnrel(-exp(dcmplx(0.,2.*pi)*z)) | 58 | 59 | c tmp_arg = -exp(cmplx(0.,2.*pi,kind(1.d0))*z) | |
| 60 | corr = -cmplx(0.0,pi,kind(1.d0)) + | |||
| 61 | 1 zlnrel(-exp(cmplx(0.,2.*pi,kind(1.d0))*z)) | |||
| if (y.lt.0.0) corr = conjg(corr) | 59 | 62 | if (y.lt.0.0) corr = conjg(corr) | |
| c | 60 | 63 | c | |
| z8lgmc = corr - z8lgmc | 61 | 64 | z8lgmc = corr - z8lgmc | |
| c | 62 | 65 | c | |
| call erroff | 63 | 66 | call erroff | |
| call entsrc (ir, irold) | 64 | 67 | call entsrc (ir, irold) | |
| if (test.lt.sqeps) call seteru ( 63hc8lgmc answer lt half precisi | 65 | 68 | if (test.lt.sqeps) call seteru ( 63hc8lgmc answer lt half precisi | |
| 1on, z too near a negative integer, 63, 1, 1) | 66 | 69 | 1on, z too near a negative integer, 63, 1, 1) | |
| c | 67 | 70 | c | |
| return | 68 | 71 | return | |
| end | 69 | 72 | end | |
| 70 | 73 | |||
fvn_fnlib/z9lgmc.f
| complex(8) 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(8) zin, z, z2inv | 10 | 10 | complex(kind(1.d0)) zin, z, z2inv | |
| real(8) 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 3333333e0 / | |
| data bern( 2) / -.002777777777 7777778e0 / | 17 | 17 | data bern( 2) / -.002777777777 7777778e0 / | |
| data bern( 3) / .0007936507936 5079365e0 / | 18 | 18 | data bern( 3) / .0007936507936 5079365e0 / | |
| data bern( 4) / -.0005952380952 3809524e0 / | 19 | 19 | data bern( 4) / -.0005952380952 3809524e0 / | |
| data bern( 5) / .0008417508417 5084175e0 / | 20 | 20 | data bern( 5) / .0008417508417 5084175e0 / | |
| data bern( 6) / -.001917526917 5269175e0 / | 21 | 21 | data bern( 6) / -.001917526917 5269175e0 / | |
| data bern( 7) / .006410256410 2564103e0 / | 22 | 22 | data bern( 7) / .006410256410 2564103e0 / | |
| data bern( 8) / -.02955065359 4771242e0 / | 23 | 23 | data bern( 8) / -.02955065359 4771242e0 / | |
| data bern( 9) / .1796443723 6883057e0 / | 24 | 24 | data bern( 9) / .1796443723 6883057e0 / | |
| data bern(10) / -1.392432216 9059011e0 / | 25 | 25 | data bern(10) / -1.392432216 9059011e0 / | |
| data bern(11) / 13.40286404 4168392e0 / | 26 | 26 | data bern(11) / 13.40286404 4168392e0 / | |
| c | 27 | 27 | c | |
| data nterm, bound, xbig, xmax / 0, 3*0.0 / | 28 | 28 | data nterm, bound, xbig, xmax / 0, 3*0.0 / | |
| 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 | |||
fvn_fnlib/z9ln2r.f
| complex(8) function z9ln2r (z) | 1 | 1 | complex(kind(1.d0)) function z9ln2r (z) | |
| 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 evaluate clog(1+z) from 2-nd order with relative error accuracy so | 5 | 5 | c evaluate clog(1+z) from 2-nd order with relative error accuracy so | |
| c that clog(1+z) = z - z**2/2 + z**3*c9ln2r(z). | 6 | 6 | c that clog(1+z) = z - z**2/2 + z**3*c9ln2r(z). | |
| c | 7 | 7 | c | |
| c now clog(1+z) = 0.5*alog(1+2*x+cabs(z)**2) + i*carg(1+z), | 8 | 8 | c now clog(1+z) = 0.5*alog(1+2*x+cabs(z)**2) + i*carg(1+z), | |
| c where x = real(z) and y = aimag(z). | 9 | 9 | c where x = real(z) and y = aimag(z). | |
| c we find | 10 | 10 | c we find | |
| c z**3 * c9ln2r(z) = -x*cabs(z)**2 - 0.25*cabs(z)**4 | 11 | 11 | c z**3 * c9ln2r(z) = -x*cabs(z)**2 - 0.25*cabs(z)**4 | |
| c + (2*x+cabs(z)**2)**3 * r9ln2r(2*x+cabs(z)**2) | 12 | 12 | c + (2*x+cabs(z)**2)**3 * r9ln2r(2*x+cabs(z)**2) | |
| c + i * (carg(1+z) + (x-1)*y) | 13 | 13 | c + i * (carg(1+z) + (x-1)*y) | |
| c the imaginary part must be evaluated carefully as | 14 | 14 | c the imaginary part must be evaluated carefully as | |
| c (atan(y/(1+x)) - y/(1+x)) + y/(1+x) - (1-x)*y | 15 | 15 | c (atan(y/(1+x)) - y/(1+x)) + y/(1+x) - (1-x)*y | |
| c = (y/(1+x))**3 * r9atn1(y/(1+x)) + x**2*y/(1+x) | 16 | 16 | c = (y/(1+x))**3 * r9atn1(y/(1+x)) + x**2*y/(1+x) | |
| c | 17 | 17 | c | |
| c now we divide through by z**3 carefully. write | 18 | 18 | c now we divide through by z**3 carefully. write | |
| c 1/z**3 = (x-i*y)/cabs(z)**3 * (1/cabs(z)**3) | 19 | 19 | c 1/z**3 = (x-i*y)/cabs(z)**3 * (1/cabs(z)**3) | |
| c then c9ln2r(z) = ((x-i*y)/cabs(z))**3 * (-x/cabs(z) - cabs(z)/4 | 20 | 20 | c then c9ln2r(z) = ((x-i*y)/cabs(z))**3 * (-x/cabs(z) - cabs(z)/4 | |
| c + 0.5*((2*x+cabs(z)**2)/cabs(z))**3 * r9ln2r(2*x+cabs(z)**2) | 21 | 21 | c + 0.5*((2*x+cabs(z)**2)/cabs(z))**3 * r9ln2r(2*x+cabs(z)**2) | |
| c + i*y/(cabs(z)*(1+x)) * ((x/cabs(z))**2 + | 22 | 22 | c + i*y/(cabs(z)*(1+x)) * ((x/cabs(z))**2 + | |
| c + (y/(cabs(z)*(1+x)))**2 * r9atn1(y/(1+x)) ) ) | 23 | 23 | c + (y/(cabs(z)*(1+x)))**2 * r9atn1(y/(1+x)) ) ) | |
| c | 24 | 24 | c | |
| c if we let xz = x/cabs(z) and yz = y/cabs(z) we may write | 25 | 25 | c if we let xz = x/cabs(z) and yz = y/cabs(z) we may write | |
| c c9ln2r(z) = (xz-i*yz)**3 * (-xz - cabs(z)/4 | 26 | 26 | c c9ln2r(z) = (xz-i*yz)**3 * (-xz - cabs(z)/4 | |
| c + 0.5*(2*xz+cabs(z))**3 * r9ln2r(2*x+cabs(z)**2) | 27 | 27 | c + 0.5*(2*xz+cabs(z))**3 * r9ln2r(2*x+cabs(z)**2) | |
| c + i*yz/(1+x) * (xz**2 + (yz/(1+x))**2*r9atn1(y/(1+x)) )) | 28 | 28 | c + i*yz/(1+x) * (xz**2 + (yz/(1+x))**2*r9atn1(y/(1+x)) )) | |
| c | 29 | 29 | c | |
| complex(8) z | 30 | 30 | complex(kind(1.d0)) z | |
| real(8) x,y,xz,yz,cabsz,d9ln2r,d9atn1,arg,rpart,aipart,y1x | 31 | 31 | real(kind(1.d0)) x,y,xz,yz,cabsz,d9ln2r,d9atn1,arg,rpart,aipart, | |
| 32 | 1 y1x | |||
| 32 | 33 | |||
| external d9atn1, d9ln2r | 33 | 34 | external d9atn1, d9ln2r | |
| c | 34 | 35 | c | |
| x = real (z) | 35 | 36 | x = real (z) | |
| y = aimag (z) | 36 | 37 | y = aimag (z) | |
| c | 37 | 38 | c | |
| cabsz = abs(z) | 38 | 39 | cabsz = abs(z) | |
| if (cabsz.gt.0.8125) go to 20 | 39 | 40 | if (cabsz.gt.0.8125) go to 20 | |
| c | 40 | 41 | c | |
| z9ln2r = dcmplx (1.0/3.0, 0.0) | 41 | 42 | z9ln2r = cmplx (1.0/3.0, 0.0, kind(1.d0)) | |
| if (cabsz.eq.0.0) return | 42 | 43 | if (cabsz.eq.0.0) return | |
| c | 43 | 44 | c | |
| xz = x/cabsz | 44 | 45 | xz = x/cabsz | |
| yz = y/cabsz | 45 | 46 | yz = y/cabsz | |
| c | 46 | 47 | c | |
| arg = 2.0*xz + cabsz | 47 | 48 | arg = 2.0*xz + cabsz | |
| rpart = 0.5*arg**3*d9ln2r(cabsz*arg) - xz - 0.25*cabsz | 48 | 49 | rpart = 0.5*arg**3*d9ln2r(cabsz*arg) - xz - 0.25*cabsz | |
| y1x = yz/(1.0+x) | 49 | 50 | y1x = yz/(1.0+x) | |
| aipart = y1x * (xz**2 + y1x**2*d9atn1(cabsz*y1x) ) | 50 | 51 | aipart = y1x * (xz**2 + y1x**2*d9atn1(cabsz*y1x) ) | |
| c | 51 | 52 | c | |
| z9ln2r = dcmplx(xz,-yz)**3 * dcmplx(rpart,aipart) | 52 | 53 | z9ln2r = cmplx(xz,-yz,kind(1.d0))**3 * | |
| 54 | 1 cmplx(rpart,aipart,kind(1.d0)) | |||
| return | 53 | 55 | return | |
| c | 54 | 56 | c | |
| 20 z9ln2r = (log(1.0+z) - z*(1.0-0.5*z)) / z**3 | 55 | 57 | 20 z9ln2r = (log(1.0+z) - z*(1.0-0.5*z)) / z**3 | |
| return | 56 | 58 | return | |
| c | 57 | 59 | c | |
| end | 58 | 60 | end | |
| 59 | 61 | |||
fvn_fnlib/zacos.f
| complex(8) function zacos (z) | 1 | 1 | complex(kind(1.d0)) function zacos (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(8) z, zasin | 4 | 4 | complex(kind(1.d0)) z, zasin | |
| real(8) pi2 | 5 | 5 | real(kind(1.d0)) pi2 | |
| external zasin | 6 | 6 | external zasin | |
| data pi2 /1.5707963267 9489661923d0/ | 7 | 7 | data pi2 /1.5707963267 9489661923d0/ | |
| c | 8 | 8 | c | |
| zacos = pi2 - zasin (z) | 9 | 9 | zacos = pi2 - zasin (z) | |
| c | 10 | 10 | c | |
| return | 11 | 11 | return | |
| end | 12 | 12 | end | |
| 13 | 13 | |||
fvn_fnlib/zacosh.f
| complex(8) function zacosh (z) | 1 | 1 | complex(kind(1.d0)) function zacosh (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(8) z, ci, zacos | 4 | 4 | complex(kind(1.d0)) z, ci, zacos | |
| external zacos | 5 | 5 | external zacos | |
| data ci /(0.d0,1.d0)/ | 6 | 6 | data ci /(0.d0,1.d0)/ | |
| c | 7 | 7 | c | |
| zacosh = ci*zacos(z) | 8 | 8 | zacosh = ci*zacos(z) | |
| c | 9 | 9 | c | |
| return | 10 | 10 | return | |
| end | 11 | 11 | end | |
| 12 | 12 | |||
fvn_fnlib/zarg.f
| real(8) function zarg (z) | 1 | 1 | real(kind(1.d0)) function zarg (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(8) z | 4 | 4 | complex(kind(1.d0)) z | |
| 5 | 5 | |||
| c | 6 | 6 | c | |
| zarg = 0.0 | 7 | 7 | zarg = 0.0 | |
| if (real(z).ne.0. .or. aimag(z).ne.0.) zarg = | 8 | 8 | if (real(z).ne.0. .or. aimag(z).ne.0.) zarg = | |
| 1 atan2 (aimag(z), real(z)) | 9 | 9 | 1 atan2 (aimag(z), real(z)) | |
| c | 10 | 10 | c | |
| return | 11 | 11 | return | |
| end | 12 | 12 | end | |
| 13 | 13 |
fvn_fnlib/zasin.f
| complex(8) 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(8) zinp, z, z2, sqzp1, ci | 8 | 8 | complex(kind(1.d0)) zinp, z, z2, sqzp1, ci | |
| real(8) 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.,1.)/ | |
| data nterms, rmin / 0, 0.0 / | 16 | 16 | data nterms, rmin / 0, 0.0 / | |
| 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 | |||
fvn_fnlib/zasinh.f
| complex(8) function zasinh (z) | 1 | 1 | complex(kind(1.d0)) function zasinh (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 z, ci, zasin | 4 | 4 | complex z, ci, zasin | |
| external zasin | 5 | 5 | external zasin | |
| data ci /(0.,1.)/ | 6 | 6 | data ci /(0.,1.)/ | |
| c | 7 | 7 | c | |
| zasinh = -ci*zasin (ci*z) | 8 | 8 | zasinh = -ci*zasin (ci*z) | |
| c | 9 | 9 | c | |
| return | 10 | 10 | return | |
| end | 11 | 11 | end | |
| 12 | 12 |
fvn_fnlib/zatan.f
| complex(8) function zatan (z) | 1 | 1 | complex(kind(1.d0)) function zatan (z) | |
| implicit none | 2 | 2 | implicit none | |
| c jan 1984 edition. w. fullerton, los alamos scientific lab. | 3 | 3 | c jan 1984 edition. w. fullerton, los alamos scientific lab. | |
| c | 4 | 4 | c | |
| complex(8) z, z2 | 5 | 5 | complex(kind(1.d0)) z, z2 | |
| real(8) pi2,sqeps,rmin,rmax,r2sml,r,d1mach,x,y,r2,xans,yans | 6 | 6 | real(kind(1.d0)) pi2,sqeps,rmin,rmax,r2sml,r,x,y,r2,xans,yans | |
| external d1mach | 7 | 7 | real(kind(1.d0)) d1mach | |
| integer nterms,i,twoi | 8 | 8 | integer nterms,i,twoi | |
| data pi2 /1.5707963267 9489661923d0 / | 9 | 9 | data pi2 /1.5707963267 9489661923d0 / | |
| data nterms, sqeps, rmin, rmax, r2sml / 0, 4*0.0 / | 10 | 10 | data nterms, sqeps, rmin, rmax, r2sml / 0, 4*0.0d0 / | |
| 11 | external d1mach | |||
| 12 | ||||
| c | 11 | 13 | c | |
| if (nterms.ne.0) go to 10 | 12 | 14 | if (nterms.ne.0) go to 10 | |
| c nterms = alog(eps)/alog(rbnd) where rbnd = 0.1 | 13 | 15 | c nterms = alog(eps)/alog(rbnd) where rbnd = 0.1 | |
| nterms = -0.4343*log(d1mach(3)) + 1.0 | 14 | 16 | nterms = -0.4343*log(d1mach(3)) + 1.0 | |
| sqeps = sqrt (d1mach(4)) | 15 | 17 | sqeps = sqrt (d1mach(4)) | |
| rmin = sqrt (3.0*d1mach(3)) | 16 | 18 | rmin = sqrt (3.0*d1mach(3)) | |
| rmax = 1.0/d1mach(3) | 17 | 19 | rmax = 1.0/d1mach(3) | |
| r2sml = 10.0*d1mach (4) | 18 | 20 | r2sml = 10.0*d1mach (4) | |
| c | 19 | 21 | c | |
| 10 r = abs(z) | 20 | 22 | 10 r = abs(z) | |
| if (r.gt.0.1) go to 30 | 21 | 23 | if (r.gt.0.1) go to 30 | |
| c | 22 | 24 | c | |
| zatan = z | 23 | 25 | zatan = z | |
| if (r.lt.rmin) return | 24 | 26 | if (r.lt.rmin) return | |
| c | 25 | 27 | c | |
| zatan = (0.0, 0.0) | 26 | 28 | zatan = (0.0, 0.0) | |
| z2 = z*z | 27 | 29 | z2 = z*z | |
| do 20 i=1,nterms | 28 | 30 | do 20 i=1,nterms | |
| twoi = 2*(nterms-i) + 1 | 29 | 31 | twoi = 2*(nterms-i) + 1 | |
| zatan = 1.0/twoi - z2*zatan | 30 | 32 | zatan = 1.0/twoi - z2*zatan | |
| 20 continue | 31 | 33 | 20 continue | |
| zatan = z*zatan | 32 | 34 | zatan = z*zatan | |
| return | 33 | 35 | return | |
| c | 34 | 36 | c | |
| 30 if (r.gt.rmax) go to 50 | 35 | 37 | 30 if (r.gt.rmax) go to 50 | |
| x = real(z) | 36 | 38 | x = real(z) | |
| y = aimag(z) | 37 | 39 | y = aimag(z) | |
| r2 = r*r | 38 | 40 | r2 = r*r | |
| c | 39 | 41 | c | |
| if (abs(r2-1.0).gt.sqeps) go to 40 | 40 | 42 | 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 | 41 | 43 | if (r2+2.0*y+1.0.lt.r2sml .or. r2-2.0*y+1.0.lt.r2sml) call | |
| 1 seteru ( | 42 | 44 | 1 seteru ( | |
| 2 55hzatan no precision because z is too close to +i or -i, | 43 | 45 | 2 55hzatan no precision because z is too close to +i or -i, | |
| 3 55, 2, 2) | 44 | 46 | 3 55, 2, 2) | |
| if (abs(dcmplx(1.0,0.0)+z*z).lt.sqeps) call seteru ( | 45 | 47 | if (abs(cmplx(1.0d0,0.0d0,kind(1.d0))+z*z).lt.sqeps) call seteru ( | |
| 1 50hzatan answer lt half precision, z**2 close to -1, 50, 1, 1) | 46 | 48 | 1 50hzatan answer lt half precision, z**2 close to -1, 50, 1, 1) | |
| c | 47 | 49 | c | |
| 40 xans = 0.5*atan2 (2.0*x, 1.0-r2) | 48 | 50 | 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)) | 49 | 51 | yans = 0.25*log((r2+2.0*y+1.0)/(r2-2.0*y+1.0)) | |
| zatan = dcmplx (xans, yans) | 50 | 52 | zatan = cmplx (xans, yans, kind(1.d0)) | |
| return | 51 | 53 | return | |
| c | 52 | 54 | c | |
| 50 zatan = dcmplx (pi2, 0.0) | 53 | 55 | 50 zatan = cmplx (pi2, 0.0d0, kind(1.d0)) | |
| if (real(z).lt.0.0) zatan = dcmplx (-pi2, 0.0) | 54 | 56 | if (real(z).lt.0.0) zatan = cmplx (-pi2, 0.0d0, kind(1.d0)) | |
| return | 55 | 57 | return | |
| c | 56 | 58 | c | |
| end | 57 | 59 | end | |
| 58 | 60 | |||
fvn_fnlib/zatan2.f
| complex(8) function zatan2 (csn, ccs) | 1 | 1 | complex(kind(1.d0)) function zatan2 (csn, ccs) | |
| 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(8) csn, ccs, zatan | 4 | 4 | complex(kind(1.d0)) csn, ccs, zatan | |
| real(8) pi | 5 | 5 | real(kind(1.d0)) pi | |
| external zatan | 6 | 6 | external zatan | |
| data pi / 3.1415926535 8979323846d0 / | 7 | 7 | data pi / 3.1415926535 8979323846d0 / | |
| c | 8 | 8 | c | |
| if (abs(ccs).eq.0.) go to 10 | 9 | 9 | if (abs(ccs).eq.0.) go to 10 | |
| c | 10 | 10 | c | |
| zatan2 = zatan (csn/ccs) | 11 | 11 | zatan2 = zatan (csn/ccs) | |
| if (real(ccs).lt.0.) zatan2 = zatan2 + pi | 12 | 12 | if (real(ccs).lt.0.) zatan2 = zatan2 + pi | |
| if (real(zatan2).gt.pi) zatan2 = zatan2 - 2.0*pi | 13 | 13 | if (real(zatan2).gt.pi) zatan2 = zatan2 - 2.0*pi | |
| return | 14 | 14 | return | |
| c | 15 | 15 | c | |
| 10 if (abs(csn).eq.0.) call seteru ( | 16 | 16 | 10 if (abs(csn).eq.0.) call seteru ( | |
| 1 34hzatan2 called with both args zero, 34, 1, 2) | 17 | 17 | 1 34hzatan2 called with both args zero, 34, 1, 2) | |
| c | 18 | 18 | c | |
| zatan2 = dcmplx (sign(0.5*pi,real(csn)), 0.0) | 19 | 19 | zatan2 = cmplx (sign(0.5*pi,real(csn)), 0.0, kind(1.d0)) | |
| c | 20 | 20 | c | |
| return | 21 | 21 | return | |
| end | 22 | 22 | end | |
| 23 | 23 | |||
fvn_fnlib/zatanh.f
| complex(8) 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(8) 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.,1.)/ | |
| 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 | |||
fvn_fnlib/zbeta.f
| complex(8) 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(8) a, b, zgamma, zlbeta | 4 | 4 | complex(kind(1.d0)) a, b, zgamma, zlbeta | |
| real(8) 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.0 / | |
| 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 | |||
fvn_fnlib/zcbrt.f
| complex(8) function zcbrt (z) | 1 | 1 | complex(kind(1.d0)) function zcbrt (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(8) z | 4 | 4 | complex(kind(1.d0)) z | |
| real(8) theta,zarg,r,dcbrt | 5 | 5 | real(kind(1.d0)) theta,zarg,r,dcbrt | |
| external zarg, dcbrt | 6 | 6 | external zarg, dcbrt | |
| c | 7 | 7 | c | |
| theta = zarg(z) / 3.0 | 8 | 8 | theta = zarg(z) / 3.0 | |
| r = dcbrt (abs(z)) | 9 | 9 | r = dcbrt (abs(z)) | |
| c | 10 | 10 | c | |
| zcbrt = dcmplx (r*cos(theta), r*sin(theta)) | 11 | 11 | zcbrt = cmplx (r*cos(theta), r*sin(theta), kind(1.d0)) | |
| c | 12 | 12 | c | |
| return | 13 | 13 | return | |
| end | 14 | 14 | end | |
| 15 | 15 | |||
fvn_fnlib/zcosh.f
| complex(8) 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(8) z, ci | 4 | 4 | complex(kind(1.d0)) z, ci | |
| data ci /(0.,1.)/ | 5 | 5 | data ci /(0.,1.)/ | |
| 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 | |||
fvn_fnlib/zcot.f
| complex(8) 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(8) z | 4 | 4 | complex(kind(1.d0)) z | |
| real(8) d1mach | 5 | 5 | real(kind(1.d0)) d1mach | |
| real(8) eps, xmax, ylarge, ybig, rmin, ymin | 6 | 6 | real(kind(1.d0)) eps, xmax, ylarge, ybig, rmin, ymin | |
| real(8) 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.0, 1.5 / | |
| 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 = dcmplx (sn2x/den, -sinh(y2)/den) | 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 = dcmplx (0.0, 1.0/tanh(y2)) | 49 | 49 | zcot = cmplx (0.0, 1.0/tanh(y2), kind(1.d0)) | |
| return | 50 | 50 | return | |
| c | 51 | 51 | c | |
| 30 zcot = dcmplx (0.0, sign (1.0d0, y)) | 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 | |||
fvn_fnlib/zexprl.f
| complex(8) 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(8) z | 11 | 11 | complex(kind(1.d0)) z | |
| real(8) 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 / 0, 2*0.0 / | |
| 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 | |||
fvn_fnlib/zgamma.f
| complex(8) function zgamma (z) | 1 | 1 | complex(kind(1.d0)) function zgamma (z) | |
| 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. | |
| c a preliminary version that is portable, but not accurate enough. | 4 | 4 | c a preliminary version that is portable, but not accurate enough. | |
| complex(8) z, zlngam | 5 | 5 | complex(kind(1.d0)) z, zlngam | |
| external zlngam | 6 | 6 | external zlngam | |
| c | 7 | 7 | c | |
| zgamma = exp (zlngam(z)) | 8 | 8 | zgamma = exp (zlngam(z)) | |
| c | 9 | 9 | c | |
| return | 10 | 10 | return | |
| end | 11 | 11 | end | |
| 12 | 12 | |||
fvn_fnlib/zgamr.f
| complex(8) function zgamr (z) | 1 | 1 | complex(kind(1.d0)) function zgamr (z) | |
| 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. | |
| c this version is an inaccurate preliminary one. eventually this | 4 | 4 | c this version is an inaccurate preliminary one. eventually this | |
| c routine should be a fundamental routine with no dependence on cgamma. | 5 | 5 | c routine should be a fundamental routine with no dependence on cgamma. | |
| c | 6 | 6 | c | |
| complex(8) z, zlngam | 7 | 7 | complex(kind(1.d0)) z, zlngam | |
| real(8) x | 8 | 8 | real(kind(1.d0)) x | |
| integer irold,ir | 9 | 9 | integer irold,ir | |
| external zlngam | 10 | 10 | external zlngam | |
| c | 11 | 11 | c | |
| zgamr = (0.0, 0.0) | 12 | 12 | zgamr = (0.0, 0.0) | |
| x = real (z) | 13 | 13 | x = real (z) | |
| if (x.le.0.0 .and. aint(x).eq.x .and. 0.0.eq.aimag(z)) return | 14 | 14 | if (x.le.0.0 .and. aint(x).eq.x .and. 0.0.eq.aimag(z)) return | |
| c | 15 | 15 | c | |
| call entsrc (irold, 1) | 16 | 16 | call entsrc (irold, 1) | |
| zgamr = zlngam(z) | 17 | 17 | zgamr = zlngam(z) | |
| call erroff | 18 | 18 | call erroff | |
| call entsrc (ir, irold) | 19 | 19 | call entsrc (ir, irold) | |
| zgamr = exp (-zgamr) | 20 | 20 | zgamr = exp (-zgamr) | |
| c | 21 | 21 | c | |
| return | 22 | 22 | return | |
| end | 23 | 23 | end | |
| 24 | 24 | |||
fvn_fnlib/zlbeta.f
| complex(8) function zlbeta (a, b) | 1 | 1 | complex(kind(1.d0)) function zlbeta (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. | |
| c a preliminary version that is portable, but not accurate enough. | 4 | 4 | c a preliminary version that is portable, but not accurate enough. | |
| complex(8) a, b, zlngam | 5 | 5 | complex(kind(1.d0)) a, b, zlngam | |
| external zlngam | 6 | 6 | external zlngam | |
| c | 7 | 7 | c | |
| if (real(a).le.0.0 .or. real(b).le.0.0) call seteru ( | 8 | 8 | if (real(a).le.0.0 .or. real(b).le.0.0) call seteru ( | |
| 1 48hclbeta real part of both arguments must be gt 0, 48, 1, 2) | 9 | 9 | 1 48hclbeta real part of both arguments must be gt 0, 48, 1, 2) | |
| c | 10 | 10 | c | |
| zlbeta = zlngam(a) + zlngam(b) - zlngam(a+b) | 11 | 11 | zlbeta = zlngam(a) + zlngam(b) - zlngam(a+b) | |
| c | 12 | 12 | c | |
| return | 13 | 13 | return | |
| end | 14 | 14 | end | |
| 15 | 15 | |||
fvn_fnlib/zlngam.f
| complex(8) 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(8) zin, z, corr, zlnrel, z9lgmc | 6 | 6 | complex(kind(1.d0)) zin, z, corr, zlnrel, z9lgmc | |
| real(8) d1mach,pi,sq2pil,bound,dxrel,rmax,cabsz | 7 | 7 | real(kind(1.d0)) d1mach,pi,sq2pil,bound,dxrel,rmax,cabsz | |
| real(8) 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 8979324e0 / | 11 | 11 | data pi / 3.1415926535 8979324d0 / | |
| data sq2pil / 0.9189385332 0467274e0 / | 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.0 / | |
| 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 (-dcmplx(0.0,2.0*pi)*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 - dcmplx(0.0,pi)*(z-0.5) - zlnrel(-corr) | 45 | 45 | zlngam = sq2pil + 1.0 - cmplx(0.0,pi,kind(1.d0))*(z-0.5) - | |
| 1 + (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 = -dcmplx (log(abs(corr)), argsum) | 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 | |||
fvn_fnlib/zlnrel.f
| complex(8) 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(8) z | 13 | 13 | complex(kind(1.d0)) z | |
| real(8) 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.0/ | |
| 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 = dcmplx (0.5*dlnrel(2.*x+rho**2), zarg(1.0+z)) | 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 | |||
fvn_fnlib/zlog10.f
| complex(8) function zlog10 (z) | 1 | 1 | complex(kind(1.d0)) function zlog10 (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(8) z | 4 | 4 | complex(kind(1.d0)) z | |
| real(8) aloge | 5 | 5 | real(kind(1.d0)) aloge | |
| data aloge / 0.4342944819 0325182765d0 / | 6 | 6 | data aloge / 0.4342944819 0325182765d0 / | |
| c | 7 | 7 | c | |
| zlog10 = aloge * log(z) | 8 | 8 | zlog10 = aloge * log(z) | |
| c | 9 | 9 | c | |
| return | 10 | 10 | return | |
| end | 11 | 11 | end | |
| 12 | 12 | |||
fvn_fnlib/zpsi.f
| complex(8) 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(8) 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(8) bern,d1mach,pi,bound, dxrel, rmin, rbig | 6 | 6 | real(kind(1.d0)) bern,d1mach,pi,bound, dxrel, rmin, rbig | |
| real(8) 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 e-1 / | |
| data bern( 2) / -.8333333333 3333333 e-2 / | 12 | 12 | data bern( 2) / -.8333333333 3333333 e-2 / | |
| data bern( 3) / .3968253968 2539683 e-2 / | 13 | 13 | data bern( 3) / .3968253968 2539683 e-2 / | |
| data bern( 4) / -.4166666666 6666667 e-2 / | 14 | 14 | data bern( 4) / -.4166666666 6666667 e-2 / | |
| data bern( 5) / .7575757575 7575758 e-2 / | 15 | 15 | data bern( 5) / .7575757575 7575758 e-2 / | |
| data bern( 6) / -.2109279609 2796093 e-1 / | 16 | 16 | data bern( 6) / -.2109279609 2796093 e-1 / | |
| data bern( 7) / .8333333333 3333333 e-1 / | 17 | 17 | data bern( 7) / .8333333333 3333333 e-1 / | |
| data bern( 8) / -.4432598039 2156863 e0 / | 18 | 18 | data bern( 8) / -.4432598039 2156863 e0 / | |
| data bern( 9) / .3053954330 2701197 e1 / | 19 | 19 | data bern( 9) / .3053954330 2701197 e1 / | |
| data bern(10) / -.2645621212 1212121 e2 / | 20 | 20 | data bern(10) / -.2645621212 1212121 e2 / | |
| data bern(11) / .2814601449 2753623 e3 / | 21 | 21 | data bern(11) / .2814601449 2753623 e3 / | |
| data bern(12) / -.3454885393 7728938 e4 / | 22 | 22 | data bern(12) / -.3454885393 7728938 e4 / | |
| data bern(13) / .5482758333 3333333 e5 / | 23 | 23 | data bern(13) / .5482758333 3333333 e5 / | |
| c | 24 | 24 | c | |
| data pi / 3.141592653 589793 e0 / | 25 | 25 | data pi / 3.141592653 589793 e0 / | |
| data nterm, bound, dxrel, rmin, rbig / 0, 4*0.0 / | 26 | 26 | data nterm, bound, dxrel, rmin, rbig / 0, 4*0.0 / | |
| 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 | |||
fvn_fnlib/zsinh.f
| complex(8) 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(8) z, ci | 4 | 4 | complex(kind(1.d0)) z, ci | |
| data ci /(0.,1.)/ | 5 | 5 | data ci /(0.,1.)/ | |
| 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 | |||
fvn_fnlib/ztan.f
| complex(8) 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(8) z | 4 | 4 | complex(kind(1.d0)) z | |
| real(8) d1mach,eps, xmax, ylarge, ybig, ymin,sn2x,den | 5 | 5 | real(kind(1.d0)) d1mach,eps, xmax, ylarge, ybig, ymin,sn2x,den | |
| real(8) 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.0, 1.50 / | |
| 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 = dcmplx (sn2x/den, sinh(y2)/den) | 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 = dcmplx (0.0, tanh(y2)) | 45 | 45 | ztan = cmplx (0.0, tanh(y2), kind(1.d0)) | |
| return | 46 | 46 | return | |
| c | 47 | 47 | c | |
| 30 ztan = dcmplx (0.0, -sign(1.0d0, y)) | 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 | |||
fvn_fnlib/ztanh.f
| complex(8) 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(8) 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.,1.)/ | |
| 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 | |||