subroutine jn_zeros ( n, nt, z )

!*****************************************************************************80
!
!! jn_zeros() computes the first NT zeros of Bessel function Jn(x).
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    07 June 2025
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the order of the Bessel functions.
!
!    integer NT, the number of zeros.
!
!  Output:
!
!    real ( kind = rk8 ) Z(NT), the first NT zeros of Jn(x).
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer nt

  integer n
  real ( kind = rk8 ) rj1(nt)
  real ( kind = rk8 ) ry0(nt)
  real ( kind = rk8 ) ry1(nt)
  real ( kind = rk8 ) z(nt)

  call jyzo ( n, nt, z, rj1, ry0, ry1 )

  return
end
subroutine yn_zeros ( n, nt, z )

!*****************************************************************************80
!
!! yn_zeros() computes the first NT zeros of Bessel function Yn(x).
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    07 June 2025
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the order of the Bessel functions.
!
!    integer NT, the number of zeros.
!
!  Output:
!
!    real ( kind = rk8 ) Z(NT), the first NT zeros of Yn(x).
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer nt

  integer n
  real ( kind = rk8 ) rj0(nt)
  real ( kind = rk8 ) rj1(nt)
  real ( kind = rk8 ) ry1(nt)
  real ( kind = rk8 ) z(nt)

  call jyzo ( n, nt, rj0, rj1, z, ry1 )

  return
end
subroutine jyndd ( n, x, bjn, djn, fjn, byn, dyn, fyn )

!*****************************************************************************80
!
!! jyndd(): Bessel functions Jn(x) and Yn(x), first and second 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:
!
!    Original FORTRAN77 version by Shanjie Zhang, Jianming Jin.
!    This version by John Burkardt.
!
!  Reference:
!
!    Shanjie Zhang, Jianming Jin,
!    Computation of Special Functions,
!    Wiley, 1996,
!    ISBN: 0-471-11963-6,
!    LC: QA351.C45.
!
!  Input:
!
!    integer N, the order.
!
!    real ( kind = rk8 ) X, the argument.
!
!  Output:
!
!    real ( kind = rk8 ) BJN, DJN, FJN, BYN, DYN, FYN, the values of
!    Jn(x), Jn'(x), Jn"(x), Yn(x), Yn'(x), Yn"(x).
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ) bj(102)
  real ( kind = rk8 ) bjn
  real ( kind = rk8 ) byn
  real ( kind = rk8 ) bs
  real ( kind = rk8 ) by(102)
  real ( kind = rk8 ) djn
  real ( kind = rk8 ) dyn
  real ( kind = rk8 ) e0
  real ( kind = rk8 ) ec
  real ( kind = rk8 ) f
  real ( kind = rk8 ) f0
  real ( kind = rk8 ) f1
  real ( kind = rk8 ) fjn
  real ( kind = rk8 ) fyn
  integer k
  integer m
  integer mt
  integer n
  integer nt
  real ( kind = rk8 ) s1
  real ( kind = rk8 ) su
  real ( kind = rk8 ) x

  do nt = 1, 900
    mt = int ( 0.5D+00 * log10 ( 6.28D+00 * nt ) &
      - nt * log10 ( 1.36D+00 * abs ( x ) / nt ) )
    if ( 20 < mt ) then
      exit
    end if
  end do

  m = nt
  bs = 0.0D+00
  f0 = 0.0D+00
  f1 = 1.0D-35
  su = 0.0D+00
  do k = m, 0, -1
    f = 2.0D+00 * ( k + 1.0D+00 ) * f1 / x - f0
    if ( k <= n + 1 ) then
      bj(k+1) = f
    end if
    if ( k == 2 * int ( k / 2 ) ) then
      bs = bs + 2.0D+00 * f
      if ( k /= 0 ) then
        su = su + ( -1.0D+00 ) ** ( k / 2 ) * f / k
      end if
    end if
    f0 = f1
    f1 = f
  end do

  do k = 0, n + 1
    bj(k+1) = bj(k+1) / ( bs - f )
  end do

  bjn = bj(n+1)
  ec = 0.5772156649015329D+00
  e0 = 0.3183098861837907D+00
  s1 = 2.0D+00 * e0 * ( log ( x / 2.0D+00 ) + ec ) * bj(1)
  f0 = s1 - 8.0D+00 * e0 * su / ( bs - f )
  f1 = ( bj(2) * f0 - 2.0D+00 * e0 / x ) / bj(1)

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

  byn = by(n+1)
  djn = - bj(n+2) + n * bj(n+1) / x
  dyn = - by(n+2) + n * by(n+1) / x
  fjn = ( n * n / ( x * x ) - 1.0D+00 ) * bjn - djn / x
  fyn = ( n * n / ( x * x ) - 1.0D+00 ) * byn - dyn / x

  return
end
subroutine jyzo ( n, nt, rj0, rj1, ry0, ry1 )

!*****************************************************************************80
!
!! jyzo() computes the zeros of Bessel functions Jn(x), 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:
!
!    28 July 2012
!
!  Author:
!
!    Original FORTRAN77 version by Shanjie Zhang, Jianming Jin.
!    This version by John Burkardt.
!
!  Reference:
!
!    Shanjie Zhang, Jianming Jin,
!    Computation of Special Functions,
!    Wiley, 1996,
!    ISBN: 0-471-11963-6,
!    LC: QA351.C45.
!
!  Input:
!
!    integer N, the order of the Bessel functions.
!
!    integer NT, the number of zeros.
!
!  Output:
!
!    real ( kind = rk8 ) RJ0(NT), RJ1(NT), RY0(NT), RY1(NT), the zeros
!    of Jn(x), Jn'(x), Yn(x), Yn'(x).
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer nt

  real ( kind = rk8 ) bjn
  real ( kind = rk8 ) byn
  real ( kind = rk8 ) djn
  real ( kind = rk8 ) dyn
  real ( kind = rk8 ) fjn
  real ( kind = rk8 ) fyn
  integer l
  integer n
  real ( kind = rk8 ) n_r8
  real ( kind = rk8 ) rj0(nt)
  real ( kind = rk8 ) rj1(nt)
  real ( kind = rk8 ) ry0(nt)
  real ( kind = rk8 ) ry1(nt)
  real ( kind = rk8 ) x
  real ( kind = rk8 ) x0

  n_r8 = real ( n, kind = rk8 )

  if ( n <= 20 ) then
    x = 2.82141D+00 + 1.15859D+00 * n_r8
  else
    x = n + 1.85576D+00 * n_r8 ** 0.33333D+00 &
      + 1.03315D+00 / n_r8 ** 0.33333D+00
  end if

  l = 0

  do

    x0 = x
    call jyndd ( n, x, bjn, djn, fjn, byn, dyn, fyn )
    x = x - bjn / djn

    if ( 1.0D-09 < abs ( x - x0 ) ) then
      cycle
    end if

    l = l + 1
    rj0(l) = x
    x = x + 3.1416D+00 + ( 0.0972D+00 + 0.0679D+00 * n_r8 &
      - 0.000354D+00 * n_r8 ** 2 ) / l

    if ( nt <= l ) then
      exit
    end if

  end do

  if ( n <= 20 ) then
    x = 0.961587D+00 + 1.07703D+00 * n_r8
  else
    x = n_r8 + 0.80861D+00 * n_r8 ** 0.33333D+00 &
      + 0.07249D+00 / n_r8 ** 0.33333D+00
  end if

  if ( n == 0 ) then
    x = 3.8317D+00
  end if

  l = 0

  do

    x0 = x
    call jyndd ( n, x, bjn, djn, fjn, byn, dyn, fyn )
    x = x - djn / fjn
    if ( 1.0D-09 < abs ( x - x0 ) ) then
      cycle
    end if
    l = l + 1
    rj1(l) = x
    x = x + 3.1416D+00 + ( 0.4955D+00 + 0.0915D+00 * n_r8 &
      - 0.000435D+00 * n_r8 ** 2 ) / l

    if ( nt <= l ) then
      exit
    end if

  end do

  if ( n <= 20 ) then
    x = 1.19477D+00 + 1.08933D+00 * n_r8
  else
    x = n_r8 + 0.93158D+00 * n_r8 ** 0.33333D+00 &
      + 0.26035D+00 / n_r8 ** 0.33333D+00
  end if

  l = 0

  do

    x0 = x
    call jyndd ( n, x, bjn, djn, fjn, byn, dyn, fyn )
    x = x - byn / dyn

    if ( 1.0D-09 < abs ( x - x0 ) ) then
      cycle
    end if

    l = l + 1
    ry0(l) = x
    x = x + 3.1416D+00 + ( 0.312D+00 + 0.0852D+00 * n_r8 &
      - 0.000403D+00 * n_r8 ** 2 ) / l

    if ( nt <= l ) then
      exit
    end if

  end do

  if ( n <= 20 ) then
    x = 2.67257D+00 + 1.16099D+00 * n_r8
  else
    x = n_r8 + 1.8211D+00 * n_r8 ** 0.33333D+00 &
      + 0.94001D+00 / n_r8 ** 0.33333D+00
  end if

  l = 0

  do

    x0 = x
    call jyndd ( n, x, bjn, djn, fjn, byn, dyn, fyn )
    x = x - dyn / fyn

    if ( 1.0D-09 < abs ( x - x0 ) ) then
      cycle
    end if

    l = l + 1
    ry1(l) = x
    x = x + 3.1416D+00 + ( 0.197D+00 + 0.0643D+00 * n_r8 &
      -0.000286D+00 * n_r8 ** 2 ) / l

    if ( nt <= l ) then
      exit
    end if

  end do

  return
end