ciklv Subroutine

subroutine ciklv ( v, z, cbiv, cdiv, cbkv, cdkv )

  !*****************************************************************************80
  !
  !! CIKLV: modified Bessel functions Iv(z), Kv(z), complex argument, large order.
  !
  !  Discussion:
  !
  !    This procedure computes modified Bessel functions Iv(z) and
  !    Kv(z) and their derivatives with a complex argument and a large order.
  !
  !  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 Iv(z) and Kv(z).
  !
  !    Input, complex ( kind = 8 ) Z, the argument.
  !
  !    Output, real ( kind = 8 ) CBIV, CDIV, CBKV, CDKV, the values of
  !    Iv(z), Iv'(z), Kv(z), Kv'(z).
  !
  implicit none

  real ( kind = 8 ) a(91)
  complex ( kind = 8 ) cbiv
  complex ( kind = 8 ) cbkv
  complex ( kind = 8 ) cdiv
  complex ( kind = 8 ) cdkv
  complex ( kind = 8 ) ceta
  complex ( kind = 8 ) cf(12)
  complex ( kind = 8 ) cfi
  complex ( kind = 8 ) cfk
  complex ( kind = 8 ) csi
  complex ( kind = 8 ) csk
  complex ( kind = 8 ) ct
  complex ( kind = 8 ) ct2
  complex ( kind = 8 ) cws
  integer ( kind = 4 ) i
  integer ( kind = 4 ) k
  integer ( kind = 4 ) km
  integer ( kind = 4 ) l
  integer ( kind = 4 ) l0
  integer ( kind = 4 ) lf
  real ( kind = 8 ) pi
  real ( kind = 8 ) v
  real ( kind = 8 ) v0
  real ( kind = 8 ) vr
  complex ( kind = 8 ) z

  pi = 3.141592653589793D+00
  km = 12
  call cjk ( km, a )

  do l = 1, 0, -1

     v0 = v - l
     cws = sqrt ( 1.0D+00 + ( z / v0 ) * ( z / v0 ) )
     ceta = cws + log ( z / v0 / ( 1.0D+00 + cws ) )
     ct = 1.0D+00 / cws
     ct2 = ct * ct
     do k = 1, km
        l0 = k * ( k + 1 ) / 2 + 1
        lf = l0 + k
        cf(k) = a(lf)
        do i = lf - 1, l0, -1
           cf(k) = cf(k) * ct2 + a(i)
        end do
        cf(k) = cf(k) * ct ** k
     end do
     vr = 1.0D+00 / v0
     csi = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     do k = 1, km
        csi = csi + cf(k) * vr ** k
     end do
     cbiv = sqrt ( ct / ( 2.0D+00 * pi * v0 ) ) * exp ( v0 * ceta ) * csi
     if ( l == 1 ) then
        cfi = cbiv
     end if
     csk = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     do k = 1, km
        csk = csk + ( - 1 ) ** k * cf(k) * vr ** k
     end do
     cbkv = sqrt ( pi * ct / ( 2.0D+00 * v0 ) ) * exp ( - v0 * ceta ) * csk

     if ( l == 1 ) then
        cfk = cbkv
     end if

  end do

  cdiv =   cfi - v / z * cbiv
  cdkv = - cfk - v / z * cbkv

  return
end subroutine ciklv