1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
|
! Copyright (C) 2002-2005 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
! This file is distributed under the terms of the GNU Lesser General Public
! License. See the file COPYING for license details.
!BOP
! !ROUTINE: sbesseldm
! !INTERFACE:
subroutine sbesseldm(m,lmax,x,djl)
! !INPUT/OUTPUT PARAMETERS:
! m : order of derivatve (in,integer)
! lmax : maximum order of Bessel function (in,integer)
! x : real argument (in,real)
! djl : array of returned values (out,real(0:lmax))
! !DESCRIPTION:
! Computes the $m$th derivative of the spherical Bessel function of the first
! kind, $j_l(x)$, for argument $x$ and $l=0,1,\ldots,l_{\rm max}$. For
! $x\ge 1$ this is done by repeatedly using the relations
! \begin{align*}
! \frac{d}{dx}j_l(x)&=\frac{l}{x}j_l(x)-j_{l+1}(x) \\
! j_{l+1}(x)&=\frac{2l+1}{x}j_l(x)-j_{l-1}(x).
! \end{align*}
! While for $x<1$ the series expansion of the Bessel function is used
! $$ \frac{d^m}{dx^m}j_l(x)=\sum_{i=0}^{\infty}
! \frac{(2i+l)!}{(-2)^ii!(2i+l-m)!(2i+2l+1)!!}x^{2i+l-m}. $$
! This procedure is numerically stable and accurate to near machine precision
! for $l\le 30$ and $m\le 6$.
!
! !REVISION HISTORY:
! Created March 2003 (JKD)
! Modified to return an array of values, October 2004 (JKD)
!EOP
!BOC
implicit none
! arguments
integer, intent(in) :: m,lmax
real(8), intent(in) :: x
real(8), intent(out) :: djl(0:lmax)
! local variables
integer, parameter :: maxm=6,maxns=20
integer i,j,l,i0
real(8) t1,sm,x2
! automatic arrays
integer a(0:maxm+1),a1(0:maxm+1)
integer b(0:maxm+1),b1(0:maxm+1)
real(8) jl(0:lmax+1)
! external functions
real(8), external :: factn,factn2,factr
if (m == 0) then
call sbessel(lmax,x,djl)
return
end if
if ((m < 0).or.(m > maxm)) then
write(*,*)
write(*,'("Error(sbesseldm): m out of range : ",I0)') m
write(*,*)
stop
end if
if ((lmax < 0).or.(lmax > 30)) then
write(*,*)
write(*,'("Error(sbesseldm): lmax out of range : ",I0)') lmax
write(*,*)
stop
end if
if ((x < 0.d0).or.(x > 1.d5)) then
write(*,*)
write(*,'("Error(sbesseldm): x out of range : ",G18.10)') x
write(*,*)
stop
end if
if (x > 1.d0) then
call sbessel(lmax+1,x,jl)
do l=0,lmax
a(1:m+1)=0
a(0)=1
a1(0:m+1)=0
do i=1,m
b(0)=0
b1(0)=0
do j=0,i
b(j+1)=a(j)*(l-j)
b1(j+1)=-a1(j)*(j+l+2)
end do
do j=0,i
b1(j)=b1(j)-a(j)
b(j)=b(j)+a1(j)
end do
a(0:i+1)=b(0:i+1)
a1(0:i+1)=b1(0:i+1)
end do
t1=1.d0
sm=dble(a(0))*jl(l)+dble(a1(0))*jl(l+1)
do i=1,m+1
t1=t1*x
sm=sm+(dble(a(i))*jl(l)+dble(a1(i))*jl(l+1))/t1
end do
djl(l)=sm
end do
else
x2=x**2
do l=0,lmax
i0=max((m-l+1)/2,0)
j=2*i0+l-m
if (j == 0) then
t1=1.d0
else
t1=x**j
end if
t1=factr(j+m,j)*t1/(factn(i0)*factn2(j+l+m+1)*dble((-2)**i0))
sm=t1
do i=i0+1,maxns
j=2*i+l
t1=-t1*dble((j-1)*j)*x2/dble((j-l)*(j-m-1)*(j-m)*(j+l+1))
if (abs(t1) < 1.d-40) exit
sm=sm+t1
end do
djl(l)=sm
end do
end if
end subroutine
!EOC
|
