!--------------------------------------------------------------------------------
! 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_modcylk
      use m_juDFT
c     generates Makdonald's function K_m(x) for a given
c     order m and point x
c                         Y.Mokrousov
      CONTAINS
      SUBROUTINE modcylk(
     >                   m,x,
     <                   kJ) 
      

      implicit none
                             
      integer, intent  (in) :: m
      real,    intent  (in) :: x
      real,    intent (out) :: kJ
                               
      real,    parameter :: zero = 0.0
      real,    parameter :: gamma = 
     =     0.5772156649015328608
                               
      integer :: mm,i,mass,fact,j,m1
      real :: quot,t,k0,k1,kJ1,kJ2,psi,a
c     real, allocatable :: aux(:)
                                
    
      if (x.le.zero)  CALL juDFT_error("x.le.zero",calledby="modcylk")

c     K_0 and K_1 part
      
      if (x.gt.1.0) then
         t = 1./x
         k0 = exp(-x)*sqrt(t)*(1.2533141373 -
     -    0.1566641816*t      + 0.0881112782*(t**2)  -
     -    0.0913909546*(t**3) + 0.1344569228*(t**4)  -
     -    0.2299850328*(t**5) + 0.3792409730*(t**6)  -
     -    0.5247277331*(t**7) + 0.5575368367*(t**8)  -
     -    0.4262632912*(t**9) + 0.2184518096*(t**10) - 
     -    0.0668097672*(t**11)+ 0.0091893830*(t**12))
         k1 = exp(-x)*sqrt(t)*(1.2533141373 +
     +    0.4699927013*t      - 0.1468582957*(t**2)  +
     +    0.1280426636*(t**3) - 0.1736431637*(t**4)  +
     +    0.2847618149*(t**5) - 0.4594342117*(t**6)  +
     +    0.6283380681*(t**7) - 0.6632295430*(t**8)  +
     +    0.5050238576*(t**9) - 0.2581303765*(t**10) + 
     +    0.0788000118*(t**11)- 0.0108241775*(t**12))
      end if
      if (x.le.1.0) then
         t = 1./x
         k0 = - gamma - log(x/2.)
         k1 = t
         do i = 1,8
            psi = -gamma
            fact = 1
            do j = 1,i
               psi = psi + 1./j
               fact = fact*j
            end do
            if (i.le.6) then
            k0 = k0 + ((x/2.)**(2*i))*(psi-log(x/2.))/
     /           (fact*fact)
            end if
            k1 = k1 + ((x/2.)**(2*i-1))*(0.5 - i*
     *           (psi - log(x/2.)))/(fact*fact)
         end do
      end if
      if (m.eq.0) then
         kJ = k0
         return
      end if
      if (m.eq.1 .or. m.eq.-1) then 
         kJ = k1
         return
      end if

c     forward recursion

      kJ1 = k0
      kJ2 = k1
      m1 = int(abs(m))
      do mm = 2,m1
         a = kJ2
         kJ2 = 2*(mm-1)*t*kJ2 + kJ1
         kJ1 = a
      end do
      kJ = kJ2
     
      END SUBROUTINE modcylk
      END MODULE m_modcylk