cjyva Subroutine

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

  !*****************************************************************************80
  !
  !! CJYVA: 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 ), 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 ) cf0
  complex ( kind = 8 ) cf1
  complex ( kind = 8 ) cf2
  complex ( kind = 8 ) cfac0
  complex ( kind = 8 ) cfac1
  complex ( kind = 8 ) cg0
  complex ( kind = 8 ) cg1
  complex ( kind = 8 ) ch0
  complex ( kind = 8 ) ch1
  complex ( kind = 8 ) ch2
  complex ( kind = 8 ) ci
  complex ( kind = 8 ) cju0
  complex ( kind = 8 ) cju1
  complex ( kind = 8 ) cjv0
  complex ( kind = 8 ) cjv1
  complex ( kind = 8 ) cjvl
  complex ( kind = 8 ) cp11
  complex ( kind = 8 ) cp12
  complex ( kind = 8 ) cp21
  complex ( kind = 8 ) cp22
  complex ( kind = 8 ) cpz
  complex ( kind = 8 ) cqz
  complex ( kind = 8 ) cr
  complex ( kind = 8 ) cr0
  complex ( kind = 8 ) cr1
  complex ( kind = 8 ) crp
  complex ( kind = 8 ) crq
  complex ( kind = 8 ) cs
  complex ( kind = 8 ) cs0
  complex ( kind = 8 ) cs1
  complex ( kind = 8 ) csk
  complex ( kind = 8 ) cyk
  complex ( kind = 8 ) cyl1
  complex ( kind = 8 ) cyl2
  complex ( kind = 8 ) cylk
  complex ( kind = 8 ) cyv0
  complex ( kind = 8 ) cyv1
  real ( kind = 8 ) ga
  real ( kind = 8 ) gb
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) k0
  integer ( kind = 4 ) l
  integer ( kind = 4 ) lb
  integer ( kind = 4 ) lb0
  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 ) pv1
  real ( kind = 8 ) rp2
  real ( kind = 8 ) v
  real ( kind = 8 ) v0
  real ( kind = 8 ) vg
  real ( kind = 8 ) vl
  real ( kind = 8 ) vm
  real ( kind = 8 ) vv
  real ( kind = 8 ) w0
  real ( kind = 8 ) w1
  real ( kind = 8 ) wa
  real ( kind = 8 ) ya0
  real ( kind = 8 ) ya1
  real ( kind = 8 ) yak
  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
  pv1 = pi * ( 1.0D+00 + 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

     do l = 0, 1
        vl = v0 + l
        cjvl = 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 + vl ) )
           cjvl = cjvl + cr
           if ( abs ( cr ) < abs ( cjvl ) * 1.0D-15 ) then
              exit
           end if
        end do

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

        if ( l == 0 ) then
           cjv0 = cjvl * ca
        else
           cjv1 = cjvl * ca
        end if

     end do

  else

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

     do j = 0, 1
        vv = 4.0D+00 * ( j + v0 ) * ( j + 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 * ( j + v0 ) + 0.25D+00 ) * pi
        ca0 = sqrt ( rp2 / z1 )
        cck = cos ( zk )
        csk = sin ( zk )
        if ( j == 0 ) then
           cjv0 = ca0 * ( cpz * cck - cqz * csk )
           cyv0 = ca0 * ( cpz * csk + cqz * cck )
        else if ( j == 1 ) then
           cjv1 = ca0 * ( cpz * cck - cqz * csk )
           cyv1 = ca0 * ( cpz * csk + cqz * cck )
        end if
     end do

  end if

  if ( a0 <= 12.0D+00 ) then

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

        do l = 0, 1
           vl = v0 + l
           cjvl = 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 - vl ) )
              cjvl = cjvl + cr
              if ( abs ( cr ) < abs ( cjvl ) * 1.0D-15 ) then
                 exit
              end if
           end do

           vg = 1.0D+00 - vl
           call gammaf ( vg, gb )
           cb = ( 2.0D+00 / z1 ) ** vl / gb
           if ( l == 0 ) then
              cju0 = cjvl * cb
           else
              cju1 = cjvl * cb
           end if
        end do
        cyv0 = ( cjv0 * cos ( pv0 ) - cju0 ) / sin ( pv0 )
        cyv1 = ( cjv1 * cos ( pv1 ) - cju1 ) / sin ( pv1 )

     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 )
        cs1 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
        w1 = 0.0D+00
        cr1 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
        do k = 1, 30
           w1 = w1 + 1.0D+00 / k
           cr1 = -0.25D+00 * cr1 / ( k * ( k + 1 ) ) * z2
           cs1 = cs1 + cr1 * ( 2.0D+00 * w1 + 1.0D+00 / ( k + 1.0D+00 ) )
        end do
        cyv1 = rp2 * ( cec * cjv1 - 1.0D+00 / z1 - 0.25D+00 * z1 * cs1 )

     end if

  end if

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

     cfac0 = exp ( pv0 * ci )
     cfac1 = exp ( pv1 * ci )

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

  end if

  cbj(0) = cjv0
  cbj(1) = cjv1

  if ( 2 <= n .and. n <= int ( 0.25D+00 * a0 ) ) then

     cf0 = cjv0
     cf1 = cjv1
     do k = 2, n
        cf = 2.0D+00 * ( k + v0 - 1.0D+00 ) / z * cf1 - cf0
        cbj(k) = cf
        cf0 = cf1
        cf1 = cf
     end do

  else if ( 2 <= n ) then

     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 ) / z * cf1 - cf2
        if ( k <= n ) then
           cbj(k) = cf
        end if
        cf2 = cf1
        cf1 = cf
     end do
     if ( abs ( cjv1 ) < abs ( cjv0 ) ) then
        cs = cjv0 / cf
     else
        cs = cjv1 / cf2
     end if

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

  end if

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

  cby(0) = cyv0
  cby(1) = cyv1
  ya0 = abs ( cyv0 )
  lb = 0
  cg0 = cyv0
  cg1 = cyv1
  do k = 2, n
     cyk = 2.0D+00 * ( v0 + k - 1.0D+00 ) / z * cg1 - cg0
     if ( abs ( cyk ) <= 1.0D+290 ) then
        yak = abs ( cyk )
        ya1 = abs ( cg0 )
        if ( yak < ya0 .and. yak < ya1 ) then
           lb = k
        end if
        cby(k) = cyk
        cg0 = cg1
        cg1 = cyk
     end if
  end do

  if ( 4 < lb .and. imag ( z ) /= 0.0D+00 ) then

     do

        if ( lb == lb0 ) then
           exit
        end if

        ch2 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
        ch1 = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
        lb0 = lb
        do k = lb, 1, -1
           ch0 = 2.0D+00 * ( k + v0 ) / z * ch1 - ch2
           ch2 = ch1
           ch1 = ch0
        end do
        cp12 = ch0
        cp22 = ch2
        ch2 = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
        ch1 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
        do k = lb, 1, -1
           ch0 = 2.0D+00 * ( k + v0 ) / z * ch1 - ch2
           ch2 = ch1
           ch1 = ch0
        end do
        cp11 = ch0
        cp21 = ch2

        if ( lb == n ) then
           cbj(lb+1) = 2.0D+00 * ( lb + v0 ) / z * cbj(lb) - cbj(lb-1)
        end if

        if ( abs ( cbj(1) ) < abs ( cbj(0) ) ) then
           cby(lb+1) = ( cbj(lb+1) * cyv0 - 2.0D+00 * cp11 / ( pi * z ) ) &
                / cbj(0)
           cby(lb) = ( cbj(lb) * cyv0 + 2.0D+00 * cp12 / ( pi * z ) ) / cbj(0)
        else
           cby(lb+1) = ( cbj(lb+1) * cyv1 - 2.0D+00 * cp21 / ( pi * z ) ) &
                / cbj(1)
           cby(lb) = ( cbj(lb) * cyv1 + 2.0D+00 * cp22 / ( pi * z ) ) / cbj(1)
        end if

        cyl2 = cby(lb+1)
        cyl1 = cby(lb)
        do k = lb - 1, 0, -1
           cylk = 2.0D+00 * ( k + v0 + 1.0D+00 ) / z * cyl1 - cyl2
           cby(k) = cylk
           cyl2 = cyl1
           cyl1 = cylk
        end do

        cyl1 = cby(lb)
        cyl2 = cby(lb+1)
        do k = lb + 1, n - 1
           cylk = 2.0D+00 * ( k + v0 ) / z * cyl2 - cyl1
           cby(k+1) = cylk
           cyl1 = cyl2
           cyl2 = cylk
        end do

        do k = 2, n
           wa = abs ( cby(k) )
           if ( wa < abs ( cby(k-1) ) ) then
              lb = k
           end if
        end do

     end do

  end if

  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 cjyva