cikva Subroutine

subroutine cikva ( v, z, vm, cbi, cdi, cbk, cdk )

  !*****************************************************************************80
  !
  !! CIKVA: modified Bessel functions Iv(z), Kv(z), arbitrary order, complex.
  !
  !  Discussion:
  !
  !    Compute the modified Bessel functions Iv(z), Kv(z)
  !    and their derivatives for an arbitrary order and
  !    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:
  !
  !    31 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, real ( kind = 8 ) V, the order of the functions.
  !
  !    Input, complex ( kind = 8 ) Z, the argument.
  !
  !    Output, real ( kind = 8 ) VM, the highest order computed.
  !
  !    Output, real ( kind = 8 ) CBI(0:N), CDI(0:N), CBK(0:N), CDK(0:N),
  !    the values of In+v0(z), In+v0'(z), Kn+v0(z), Kn+v0'(z).
  !
  implicit none

  real ( kind = 8 ) a0
  complex ( kind = 8 ) ca
  complex ( kind = 8 ) ca1
  complex ( kind = 8 ) ca2
  complex ( kind = 8 ) cb
  complex ( kind = 8 ) cbi(0:*)
  complex ( kind = 8 ) cbi0
  complex ( kind = 8 ) cdi(0:*)
  complex ( kind = 8 ) cbk(0:*)
  complex ( kind = 8 ) cbk0
  complex ( kind = 8 ) cbk1
  complex ( kind = 8 ) cdk(0:*)
  complex ( kind = 8 ) cf
  complex ( kind = 8 ) cf1
  complex ( kind = 8 ) cf2
  complex ( kind = 8 ) cg0
  complex ( kind = 8 ) cg1
  complex ( kind = 8 ) cgk
  complex ( kind = 8 ) ci
  complex ( kind = 8 ) ci0
  complex ( kind = 8 ) cp
  complex ( kind = 8 ) cr
  complex ( kind = 8 ) cr1
  complex ( kind = 8 ) cr2
  complex ( kind = 8 ) cs
  complex ( kind = 8 ) csu
  complex ( kind = 8 ) ct
  complex ( kind = 8 ) cvk
  real ( kind = 8 ) gan
  real ( kind = 8 ) gap
  integer ( kind = 4 ) k
  integer ( kind = 4 ) k0
  integer ( kind = 4 ) m
    ! integer ( kind = 4 ) msta1
  ! integer ( kind = 4 ) msta2
  integer ( kind = 4 ) n
  real ( kind = 8 ) pi
  real ( kind = 8 ) piv
  real ( kind = 8 ) v
  real ( kind = 8 ) v0
  real ( kind = 8 ) v0n
  real ( kind = 8 ) v0p
  real ( kind = 8 ) vm
  real ( kind = 8 ) vt
  real ( kind = 8 ) w0
  real ( kind = 8 ) ws
  real ( kind = 8 ) ws0
  complex ( kind = 8 ) z
  complex ( kind = 8 ) z1
  complex ( kind = 8 ) z2

  pi = 3.141592653589793D+00
  ci = cmplx ( 0.0D+00, 1.0D+00, kind = 8 )
  a0 = abs ( z )
  z1 = z
  z2 = z * z
  n = int ( v )
  v0 = v - n
  piv = pi * v0
  vt = 4.0D+00 * v0 * v0

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

  if ( a0 < 1.0D-100 ) then

     do k = 0, n
        cbi(k) = 0.0D+00
        cdi(k) = 0.0D+00
        cbk(k) = -1.0D+300
        cdk(k) = 1.0D+300
     end do

     if ( v0 == 0.0D+00 ) then
        cbi(0) = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
        cdi(1) = cmplx ( 0.5D+00, 0.0D+00, kind = 8 )
     end if

     vm = v
     return

  end if

  if ( a0 < 35.0D+00 ) then
     k0 = 14
  else if ( a0 < 50.0D+00 ) then
     k0 = 10
  else
     k0 = 8
  end if

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

  if ( a0 < 18.0D+00 ) then

     if ( v0 == 0.0D+00 ) then
        ca1 = cmplx (1.0D+00, 0.0D+00, kind = 8 )
     else
        v0p = 1.0D+00 + v0
        call gammaf ( v0p, gap )
        ca1 = ( 0.5D+00 * z1 ) ** v0 / gap
     end if

     ci0 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     cr = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     do k = 1, 50
        cr = 0.25D+00 * cr * z2 / ( k * ( k + v0 ) )
        ci0 = ci0 + cr
        if ( abs ( cr ) < abs ( ci0 ) * 1.0D-15 ) then
           exit
        end if
     end do

     cbi0 = ci0 * ca1

  else

     ca = exp ( z1 ) / sqrt ( 2.0D+00 * pi * z1 )
     cs = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     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 )
        cs = cs + cr
     end do
     cbi0 = ca * cs

  end if

  m = msta1 ( a0, 200 )

  if ( m < n ) then
     n = 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 * ( v0 + k + 1.0D+00 ) / z1 * cf1 + cf2
     if ( k <= n ) then
        cbi(k) = cf
     end if
     cf2 = cf1
     cf1 = cf
  end do

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

  if ( a0 <= 9.0D+00 ) then

     if ( v0 == 0.0D+00 ) then
        ct = - log ( 0.5D+00 * z1 ) - 0.5772156649015329D+00
        cs = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
        w0 = 0.0D+00
        cr = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
        do k = 1, 50
           w0 = w0 + 1.0D+00 / k
           cr = 0.25D+00 * cr / ( k * k ) * z2
           cp = cr * ( w0 + ct )
           cs = cs + cp
           if ( 10 <= k .and. abs ( cp / cs ) < 1.0D-15 ) then
              exit
           end if
        end do

        cbk0 = ct + cs

     else

        v0n = 1.0D+00 - v0
        call gammaf ( v0n, gan )
        ca2 = 1.0D+00 / ( gan * ( 0.5D+00 * z1 ) ** v0 )
        ca1 = ( 0.5D+00 * z1 ) ** v0 / gap
        csu = ca2 - ca1
        cr1 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
        cr2 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
        do k = 1, 50
           cr1 = 0.25D+00 * cr1 * z2 / ( k * ( k - v0 ) )
           cr2 = 0.25D+00 * cr2 * z2 / ( k * ( k + v0 ) )
           csu = csu + ca2 * cr1 - ca1 * cr2
           ws = abs ( csu )
           if ( 10 <= k .and. abs ( ws - ws0 ) / ws < 1.0D-15 ) then
              exit
           end if
           ws0 = ws
        end do

        cbk0 = 0.5D+00 * pi * csu / sin ( piv )

     end if

  else

     cb = exp ( - z1 ) * sqrt ( 0.5D+00 * pi / z1 )
     cs = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     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 )
        cs = cs + cr
     end do
     cbk0 = cb * cs

  end if

  cbk1 = ( 1.0D+00 / z1 - cbi(1) * cbk0 ) / cbi(0)
  cbk(0) = cbk0
  cbk(1) = cbk1
  cg0 = cbk0
  cg1 = cbk1

  do k = 2, n
     cgk = 2.0D+00 * ( v0 + k - 1.0D+00 ) / z1 * cg1 + cg0
     cbk(k) = cgk
     cg0 = cg1
     cg1 = cgk
  end do

  if ( real ( z, kind = 8 ) < 0.0D+00 ) then
     do k = 0, n
        cvk = exp ( ( k + v0 ) * pi * ci )
        if ( imag ( z ) < 0.0D+00 ) then
           cbk(k) = cvk * cbk(k) + pi * ci * cbi(k)
           cbi(k) = cbi(k) / cvk
        else if ( 0.0D+00 < imag ( z ) ) then
           cbk(k) = cbk(k) / cvk - pi * ci * cbi(k)
           cbi(k) = cvk * cbi(k)
        end if
     end do
  end if

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

  vm = n + v0

  return
end subroutine cikva