LCOV - code coverage report
Current view: top level - math - sphbes.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 33 33 100.0 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

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

Generated by: LCOV version 1.13