# include <math.h>
# include <stdbool.h>
# include <stdlib.h>

# include "besselj_zero.h"

/******************************************************************************/

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

/******************************************************************************/
/*
  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:

    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.

    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;

  bj = ( double * ) malloc ( ( n + 2 ) * sizeof ( double ) );
 
  for ( nt = 1; nt <= 900; nt++ )
  {
    mt = ( int ) ( 0.5 * log10 ( 6.28 * nt ) 
      - nt * log10 ( 1.36 * abs ( 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;

  free ( bj );

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

double *jn_zeros ( int n, int nt )

/******************************************************************************/
/*
  Purpose:

    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:

    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 = ( double * ) malloc ( nt * sizeof ( double ) );

  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;
}