# include <math.h>
# include <stdio.h>
# include <stdlib.h>
# include <time.h>

# include "besselj.h"

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

void bessel_jx_values ( int *n_data, double *nu, double *x, double *fx )

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

    BESSEL_JX_VALUES returns some values of the Jx Bessel function.

  Discussion:

    This set of data considers the less common case in which the
    index of the Bessel function Jn is actually not an integer.
    We may suggest this case by occasionally replacing the symbol
    "Jn" by "Jx".

    In Mathematica, the function can be evaluated by:

      BesselJ[n,x]

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    01 April 2007

  Author:

    John Burkardt

  Reference:

    Milton Abramowitz, Irene Stegun,
    Handbook of Mathematical Functions,
    National Bureau of Standards, 1964,
    ISBN: 0-486-61272-4,
    LC: QA47.A34.

    Stephen Wolfram,
    The Mathematica Book,
    Fourth Edition,
    Cambridge University Press, 1999,
    ISBN: 0-521-64314-7,
    LC: QA76.95.W65.

  Parameters:

    Input/output, int *N_DATA.  The user sets N_DATA to 0 before the
    first call.  On each call, the routine increments N_DATA by 1, and
    returns the corresponding data; when there is no more data, the
    output value of N_DATA will be 0 again.

    Output, double *NU, the order of the function.

    Output, double *X, the argument of the function.

    Output, double *FX, the value of the function.
*/
{
# define N_MAX 28

  static double fx_vec[N_MAX] = {
       0.3544507442114011E+00,
        0.6713967071418031E+00,
        0.5130161365618278E+00,
        0.3020049060623657E+00,
        0.06500818287737578E+00,
       -0.3421679847981618E+00,
       -0.1372637357550505E+00,
        0.1628807638550299E+00,
        0.2402978391234270E+00,
        0.4912937786871623E+00,
       -0.1696513061447408E+00,
        0.1979824927558931E+00,
       -0.1094768729883180E+00,
        0.04949681022847794E+00,
        0.2239245314689158E+00,
        0.2403772011113174E+00,
        0.1966584835818184E+00,
        0.02303721950962553E+00,
        0.3314145508558904E+00,
        0.5461734240402840E+00,
       -0.2616584152094124E+00,
        0.1296035513791289E+00,
       -0.1117432171933552E+00,
        0.03142623570527935E+00,
        0.1717922192746527E+00,
        0.3126634069544786E+00,
        0.1340289119304364E+00,
        0.06235967135106445E+00 };
  static double nu_vec[N_MAX] = {
    0.50E+00,
    0.50E+00,
    0.50E+00,
    0.50E+00,
    0.50E+00,
    0.50E+00,
    0.50E+00,
    0.50E+00,
    1.50E+00,
    1.50E+00,
    1.50E+00,
    1.50E+00,
    1.50E+00,
    2.50E+00,
    2.50E+00,
    2.50E+00,
    2.50E+00,
    2.50E+00,
    1.25E+00,
    1.25E+00,
    1.25E+00,
    1.25E+00,
    1.25E+00,
    2.75E+00,
    2.75E+00,
    2.75E+00,
    2.75E+00,
    2.75E+00 };

  static double x_vec[N_MAX] = {
      0.2E+00,
      1.0E+00,
      2.0E+00,
      2.5E+00,
      3.0E+00,
      5.0E+00,
     10.0E+00,
     20.0E+00,
      1.0E+00,
      2.0E+00,
      5.0E+00,
     10.0E+00,
     50.0E+00,
      1.0E+00,
      2.0E+00,
      5.0E+00,
     10.0E+00,
     50.0E+00,
      1.0E+00,
      2.0E+00,
      5.0E+00,
     10.0E+00,
     50.0E+00,
      1.0E+00,
      2.0E+00,
      5.0E+00,
     10.0E+00,
     50.0E+00 };

  if ( *n_data < 0 )
  {
    *n_data = 0;
  }

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *nu = 0.0;
    *x = 0.0;
    *fx = 0.0;
  }
  else
  {
    *nu = nu_vec[*n_data-1];
    *x = x_vec[*n_data-1];
    *fx = fx_vec[*n_data-1];
  }

  return;
# undef N_MAX
}
/******************************************************************************/

