csphik Subroutine

subroutine csphik ( n, z, nm, csi, cdi, csk, cdk )

  !*****************************************************************************80
  !
  !! CSPHIK: complex modified spherical Bessel functions and 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:
  !
  !    29 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 of in(z) and kn(z).
  !
  !    Input, complex ( kind = 8 ) Z, the argument.
  !
  !    Output, integer ( kind = 4 ) NM, the highest order computed.
  !
  !    Output, complex ( kind = 8 ) CSI(0:N), CDI(0:N), CSK(0:N), CDK(0:N),
  !    the values of in(z), in'(z), kn(z), kn'(z).
  !
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) a0
  complex ( kind = 8 ) ccosh1
  complex ( kind = 8 ) cdi(0:n)
  complex ( kind = 8 ) cdk(0:n)
  complex ( kind = 8 ) cf
  complex ( kind = 8 ) cf0
  complex ( kind = 8 ) cf1
  complex ( kind = 8 ) ci
  complex ( kind = 8 ) cs
  complex ( kind = 8 ) csi(0:n)
  complex ( kind = 8 ) csi0
  complex ( kind = 8 ) csi1
  complex ( kind = 8 ) csinh1
  complex ( kind = 8 ) csk(0:n)
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
    ! integer ( kind = 4 ) msta1
  ! integer ( kind = 4 ) msta2
  integer ( kind = 4 ) nm
  real ( kind = 8 ) pi
  complex ( kind = 8 ) z

  pi = 3.141592653589793D+00
  a0 = abs ( z )    
  nm = n

  if ( a0 < 1.0D-60 ) then
     do k = 0, n
        csi(k) = 0.0D+00
        cdi(k) = 0.0D+00
        csk(k) = 1.0D+300
        cdk(k) = -1.0D+300
     end do
     csi(0) = 1.0D+00
     cdi(1) = 0.3333333333333333D+00
     return
  end if

  ci = cmplx ( 0.0D+00, 1.0D+00, kind = 8 )
  csinh1 = sin ( ci * z ) / ci
  ccosh1 = cos ( ci * z )
  csi0 = csinh1 / z
  csi1 = ( - csinh1 / z + ccosh1 ) / z
  csi(0) = csi0
  csi(1) = csi1

  if ( 2 <= n ) then

     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
           csi(k) = cf
        end if
        cf0 = cf1
        cf1 = cf
     end do

     if ( abs ( csi0 ) <= abs ( csi1 ) ) then
        cs = csi1 / cf0
     else
        cs = csi0 / cf
     end if

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

  end if

  cdi(0) = csi(1)
  do k = 1, nm
     cdi(k) = csi(k-1) - ( k + 1.0D+00 ) * csi(k) / z
  end do

  csk(0) = 0.5D+00 * pi / z * exp ( - z )
  csk(1) = csk(0) * ( 1.0D+00 + 1.0D+00 / z )
  do k = 2, nm
     if ( abs ( csi(k-2) ) < abs ( csi(k-1) ) ) then
        csk(k) = ( 0.5D+00 * pi / ( z * z ) - csi(k) * csk(k-1) ) / csi(k-1)
     else
        csk(k) = ( csi(k) * csk(k-2) + ( k - 0.5D+00 ) * pi / z ** 3 ) / csi(k-2)
     end if
  end do

  cdk(0) = -csk(1)
  do k = 1, nm
     cdk(k) = - csk(k-1) - ( k + 1.0D+00 ) * csk(k) / z
  end do

  return
end subroutine csphik