Line data    Source code

       1             :       MODULE m_sphbes
       2             :       use m_juDFT
       3             :       IMPLICIT NONE
       4             : c********************************************************************
       5             : c calculate spherical Bessel functions and derivatives of sqrt(e)*r
       6             : c                                   P. Marksteiner and E. Badralexe
       7             : c********************************************************************
       8             :       CONTAINS
       9   757007527 :       SUBROUTINE sphbes(
      10             :      >                  lmax,x,
      11   757007527 :      <                  fj)
      12             : !$acc routine
      13             : 
      14             : !     ..
      15             : !     .. Arguments ..
      16             :       INTEGER, INTENT  (IN) :: lmax
      17             :       REAL,    INTENT  (IN) :: x
      18             :       REAL,    INTENT (OUT) :: fj(0:lmax)
      19             : !
      20             : !     .. Parameters ..
      21             :       REAL,    PARAMETER :: small = 1.0e-03 , zero = 0.0
      22             : !     ..
      23             : !     .. Locals ..
      24             :       INTEGER i,l,min,n
      25             :       REAL fac,quot,xinv,xx
      26   757007527 :       REAL :: aux(0:int(lmax+10+x))
      27             :       !REAL, ALLOCATABLE :: aux(:)
      28             : !     ..
      29             : !     .. Intrinsic Functions ..
      30             :       INTRINSIC abs,cos,sin
      31             : 
      32   757007527 :       IF (lmax==0) THEN
      33   693165539 :          IF (x.GE.small) THEN
      34   557011205 :             fj(0) = sin(x)/x
      35             :          ELSE
      36   136154334 :             fj(0) = 1 - x*x/6.* (1.-x*x/20.* (1.-x*x/42.))
      37             :          END IF
      38   693165539 :          RETURN
      39             :       ENDIF
      40             : !#ifdef _OPENACC
      41             : !      IF (x.LT.zero)  stop "sphbes2"
      42             : !#else
      43             : !      IF (x.LT.zero)  CALL juDFT_error("sphbes2",calledby="sphbes")
      44             : !#endif
      45    63841988 :       xx = x*x
      46    63841988 :       IF (x.GE.small) THEN
      47    53694936 :          xinv = 1./x
      48    53694936 :          fj(0) = sin(x)*xinv
      49    53694936 :          fj(1) = (fj(0)-cos(x))*xinv
      50             :       ELSE
      51    10147052 :          fj(0) = 1 - xx/6.* (1.-xx/20.* (1.-xx/42.))
      52    10147052 :          fj(1) = (x/3.)* (1.-xx/10.* (1.-xx/28.))
      53    10147052 :          fac = xx/15.
      54    87669342 :          DO l=2,lmax
      55    77522290 :            fj(l) = fac * ( 1. - xx / (4*l+6) )
      56    87669342 :            fac = x * fac / (2*l+3)
      57             :          ENDDO
      58             :          RETURN
      59             :       END IF
      60    53694936 :       IF (lmax.LT.x) THEN
      61             : 
      62    49109337 :          DO l = 2,lmax
      63    49109337 :             fj(l) = (2*l-1)*xinv*fj(l-1) - fj(l-2)
      64             :          ENDDO
      65             : 
      66    48432502 :       ELSE IF (lmax.GE.2) THEN
      67    48432502 :          n = INT( lmax + 10 + x )
      68             :        !  ALLOCATE( aux(0:n) )
      69             : !
      70             : ! downward recursion from arbitrary starting values
      71             : !
      72    48432502 :          aux(n) = 0.
      73    48432502 :          aux(n-1) = 1.
      74   954368933 :          DO i = n - 1,1,-1
      75   954368933 :             aux(i-1) = (2*i+1)*xinv*aux(i) - aux(i+1)
      76             :          ENDDO
      77             : !
      78             : ! normalize with j0 or j1, whichever is larger
      79             : !
      80    48432502 :          min = 0
      81    48432502 :          IF (abs(fj(0)).LT.abs(fj(1))) min = 1
      82    48432502 :          quot = fj(min)/aux(min)
      83   412571346 :          DO l = 2,lmax
      84   412571346 :             fj(l) = aux(l)*quot
      85             :          ENDDO
      86             :         ! DEALLOCATE( aux )
      87             :       END IF
      88             : 
      89             :       RETURN
      90             :       END SUBROUTINE sphbes
      91             :       END MODULE m_sphbes