dsphbs.f 1.33 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
      MODULE m_dsphbs
!********************************************************************
!     calculates the derivative of the spherical bessel functions
!     dfj(l) = d jl(x)/dx
!     for l=0,lmax and argument x
!     note that the spherical bessel functions fj(l), l=0,lmax are
!     needed (call sphbes to generate these)
!                                            m. weinert
!********************************************************************
      CONTAINS
      SUBROUTINE dsphbs(
     >                  lmax,x,fj,
     <                  dfj)

      IMPLICIT NONE
!     ..
!     .. Arguments ..
      INTEGER, INTENT  (IN) :: lmax
      REAL,    INTENT  (IN) :: x
      REAL,    INTENT  (IN) ::  fj(0:lmax)
      REAL,    INTENT (OUT) :: dfj(0:lmax)
!
!     .. Parameters ..
      REAL,    PARAMETER :: xlim = 1.0e-04
!     ..
!     .. Local Scalars ..
      REAL fac,x2
      INTEGER l

      dfj(0) = -fj(1)
!--->    small x limit
      IF (x.LT.xlim) THEN
         x2 = 0.5*x*x
         fac = 1./3.
         DO l = 1,lmax
            dfj(l) = fac* (l-x2* (l+2)/ (2*l+3))
            fac = x*fac/ (2*l+3)
         ENDDO
      ELSE
!--->    obtain dfj using recurrence relationship
         DO l = 1,lmax
            dfj(l) = fj(l-1) - (l+1)*fj(l)/x
         ENDDO
      END IF

      RETURN
      END SUBROUTINE dsphbs
      END MODULE m_dsphbs