cylbes.f 1.49 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
      MODULE m_cylbes     
      use m_juDFT
c********************************************************************
c generates cylindrical Bessel functions for a given x and for the orders 
c from mmax to -mmax
c                                  Y. Mokrousov
c********************************************************************
      CONTAINS
      SUBROUTINE cylbes(
     >                  mmax,x,
     <                  fJ) 

      IMPLICIT NONE
!     ..
!     ..Arguments ..
      INTEGER, INTENT  (IN) :: mmax
      REAL,    INTENT  (IN) :: x
      REAL,    INTENT (OUT) :: fJ(-mmax:mmax)
!
!     .. Parameters ..
      REAL,    PARAMETER :: zero = 0.0
!     ..Locals ..
      INTEGER :: m,i,mass
      REAL :: quot
      REAL, ALLOCATABLE :: aux(:)
!     ..

      IF (x.LT.zero)  CALL juDFT_error("cylbes2",calledby="cylbes")

      IF (x.EQ.zero) THEN
        fJ(0) = 1.

        DO m=1,mmax
         fJ(m) = 0.
         fJ(-m) = 0.
        END DO
        RETURN
      END IF 

      mass = INT( mmax + 50 + x )
      ALLOCATE ( aux(0:mass) )       
      aux(mass) = 0.0
      aux(mass-1) = 1.0e-22   
     
      DO i=mass-2,0,-1
        aux(i) = 2*(i+1)*aux(i+1)/x - aux(i+2)
      END DO

      quot = aux(0)

      DO i=1,INT( mass/2. )
        quot = quot + 2*aux(2*i)
      END DO  

      fJ(0) = aux(0)/quot

      DO m=1,mmax
        fJ(m) = aux(m)/quot
        fJ(-m) = ((-1)**m)*fJ(m)
      END DO

      DEALLOCATE ( aux )

      RETURN
      END SUBROUTINE cylbes
      END MODULE m_cylbes