MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
specfun.f90
Go to the documentation of this file.
1 
7 module specfun
8  use constants, only: dp,pi,iu,eulergamma
9  implicit none
10  private
11  public factorial,digamma,bessel,hankel
12 contains
13 !********************************BEGINNING*************************************
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 !maximum order calculated
22  complex(dp), intent(in) :: x !argument
23  real(dp), intent(in) :: tol !relative tolerance, currently not used
24  complex(dp), dimension(0:n), intent(out) :: &
25  j,& !bessel function of the first kind
26  dj,& !derivative of bessel function of the first kind
27  y,& !bessel function of the second kind
28  dy !derivative of bessel function of the second kind
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)
32  if (nz>0) then
33  write(*,*) 'underflow encountered'
34  write(*,*) 'bessel function J not calculated',n,x
35  stop
36  endif
37  if (ierr /= 0) then
38  write(*,*) 'did not converge'
39  write(*,*) 'bessel function J not calculated',n,x
40  stop
41  endif
42  do i=0,n
43  j(i)=cmplx(cyr(i+1),cyi(i+1),kind=dp)
44  enddo
45 ! write(*,*) x
46 ! write(*,*) j
47  call zbesy(realpart(x),imagpart(x),0.e0_dp,1,n+1,cyr,cyi,nz,cwr,cwi,ierr)
48  if (nz>0) then
49  write(*,*) 'underflow encountered'
50  write(*,*) 'bessel function Y not calculated',n,x
51  stop
52  endif
53  if (ierr /= 0) then
54  write(*,*) 'did not converge'
55  write(*,*) 'bessel function Y not calculated',n,x
56  stop
57  endif
58  do i=0,n
59  y(i)=cmplx(cyr(i+1),cyi(i+1),kind=dp)
60  enddo
61 ! write(*,*) x
62 ! write(*,*) y
63 ! stop
64  !http://keisan.casio.com/exec/system/1180573474
65  dj(0)=-j(1)
66  dy(0)=-y(1)
67  do m=1,n
68  dj(m)=j(m-1)-m/x*j(m)
69  dy(m)=y(m-1)-m/x*y(m)
70  enddo
71 end subroutine bessel
72 !***********************************END****************************************
73 
74 
75 !********************************BEGINNING*************************************
80 subroutine hankel(n,x,tol,h1,dh1,h2,dh2)
81  integer, intent(in) :: n !maximum order calculated
82  complex(dp), intent(in) :: x !argument
83  real(dp), intent(in) :: tol !relative tolerance
84  complex(dp), dimension(0:n), intent(out) :: &
85  h1,& !hankel function of the first kind
86  dh1,& !derivative of hankel function of the first kind
87  h2,& !hankel function of the second kind
88  dh2 !derivative of hankel function of the second kind
89  integer :: i
90  complex(dp), dimension(0:n) :: j,dj,y,dy
91  !http://keisan.casio.com/has10/SpecExec.cgi?id=system/2006/1222514923
92  call bessel(n,x,tol,j,dj,y,dy)
93  h1=j+iu*y
94  h2=j-iu*y
95  dh1(0)=-h1(1)
96  dh2(0)=-h2(1)
97  do i=1,n
98  dh1(i)=h1(i-1)-i/x*h1(i)
99  dh2(i)=h2(i-1)-i/x*h2(i)
100  enddo
101 end subroutine hankel
102 !***********************************END****************************************
103 
104 
105 !********************************BEGINNING*************************************
107 recursive function factorial(n) result(fact)
108  integer, intent(in) :: n !argument
109  real(dp) :: fact
110  if (n<0) then
111  stop 'factorial defined only for nonnegative integers'
112  elseif (n==0) then
113  fact=1
114  else
115  fact=n*factorial(n-1)
116  endif
117 end function factorial
118 !***********************************END****************************************
119 
120 
121 !********************************BEGINNING*************************************
125 real(dp) function digamma(n)
126  integer, intent(in) :: n !argument
127  integer :: i
128  digamma=-eulergamma
129  do i=1,n-1
130  digamma=digamma+1.0_dp/i
131  enddo
132 end function digamma
133 !***********************************END****************************************
134 
135 
136 !SUBROUTINE JNCOMP(Z,NS,JN)
137 ! !taken from Daniel Mackowski code for scattering of radiation by long cylinders
138 ! !it might be more effective, I did not test it thoroughly
139 ! IMPLICIT none
140 ! COMPLEX*16 JN(0:*),A,Z
141 ! integer :: n,nd,ns
142 !
143 !! CALCULATES THE INTEGER-ORDER BESSEL FUNCTION JN(Z)
144 !! N=0,1,...NS, FOR COMPLEX ARGUMENT Z. REFER TO BOHREN & HUFFMAN
145 !! FOR DETAILS. CODED 9/20/90
146 !
147 ! ND=NINT((101+CDABS(Z))**.499+NS)
148 ! ND=2*NINT(ND/2.)
149 ! JN(ND)=0.
150 ! JN(ND-1)=1.D-32
151 ! A=0.
152 ! DO N=ND-1,3,-2
153 ! JN(N-1)=2.*N*JN(N)/Z-JN(N+1)
154 ! JN(N-2)=2.*(N-1)*JN(N-1)/Z-JN(N)
155 ! A=A+JN(N-1)
156 ! enddo
157 ! JN(0)=2.*JN(1)/Z-JN(2)
158 ! A=JN(0)+2.*A
159 ! DO N=0,NS
160 ! JN(N)=JN(N)/A
161 ! enddo
162 !RETURN
163 !END subroutine
164 
165 
166 !!********************************BEGINNING*************************************
167 !!calculates Bessel functions of the first and second kind and their derivatives for complex argument
168 !!doesn't work properly for high orders
169 !subroutine bessel(n,x,tol,j,dj,y,dy)
170 ! integer, intent(in) :: n !maximum order calculated
171 ! complex(dp), intent(in) :: x !argument
172 ! real(dp), intent(in) :: tol !relative tolerance
173 ! complex(dp), dimension(0:n), intent(out) :: j !bessel function of the first kind
174 ! complex(dp), dimension(0:n), intent(out) :: dj !derivative of bessel function of the first kind
175 ! complex(dp), dimension(0:n), intent(out) :: y !bessel function of the second kind
176 ! complex(dp), dimension(0:n), intent(out) :: dy !derivative of bessel function of the second kind
177 ! integer :: l,m,lmax=200
178 ! complex(dp) :: jold,yold
179 ! jold=0;yold=0
180 ! j=0;dj=0;y=0;dy=0
181 !! write(*,*) x
182 ! do m=0,n
183 ! do l=0,lmax
184 !! write(*,*) l,m,j(m)
185 !! j(m)=j(m)+(-1)**l/(2**(2*l+m)*factorial(l)*factorial(m+l)+0.0_dp)*x**(2*l+m) !http://mathworld.wolfram.com/BesselFunctionoftheFirstKind.html (it works just as well)
186 ! j(m)=j(m)+(-1)**l/(factorial(l)*factorial(m+l)+0.0_dp)*(x/2)**(2*l+m) !http://keisan.casio.com/exec/system/1180573474
187 ! if (abs(j(m)-jold)/abs(j(m))<tol) then
188 !! write(*,*) m,l
189 ! exit
190 ! elseif (l==lmax) then
191 ! write(*,*) 'bessel function J not calculated',n,x
192 ! stop
193 ! endif
194 ! jold=j(m)
195 ! enddo
196 ! write(*,*) j(m)
197 ! !http://mathworld.wolfram.com/BesselFunctionoftheSecondKind.html
198 ! y(m)=2/pi*log(x/2)*j(m) !second term
199 ! do l=0,m-1
200 ! y(m)=y(m)-(x/2)**(-m)/pi*factorial(m-l-1)/factorial(l)*(x**2/4)**l !first term
201 ! enddo
202 ! do l=0,lmax
203 !! write(*,*) l,m,y(m)
204 ! y(m)=y(m)-(x/2)**m/pi*(digamma(l+1)+digamma(m+l+1))*(-x**2/4)**l/(factorial(l)*factorial(m+l)) !third term
205 !! write(*,*) l
206 ! if (abs(y(m)-yold)/abs(y(m))<tol) then
207 !! write(*,*) l
208 ! exit
209 ! elseif (l==lmax) then
210 ! write(*,*) 'bessel function Y not calculated',n,x
211 ! stop
212 ! endif
213 ! yold=y(m)
214 ! enddo
215 ! enddo
216 ! !http://keisan.casio.com/exec/system/1180573474
217 ! dj(0)=-j(1)
218 ! dy(0)=-y(1)
219 ! do m=1,n
220 ! dj(m)=j(m-1)-m/x*j(m)
221 ! dy(m)=y(m-1)-m/x*y(m)
222 ! enddo
223 !end subroutine bessel
224 !!***********************************END****************************************
225 end module specfun