void rjbesl ( double x, double alpha, int nb, double b[], int *ncalc )

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

    RJBESL evaluates a sequence of Bessel J functions.

  Discussion:

    This routine calculates Bessel functions J sub(N+ALPHA) (X)
    for non-negative argument X, and non-negative order N+ALPHA.

    This program is based on a program written by David Sookne
    that computes values of the Bessel functions J or I of real
    argument and integer order.  Modifications include the restriction
    of the computation to the J Bessel function of non-negative real
    argument, the extension of the computation to arbitrary positive
    order, and the elimination of most underflow.

  Modified:

    15 January 2016

  Author:

    Original FORTRAN77 version by William Cody.
    C version by John Burkardt.

  Reference: 

    F Olver, David Sookne,
    A Note on Backward Recurrence Algorithms," 
    Math. Comp.,
    Volume 26, 1972, pages 941-947.

    David Sookne,
    Bessel Functions of Real Argument and Integer Order,
    NBS Journal of Res. B,
    Volume 77B, 1973, pages 125-132.

  Parameters:

    Input, double X, non-negative real argument for which
    J's are to be calculated.

    Input, double ALPHA, fractional part of order for which
    J's or exponentially scaled J'r (J*exp(X)) are
    to be calculated.  0 <= ALPHA < 1.0.

    Input, int NB, number of functions to be calculated, 
    NB > 0.  The first function calculated is of order ALPHA, and the
    last is of order (NB - 1 + ALPHA).

    Output, double B(NB).  If RJBESL
    terminates normally (NCALC=NB), the vector B contains the
    functions J/ALPHA/(X) through J/NB-1+ALPHA/(X), or the
    corresponding exponentially scaled functions.

    Output, int *NCALC, indicates possible errors.
    Before using the vector B, the user should check that
    NCALC=NB, i.e., all orders have been calculated to
    the desired accuracy.  See Error Returns below.

  Internal Parameters:

    IT = Number of bits in the mantissa of a working precision variable

    NSIG   = Decimal significance desired.  Should be set to
    INT(LOG10(2)*it+1).  Setting NSIG lower will result
    in decreased accuracy while setting NSIG higher will
    increase CPU time without increasing accuracy.  The
    truncation error is limited to a relative error of
    T=.5*10**(-NSIG).

    Then the following machine-dependent constants must be declared
    in DATA statements.  IEEE values are provided as a default.

    ENTEN  = 10.0 ** K, where K is the largest integer such that
    ENTEN is machine-representable in working precision.

    ENSIG  = 10.0 ** NSIG

    RTNSIG = 10.0 ** (-K) for the smallest integer K such that K >= NSIG/4

    ENMTEN = Smallest ABS(X) such that X/4 does not underflow

    XLARGE = Upper limit on the magnitude of X.  If ABS(X)=N,
    then at least N iterations of the backward recursion
    will be executed.  The value of 10.0 ** 4 is used on
    every machine.

  Error returns:

    In case of an error,  NCALC != NB, and not all J's are
    calculated to the desired accuracy.

    NCALC < 0:  An argument is out of range. For example,
    NBES <= 0, ALPHA < 0 or > 1, or X is too large.
    In this case, B(1) is set to zero, the remainder of the
    B-vector is not calculated, and NCALC is set to
    MIN(NB,0)-1 so that NCALC != NB.

    NB > NCALC > 0: Not all requested function values could
    be calculated accurately.  This usually occurs because NB is
    much larger than ABS(X).  In this case, B(N) is calculated
    to the desired accuracy for N <= NCALC, but precision
    is lost for NCALC < N <= NB.  If B(N) does not vanish
    for N > NCALC (because it is too small to be represented),
    and B(N)/B(NCALC) = 10**(-K), then only the first NSIG-K
    significant figures of B(N) can be trusted.
