Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : module m_SphBessel
8 :
9 : !-------------------------------------------------------------------------
10 : ! SphBessel calculates spherical Bessel functions of the first,
11 : ! second and third kind (Bessel, Neumann and Hankel functions).
12 : ! ModSphBessel calculates modified spherical Bessel functions
13 : ! of the first and second kind.
14 : !
15 : ! jl : spherical Bessel function of the first kind (Bessel)
16 : ! nl : spherical Bessel function of the second kind (Neumann)
17 : ! hl : spherical Bessel function of the third kind (Hankel)
18 : ! il : modified spherical Bessel function of the first kind
19 : ! kl : modified spherical Bessel function of the second kind
20 : !
21 : ! z : Bessel functions are calculated for this value
22 : ! lmax: Bessel functions are calculated for all the indices l
23 : ! from 0 to lmax
24 : !
25 : ! intent(in):
26 : ! z : complex or real scalar
27 : ! lmax: integer
28 : !
29 : ! intent(out):
30 : ! * SphBessel( jl, nl, hl, z, lmax )
31 : ! jl: complex or real, dimension(0:lmax)
32 : ! nl: complex or real, dimension(0:lmax)
33 : ! hl: complex, dimension(0:lmax)
34 : ! * ModSphBessel( il, kl, z, lmax )
35 : ! il: complex or real, dimension(0:lmax)
36 : ! kl: complex or real, dimension(0:lmax)
37 : !
38 : ! All subroutines are pure and therefore can be called for a range of
39 : ! z-values concurrently, f.e. this way:
40 : ! allocate( il(0:lmax, size(z)), kl(0:lmax, size(z)) )
41 : ! do concurrent (i = 1: size(z))
42 : ! call ModSphBessel( il(:,i), kl(:,i), z(i), lmax )
43 : ! end do
44 : !
45 : ! details on implementation:
46 : ! For |z| <= 1 the taylor expansions of jl and nl are used.
47 : ! For |z| > 1 the explicit expressions for hl(+), hl(-) are used.
48 : ! For modified spherical Bessel functions il and kl the relations
49 : ! il(z) = I^{-l} * jl(I*z)
50 : ! kl(z) = -I^{l} * hl(I*z)
51 : ! are used.
52 : !
53 : ! authors:
54 : ! originally written by R. Zeller (1990)
55 : ! modernised and extended by M. Hinzen (2016)
56 : !-------------------------------------------------------------------------
57 :
58 : use m_constants, only: ImagUnit
59 : implicit none
60 :
61 :
62 : interface SphBessel
63 : module procedure SphBesselComplex, SphBesselReal
64 : end interface
65 :
66 : interface ModSphBessel
67 : ! variant Complex2 takes workspace as an argument.
68 : ! this is not possible for the subroutine working on reals.
69 : module procedure ModSphBesselComplex, ModSphBesselReal, ModSphBesselComplex2
70 : end interface
71 :
72 :
73 : contains
74 :
75 :
76 41925 : pure subroutine SphBesselComplex ( jl, nl, hl, z, lmax )
77 :
78 : complex, intent(in) :: z
79 : integer, intent(in) :: lmax
80 : complex, dimension(0:lmax), intent(out) :: jl, nl, hl
81 :
82 : complex :: termj, termn, z2, zj, zn
83 : real :: rl, rn
84 41925 : complex, dimension(0:lmax) :: rnm
85 : integer :: l, m, n
86 :
87 :
88 41925 : zj = 1.0
89 41925 : zn = 1.0 / z
90 41925 : z2 = z * z
91 438615 : jl(:) = 1.0
92 438615 : nl(:) = 1.0
93 41925 : if ( abs( z ) < lmax + 1.0 ) then
94 438570 : SERIAL_L_LOOP: do l = 0, lmax
95 396660 : rl = l + l
96 396660 : termj = 1.0
97 396660 : termn = 1.0
98 10313160 : EXPANSION: do n = 1, 25
99 9916500 : rn = n + n
100 9916500 : termj = -termj / ( rl + rn + 1.0 ) / rn * z2
101 9916500 : termn = termn / ( rl - rn + 1.0 ) / rn * z2
102 9916500 : jl(l) = jl(l) + termj
103 10313160 : nl(l) = nl(l) + termn
104 : end do EXPANSION
105 396660 : jl(l) = jl(l) * zj
106 396660 : nl(l) = -nl(l) * zn
107 396660 : hl(l) = jl(l) + nl(l) * ImagUnit
108 396660 : zj = zj * z / ( rl + 3.0 )
109 438570 : zn = zn / z * ( rl + 1.0 )
110 : end do SERIAL_L_LOOP
111 : end if
112 :
113 438615 : rnm(:) = 1.0
114 : !PARALLEL_L_LOOP: do concurrent (l = 0: lmax)
115 438615 : PARALLEL_L_LOOP: do l = 0,lmax
116 438615 : if ( abs( z ) >= l + 1.0 ) then
117 40950 : hl(l) = 0.0
118 40950 : nl(l) = 0.0
119 101715 : SERIAL_M_LOOP: do m = 0, l
120 60765 : hl(l) = hl(l) + (-1) ** m * rnm(l)
121 60765 : nl(l) = nl(l) + rnm(l)
122 101715 : rnm(l) = rnm(l) / ( m + 1.0 ) * ( l * ( l + 1 ) - m * ( m + 1 ) ) / ( ImagUnit * ( z + z ) )
123 : end do SERIAL_M_LOOP
124 40950 : hl(l) = hl(l) * (-ImagUnit) ** l * exp( ImagUnit * z ) / ( ImagUnit * z )
125 40950 : nl(l) = nl(l) * ImagUnit ** l * exp( -ImagUnit * z ) / ( -ImagUnit * z )
126 40950 : jl(l) = ( hl(l) + nl(l) ) / 2.0
127 40950 : nl(l) = ( hl(l) - jl(l) ) * (-ImagUnit)
128 : end if
129 : end do PARALLEL_L_LOOP
130 :
131 41925 : end subroutine SphBesselComplex
132 :
133 :
134 :
135 0 : pure subroutine SphBesselReal ( jl, nl, hl, x, lmax )
136 :
137 : real, intent(in) :: x
138 : integer, intent(in) :: lmax
139 : real, dimension(0:lmax), intent(out) :: jl, nl
140 : complex, dimension(0:lmax), intent(out) :: hl
141 :
142 0 : complex, dimension(0:lmax) :: jl_complex, nl_complex
143 :
144 0 : call SphBesselComplex( jl_complex, nl_complex, hl, cmplx(x,0.0), lmax )
145 :
146 0 : jl = real(jl_complex)
147 0 : nl = real(nl_complex)
148 :
149 0 : end subroutine SphBesselReal
150 :
151 :
152 :
153 0 : pure subroutine ModSphBesselComplex ( il, kl, z, lmax )
154 :
155 : complex, intent(in) :: z
156 : integer, intent(in) :: lmax
157 : complex, dimension(0:lmax), intent(out) :: il, kl
158 :
159 0 : complex, dimension(0:lmax) :: nl
160 : integer :: l
161 :
162 :
163 0 : call SphBesselComplex( il, nl, kl, ImagUnit * z, lmax )
164 :
165 0 : do l = 0, lmax
166 0 : il(l) = (-ImagUnit) ** l * il(l)
167 0 : kl(l) = - ImagUnit ** l * kl(l)
168 : end do
169 :
170 0 : end subroutine ModSphBesselComplex
171 :
172 :
173 : !another implementation of ModSphBesselComplex, where nl is allocated outside for performance reasons
174 0 : pure subroutine ModSphBesselComplex2 ( il, kl, nl, z, lmax )
175 :
176 : complex, intent(in) :: z
177 : integer, intent(in) :: lmax
178 : complex, dimension(0:lmax), intent(out) :: il, kl, nl
179 :
180 : integer :: l
181 :
182 :
183 0 : call SphBesselComplex( il, nl, kl, ImagUnit * z, lmax )
184 :
185 0 : do l = 0, lmax
186 0 : il(l) = (-ImagUnit) ** l * il(l)
187 0 : kl(l) = - ImagUnit ** l * kl(l)
188 : end do
189 :
190 0 : end subroutine ModSphBesselComplex2
191 :
192 :
193 :
194 41925 : pure subroutine ModSphBesselReal ( il, kl, x, lmax )
195 :
196 : real, intent(in) :: x
197 : integer, intent(in) :: lmax
198 : real, dimension(0:lmax), intent(out) :: il, kl
199 :
200 41925 : complex, dimension(0:lmax) :: jl, nl, hl
201 : integer :: l
202 : complex :: z
203 :
204 :
205 41925 : z = ImagUnit * x
206 41925 : call SphBesselComplex( jl, nl, hl, z, lmax )
207 :
208 438615 : do l = 0, lmax
209 396690 : il(l) = real((-ImagUnit) ** l * jl(l))
210 438615 : kl(l) = real(- ImagUnit ** l * hl(l))
211 : end do
212 :
213 41925 : end subroutine ModSphBesselReal
214 :
215 :
216 : end module m_SphBessel
|
