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_cylbes 8 : use m_juDFT 9 : c******************************************************************** 10 : c generates cylindrical Bessel functions for a given x and for the orders 11 : c from mmax to -mmax 12 : c Y. Mokrousov 13 : c******************************************************************** 14 : CONTAINS 15 0 : SUBROUTINE cylbes( 16 : > mmax,x, 17 0 : < fJ) 18 : 19 : IMPLICIT NONE 20 : ! .. 21 : ! ..Arguments .. 22 : INTEGER, INTENT (IN) :: mmax 23 : REAL, INTENT (IN) :: x 24 : REAL, INTENT (OUT) :: fJ(-mmax:mmax) 25 : ! 26 : ! .. Parameters .. 27 : REAL, PARAMETER :: zero = 0.0 28 : ! ..Locals .. 29 : INTEGER :: m,i,mass 30 : REAL :: quot 31 0 : REAL, ALLOCATABLE :: aux(:) 32 : ! .. 33 : 34 0 : IF (x.LT.zero) CALL juDFT_error("cylbes2",calledby="cylbes") 35 : 36 0 : IF (x.EQ.zero) THEN 37 0 : fJ(0) = 1. 38 : 39 0 : DO m=1,mmax 40 0 : fJ(m) = 0. 41 0 : fJ(-m) = 0. 42 : END DO 43 : RETURN 44 : END IF 45 : 46 0 : mass = INT( mmax + 50 + x ) 47 0 : ALLOCATE ( aux(0:mass) ) 48 0 : aux(mass) = 0.0 49 0 : aux(mass-1) = 1.0e-22 50 : 51 0 : DO i=mass-2,0,-1 52 0 : aux(i) = 2*(i+1)*aux(i+1)/x - aux(i+2) 53 : END DO 54 : 55 0 : quot = aux(0) 56 : 57 0 : DO i=1,INT( mass/2. ) 58 0 : quot = quot + 2*aux(2*i) 59 : END DO 60 : 61 0 : fJ(0) = aux(0)/quot 62 : 63 0 : DO m=1,mmax 64 0 : fJ(m) = aux(m)/quot 65 0 : fJ(-m) = ((-1)**m)*fJ(m) 66 : END DO 67 : 68 0 : DEALLOCATE ( aux ) 69 : 70 0 : RETURN 71 : END SUBROUTINE cylbes 72 : END MODULE m_cylbes 73 : 74 : 75 : 76 : 77 : 78 :