*/
{
  double alpem;
  double alp2em;
  double capp;
  double capq;
  const double eighth = 0.125E+00;
  double em;
  double en;
  const double enmten = 8.90E-308;
  const double ensig = 1.0E+16;
  const double enten = 1.0E+308;
  const double fact[25] = {
    1.0E+00, 
    1.0E+00, 
    2.0E+00, 
    6.0E+00, 
    24.0E+00, 
    1.2E+02, 
    7.2E+02, 
    5.04E+03, 
    4.032E+04, 
    3.6288E+05, 
    3.6288E+06, 
    3.99168E+07, 
    4.790016E+08, 
    6.2270208E+09, 
    8.71782912E+10, 
    1.307674368E+12, 
    2.0922789888E+13, 
    3.55687428096E+14, 
    6.402373705728E+15, 
    1.21645100408832E+17, 
    2.43290200817664E+18, 
    5.109094217170944E+19, 
    1.1240007277776076E+21, 
    2.585201673888497664E+22, 
    6.2044840173323943936E+23 };
  const double four = 4.0E+00;
  double gnu;
  const double half = 0.5E+00;
  double halfx;
  int i;
  int j;
  int jump;
  int k;
  int l;
  int m;
  int magx;
  int n;
  int nbmx;
  int nend;
  int nstart;
  const double one = 1.0E+00;
  const double one30 = 130.0E+00;
  double p;
  const double pi2 = 0.636619772367581343075535E+00;
  double plast;
  double pold;
  double psave;
  double psavel;
  const double rtnsig = 1.0E-04;
  double s;
  double sum;
  double t;
  double t1;
  double tempa;
  double tempb;
  double tempc;
  double test;
  const double three = 3.0E+00;
  const double three5 = 35.0E+00;
  double tover;
  const double two = 2.0E+00;
  const double twofiv = 25.0E+00;
  const double twopi1 = 6.28125E+00;
  const double twopi2 = 1.935307179586476925286767E-03;
  double xc;
  double xin;
  double xk;
  const double xlarge = 1.0E+04;
  double xm;
  double vcos;
  double vsin;
  double z;
  const double zero = 0.0E+00;

  jump = 0;
/*
  Check for out of range arguments.
*/
  magx = ( int ) ( x );

  if ( 
    ( 0 < nb ) &&
    ( zero <= x ) &&
    ( x <= xlarge ) &&
    ( zero <= alpha ) &&
    ( alpha < one ) )
  {
/*
  Initialize result array to zero.
*/
    *ncalc = nb;
    for ( i = 1; i <= nb; i++ )
    {
      b[i-1] = zero;
    }
/*
  Branch to use 2-term ascending series for small X and asymptotic
  form for large X when NB is not too large.
*/
    if ( x < rtnsig )
    {
/*
  Two-term ascending series for small X.
*/
      tempa = one;
      alpem = one + alpha;

      if ( enmten < x )
      {
        halfx = half * x;
      }
      else
      {
        halfx = zero;
      }

      if ( alpha != zero )
      {
        tempa = pow ( halfx, alpha ) / ( alpha * tgamma ( alpha ) );
      }

      if ( one < ( x + one ) )
      {
        tempb = - halfx * halfx;
      }
      else
      {
        tempb = zero;
      }

      b[0] = tempa + tempa * tempb / alpem;

      if ( ( x != zero ) && ( b[0] == zero ) )
      {
        *ncalc = 0;
      }

      if ( nb != 1 )
      {
        if ( x <= zero )
        {
          for ( n = 2; n <= nb; n++ )
          {
            b[n-1] = zero;
          }
        }
/*
  Calculate higher order functions.
*/
        else
        {
          tempc = halfx;

          if ( tempb != zero )
          {
            tover = enmten / tempb;
          }
          else
          {
            tover = ( enmten + enmten ) / x;
          }

          for ( n = 2; n <= nb; n++ )
          {
            tempa = tempa / alpem;
            alpem = alpem + one;

            tempa = tempa * tempc;
            if ( tempa <= tover * alpem )
            {
              tempa = zero;
            }

            b[n-1] = tempa + tempa * tempb / alpem;

            if ( ( b[n-1] == zero ) && ( n < *ncalc ) )
            {
              *ncalc = n - 1;
            }
          }
        }
      }
    }
/*
  Asymptotic series for 21.0 < X.
*/
    else if ( ( twofiv < x ) && ( nb <= magx + 1 ) )
    {
      xc = sqrt ( pi2 / x );
      xin = pow ( eighth / x, 2 );

      if ( x < three5 )
      {
        m = 11;
      }
      else if ( x < one30 )
      {
        m = 8;
      }
      else
      {
        m = 4;
      }

      xm = four * ( double ) ( m );
/*
  Argument reduction for SIN and COS routines.
*/
      t = trunc ( x / ( twopi1 + twopi2 ) + half );
      z = ( ( x - t * twopi1 ) - t * twopi2 ) - ( alpha + half ) / pi2;
      vsin = sin ( z );
      vcos = cos ( z );
      gnu = alpha + alpha;

      for ( i = 1; i <= 2; i++ )
      {
        s = ( ( xm - one ) - gnu ) * ( ( xm - one ) + gnu ) * xin * half;
        t = ( gnu - ( xm - three ) ) * ( gnu + ( xm - three ) );
        capp = s * t / fact[2*m];
        t1 = ( gnu - ( xm + one ) ) * ( gnu + ( xm + one ) );
        capq = s * t1 / fact[2*m+1];
        xk = xm;
        k = m + m;
        t1 = t;

        for ( j = 2; j <= m; j++ )
        {
          xk = xk - four;
          s = ( ( xk - one ) - gnu ) * ( ( xk - one ) + gnu );
          t = ( gnu - ( xk - three ) ) * ( gnu + ( xk - three ) );
          capp = ( capp + one / fact[k-2] ) * s * t * xin;
          capq = ( capq + one / fact[k-1] ) * s * t1 * xin;
          k = k - 2;
          t1 = t;
        }

        capp = capp + one;
        capq = ( capq + one ) * ( gnu * gnu - one ) * ( eighth / x );
        b[i-1] = xc * ( capp * vcos - capq * vsin );

        if ( nb == 1 )
        {
          return;
        }

        t = vsin;
        vsin = - vcos;
        vcos = t;
        gnu = gnu + two;
      }
/*
  If 2 < NB, compute J(X,ORDER+I)  I = 2, NB-1
*/
      if ( 2 < nb )
      {
        gnu = alpha + alpha + two;
        for ( j = 3; j <= nb; j++ )
        {
          b[j-1] = gnu * b[j-2] / x - b[j-3];
          gnu = gnu + two;
        }
      }
    }
/*
  Use recurrence to generate results.  First initialize the calculation of P*S.
*/
    else
    {
      nbmx = nb - magx;
      n = magx + 1;
      en = ( double ) ( n + n ) + ( alpha + alpha );
      plast = one;
      p = en / x;
/*
  Calculate general significance test.
*/
      test = ensig + ensig;
/*
  Calculate P*S until N = NB-1.  Check for possible overflow.
*/
      if ( 3 <= nbmx )
      {
        tover = enten / ensig;
        nstart = magx + 2;
        nend = nb - 1;
        en = ( double ) ( nstart + nstart ) - two + ( alpha + alpha );

        for ( k = nstart; k <= nend; k++ )
        {
          n = k;
          en = en + two;
          pold = plast;
          plast = p;
          p = en * plast / x - pold;
/*
  To avoid overflow, divide P*S by TOVER.  Calculate P*S until 1 < ABS(P).
*/
          if ( tover < p )
          {
            tover = enten;
            p = p / tover;
            plast = plast / tover;
            psave = p;
            psavel = plast;
            nstart = n + 1;

            for ( ; ; )
            {
              n = n + 1;
              en = en + two;
              pold = plast;
              plast = p;
              p = en * plast / x - pold;
              if ( one < p )
              {
                break;
              }
            }

            tempb = en / x;
/*
  Calculate backward test and find NCALC, the highest N such that
  the test is passed.
*/
            test = pold * plast * ( half - half / ( tempb * tempb ) );
            test = test / ensig;
            p = plast * tover;
            n = n - 1;
            en = en - two;
            if ( nb < n )
            {
              nend = nb;
            }
            else
            {
              nend = n;
            }

            for ( l = nstart; l <= nend; l++ )
            {
              pold = psavel;
              psavel = psave;
              psave = en * psavel / x - pold;
              if ( test < psave * psavel )
              {
                *ncalc = l - 1;
                jump = 1;
                break;
              }
            }

            if ( jump )
            {
              break;
            }

            *ncalc = nend;
            jump = 1;
            break;
          }
        }

        if ( ! jump )
        {
          n = nend;
          en = ( double ) ( n + n ) + ( alpha + alpha );
/*
  Calculate special significance test for 2 < NBMX.
*/
          if ( test < sqrt ( plast * ensig ) * sqrt ( p + p ) )
          {
            test = sqrt ( plast * ensig ) * sqrt ( p + p );
          }
        }
      }
/*
  Calculate P*S until significance test passes.
*/
      if ( ! jump )
      {
        for ( ; ; )
        {
          n = n + 1;
          en = en + two;
          pold = plast;
          plast = p;
          p = en * plast / x - pold;

          if ( test <= p )
          {
            break;
          }
        }
      }
/*
  Initialize the backward recursion and the normalization sum.
*/
      n = n + 1;
      en = en + two;
      tempb = zero;
      tempa = one / p;
      m = 2 * n - 4 * ( n / 2 );
      sum = zero;
      em = ( double ) ( n / 2 );
      alpem = ( em - one ) + alpha;
      alp2em = ( em + em ) + alpha;
      if ( m != 0 )
      {
        sum = tempa * alpem * alp2em / em;
      }
      nend = n - nb;
/*
  Recur backward via difference equation, calculating (but not
  storing) B, until N = NB.
*/
      if ( 0 < nend )
      {
        for ( l = 1; l <= nend; l++ )
        {
          n = n - 1;
          en = en - two;
          tempc = tempb;
          tempb = tempa;
          tempa = ( en * tempb ) / x - tempc;
          m = 2 - m;

          if ( m != 0 )
          {
            em = em - one;
            alp2em = ( em + em ) + alpha;
            if ( n == 1 )
            {
              break;
            }
            alpem = ( em - one ) + alpha;
            if ( alpem == zero )
            {
              alpem = one;
            }
            sum = ( sum + tempa * alp2em ) * alpem / em;
          }
        }
      }
/*
  Store B[NB-1].
*/
      b[n-1] = tempa;

      if ( 0 <= nend )
      {
        if ( nb <= 1 )
        {
          alp2em = alpha;
          if ( ( alpha + one ) == one )
          {
            alp2em = one;
          }
          sum = sum + b[0] * alp2em;

          if ( ( alpha + one ) != one )
          {
            sum = sum * tgamma ( alpha ) * pow ( x * half, - alpha );
          }

          tempa = enmten;

          if ( one < sum )
          {
            tempa = tempa * sum;
          }

          for ( n = 1; n <= nb; n++ )
          {
            if ( fabs ( b[n-1] ) < tempa )
            {
              b[n-1] = zero;
            }
            b[n-1] = b[n-1] / sum;
          }

          return;
        }
/*
  Calculate and store B[NB-2].
*/
        else
        {
          n = n - 1;
          en = en - two;
          b[n-1] = ( en * tempa ) / x - tempb;

          if ( n == 1 )
          {
            em = em - one;
            alp2em = ( em + em ) + alpha;
            if ( alp2em == zero )
            {
              alp2em = one;
            }
            sum = sum + b[0] * alp2em;
/*
  Normalize.  Divide all B by sum.
*/
            if ( ( alpha + one ) != one )
            {
              sum = sum * tgamma ( alpha ) * pow ( x * half, - alpha );
            }

            tempa = enmten;

            if ( one < sum )
            {
              tempa = tempa * sum;
            }
  
            for ( n = 1; n <= nb; n++ )
            {
              if ( fabs ( b[n-1] ) < tempa )
              {
                b[n-1] = zero;
              }
              b[n-1] = b[n-1] / sum;
            }
            return;
          }

          m = 2 - m;

          if ( m != 0 )
          {
            em = em - one;
            alp2em = ( em + em ) + alpha;
            alpem = ( em - one ) + alpha;
            if ( alpem == zero )
            {
              alpem = one;
            }
            sum = ( sum + b[n-1] * alp2em ) * alpem / em;
          }
        }
      }

      nend = n - 2;
/*
  Calculate via difference equation and store B, until N = 2.
*/
      if ( nend != 0 )
      {
        for ( l = 1; l <= nend; l++ )
        {
          n = n - 1;
          en = en - two;
          b[n-1] = ( en * b[n] ) / x - b[n+1];
          m = 2 - m;

          if ( m != 0 )
          {
            em = em - one;
            alp2em = ( em + em ) + alpha;
            alpem = ( em - one ) + alpha;
            if ( alpem == zero )
            {
              alpem = one;
            }
            sum = ( sum + b[n-1] * alp2em ) * alpem / em;
          }
        }
      }
/*
  Calculate B[0].
*/
      b[0] = two * ( alpha + one ) * b[1] / x - b[2];

      em = em - one;
      alp2em = ( em + em ) + alpha;
      if ( alp2em == zero )
      {
        alp2em = one;
      }
      sum = sum + b[0] * alp2em;
/*
  Normalize.  Divide all B by sum.
*/
      if ( ( alpha + one ) != one )
      {
        sum = sum * tgamma ( alpha ) * pow ( x * half, - alpha );
      }

      tempa = enmten;

      if ( one < sum )
      {
        tempa = tempa * sum;
      }

      for ( n = 1; n <= nb; n++ )
      {
        if ( fabs ( b[n-1] ) < tempa )
        {
          b[n-1] = zero;
        }
        b[n-1] = b[n-1] / sum;
      }
    }
  }
/*
  Error return.
*/
  else
  {
    b[0] = zero;
    if ( nb < 0 )
    {
      *ncalc = nb - 1;
    }
    else
    {
      *ncalc = -1;
    }
  }

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

void timestamp ( )

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

    TIMESTAMP prints the current YMDHMS date as a time stamp.

  Example:

    17 June 2014 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 June 2014

  Author:

    John Burkardt

  Parameters:

    None
*/
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  printf ( "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}