SphBessel.f90 6.89 KB
Newer Older
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
!--------------------------------------------------------------------------------
! 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
115 116
    !PARALLEL_L_LOOP: do concurrent (l = 0: lmax)
    PARALLEL_L_LOOP: do l = 0,lmax
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