sphbes.f 2.12 KB
 Markus Betzinger committed Apr 26, 2016 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 73 74 75 76 77 78 79 `````` MODULE m_sphbes use m_juDFT c******************************************************************** c calculate spherical Bessel functions and derivatives of sqrt(e)*r c P. Marksteiner and E. Badralexe c******************************************************************** CONTAINS SUBROUTINE sphbes( > lmax,x, < fj) IMPLICIT NONE ! .. ! .. Arguments .. INTEGER, INTENT (IN) :: lmax REAL, INTENT (IN) :: x REAL, INTENT (OUT) :: fj(0:lmax) ! ! .. Parameters .. REAL, PARAMETER :: small = 1.0e-03 , zero = 0.0 ! .. ! .. Locals .. INTEGER i,l,min,n REAL fac,quot,xinv,xx REAL :: aux(0:int(lmax+10+x)) !REAL, ALLOCATABLE :: aux(:) ! .. ! .. Intrinsic Functions .. INTRINSIC abs,cos,sin IF (lmax.LT.1) CALL juDFT_error("sphbes1",calledby="sphbes") IF (x.LT.zero) CALL juDFT_error("sphbes2",calledby="sphbes") xx = x*x IF (x.GE.small) THEN xinv = 1./x fj(0) = sin(x)*xinv fj(1) = (fj(0)-cos(x))*xinv ELSE fj(0) = 1 - xx/6.* (1.-xx/20.* (1.-xx/42.)) fj(1) = (x/3.)* (1.-xx/10.* (1.-xx/28.)) fac = xx/15. DO l=2,lmax fj(l) = fac * ( 1. - xx / (4*l+6) ) fac = x * fac / (2*l+3) ENDDO RETURN END IF IF (lmax.LT.x) THEN DO l = 2,lmax fj(l) = (2*l-1)*xinv*fj(l-1) - fj(l-2) ENDDO ELSE IF (lmax.GE.2) THEN n = INT( lmax + 10 + x ) ! ALLOCATE( aux(0:n) ) ! ! downward recursion from arbitrary starting values ! aux(n) = 0. aux(n-1) = 1. DO i = n - 1,1,-1 aux(i-1) = (2*i+1)*xinv*aux(i) - aux(i+1) ENDDO ! ! normalize with j0 or j1, whichever is larger ! min = 0 IF (abs(fj(0)).LT.abs(fj(1))) min = 1 quot = fj(min)/aux(min) DO l = 2,lmax fj(l) = aux(l)*quot ENDDO ! DEALLOCATE( aux ) END IF RETURN END SUBROUTINE sphbes END MODULE m_sphbes``````