cjylv Subroutine

subroutine cjylv ( v, z, cbjv, cdjv, cbyv, cdyv )

  !*****************************************************************************80
  !
  !! CJYLV: Bessel functions Jv(z), Yv(z) of complex argument and large order v.
  !
  !  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:
  !
  !    25 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 Jv(z) and Yv(z).
  !
  !    Input, complex ( kind = 8 ) Z, the argument.
  !
  !    Output, complex ( kind = 8 ) CBJV, CDJV, CBYV, CDYV, the values of Jv(z), 
  !    Jv'(z), Yv(z), Yv'(z).
  !
  implicit none

  real ( kind = 8 ) a(91)
  complex ( kind = 8 ) cbjv
  complex ( kind = 8 ) cbyv
  complex ( kind = 8 ) cdjv
  complex ( kind = 8 ) cdyv
  complex ( kind = 8 ) ceta
  complex ( kind = 8 ) cf(12)
  complex ( kind = 8 ) cfj
  complex ( kind = 8 ) cfy
  complex ( kind = 8 ) csj
  complex ( kind = 8 ) csy
  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

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

  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
     csj = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     do k = 1, km
        csj = csj + cf(k) * vr ** k
     end do
     cbjv = sqrt ( ct / ( 2.0D+00 * pi * v0 ) ) * exp ( v0 * ceta ) * csj
     if ( l == 1 ) then
        cfj = cbjv
     end if
     csy = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     do k = 1, km
        csy = csy + ( -1.0D+00 ) ** k * cf(k) * vr ** k
     end do
     cbyv = - sqrt ( 2.0D+00 * ct / ( pi * v0 ) ) * exp ( - v0 * ceta ) * csy
     if ( l == 1 ) then
        cfy = cbyv
     end if

  end do

  cdjv = - v / z * cbjv + cfj
  cdyv = - v / z * cbyv + cfy

  return
end subroutine cjylv