sphj Subroutine

subroutine sphj ( n, x, nm, sj, dj )

  !*****************************************************************************80
  !
  !! SPHJ computes spherical Bessel functions jn(x) and their derivatives.
  !
  !  Licensing:
  !
  !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
  !    they give permission to incorporate this routine into a user program 
  !    provided that the copyright is acknowledged.
  !
  !  Modified:
  !
  !    15 July 2012
  !
  !  Author:
  !
  !    Shanjie Zhang, Jianming Jin
  !
  !  Reference:
  !
  !    Shanjie Zhang, Jianming Jin,
  !    Computation of Special Functions,
  !    Wiley, 1996,
  !    ISBN: 0-471-11963-6,
  !    LC: QA351.C45.
  !
  !  Parameters:
  !
  !    Input, integer ( kind = 4 ) N, the order.
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, integer ( kind = 4 ) NM, the highest order computed.
  !
  !    Output, real ( kind = 8 ) SJ(0:N), the values of jn(x).
  !
  !    Output, real ( kind = 8 ) DJ(0:N), the values of jn'(x).
  !
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) cs
  real ( kind = 8 ) dj(0:n)
  real ( kind = 8 ) f
  real ( kind = 8 ) f0
  real ( kind = 8 ) f1
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  ! integer ( kind = 4 ) msta1
  ! integer ( kind = 4 ) msta2
  integer ( kind = 4 ) nm
  real ( kind = 8 ) sa
  real ( kind = 8 ) sb
  real ( kind = 8 ) sj(0:n)
  real ( kind = 8 ) x

  nm = n

  if ( abs ( x ) <= 1.0D-100 ) then
     do k = 0, n
        sj(k) = 0.0D+00
        dj(k) = 0.0D+00
     end do
     sj(0) = 1.0D+00
     dj(1) = 0.3333333333333333D+00
     return
  end if

  sj(0) = sin ( x ) / x
  sj(1) = ( sj(0) - cos ( x ) ) / x

  if ( 2 <= n ) then

     sa = sj(0)
     sb = sj(1)
     m = msta1 ( x, 200 )
     if ( m < n ) then
        nm = m
     else
        m = msta2 ( x, n, 15 )
     end if

     f0 = 0.0D+00
     f1 = 1.0D+00-100
     do k = m, 0, -1
        f = ( 2.0D+00 * k + 3.0D+00 ) * f1 / x - f0
        if ( k <= nm ) then
           sj(k) = f
        end if
        f0 = f1
        f1 = f
     end do

     if ( abs ( sa ) <= abs ( sb ) ) then
        cs = sb / f0
     else
        cs = sa / f
     end if

     do k = 0, nm
        sj(k) = cs * sj(k)
     end do

  end if

  dj(0) = ( cos(x) - sin(x) / x ) / x
  do k = 1, nm
     dj(k) = sj(k-1) - ( k + 1.0D+00 ) * sj(k) / x
  end do

  return
end subroutine sphj