SphBessel.f90 6.98 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
!--------------------------------------------------------------------------------
! 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)
  !-------------------------------------------------------------------------

58
  use m_constants, only: ImagUnit
59 60 61 62
  implicit none


  interface SphBessel
63
    module procedure  SphBesselComplex, SphBesselReal
64 65 66 67 68
  end interface

  interface ModSphBessel
    ! variant Complex2 takes workspace as an argument.
    ! this is not possible for the subroutine working on reals.
69
    module procedure  ModSphBesselComplex, ModSphBesselReal, ModSphBesselComplex2
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
  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
107
        hl(l) =  jl(l) + nl(l) * ImagUnit
108 109 110 111 112 113
        zj = zj * z / ( rl + 3.0 )
        zn = zn / z * ( rl + 1.0 )
      end do SERIAL_L_LOOP
    end if

    rnm(:) = 1.0
114 115
    !PARALLEL_L_LOOP: do concurrent (l = 0: lmax)
    PARALLEL_L_LOOP: do l = 0,lmax
116 117 118 119 120 121
      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)
122
          rnm(l) = rnm(l) / ( m + 1.0 ) * ( l * ( l + 1 ) - m * ( m + 1 ) ) / ( ImagUnit * ( z + z ) )
123
        end do SERIAL_M_LOOP
124 125
        hl(l) = hl(l) * (-ImagUnit) ** l * exp(  ImagUnit * z ) / (  ImagUnit * z )
        nl(l) = nl(l) *   ImagUnit  ** l * exp( -ImagUnit * z ) / ( -ImagUnit * z )
126
        jl(l) = ( hl(l) + nl(l) ) / 2.0
127
        nl(l) = ( hl(l) - jl(l) ) * (-ImagUnit)
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
      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 


167
    call SphBesselComplex( il, nl, kl, ImagUnit * z, lmax )
168 169

    do l = 0, lmax
170 171
      il(l) = (-ImagUnit) ** l * il(l)
      kl(l) = - ImagUnit  ** l * kl(l)
172 173 174 175 176 177 178 179 180 181 182 183 184 185 186
    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 


187
    call SphBesselComplex( il, nl, kl, ImagUnit * z, lmax )
188 189

    do l = 0, lmax
190 191
      il(l) = (-ImagUnit) ** l * il(l)
      kl(l) = - ImagUnit  ** l * kl(l)
192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208
    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


209
    z = ImagUnit * x
210 211 212
    call SphBesselComplex( jl, nl, hl, z, lmax )

    do l = 0, lmax
213 214
      il(l) = (-ImagUnit) ** l * jl(l)
      kl(l) = - ImagUnit  ** l * hl(l)
215 216 217 218 219 220
    end do     

  end subroutine ModSphBesselReal


end module m_SphBessel