csphjy Subroutine

subroutine csphjy ( n, z, nm, csj, cdj, csy, cdy )

  !*****************************************************************************80
  !
  !! CSPHJY: spherical Bessel functions jn(z) and yn(z) for complex argument.
  !
  !  Discussion:
  !
  !    This procedure computes spherical Bessel functions jn(z) and yn(z)
  !    and their derivatives for a complex argument.
  !
  !  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:
  !
  !    01 August 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 of jn(z) and yn(z).
  !
  !    Input, complex ( kind = 8 ) Z, the argument.
  !
  !    Output, integer ( kind = 4 ) NM, the highest order computed.
  !
  !    Output, complex ( kind = 8 ) CSJ(0:N0, CDJ(0:N), CSY(0:N), CDY(0:N),
  !    the values of jn(z), jn'(z), yn(z), yn'(z).
  !
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) a0
  complex ( kind = 8 ) csj(0:n)
  complex ( kind = 8 ) cdj(0:n)
  complex ( kind = 8 ) csy(0:n)
  complex ( kind = 8 ) cdy(0:n)
  complex ( kind = 8 ) cf
  complex ( kind = 8 ) cf0
  complex ( kind = 8 ) cf1
  complex ( kind = 8 ) cs
  complex ( kind = 8 ) csa
  complex ( kind = 8 ) csb
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
    ! integer ( kind = 4 ) msta1
  ! integer ( kind = 4 ) msta2
  integer ( kind = 4 ) nm
  complex ( kind = 8 ) z

  a0 = abs ( z )
  nm = n

  if ( a0 < 1.0D-60 ) then
     do k = 0, n
        csj(k) = 0.0D+00
        cdj(k) = 0.0D+00
        csy(k) = -1.0D+300
        cdy(k) = 1.0D+300
     end do
     csj(0) = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     cdj(1) = cmplx ( 0.333333333333333D+00, 0.0D+00, kind = 8 )
     return
  end if

  csj(0) = sin ( z ) / z
  csj(1) = ( csj(0) - cos ( z ) ) / z

  if ( 2 <= n ) then
     csa = csj(0)
     csb = csj(1)
     m = msta1 ( a0, 200 )
     if ( m < n ) then
        nm = m
     else
        m = msta2 ( a0, n, 15 )
     end if
     cf0 = 0.0D+00
     cf1 = 1.0D+00-100
     do k = m, 0, -1
        cf = ( 2.0D+00 * k + 3.0D+00 ) * cf1 / z - cf0
        if ( k <= nm ) then
           csj(k) = cf
        end if
        cf0 = cf1
        cf1 = cf
     end do

     if ( abs ( csa ) <= abs ( csb ) ) then
        cs = csb / cf0
     else
        cs = csa / cf
     end if

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

  end if

  cdj(0) = ( cos ( z ) - sin ( z ) / z ) / z
  do k = 1, nm
     cdj(k) = csj(k-1) - ( k + 1.0D+00 ) * csj(k) / z
  end do
  csy(0) = - cos ( z ) / z
  csy(1) = ( csy(0) - sin ( z ) ) / z
  cdy(0) = ( sin ( z ) + cos ( z ) / z ) / z
  cdy(1) = ( 2.0D+00 * cdy(0) - cos ( z ) )  / z

  do k = 2, nm
     if ( abs ( csj(k-2) ) < abs ( csj(k-1) ) ) then 
        csy(k) = ( csj(k) * csy(k-1) - 1.0D+00 / ( z * z ) ) / csj(k-1)
     else
        csy(k) = ( csj(k) * csy(k-2) &
             - ( 2.0D+00 * k - 1.0D+00 ) / z ** 3 ) / csj(k-2)
     end if
  end do

  do k = 2, nm
     cdy(k) = csy(k-1) - ( k + 1.0D+00 ) * csy(k) / z
  end do

  return
end subroutine csphjy