LCOV - code coverage report
Current view: top level - math - SphBessel.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 45 65 69.2 %
Date: 2019-09-08 04:53:50 Functions: 2 5 40.0 %

          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       83850 :     real,    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      835230 :       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      835305 :     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             :     complex                                 :: z
     144             : 
     145             : 
     146           0 :     z = x           ! internal conversion from real to complex
     147             : 
     148           0 :     call SphBesselComplex( jl_complex, nl_complex, hl, z, lmax )
     149             : 
     150           0 :     jl = jl_complex ! internal conversion from complex to real
     151           0 :     nl = nl_complex ! internal conversion from complex to real
     152             : 
     153           0 :   end subroutine SphBesselReal
     154             : 
     155             : 
     156             : 
     157           0 :   pure subroutine ModSphBesselComplex ( il, kl, z, lmax )
     158             : 
     159             :     complex,                    intent(in)  :: z
     160             :     integer,                    intent(in)  :: lmax
     161             :     complex, dimension(0:lmax), intent(out) :: il, kl
     162             : 
     163           0 :     complex, dimension(0:lmax)              :: nl
     164             :     integer                                 :: l 
     165             : 
     166             : 
     167           0 :     call SphBesselComplex( il, nl, kl, ImagUnit * z, lmax )
     168             : 
     169           0 :     do l = 0, lmax
     170           0 :       il(l) = (-ImagUnit) ** l * il(l)
     171           0 :       kl(l) = - ImagUnit  ** l * kl(l)
     172             :     end do     
     173             : 
     174           0 :   end subroutine ModSphBesselComplex
     175             : 
     176             : 
     177             :   !another implementation of ModSphBesselComplex, where nl is allocated outside for performance reasons
     178           0 :   pure subroutine ModSphBesselComplex2 ( il, kl, nl, z, lmax )
     179             : 
     180             :     complex,                    intent(in)  :: z
     181             :     integer,                    intent(in)  :: lmax
     182             :     complex, dimension(0:lmax), intent(out) :: il, kl, nl
     183             : 
     184             :     integer                                 :: l 
     185             : 
     186             : 
     187           0 :     call SphBesselComplex( il, nl, kl, ImagUnit * z, lmax )
     188             : 
     189           0 :     do l = 0, lmax
     190           0 :       il(l) = (-ImagUnit) ** l * il(l)
     191           0 :       kl(l) = - ImagUnit  ** l * kl(l)
     192             :     end do     
     193             : 
     194           0 :   end subroutine ModSphBesselComplex2
     195             : 
     196             : 
     197             : 
     198       41925 :   pure subroutine ModSphBesselReal ( il, kl, x, lmax )
     199             : 
     200             :     real,                       intent(in)  :: x
     201             :     integer,                    intent(in)  :: lmax
     202             :     real,    dimension(0:lmax), intent(out) :: il, kl
     203             : 
     204       83850 :     complex, dimension(0:lmax)              :: jl, nl, hl
     205             :     integer                                 :: l 
     206             :     complex                                 :: z
     207             : 
     208             : 
     209       41925 :     z = ImagUnit * x
     210       41925 :     call SphBesselComplex( jl, nl, hl, z, lmax )
     211             : 
     212      438615 :     do l = 0, lmax
     213      396690 :       il(l) = (-ImagUnit) ** l * jl(l)
     214      438615 :       kl(l) = - ImagUnit  ** l * hl(l)
     215             :     end do     
     216             : 
     217       41925 :   end subroutine ModSphBesselReal
     218             : 
     219             : 
     220             : end module m_SphBessel

Generated by: LCOV version 1.13