jyv Subroutine

subroutine jyv ( v, x, vm, bj, dj, by, dy )

  !*****************************************************************************80
  !
  !! JYV computes Bessel functions Jv(x) and Yv(x) and their derivatives.
  !
  !  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:
  !
  !    02 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(x) and Yv(x).
  !
  !    Input, real ( kind = 8 ) X, the argument of Jv(x) and Yv(x).
  !
  !    Output, real ( kind = 8 ) VM, the highest order computed.
  !
  !    Output, real ( kind = 8 ) BJ(0:N), DJ(0:N), BY(0:N), DY(0:N),
  !    the values of Jn+v0(x), Jn+v0'(x), Yn+v0(x), Yn+v0'(x).
  !
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) a0
  real ( kind = 8 ) b
  real ( kind = 8 ) bj(0:*)
  real ( kind = 8 ) bju0
  real ( kind = 8 ) bju1
  real ( kind = 8 ) bjv0
  real ( kind = 8 ) bjv1
  real ( kind = 8 ) bjvl
  real ( kind = 8 ) by(0:*)
  real ( kind = 8 ) byv0
  real ( kind = 8 ) byv1
  real ( kind = 8 ) byvk
  real ( kind = 8 ) ck
  real ( kind = 8 ) cs
  real ( kind = 8 ) cs0
  real ( kind = 8 ) cs1
  real ( kind = 8 ) dj(0:*)
  real ( kind = 8 ) dy(0:*)
  real ( kind = 8 ) ec
  real ( kind = 8 ) el
  real ( kind = 8 ) f
  real ( kind = 8 ) f0
  real ( kind = 8 ) f1
  real ( kind = 8 ) f2
  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 ) 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 ) px
  real ( kind = 8 ) qx
  real ( kind = 8 ) r
  real ( kind = 8 ) r0
  real ( kind = 8 ) r1
  real ( kind = 8 ) rp
  real ( kind = 8 ) rp2
  real ( kind = 8 ) rq
  real ( kind = 8 ) sk
  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 ) x
  real ( kind = 8 ) x2
  real ( kind = 8 ) xk

  el = 0.5772156649015329D+00
  pi = 3.141592653589793D+00
  rp2 = 0.63661977236758D+00
  x2 = x * x
  n = int ( v )
  v0 = v - n

  if ( x < 1.0D-100 ) then

     do k = 0, n
        bj(k) = 0.0D+00
        dj(k) = 0.0D+00
        by(k) = -1.0D+300
        dy(k) = 1.0D+300
     end do

     if ( v0 == 0.0D+00 ) then
        bj(0) = 1.0D+00
        dj(1) = 0.5D+00
     else
        dj(0) = 1.0D+300
     end if
     vm = v  
     return

  end if

  if ( x <= 12.0D+00 ) then

     do l = 0, 1
        vl = v0 + l
        bjvl = 1.0D+00
        r = 1.0D+00
        do k = 1, 40
           r = -0.25D+00 * r * x2 / ( k * ( k + vl ) )
           bjvl = bjvl + r
           if ( abs ( r ) < abs ( bjvl ) * 1.0D-15 ) then
              exit
           end if
        end do

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

        if ( l == 0 ) then
           bjv0 = bjvl * a
        else
           bjv1 = bjvl * a
        end if

     end do

  else

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

     do j = 0, 1

        vv = 4.0D+00 * ( j + v0 ) * ( j + v0 )
        px = 1.0D+00
        rp = 1.0D+00
        do k = 1, k0
           rp = -0.78125D-02 * rp &
                * ( 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 ) * x2 )
           px = px + rp
        end do
        qx = 1.0D+00
        rq = 1.0D+00
        do k = 1, k0
           rq = -0.78125D-02 * rq &
                * ( 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 ) * x2 )
           qx = qx + rq
        end do
        qx = 0.125D+00 * ( vv - 1.0D+00 ) * qx / x
        xk = x - ( 0.5D+00 * ( j + v0 ) + 0.25D+00 ) * pi
        a0 = sqrt ( rp2 / x )
        ck = cos ( xk )
        sk = sin ( xk )
        if ( j == 0 ) then
           bjv0 = a0 * ( px * ck - qx * sk )
           byv0 = a0 * ( px * sk + qx * ck )
        else if ( j == 1 ) then
           bjv1 = a0 * ( px * ck - qx * sk )
           byv1 = a0 * ( px * sk + qx * ck )
        end if

     end do

  end if

  bj(0) = bjv0
  bj(1) = bjv1
  dj(0) = v0 / x * bj(0) - bj(1)
  dj(1) = - ( 1.0D+00 + v0 ) / x * bj(1) + bj(0)

  if ( 2 <= n .and. n <= int ( 0.9D+00 * x ) ) then
     f0 = bjv0
     f1 = bjv1
     do k = 2, n
        f = 2.0D+00 * ( k + v0 - 1.0D+00 ) / x * f1 - f0
        bj(k) = f
        f0 = f1
        f1 = f
     end do
  else if ( 2 <= n ) then
     m = msta1 ( x, 200 )
     if ( m < n ) then
        n = m
     else
        m = msta2 ( x, n, 15 )
     end if
     f2 = 0.0D+00
     f1 = 1.0D-100
     do k = m, 0, -1
        f = 2.0D+00 * ( v0 + k + 1.0D+00 ) / x * f1 - f2
        if ( k <= n ) then
           bj(k) = f
        end if
        f2 = f1
        f1 = f
     end do

     if ( abs ( bjv1 ) < abs ( bjv0 ) ) then
        cs = bjv0 / f
     else
        cs = bjv1 / f2
     end if
     do k = 0, n
        bj(k) = cs * bj(k)
     end do
  end if

  do k = 2, n
     dj(k) = - ( k + v0 ) / x * bj(k) + bj(k-1)
  end do

  if ( x <= 12.0D+00 ) then

     if ( v0 /= 0.0D+00 ) then

        do l = 0, 1

           vl = v0 + l
           bjvl = 1.0D+00
           r = 1.0D+00
           do k = 1, 40
              r = -0.25D+00 * r * x2 / ( k * ( k - vl ) )
              bjvl = bjvl + r
              if ( abs ( r ) < abs ( bjvl ) * 1.0D-15 ) then
                 exit
              end if
           end do

           vg = 1.0D+00 - vl
           call gammaf ( vg, gb )
           b = ( 2.0D+00 / x ) ** vl / gb

           if ( l == 0 ) then
              bju0 = bjvl * b
           else
              bju1 = bjvl * b
           end if

        end do

        pv0 = pi * v0
        pv1 = pi * ( 1.0D+00 + v0 )
        byv0 = ( bjv0 * cos ( pv0 ) - bju0 ) / sin ( pv0 )
        byv1 = ( bjv1 * cos ( pv1 ) - bju1 ) / sin ( pv1 )

     else

        ec = log ( x / 2.0D+00 ) + el
        cs0 = 0.0D+00
        w0 = 0.0D+00
        r0 = 1.0D+00
        do k = 1, 30
           w0 = w0 + 1.0D+00 / k
           r0 = -0.25D+00 * r0 / ( k * k ) * x2
           cs0 = cs0 + r0 * w0
        end do
        byv0 = rp2 * ( ec * bjv0 - cs0 )
        cs1 = 1.0D+00
        w1 = 0.0D+00
        r1 = 1.0D+00
        do k = 1, 30
           w1 = w1 + 1.0D+00 / k
           r1 = -0.25D+00 * r1 / ( k * ( k + 1 ) ) * x2
           cs1 = cs1 + r1 * ( 2.0D+00 * w1 + 1.0D+00 / ( k + 1.0D+00 ) )
        end do
        byv1 = rp2 * ( ec * bjv1 - 1.0D+00 / x - 0.25D+00 * x * cs1 )

     end if

  end if

  by(0) = byv0
  by(1) = byv1
  do k = 2, n
     byvk = 2.0D+00 * ( v0 + k - 1.0D+00 ) / x * byv1 - byv0
     by(k) = byvk
     byv0 = byv1
     byv1 = byvk
  end do

  dy(0) = v0 / x * by(0) - by(1)
  do k = 1, n
     dy(k) = - ( k + v0 ) / x * by(k) + by(k-1)
  end do

  vm = n + v0

  return
end subroutine jyv