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_modcyli
       8             :       use m_juDFT
       9             : 
      10             : c     generates modified bessel functions I_m for a given x and
      11             : c     order m
      12             : c                Y.Mokrousov, 21 oct 2002
      13             :       CONTAINS
      14           0 :       SUBROUTINE modcyli(
      15             :      >                   m1,x,
      16             :      <                   iJ)
      17             : 
      18             : 
      19             :       USE m_constants, ONLY : pimach
      20             :       IMPLICIT NONE
      21             : 
      22             :       integer, intent  (in) :: m1
      23             :       real,    intent  (in) :: x
      24             :       real,    intent (out) :: iJ
      25             : 
      26             :       real, parameter :: zero = 0.0
      27             :       real, parameter :: twelve = 12.0
      28             : 
      29             :       integer :: mm,i,mass,m,j
      30             :       real    :: quot,sm,qm,iJ0,fact1,fact2,iiJ,nk,fact
      31           0 :       real, allocatable :: aux(:)
      32             : 
      33             :       intrinsic abs,int,sqrt
      34             : 
      35           0 :       m = int(abs(m1))
      36             : 
      37           0 :       iJ = 0.
      38           0 :       if (x.lt.zero)  CALL juDFT_error("x.lt.zero",calledby="modcyli")
      39           0 :       if (x.eq.zero) then
      40           0 :         if (m.eq.0) then
      41           0 :            iJ = 1.0
      42           0 :            return
      43             :         else
      44           0 :            iJ = 0.0
      45           0 :            return
      46             :         end if
      47             :       else
      48           0 :          sm = sqrt(real(m))
      49           0 :          qm = m*m
      50           0 :          if (x.le.twelve .and. x.le.4*sm) then
      51             : c            write (*,*) '<12,<4sqrt(n)'
      52           0 :             fact1 = 1.
      53             :             fact = 1.
      54           0 :             do i = 1,m
      55           0 :                fact = fact*i
      56             :             end do
      57           0 :             iJ = ((x/2.)**m)/fact
      58           0 :             fact2 = 1.
      59           0 :             do i=1,31
      60           0 :                fact1 = fact1*i
      61           0 :                fact2 = fact
      62           0 :                do j=1,i
      63           0 :                   fact2 = fact2*(m+j)
      64             :                end do
      65             :                iJ = iJ + ((x/2.)**m)*(((x*x)/4.)**i)/
      66           0 :      /              (fact1*fact2)
      67             :             end do
      68             :             return
      69           0 :          elseif (x.gt.twelve .and. x.ge.6*qm) then
      70             : c           write (*,*) '>12,>6n2'
      71           0 :             iJ = exp(x)/sqrt(2*pimach()*x)
      72           0 :             fact1 = 1.
      73           0 :             iiJ = 1.
      74           0 :             nk = 1.
      75           0 :             do i = 1,31
      76           0 :                fact1 = fact1*i
      77           0 :                nk = nk*(4*(m*m)-(2*i-1)*(2*i-1))
      78           0 :                iiJ = iiJ + ((-1)**i)*nk/(fact1*((8*x)**i))
      79             :             end do
      80           0 :             iJ = iJ*iiJ
      81           0 :             return
      82             :          else
      83           0 :             mass = int( m + 10 + x )
      84           0 :             allocate ( aux(0:mass) )
      85           0 :             aux(mass) = 0.0
      86           0 :             aux(mass-1) = 1.0e-22
      87           0 :             do i=mass-2,0,-1
      88           0 :                aux(i) = 2*(i+1)*aux(i+1)/x + aux(i+2)
      89             : c               if (i.lt.m .and. x.gt.6*qm .and. i.ne.0) then
      90             : c                  quot = aux(i)
      91             : c                  iJ = aux(m)/quot
      92             : c                  return
      93             : c               end if
      94             :             end do
      95           0 :             quot = aux(0)
      96           0 :             if (x.le.twelve) then
      97             : c               write (*,*) 'kinda here'
      98             :                fact1 = 1.
      99             :                iJ0 = 1.
     100           0 :                do i=1,31
     101           0 :                   fact1 = fact1*i
     102           0 :                   iJ0 = iJ0 + (((x*x)/4.)**i)/(fact1*fact1)
     103             :                end do
     104           0 :                iJ = aux(m)*iJ0/quot
     105           0 :                deallocate (aux)
     106           0 :                return
     107             :             elseif (x.gt.twelve) then
     108           0 :                iJ0 = exp(x)/sqrt(2*pimach()*x)
     109           0 :                fact1 = 1.
     110           0 :                iiJ = 1.
     111           0 :                nk = 1.
     112           0 :                do i = 1,31
     113           0 :                   fact1 = fact1*i
     114           0 :                   nk = nk*(-(2*i-1)*(2*i-1))
     115           0 :                   iiJ = iiJ + ((-1)**i)*nk/(fact1*((8*x)**i))
     116             :                end do
     117           0 :                iJ0 = iJ0*iiJ
     118           0 :                iJ = aux(m)*iJ0/quot
     119           0 :                deallocate (aux)
     120           0 :                return
     121             :             end if
     122             :          end if
     123             :       end if
     124             : 
     125             :       END SUBROUTINE modcyli
     126             :       END MODULE m_modcyli