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

%*****************************************************************************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 A
%    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:
%
%    07 July 2024
%
%  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 zt 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,1), y(nxy,1): the Cartestian coordinates at which Z should
%    be sampled.
%
%  Output:
%
%    real Z(nxy,1): the solution value at (x,y).
%

%
%  Convert Cartesian data to polar form.
%
  theta = atan2 ( y, x );
  r = hypot ( x, y );
%
%  Evaluate angular component T(theta).
%
  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.
%
  if ( m == 0 )
    if ( n == 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'helmholtz_exact(): Fatal error!\n' );
      fprintf ( 1, '  m = 0 requests special first zero at 0.\n' );
      fprintf ( 1, '  This is illegal, because you specified n = 0 as well.\n' );
      error ( 'helmholtz_exact(): Fatal error!\n' );
    end
    rho = 0.0
  else
    kind = 1;
    rho_vec = besselzero ( n, m, kind );
    rho = rho_vec(m) * r / a;
  end
%
%  Evaluate the radial component R(r).
%
  R = gamma * besselj ( n, rho );
%
%  Z is the elementwise product.
%
  Z = R .* T;

  return 
end

