# include <cmath>
# include <cstdlib>
# include <iomanip>
# include <iostream>

using namespace std;

# include "helmholtz_exact.hpp"

//****************************************************************************80

double *helmholtz_exact ( double a, int m, int n, double alpha, double beta, 
  double gamma, int nxy, double x[], double y[] )

//****************************************************************************80
//
//  Purpose:
//
//    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:
//
//    17 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).
//
{
  double bjn;
  double djn;
  int i;
  double *jz;
  double *R;
  double *rad;
  double *rho;
  double *T;
  double *theta;
  double *Z;
//
//  Convert Cartesian data to polar form.
//
  theta = new double[nxy];
  for ( i = 0; i < nxy; i++ )
  {
    theta[i] = atan2 ( y[i], x[i] );
  }

  rad = new double[nxy];
  for ( i = 0; i < nxy; i++ )
  {
    rad[i] = sqrt ( x[i] * x[i] + y[i] * y[i] );
  }
//
//  Evaluate angular component T(theta).
//
  T = new double[nxy];
  for ( i = 0; i < nxy; i++ )
  {
    T[i] = alpha * cos ( n * theta[i] ) + beta * sin ( n * theta[i] );
  }
//
//  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.
//
  rho = new double[nxy];

  if ( m == 0 )
  {
    if ( n == 0 )
    {
      cout << "\n";
      cout << "helmholtz_exact(): Fatal error!\n";
      cout << "  m = 0 requests special first zero at 0.\n";
      cout << "  This is illegal, because n = 0 as well.\n";
      exit ( 1 );
    }
    for ( i = 0; i < nxy; i++ )
    {
      rho[i] = 0.0;
    }
  }
  else
  {
    jz = jn_zeros ( n, m );
    for ( i = 0; i < nxy; i++ )
    {
      rho[i] = jz[m-1] * rad[i] / a;
    }
    delete [] jz;
  }
//
//  Evaluate the radial component R(r), 
//  the n-th Bessel J function evaluated at rho().
//
  R = new double[nxy];
  for ( i = 0; i < nxy; i++ )
  {
    jyndd ( n, rho[i], bjn, djn );
    R[i] = gamma * bjn;
  }
//
//  Z is the elementwise product.
//
  Z = new double[nxy];
  for ( i = 0; i < nxy; i++ )
  {
    Z[i] = R[i] * T[i];
  }

  delete [] R;
  delete [] rad;
  delete [] rho;
  delete [] T;
  delete [] theta;

  return Z;
}
//****************************************************************************80

void jyndd ( int n, double x, double &bjn, double &djn )

//****************************************************************************80
//
//  Purpose:
//
//    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:
//
//    13 June 2025
//
//  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:
//
//    int N, the order.
//
//    double X, the argument.
//
//  Output:
//
//    double &BJN, &DJN, the values of Jn(x), Jn'(x).
//
{
  double *bj;
  double bs;
  double f;
  double f0;
  double f1;
  int k;
  int m;
  int mt;
  int nt;
  double su;

  if ( x <= 0.0 )
  {
    bjn = 0.0;
    djn = 0.0;
    return;
  }

  bj = new double[n+2];
 
  for ( nt = 1; nt <= 900; nt++ )
  {
    mt = ( int ) ( 0.5 * log10 ( 6.28 * nt ) 
      - nt * log10 ( 1.36 * fabs ( x ) / nt ) );

    if ( 20 < mt )
    {
      break;
    }

  }

  m = nt;
  bs = 0.0;
  f0 = 0.0;
  f1 = 1.0E-35;
  su = 0.0;

  for ( k = m; 0 <= k; k-- )
  {

    f = 2.0 * ( k + 1.0 ) * f1 / x - f0;

    if ( k <= n + 1 )
    {
      bj[k] = f;
    }

    if ( k == 2 * ( int ) ( k / 2 ) )
    {
      bs = bs + 2.0 * f;
      if ( k != 0 )
      {
        su = su + pow ( -1.0, k / 2 ) * f / k;
      }
    }
    f0 = f1;
    f1 = f;
  }

  for ( k = 0; k <= n + 1; k++ )
  {
    bj[k] = bj[k] / ( bs - f );
  }

  bjn = bj[n];
  djn = - bj[n+1] + n * bj[n] / x;

  delete [] bj;

  return;
}
//****************************************************************************80

double *jn_zeros ( int n, int nt )

//****************************************************************************80
//
//  Purpose:
//
//    jn_zeros() 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:
//
//    12 June 2025
//
//  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:
//
//    int N: the order of the Bessel functions.
//
//    int NT: the number of zeros.
//
//  Output:
//
//    double RJ0(NT): the first NT zeros of Jn(x).
//
{
  double bjn;
  double djn;
  int l;
  double n_r8;
  double *rj0;
  double x;
  double x0;

  rj0 = new double[nt];

  n_r8 = ( double ) ( n );

  if ( n <= 20 )
  {
    x = 2.82141 + 1.15859 * n_r8;
  }
  else
  {
    x = n_r8 + 1.85576 * pow ( n_r8, 0.33333 ) 
      + 1.03315 / pow ( n_r8, 0.33333 );
  }

  l = 0;

  while ( true )
  {
    x0 = x;
    jyndd ( n, x, bjn, djn );
    x = x - bjn / djn;

    if ( 1.0E-09 < fabs ( x - x0 ) )
    {
      continue;
    }

    rj0[l] = x;
    l = l + 1;
    x = x + 3.1416 
      + ( 0.0972 + 0.0679 * n_r8 - 0.000354 * pow ( n_r8, 2 ) ) / l;

    if ( nt <= l )
    {
      break;
    }

  }

  return rj0;
}