jynb Subroutine

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

  !*****************************************************************************80
  !
  !! JYNB 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:
  !
  !    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, 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 ), save, dimension ( 4 ) :: a = (/ &
       -0.7031250000000000D-01, 0.1121520996093750D+00, &
       -0.5725014209747314D+00, 0.6074042001273483D+01 /)
  real ( kind = 8 ), save, dimension ( 4 ) :: a1 = (/ &
       0.1171875000000000D+00, -0.1441955566406250D+00, &
       0.6765925884246826D+00, -0.6883914268109947D+01 /)
  real ( kind = 8 ), save, dimension ( 4 ) :: b = (/ &
       0.7324218750000000D-01, -0.2271080017089844D+00, &
       0.1727727502584457D+01, -0.2438052969955606D+02 /)
  real ( kind = 8 ), save, dimension ( 4 ) :: b1 = (/ &
       -0.1025390625000000D+00, 0.2775764465332031D+00, &
       -0.1993531733751297D+01, 0.2724882731126854D+02 /)
  real ( kind = 8 ) bj(0:n)
  real ( kind = 8 ) bj0
  real ( kind = 8 ) bj1
  real ( kind = 8 ) bjk
  real ( kind = 8 ) bs
  real ( kind = 8 ) by(0:n)
  real ( kind = 8 ) by0
  real ( kind = 8 ) by1
  real ( kind = 8 ) byk
  real ( kind = 8 ) cu
  real ( kind = 8 ) dj(0:n)
  real ( kind = 8 ) dy(0:n)
  real ( kind = 8 ) ec
  real ( kind = 8 ) f
  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 ) p0
  real ( kind = 8 ) p1
  real ( kind = 8 ) pi
  real ( kind = 8 ) q0
  real ( kind = 8 ) q1
  real ( kind = 8 ) r2p
  real ( kind = 8 ) s0
  real ( kind = 8 ) su
  real ( kind = 8 ) sv
  real ( kind = 8 ) t1
  real ( kind = 8 ) t2
  real ( kind = 8 ) x

  pi = 3.141592653589793D+00
  r2p = 0.63661977236758D+00
  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

  if ( x <= 300.0D+00 .or. int ( 0.9D+00 * x ) < n ) then

     if ( n == 0 ) then
        nm = 1
     end if

     m = msta1 ( x, 200 )

     if ( m < nm ) then
        nm = m
     else
        m = msta2 ( x, nm, 15 )
     end if

     bs = 0.0D+00
     su = 0.0D+00
     sv = 0.0D+00
     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
        if ( k == 2 * int ( k / 2 ) .and. k /= 0 ) then
           bs = bs + 2.0D+00 * f
           su = su + ( -1.0D+00 ) ** ( k / 2 ) * f / k
        else if ( 1 < k ) then
           sv = sv + ( -1.0D+00 ) ** ( k / 2 ) * k / ( k * k - 1.0D+00 ) * f
        end if
        f2 = f1
        f1 = f
     end do

     s0 = bs + f
     do k = 0, nm
        bj(k) = bj(k) / s0
     end do

     ec = log ( x / 2.0D+00 ) + 0.5772156649015329D+00
     by0 = r2p * ( ec * bj(0) - 4.0D+00 * su / s0 )
     by(0) = by0
     by1 = r2p * ( ( ec - 1.0D+00 ) * bj(1) - bj(0) / x - 4.0D+00 * sv / s0 )
     by(1) = by1

  else

     t1 = x - 0.25D+00 * pi
     p0 = 1.0D+00
     q0 = -0.125D+00 / x
     do k = 1, 4
        p0 = p0 + a(k) * x ** ( - 2 * k )
        q0 = q0 + b(k) * x ** ( - 2 * k - 1 )
     end do
     cu = sqrt ( r2p / x )
     bj0 = cu * ( p0 * cos ( t1 ) - q0 * sin ( t1 ) )
     by0 = cu * ( p0 * sin ( t1 ) + q0 * cos ( t1 ) )
     bj(0) = bj0
     by(0) = by0
     t2 = x - 0.75D+00 * pi
     p1 = 1.0D+00
     q1 = 0.375D+00 / x
     do k = 1, 4
        p1 = p1 + a1(k) * x ** ( - 2 * k )
        q1 = q1 + b1(k) * x ** ( - 2 * k - 1 )
     end do
     bj1 = cu * ( p1 * cos ( t2 ) - q1 * sin ( t2 ) )
     by1 = cu * ( p1 * sin ( t2 ) + q1 * cos ( t2 ) )
     bj(1) = bj1
     by(1) = by1
     do k = 2, nm
        bjk = 2.0D+00 * ( k - 1.0D+00 ) / x * bj1 - bj0
        bj(k) = bjk
        bj0 = bj1
        bj1 = bjk
     end do
  end if

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

  do k = 2, nm
     byk = 2.0D+00 * ( k - 1.0D+00 ) * by1 / x - by0
     by(k) = byk
     by0 = by1
     by1 = byk
  end do

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

  return
end subroutine jynb