SphBessel.f90 6.89 KB
 Miriam Hinzen committed May 30, 2018 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 ``````!-------------------------------------------------------------------------------- ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany ! This file is part of FLEUR and available as free software under the conditions ! of the MIT license as expressed in the LICENSE file in more detail. !-------------------------------------------------------------------------------- module m_SphBessel !------------------------------------------------------------------------- ! SphBessel calculates spherical Bessel functions of the first, ! second and third kind (Bessel, Neumann and Hankel functions). ! ModSphBessel calculates modified spherical Bessel functions ! of the first and second kind. ! ! jl : spherical Bessel function of the first kind (Bessel) ! nl : spherical Bessel function of the second kind (Neumann) ! hl : spherical Bessel function of the third kind (Hankel) ! il : modified spherical Bessel function of the first kind ! kl : modified spherical Bessel function of the second kind ! ! z : Bessel functions are calculated for this value ! lmax: Bessel functions are calculated for all the indices l ! from 0 to lmax ! ! intent(in): ! z : complex or real scalar ! lmax: integer ! ! intent(out): ! * SphBessel( jl, nl, hl, z, lmax ) ! jl: complex or real, dimension(0:lmax) ! nl: complex or real, dimension(0:lmax) ! hl: complex, dimension(0:lmax) ! * ModSphBessel( il, kl, z, lmax ) ! il: complex or real, dimension(0:lmax) ! kl: complex or real, dimension(0:lmax) ! ! All subroutines are pure and therefore can be called for a range of ! z-values concurrently, f.e. this way: ! allocate( il(0:lmax, size(z)), kl(0:lmax, size(z)) ) ! do concurrent (i = 1: size(z)) ! call ModSphBessel( il(:,i), kl(:,i), z(i), lmax ) ! end do ! ! details on implementation: ! For |z| <= 1 the taylor expansions of jl and nl are used. ! For |z| > 1 the explicit expressions for hl(+), hl(-) are used. ! For modified spherical Bessel functions il and kl the relations ! il(z) = I^{-l} * jl(I*z) ! kl(z) = -I^{l} * hl(I*z) ! are used. ! ! authors: ! originally written by R. Zeller (1990) ! modernised and extended by M. Hinzen (2016) !------------------------------------------------------------------------- implicit none complex, parameter :: CI = (0.0, 1.0) interface SphBessel `````` Daniel Wortmann committed Aug 15, 2018 64 `````` module procedure SphBesselComplex, SphBesselReal `````` Miriam Hinzen committed May 30, 2018 65 66 67 68 69 `````` end interface interface ModSphBessel ! variant Complex2 takes workspace as an argument. ! this is not possible for the subroutine working on reals. `````` Daniel Wortmann committed Aug 15, 2018 70 `````` module procedure ModSphBesselComplex, ModSphBesselReal, ModSphBesselComplex2 `````` Miriam Hinzen committed May 30, 2018 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 `````` end interface contains pure subroutine SphBesselComplex ( jl, nl, hl, z, lmax ) complex, intent(in) :: z integer, intent(in) :: lmax complex, dimension(0:lmax), intent(out) :: jl, nl, hl complex :: termj, termn, z2, zj, zn real :: rl, rn real, dimension(0:lmax) :: rnm integer :: l, m, n zj = 1.0 zn = 1.0 / z z2 = z * z jl(:) = 1.0 nl(:) = 1.0 if ( abs( z ) < lmax + 1.0 ) then SERIAL_L_LOOP: do l = 0, lmax rl = l + l termj = 1.0 termn = 1.0 EXPANSION: do n = 1, 25 rn = n + n termj = -termj / ( rl + rn + 1.0 ) / rn * z2 termn = termn / ( rl - rn + 1.0 ) / rn * z2 jl(l) = jl(l) + termj nl(l) = nl(l) + termn end do EXPANSION jl(l) = jl(l) * zj nl(l) = -nl(l) * zn hl(l) = jl(l) + nl(l) * CI zj = zj * z / ( rl + 3.0 ) zn = zn / z * ( rl + 1.0 ) end do SERIAL_L_LOOP end if rnm(:) = 1.0 `````` Uliana Alekseeva committed Jul 24, 2018 115 116 `````` !PARALLEL_L_LOOP: do concurrent (l = 0: lmax) PARALLEL_L_LOOP: do l = 0,lmax `````` Miriam Hinzen committed May 30, 2018 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 `````` if ( abs( z ) >= l + 1.0 ) then hl(l) = 0.0 nl(l) = 0.0 SERIAL_M_LOOP: do m = 0, l hl(l) = hl(l) + (-1) ** m * rnm(l) nl(l) = nl(l) + rnm(l) rnm(l) = rnm(l) / ( m + 1.0 ) * ( l * ( l + 1 ) - m * ( m + 1 ) ) / ( CI * ( z + z ) ) end do SERIAL_M_LOOP hl(l) = hl(l) * (-CI) ** l * exp( CI * z ) / ( CI * z ) nl(l) = nl(l) * CI ** l * exp( -CI * z ) / ( -CI * z ) jl(l) = ( hl(l) + nl(l) ) / 2.0 nl(l) = ( hl(l) - jl(l) ) * (-CI) end if end do PARALLEL_L_LOOP end subroutine SphBesselComplex pure subroutine SphBesselReal ( jl, nl, hl, x, lmax ) real, intent(in) :: x integer, intent(in) :: lmax real, dimension(0:lmax), intent(out) :: jl, nl complex, dimension(0:lmax), intent(out) :: hl complex, dimension(0:lmax) :: jl_complex, nl_complex complex :: z z = x ! internal conversion from real to complex call SphBesselComplex( jl_complex, nl_complex, hl, z, lmax ) jl = jl_complex ! internal conversion from complex to real nl = nl_complex ! internal conversion from complex to real end subroutine SphBesselReal pure subroutine ModSphBesselComplex ( il, kl, z, lmax ) complex, intent(in) :: z integer, intent(in) :: lmax complex, dimension(0:lmax), intent(out) :: il, kl complex, dimension(0:lmax) :: nl integer :: l call SphBesselComplex( il, nl, kl, CI * z, lmax ) do l = 0, lmax il(l) = (-CI) ** l * il(l) kl(l) = - CI ** l * kl(l) end do end subroutine ModSphBesselComplex !another implementation of ModSphBesselComplex, where nl is allocated outside for performance reasons pure subroutine ModSphBesselComplex2 ( il, kl, nl, z, lmax ) complex, intent(in) :: z integer, intent(in) :: lmax complex, dimension(0:lmax), intent(out) :: il, kl, nl integer :: l call SphBesselComplex( il, nl, kl, CI * z, lmax ) do l = 0, lmax il(l) = (-CI) ** l * il(l) kl(l) = - CI ** l * kl(l) end do end subroutine ModSphBesselComplex2 pure subroutine ModSphBesselReal ( il, kl, x, lmax ) real, intent(in) :: x integer, intent(in) :: lmax real, dimension(0:lmax), intent(out) :: il, kl complex, dimension(0:lmax) :: jl, nl, hl integer :: l complex :: z z = CI * x call SphBesselComplex( jl, nl, hl, z, lmax ) do l = 0, lmax il(l) = (-CI) ** l * jl(l) kl(l) = - CI ** l * hl(l) end do end subroutine ModSphBesselReal end module m_SphBessel``````