Commit 7bb515bc authored by Miriam Hinzen's avatar Miriam Hinzen

Add New Math Functions and Subroutines

parent 2ac3237f
......@@ -29,6 +29,8 @@ math/differentiate.f90
math/fft2d.F90
math/fft3d.f90
math/fft_interface.F90
math/SphBessel.f90
math/DoubleFactorial.f90
)
if (FLEUR_USE_FFTMKL)
set(fleur_F90 ${fleur_F90} math/mkl_dfti.f90)
......
module m_DoubleFactorial
implicit none
contains
real(kind=8) function DoubleFactorial( n_upper, n_lower )
! calculates ( 2 * n_upper + 1 ) !! / ( 2 * n_lower + 1 ) !! or just ( 2 * n_upper + 1 ) !!, if n_lower is not present
integer :: n_upper
integer, optional :: n_lower
integer :: i, i_lower
i_lower = 1
if( present(n_lower) ) i_lower = n_lower + 1
DoubleFactorial = 1.
do i = i_lower, n_upper
DoubleFactorial = DoubleFactorial * ( 2 * i + 1 )
end do
end function DoubleFactorial
end module m_DoubleFactorial
!--------------------------------------------------------------------------------
! 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
module procedure :: SphBesselComplex, SphBesselReal
end interface
interface ModSphBessel
! variant Complex2 takes workspace as an argument.
! this is not possible for the subroutine working on reals.
module procedure :: ModSphBesselComplex, ModSphBesselReal, ModSphBesselComplex2
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
PARALLEL_L_LOOP: do concurrent (l = 0: lmax)
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
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment