subroutine helmholtz_exact ( a, m, n, alpha, beta, gamma, nxy, x, y, Z )

!*****************************************************************************80
!
!! helmholtz_exact() evaluates an exact solution of the Helmholtz equation.
!
!  Discussion:
!
!    We consider the Helmholtz equation for Z(x,y):
!      Del^2 Z = - k^2 Z
!    where k is a parameter.  We suppose we are interested in the
!    problem when describing the vertical deflection of a vibrating circular
!    membrane, centered at the origin, of radius a.
!
!    In radial coordinates, the equation in Z(r,theta) becomes
!      Zrr + 1/r Zr + 1/r^2 Ztt + k^2 Z = 0
!    We impose the condition that Z vanishes on the circumference:
!      Z(a,theta) = 0 for all theta.
!    We try separation of variables:
!      Z(r,theta) = R(r) T(theta)
!    for which we can suppose that T is periodic
!      T'' + n^2 T = 0 for some n
!    whence
!      T(theta) = alpha * cos ( n * theta ) + beta * sin ( n * theta )
!    and
!      r^2 R'' + r R' + r^2 k^2 R - n^2 R = 0
!    for which the solution is
!      R(r) = gamma * Jn(k*r), gamma arbitrary
!
!    For each value of n, the Bessel function has infinitely many roots,
!    denoted by rho(m,n). 
!
!    To enforce the zero boundary condition at r=a, we must have that
!    the parameter k is related to some m-th zero of some n-th Bessel
!    function Jn() by:
!      k(m,n) = rho(m,n)/a
!
!    An exact solution of the Helmholtz equation, as a function of
!    r and theta, can then be written as a sum of terms of the form
!      gamma * Jn ( rho(m,n) * r / a ) 
!      * alpha * cos ( n * theta ) + beta * sin ( n * theta )
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    09 June 2025
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    This discussion relies on the Wikipedia article "Helmholtz Equation".
!
!  Input:
!
!    real a: the radius of the disk.
!
!    integer m: the index of a zero of the Jn(x) Bessel function.
!    Except when n=0, the index m=0 chooses the root at 0.
!    When n=0, there is not a root at 0, and so m=0 is illegal.
!    In any case, m=1 is the first root of Jn(x) greater than 0, 
!    2 indexes the second smallest root and so on.
!
!    integer n: the index of the J Bessel function.
!
!    real alpha, beta: the coefficients of cosine(theta) and sine(theta)
!    in the angular factor.
!
!    real gamma: the multiplier for the radial factor.
!
!    integer nxy: the number of Cartesian coordinates at which Z should
!    be sampled.
!
!    real x(nxy), y(nxy): the Cartesian coordinates at which Z should
!    be sampled.
!
!  Output:
!
!    real z(nxy): the solution value at (x,y).
!
  implicit none

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

  integer nxy

  real ( kind = rk8 ) a
  real ( kind = rk8 ) alpha
  real ( kind = rk8 ) beta
  real ( kind = rk8 ) gamma
  integer i
  real ( kind = rk8 ), allocatable :: jz(:)
  integer m
  integer n
  real ( kind = rk8 ), allocatable :: R(:)
  real ( kind = rk8 ), allocatable :: Rprime(:)
  real ( kind = rk8 ), allocatable :: rad(:)
  real ( kind = rk8 ), allocatable :: rho(:)
  real ( kind = rk8 ), allocatable :: T(:)
  real ( kind = rk8 ), allocatable :: theta(:)
  real ( kind = rk8 ) x(nxy)
  real ( kind = rk8 ) y(nxy)
  real ( kind = rk8 ) Z(nxy)
!
!  Convert Cartesian data to polar form.
!
  allocate ( rad(1:nxy) )
  allocate ( theta(1:nxy) )

  theta = atan2 ( y, x )
  rad = sqrt ( x**2 + y**2 )
!
!  Evaluate angular component T(theta).
!
  allocate ( T(1:nxy) )
  T = alpha * cos ( n * theta ) + beta * sin ( n * theta )
!
!  Evaluate the m-th zero of the n-th J bessel function.
!
!  To complicate things, Jn(x) has its first zero at zero,
!  except for n = 0.  So...an input of m=0 requests the 
!  special zero zero, which is illegal if n = 0.
!
  allocate ( rho(1:nxy) )

  if ( m == 0 ) then
    if ( n == 0 ) then
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'helmholtz_exact(): Fatal error!'
      write ( *, '(a)' ) '  m = 0 requests special first zero at 0.'
      write ( *, '(a)' ) '  This is illegal, because n = 0 as well.'
      stop ( 1 )
    end if
    rho(1:nxy) = 0.0
  else
    allocate ( jz(1:m) )
    call jn_zeros ( n, m, jz )
    rho = jz(m) * rad / a
    deallocate ( jz )
  end if
!
!  Evaluate the radial component R(r), 
!  the n-th Bessel J function evaluated at rho().
!
  allocate ( R(1:nxy) )
  allocate ( Rprime(1:nxy) )
  do i = 1, nxy
    call jyndd ( n, rho(i), R(i), Rprime(i) )
  end do
  R = gamma * R
!
!  Z is the elementwise product.
!
  Z = R * T

  deallocate ( R )
  deallocate ( Rprime )
  deallocate ( rad )
  deallocate ( rho )
  deallocate ( T )
  deallocate ( theta )

  return 
end
subroutine jyndd ( n, x, bjn, djn )

!*****************************************************************************80
!
!! jyndd() evaluates a Bessel function Jn(x) and derivative.
!
!  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, the values of Jn(x), Jn'(x).
!
  implicit none

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

  real ( kind = rk8 ) bj(102)
  real ( kind = rk8 ) bjn
  real ( kind = rk8 ) bs
  real ( kind = rk8 ) djn
  real ( kind = rk8 ) f
  real ( kind = rk8 ) f0
  real ( kind = rk8 ) f1
  integer k
  integer m
  integer mt
  integer n
  integer nt
  real ( kind = rk8 ) su
  real ( kind = rk8 ) x

  if ( x <= 0.0D+00 ) then
    bjn = 0.0D+00
    djn = 0.0D+00
    return
  end if

  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)
  djn = - bj(n+2) + n * bj(n+1) / x

  return
end
subroutine jn_zeros ( n, nt, rj0 )

!*****************************************************************************80
!
!! jyzo() computes the zeros of a Bessel function Jn(x).
!
!  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): the first NT zeros of Jn(x).
!
  implicit none

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

  integer nt

  real ( kind = rk8 ) bjn
  real ( kind = rk8 ) djn
  integer l
  integer n
  real ( kind = rk8 ) n_r8
  real ( kind = rk8 ) rj0(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_r8 + 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 )
    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

  return
end