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_modcylk
       8             :       use m_juDFT
       9             : c     generates Makdonald's function K_m(x) for a given
      10             : c     order m and point x
      11             : c                         Y.Mokrousov
      12             :       CONTAINS
      13           0 :       SUBROUTINE modcylk(
      14             :      >                   m,x,
      15             :      <                   kJ) 
      16             :       
      17             : 
      18             :       implicit none
      19             :                              
      20             :       integer, intent  (in) :: m
      21             :       real,    intent  (in) :: x
      22             :       real,    intent (out) :: kJ
      23             :                                
      24             :       real,    parameter :: zero = 0.0
      25             :       real,    parameter :: gamma = 
      26             :      =     0.5772156649015328608
      27             :                                
      28             :       integer :: mm,i,mass,fact,j,m1
      29             :       real :: quot,t,k0,k1,kJ1,kJ2,psi,a
      30             : c     real, allocatable :: aux(:)
      31             :                                 
      32             :     
      33           0 :       if (x.le.zero)  CALL juDFT_error("x.le.zero",calledby="modcylk")
      34             : 
      35             : c     K_0 and K_1 part
      36             :       
      37           0 :       if (x.gt.1.0) then
      38           0 :          t = 1./x
      39             :          k0 = exp(-x)*sqrt(t)*(1.2533141373 -
      40             :      -    0.1566641816*t      + 0.0881112782*(t**2)  -
      41             :      -    0.0913909546*(t**3) + 0.1344569228*(t**4)  -
      42             :      -    0.2299850328*(t**5) + 0.3792409730*(t**6)  -
      43             :      -    0.5247277331*(t**7) + 0.5575368367*(t**8)  -
      44             :      -    0.4262632912*(t**9) + 0.2184518096*(t**10) - 
      45           0 :      -    0.0668097672*(t**11)+ 0.0091893830*(t**12))
      46             :          k1 = exp(-x)*sqrt(t)*(1.2533141373 +
      47             :      +    0.4699927013*t      - 0.1468582957*(t**2)  +
      48             :      +    0.1280426636*(t**3) - 0.1736431637*(t**4)  +
      49             :      +    0.2847618149*(t**5) - 0.4594342117*(t**6)  +
      50             :      +    0.6283380681*(t**7) - 0.6632295430*(t**8)  +
      51             :      +    0.5050238576*(t**9) - 0.2581303765*(t**10) + 
      52           0 :      +    0.0788000118*(t**11)- 0.0108241775*(t**12))
      53             :       end if
      54           0 :       if (x.le.1.0) then
      55           0 :          t = 1./x
      56           0 :          k0 = - gamma - log(x/2.)
      57           0 :          k1 = t
      58           0 :          do i = 1,8
      59             :             psi = -gamma
      60             :             fact = 1
      61           0 :             do j = 1,i
      62           0 :                psi = psi + 1./j
      63           0 :                fact = fact*j
      64             :             end do
      65           0 :             if (i.le.6) then
      66             :             k0 = k0 + ((x/2.)**(2*i))*(psi-log(x/2.))/
      67           0 :      /           (fact*fact)
      68             :             end if
      69             :             k1 = k1 + ((x/2.)**(2*i-1))*(0.5 - i*
      70           0 :      *           (psi - log(x/2.)))/(fact*fact)
      71             :          end do
      72             :       end if
      73           0 :       if (m.eq.0) then
      74           0 :          kJ = k0
      75           0 :          return
      76             :       end if
      77           0 :       if (m.eq.1 .or. m.eq.-1) then 
      78           0 :          kJ = k1
      79           0 :          return
      80             :       end if
      81             : 
      82             : c     forward recursion
      83             : 
      84           0 :       kJ1 = k0
      85           0 :       kJ2 = k1
      86           0 :       m1 = int(abs(m))
      87           0 :       do mm = 2,m1
      88           0 :          a = kJ2
      89           0 :          kJ2 = 2*(mm-1)*t*kJ2 + kJ1
      90           0 :          kJ1 = a
      91             :       end do
      92           0 :       kJ = kJ2
      93             :      
      94             :       END SUBROUTINE modcylk
      95             :       END MODULE m_modcylk