ciknb Subroutine

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

  !*****************************************************************************80
  !
  !! CIKNB computes complex modified Bessel functions In(z) and Kn(z).
  !
  !  Discussion:
  !
  !    This procedure also evaluates the 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:
  !
  !    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 ) CB((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 ) c
  complex ( kind = 8 ) ca0
  complex ( kind = 8 ) cbi(0:n)
  complex ( kind = 8 ) cbkl
  complex ( kind = 8 ) cbs
  complex ( kind = 8 ) cdi(0:n)
  complex ( kind = 8 ) cbk(0:n)
  complex ( kind = 8 ) cdk(0:n)
  complex ( kind = 8 ) cf
  complex ( kind = 8 ) cf0
  complex ( kind = 8 ) cf1
  complex ( kind = 8 ) cg
  complex ( kind = 8 ) cg0
  complex ( kind = 8 ) cg1
  complex ( kind = 8 ) ci
  complex ( kind = 8 ) cr
  complex ( kind = 8 ) cs0
  complex ( kind = 8 ) csk0
  real ( kind = 8 ) el
  real ( kind = 8 ) fac
  integer ( kind = 4 ) k
  integer ( kind = 4 ) k0
  integer ( kind = 4 ) l
  integer ( kind = 4 ) m
    ! integer ( kind = 4 ) msta1
  ! integer ( kind = 4 ) msta2
  integer ( kind = 4 ) nm
  real ( kind = 8 ) pi
  real ( kind = 8 ) vt
  complex ( kind = 8 ) z
  complex ( kind = 8 ) z1

  pi = 3.141592653589793D+00
  el = 0.57721566490153D+00
  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 )
        cbk(k) = cmplx ( 1.0D+30, 0.0D+00, kind = 8 )
        cdi(k) = cmplx ( 0.0D+00, 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

  ci = cmplx ( 0.0D+00, 1.0D+00, kind = 8 )

  if ( real ( z, kind = 8 ) < 0.0D+00 ) then
     z1 = -z
  else
     z1 = z
  end if

  if ( n == 0 ) then
     nm = 1
  end if

  m = msta1 ( a0, 200 )

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

  cbs = 0.0D+00
  csk0 = 0.0D+00
  cf0 = 0.0D+00
  cf1 = 1.0D-100

  do k = m, 0, -1
     cf = 2.0D+00 * ( k + 1.0D+00 ) * cf1 / z1 + cf0
     if ( k <= nm ) then
        cbi(k) = cf
     end if
     if ( k /= 0 .and. k == 2 * int ( k / 2 ) ) then
        csk0 = csk0 + 4.0D+00 * cf / k
     end if
     cbs = cbs + 2.0D+00 * cf
     cf0 = cf1
     cf1 = cf
  end do

  cs0 = exp ( z1 ) / ( cbs - cf )

  do k = 0, nm
     cbi(k) = cs0 * cbi(k)
  end do

  if ( a0 <= 9.0D+00 ) then

     cbk(0) = - ( log ( 0.5D+00 * z1 ) + el ) * cbi(0) + cs0 * csk0
     cbk(1) = ( 1.0D+00 / z1 - cbi(1) * cbk(0) ) / cbi(0)

  else

     ca0 = sqrt ( pi / ( 2.0D+00 * z1 ) ) * exp ( -z1 )

     if ( a0 < 25.0D+00 ) then
        k0 = 16
     else if ( a0 < 80.0D+00 ) then
        k0 = 10
     else if ( a0 < 200.0D+00 ) then
        k0 = 8
     else
        k0 = 6
     end if

     do l = 0, 1
        cbkl = 1.0D+00
        vt = 4.0D+00 * l
        cr = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
        do k = 1, k0
           cr = 0.125D+00 * cr &
                * ( vt - ( 2.0D+00 * k - 1.0D+00 ) ** 2 ) / ( k * z1 )
           cbkl = cbkl + cr
        end do
        cbk(l) = ca0 * cbkl
     end do
  end if

  cg0 = cbk(0)
  cg1 = cbk(1)
  do k = 2, nm
     cg = 2.0D+00 * ( k - 1.0D+00 ) / z1 * cg1 + cg0
     cbk(k) = cg
     cg0 = cg1
     cg1 = cg
  end do

  if ( real ( z, kind = 8 ) < 0.0D+00 ) then
     fac = 1.0D+00
     do k = 0, nm
        if ( imag ( z ) < 0.0D+00 ) then
           cbk(k) = fac * cbk(k) + ci * pi * cbi(k)
        else
           cbk(k) = fac * cbk(k) - ci * pi * cbi(k)
        end if
        cbi(k) = fac * cbi(k)
        fac = - fac
     end do
  end if

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

  return
end subroutine ciknb