cjyvb Subroutine

subroutine cjyvb ( v, z, vm, cbj, cdj, cby, cdy )

  !*****************************************************************************80
  !
  !! CJYVB: Bessel functions and derivatives, Jv(z) and Yv(z) of 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:
  !
  !    03 August 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, real ( kind = 8 ) VM, the highest order computed.
  !
  !    Output, real ( kind = 8 ) CBJ(0:*), CDJ(0:*), CBY(0:*), CDY(0:*), 
  !    the values of Jn+v0(z), Jn+v0'(z), Yn+v0(z), Yn+v0'(z).
  !
  implicit none

  real ( kind = 8 ) a0
  complex ( kind = 8 ) ca
  complex ( kind = 8 ) ca0
  complex ( kind = 8 ) cb
  complex ( kind = 8 ) cbj(0:*)
  complex ( kind = 8 ) cby(0:*)
  complex ( kind = 8 ) cck
  complex ( kind = 8 ) cdj(0:*)
  complex ( kind = 8 ) cdy(0:*)
  complex ( kind = 8 ) cec
  complex ( kind = 8 ) cf
  complex ( kind = 8 ) cf1
  complex ( kind = 8 ) cf2
  complex ( kind = 8 ) cfac0
  complex ( kind = 8 ) ci
  complex ( kind = 8 ) cju0
  complex ( kind = 8 ) cjv0
  complex ( kind = 8 ) cjvn
  complex ( kind = 8 ) cpz
  complex ( kind = 8 ) cqz
  complex ( kind = 8 ) cr
  complex ( kind = 8 ) cr0
  complex ( kind = 8 ) crp
  complex ( kind = 8 ) crq
  complex ( kind = 8 ) cs
  complex ( kind = 8 ) cs0
  complex ( kind = 8 ) csk
  complex ( kind = 8 ) cyv0
  complex ( kind = 8 ) cyy
  real ( kind = 8 ) ga
  real ( kind = 8 ) gb
  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 ) pv0
  real ( kind = 8 ) rp2
  real ( kind = 8 ) v
  real ( kind = 8 ) v0
  real ( kind = 8 ) vg
  real ( kind = 8 ) vm
  real ( kind = 8 ) vv
  real ( kind = 8 ) w0
  complex ( kind = 8 ) z
  complex ( kind = 8 ) z1
  complex ( kind = 8 ) z2
  complex ( kind = 8 ) zk

  pi = 3.141592653589793D+00
  rp2 = 0.63661977236758D+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
  pv0 = pi * v0

  if ( a0 < 1.0D-100 ) then

     do k = 0, n
        cbj(k) = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
        cdj(k) = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
        cby(k) = - cmplx ( 1.0D+30, 0.0D+00, kind = 8 )
        cdy(k) = cmplx ( 1.0D+30, 0.0D+00, kind = 8 )
     end do

     if ( v0 == 0.0D+00 ) then
        cbj(0) = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
        cdj(1) = cmplx ( 0.5D+00, 0.0D+00, kind = 8 )
     else
        cdj(0) = cmplx ( 1.0D+30, 0.0D+00, kind = 8 )
     end if

     vm = v
     return

  end if

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

  if ( a0 <= 12.0D+00 ) then

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

     vg = 1.0D+00 + v0
     call gammaf ( vg, ga )
     ca = ( 0.5D+00 * z1 ) ** v0 / ga
     cjv0 = cjv0 * ca

  else

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

     vv = 4.0D+00 * v0 * v0
     cpz = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     crp = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     do k = 1, k0
        crp = -0.78125D-02 * crp &
             * ( vv - ( 4.0D+00 * k - 3.0D+00 ) ** 2 ) &
             * ( vv - ( 4.0D+00 * k - 1.0D+00 ) **2 ) &
             / ( k * ( 2.0D+00 * k - 1.0D+00 ) * z2 )
        cpz = cpz + crp
     end do
     cqz = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     crq = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     do k = 1, k0
        crq = -0.78125D-02 * crq &
             * ( vv - ( 4.0D+00 * k - 1.0D+00 ) ** 2 ) &
             * ( vv - ( 4.0D+00 * k + 1.0D+00 ) ** 2 ) &
             / ( k * ( 2.0D+00 * k + 1.0D+00 ) * z2 )
        cqz = cqz + crq
     end do
     cqz = 0.125D+00 * ( vv - 1.0D+00 ) * cqz / z1
     zk = z1 - ( 0.5D+00 * v0 + 0.25D+00 ) * pi
     ca0 = sqrt ( rp2 / z1 )
     cck = cos ( zk )
     csk = sin ( zk )
     cjv0 = ca0 * ( cpz * cck - cqz * csk )
     cyv0 = ca0 * ( cpz * csk + cqz * cck )

  end if

  if ( a0 <= 12.0D+00 ) then

     if ( v0 .ne. 0.0D+00 ) then

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

        vg = 1.0D+00 - v0
        call gammaf ( vg, gb )
        cb = ( 2.0D+00 / z1 ) ** v0 / gb
        cju0 = cjvn * cb
        cyv0 = ( cjv0 * cos ( pv0 ) - cju0 ) / sin ( pv0 )

     else

        cec = log ( z1 / 2.0D+00 ) + 0.5772156649015329D+00
        cs0 = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
        w0 = 0.0D+00
        cr0 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
        do k = 1, 30
           w0 = w0 + 1.0D+00 / k
           cr0 = -0.25D+00 * cr0 / ( k * k ) * z2
           cs0 = cs0 + cr0 * w0
        end do
        cyv0 = rp2 * ( cec * cjv0 - cs0 )

     end if

  end if

  if ( n .eq. 0 ) then
     n = 1
  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
        cbj(k) = cf
     end if
     cf2 = cf1
     cf1 = cf
  end do

  cs = cjv0 / cf
  do k = 0, n
     cbj(k) = cs * cbj(k)
  end do

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

     cfac0 = exp ( pv0 * ci )
     if ( imag ( z ) < 0.0D+00 ) then
        cyv0 = cfac0 * cyv0 - 2.0D+00 * ci * cos ( pv0 ) * cjv0
     else if ( 0.0D+00 < imag ( z ) ) then
        cyv0 = cyv0 / cfac0 + 2.0D+00 * ci * cos ( pv0 ) * cjv0
     end if

     do k = 0, n
        if ( imag ( z ) < 0.0D+00) then
           cbj(k) = exp ( - pi * ( k + v0 ) * ci ) * cbj(k)
        else if ( 0.0D+00 < imag ( z ) ) then
           cbj(k) = exp ( pi * ( k + v0 ) * ci ) * cbj(k)
        end if
     end do

     z1 = z1

  end if

  cby(0) = cyv0
  do k = 1, n
     cyy = ( cbj(k) * cby(k-1) - 2.0D+00 / ( pi * z ) ) / cbj(k-1)
     cby(k) = cyy
  end do

  cdj(0) = v0 / z * cbj(0) - cbj(1)
  do k = 1, n
     cdj(k) = - ( k + v0 ) / z * cbj(k) + cbj(k-1)
  end do

  cdy(0) = v0 / z * cby(0) - cby(1)
  do k = 1, n
     cdy(k) = cby(k-1) - ( k + v0 ) / z * cby(k)
  end do

  vm = n + v0

  return
end subroutine cjyvb