8 use constants
, only: dp,pi,iu,eulergamma
11 public factorial,digamma,bessel,hankel
18 subroutine bessel(n,x,tol,j,dj,y,dy)
19 use besselj
, only: zbesj
20 use bessely
, only: zbesy
21 integer,
intent(in) :: n
22 complex(dp),
intent(in) :: x
23 real(dp),
intent(in) :: tol
24 complex(dp),
dimension(0:n),
intent(out) :: &
29 real(dp),
dimension(n+1) :: cyr,cyi,cwr,cwi
30 integer :: i,m,nz,ierr
31 call zbesj(realpart(x),imagpart(x),0.e0_dp,1,n+1,cyr,cyi,nz,ierr)
33 write(*,*)
'underflow encountered' 34 write(*,*)
'bessel function J not calculated',n,x
38 write(*,*)
'did not converge' 39 write(*,*)
'bessel function J not calculated',n,x
43 j(i)=cmplx(cyr(i+1),cyi(i+1),kind=dp)
47 call zbesy(realpart(x),imagpart(x),0.e0_dp,1,n+1,cyr,cyi,nz,cwr,cwi,ierr)
49 write(*,*)
'underflow encountered' 50 write(*,*)
'bessel function Y not calculated',n,x
54 write(*,*)
'did not converge' 55 write(*,*)
'bessel function Y not calculated',n,x
59 y(i)=cmplx(cyr(i+1),cyi(i+1),kind=dp)
80 subroutine hankel(n,x,tol,h1,dh1,h2,dh2)
81 integer,
intent(in) :: n
82 complex(dp),
intent(in) :: x
83 real(dp),
intent(in) :: tol
84 complex(dp),
dimension(0:n),
intent(out) :: &
90 complex(dp),
dimension(0:n) :: j,dj,y,dy
92 call bessel(n,x,tol,j,dj,y,dy)
98 dh1(i)=h1(i-1)-i/x*h1(i)
99 dh2(i)=h2(i-1)-i/x*h2(i)
101 end subroutine hankel
107 recursive function factorial(n)
result(fact)
108 integer,
intent(in) :: n
111 stop
'factorial defined only for nonnegative integers' 115 fact=n*factorial(n-1)
117 end function factorial
125 real(dp) function digamma(n)
126 integer,
intent(in) :: n
130 digamma=digamma+1.0_dp/i