jyna Subroutine

subroutine jyna ( n, x, nm, bj, dj, by, dy )

  !*****************************************************************************80
  !
  !! JYNA computes Bessel functions Jn(x) and Yn(x) and 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:
  !
  !    29 April 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, integer ( kind = 4 ) N, the order.
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, integer ( kind = 4 ) NM, the highest order computed.
  !
  !    Output, real ( kind = 8 ) BJ(0:N), DJ(0:N), BY(0:N), DY(0:N), the values
  !    of Jn(x), Jn'(x), Yn(x), Yn'(x).
  !
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) bj(0:n)
  real ( kind = 8 ) bj0
  real ( kind = 8 ) bj1
  real ( kind = 8 ) bjk
  real ( kind = 8 ) by(0:n)
  real ( kind = 8 ) by0
  real ( kind = 8 ) by1
  real ( kind = 8 ) cs
  real ( kind = 8 ) dj(0:n)
  real ( kind = 8 ) dj0
  real ( kind = 8 ) dj1
  real ( kind = 8 ) dy(0:n)
  real ( kind = 8 ) dy0
  real ( kind = 8 ) dy1
  real ( kind = 8 ) f
  real ( kind = 8 ) f0
  real ( kind = 8 ) f1
  real ( kind = 8 ) f2
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
    ! integer ( kind = 4 ) msta1
  ! integer ( kind = 4 ) msta2
  integer ( kind = 4 ) nm
  real ( kind = 8 ) x

  nm = 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
     bj(0) = 1.0D+00
     dj(1) = 0.5D+00
     return

  end if

  call jy01b ( x, bj0, dj0, bj1, dj1, by0, dy0, by1, dy1 )
  bj(0) = bj0
  bj(1) = bj1
  by(0) = by0
  by(1) = by1
  dj(0) = dj0
  dj(1) = dj1
  dy(0) = dy0
  dy(1) = dy1

  if ( n <= 1 ) then
     return
  end if

  if ( n < int ( 0.9D+00 * x) ) then

     do k = 2, n
        bjk = 2.0D+00 * ( k - 1.0D+00 ) / x * bj1 - bj0
        bj(k) = bjk
        bj0 = bj1
        bj1 = bjk
     end do

  else

     m = msta1 ( x, 200 )

     if ( m < n ) then
        nm = 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 * ( k + 1.0D+00 ) / x * f1 - f2
        if ( k <= nm ) then
           bj(k) = f
        end if
        f2 = f1
        f1 = f
     end do

     if ( abs ( bj1 ) < abs ( bj0 ) ) then
        cs = bj0 / f
     else
        cs = bj1 / f2
     end if

     do k = 0, nm
        bj(k) = cs * bj(k)
     end do

  end if

  do k = 2, nm
     dj(k) = bj(k-1) - k / x * bj(k)
  end do

  f0 = by(0)
  f1 = by(1)
  do k = 2, nm
     f = 2.0D+00 * ( k - 1.0D+00 ) / x * f1 - f0
     by(k) = f
     f0 = f1
     f1 = f
  end do

  do k = 2, nm
     dy(k) = by(k-1) - k * by(k) / x
  end do

  return
end subroutine jyna