# include <cmath>
# include <cstdlib>

using namespace std;

# include "besselzero.hpp"

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

void jn_zeros ( int n, int nt, double z[] )

//****************************************************************************80
//
//  Purpose:
//
//    jn_zeros() computes the first NT zeros of Bessel function Jn(x).
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 June 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    integer N, the order of the Bessel functions.
//
//    integer NT, the number of zeros.
//
//  Output:
//
//    real Z(NT), the first NT zeros of Jn(x).
//
{
  double *rj1;
  double *ry0;
  double *ry1;

  rj1 = new double[nt];
  ry0 = new double[nt];
  ry1 = new double[nt];

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

  delete [] rj1;
  delete [] ry0;
  delete [] ry1;

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

void yn_zeros ( int n, int nt, double z[] )

//****************************************************************************80
//
//  Purpose:
//
//    yn_zeros() computes the first NT zeros of Bessel function Yn(x).
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 June 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    integer N, the order of the Bessel functions.
//
//    integer NT, the number of zeros.
//
//  Output:
//
//    real Z(NT), the first NT zeros of Yn(x).
//
{
  double *rj0;
  double *rj1;
  double *ry1;

  rj0 = new double[nt];
  rj1 = new double[nt];
  ry1 = new double[nt];

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

  delete [] rj0;
  delete [] rj1;
  delete [] ry1;

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

void jyndd ( int n, double x, double &bjn, double &djn, double &fjn, 
  double &byn, double &dyn, double &fyn )

//****************************************************************************80
//
//  Purpose:
//
//    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:
//
//    08 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:
//
//    integer N, the order.
//
//    real X, the argument.
//
//  Output:
//
//    real BJN, DJN, FJN, BYN, DYN, FYN, the values of
//    Jn(x), Jn'(x), Jn"(x), Yn(x), Yn'(x), Yn"(x).
//
{
  double bj[102];
  double bs;
  double by[102];
  double e0;
  double ec;
  double f;
  double f0;
  double f1;
  int k;
  int m;
  int mt;
  int nt;
  double s1;
  double su;

  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 ) == 0 )
    {
      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];
  ec = 0.5772156649015329;
  e0 = 0.3183098861837907;
  s1 = 2.0 * e0 * ( log ( x / 2.0 ) + ec ) * bj[0];
  f0 = s1 - 8.0 * e0 * su / ( bs - f );
  f1 = ( bj[1] * f0 - 2.0 * e0 / x ) / bj[0];

  by[0] = f0;
  by[1] = f1;
  for ( k = 2; k <= n + 1; k++ )
  {
    f = 2.0 * ( k - 1.0 ) * f1 / x - f0;
    by[k] = f;
    f0 = f1;
    f1 = f;
  }

  byn = by[n];
  djn = - bj[n+1] + n * bj[n] / x;
  dyn = - by[n+1] + n * by[n] / x;
  fjn = ( n * n / ( x * x ) - 1.0 ) * bjn - djn / x;
  fyn = ( n * n / ( x * x ) - 1.0 ) * byn - dyn / x;

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

void jyzo ( int n, int nt, double rj0[], double rj1[], double ry0[], 
  double ry1[] )

//****************************************************************************80
//
//  Purpose:
//
//    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:
//
//    08 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:
//
//    integer N, the order of the Bessel functions.
//
//    integer NT, the number of zeros.
//
//  Output:
//
//    real RJ0(NT), RJ1(NT), RY0(NT), RY1(NT), the zeros
//    of Jn(x), Jn'(x), Yn(x), Yn'(x).
//
{
  double bjn;
  double byn;
  double djn;
  double dyn;
  double fjn;
  double fyn;
  int l;
  double n_r8;
  double x;
  double x0;

  n_r8 = ( double ) ( n );

  if ( n <= 20 )
  {
    x = 2.82141 + 1.15859 * n_r8;
  }
  else
  {
    x = n + 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, fjn, byn, dyn, fyn );
    x = x - bjn / djn;

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

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

    if ( nt <= l )
    {
      break;
    }

  }

  if ( n <= 20 )
  {
    x = 0.961587 + 1.07703 * n_r8;
  }
  else
  {
    x = n_r8 + 0.80861 * pow ( n_r8, 0.33333 )
      + 0.07249 / pow ( n_r8, 0.33333 );
  }

  if ( n == 0 )
  {
    x = 3.8317;
  }

  l = 0;

  while ( true )
  {
    x0 = x;
    jyndd ( n, x, bjn, djn, fjn, byn, dyn, fyn );
    x = x - djn / fjn;
    if ( 1.0E-09 < fabs ( x - x0 ) )
    {
      continue;
    }
    l = l + 1;
    rj1[l-1] = x;
    x = x + 3.1416 + ( 0.4955 + 0.0915 * n_r8 
      - 0.000435 * n_r8 * n_r8 ) / l;

    if ( nt <= l )
    {
      break;
    }

  }

  if ( n <= 20 )
  {
    x = 1.19477 + 1.08933 * n_r8;
  }
  else
  {
    x = n_r8 + 0.93158 * pow ( n_r8, 0.33333 )
      + 0.26035 / pow ( n_r8, 0.33333 );
  }

  l = 0;

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

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

    l = l + 1;
    ry0[l-1] = x;
    x = x + 3.1416 + ( 0.312 + 0.0852 * n_r8 
      - 0.000403 * n_r8 * n_r8 ) / l;

    if ( nt <= l )
    {
      break;
    }

  }

  if ( n <= 20 )
  {
    x = 2.67257 + 1.16099 * n_r8;
  }
  else
  {
    x = n_r8 + 1.8211 * pow ( n_r8, 0.33333 )
      + 0.94001 / pow ( n_r8, 0.33333 );
  }

  l = 0;

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

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

    l = l + 1;
    ry1[l-1] = x;
    x = x + 3.1416 + ( 0.197 + 0.0643 * n_r8 
      -0.000286 * n_r8 * n_r8 ) / l;

    if ( nt <= l )
    {
      break;
    }

  }

  return;
}