cikna Subroutine

subroutine cikna ( n, z, nm, cbi, cdi, cbk, cdk )

  !*****************************************************************************80
  !
  !! CIKNA: modified Bessel functions In(z), Kn(z), derivatives, 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:
  !
  !    30 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 ) CBI((0:N), CDI(0:N), CBK(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 ) cbi(0:n)
  complex ( kind = 8 ) cbi0
  complex ( kind = 8 ) cbi1
  complex ( kind = 8 ) cbk(0:n)
  complex ( kind = 8 ) cbk0
  complex ( kind = 8 ) cbk1
  complex ( kind = 8 ) cdi(0:n)
  complex ( kind = 8 ) cdi0
  complex ( kind = 8 ) cdi1
  complex ( kind = 8 ) cdk(0:n)
  complex ( kind = 8 ) cdk0
  complex ( kind = 8 ) cdk1
  complex ( kind = 8 ) cf
  complex ( kind = 8 ) cf1
  complex ( kind = 8 ) cf2
  complex ( kind = 8 ) ckk
  complex ( kind = 8 ) cs
  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-100 ) then
     do k = 0, n
        cbi(k) = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
        cdi(k) = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
        cbk(k) = - cmplx ( 1.0D+30, 0.0D+00, kind = 8 )
        cdk(k) = cmplx ( 1.0D+30, 0.0D+00, kind = 8 )
     end do
     cbi(0) = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     cdi(1) = cmplx ( 0.5D+00, 0.0D+00, kind = 8 )
     return
  end if

  call cik01 ( z, cbi0, cdi0, cbi1, cdi1, cbk0, cdk0, cbk1, cdk1 )

  cbi(0) = cbi0
  cbi(1) = cbi1
  cbk(0) = cbk0
  cbk(1) = cbk1
  cdi(0) = cdi0
  cdi(1) = cdi1
  cdk(0) = cdk0
  cdk(1) = cdk1

  if ( n <= 1 ) then
     return
  end if

  m = msta1 ( a0, 200 )

  if ( m < n ) then
     nm = m
  else
     m = msta2 ( a0, n, 15 )
  end if

  cf2 = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
  cf1 = cmplx ( 1.0D-30, 0.0D+00, kind = 8 )
  do k = m, 0, -1
     cf = 2.0D+00 * ( k + 1.0D+00 ) / z * cf1 + cf2
     if ( k <= nm ) then
        cbi(k) = cf
     end if
     cf2 = cf1
     cf1 = cf
  end do

  cs = cbi0 / cf
  do k = 0, nm
     cbi(k) = cs * cbi(k)
  end do

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

  do k = 2, nm
     cdi(k) =   cbi(k-1) - k / z * cbi(k)
     cdk(k) = - cbk(k-1) - k / z * cbk(k)
  end do

  return
end subroutine cikna