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

using namespace std;

# include "test_values.H"

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

void abram0_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    ABRAM0_VALUES returns some values of the Abramowitz0 function.
//
//  Discussion:
//
//    The function is defined by:
//
//      ABRAM0(X) = integral ( 0 <= T < infinity ) exp ( -T * T - X / T ) dT
//
//    The data was reported by McLeod.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 August 2004
//
//  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.
//
//    Allan McLeod,
//    Algorithm 757:
//    MISCFUN: A software package to compute uncommon special functions,
//    ACM Transactions on Mathematical Software,
//    Volume 22, Number 3, September 1996, pages 288-301.
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 20

  double fx_vec[N_MAX] = { 
     0.87377726306985360531E+00,  
     0.84721859650456925922E+00,  
     0.77288934483988301615E+00,  
     0.59684345853450151603E+00,  
     0.29871735283675888392E+00,  
     0.15004596450516388138E+00,  
     0.11114662419157955096E+00,  
     0.83909567153151897766E-01,  
     0.56552321717943417515E-01,  
     0.49876496603033790206E-01,  
     0.44100889219762791328E-01,  
     0.19738535180254062496E-01,  
     0.86193088287161479900E-02,  
     0.40224788162540127227E-02,  
     0.19718658458164884826E-02,  
     0.10045868340133538505E-02,  
     0.15726917263304498649E-03,  
     0.10352666912350263437E-04,  
     0.91229759190956745069E-06,  
     0.25628287737952698742E-09 };
  double x_vec[N_MAX] = { 
     0.0019531250E+00,  
     0.0078125000E+00,  
     0.0312500000E+00,  
     0.1250000000E+00,  
     0.5000000000E+00,  
     1.0000000000E+00,  
     1.2500000000E+00,  
     1.5000000000E+00,  
     1.8750000000E+00,  
     2.0000000000E+00,  
     2.1250000000E+00,  
     3.0000000000E+00,  
     4.0000000000E+00,  
     5.0000000000E+00,  
     6.0000000000E+00,  
     7.0000000000E+00,  
     10.0000000000E+00,  
     15.0000000000E+00,  
     20.0000000000E+00,  
     40.0000000000E+00 };

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

  *n_data = *n_data + 1;

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

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

void abram1_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    ABRAM1_VALUES returns some values of the Abramowitz1 function.
//
//  Discussion:
//
//    The function is defined by:
//
//      ABRAM1(x) = integral ( 0 <= t < infinity ) t * exp ( -t^2 - x / t ) dt
//
//    The data was reported by McLeod.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 August 2004
//
//  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.
//
//    Allan McLeod,
//    Algorithm 757:
//    MISCFUN: A software package to compute uncommon special functions,
//    ACM Transactions on Mathematical Software,
//    Volume 22, Number 3, September 1996, pages 288-301.
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 20

  double fx_vec[N_MAX] = { 
     0.49828219848799921792E+00,  
     0.49324391773047288556E+00,  
     0.47431612784691234649E+00,  
     0.41095983258760410149E+00,  
     0.25317617388227035867E+00,  
     0.14656338138597777543E+00,  
     0.11421547056018366587E+00,  
     0.90026307383483764795E-01,  
     0.64088214170742303375E-01,  
     0.57446614314166191085E-01,  
     0.51581624564800730959E-01,  
     0.25263719555776416016E-01,  
     0.11930803330196594536E-01,  
     0.59270542280915272465E-02,  
     0.30609215358017829567E-02,  
     0.16307382136979552833E-02,  
     0.28371851916959455295E-03,  
     0.21122150121323238154E-04,  
     0.20344578892601627337E-05,  
     0.71116517236209642290E-09 };

  double x_vec[N_MAX] = { 
     0.0019531250E+00,  
     0.0078125000E+00,  
     0.0312500000E+00,  
     0.1250000000E+00,  
     0.5000000000E+00,  
     1.0000000000E+00,  
     1.2500000000E+00,  
     1.5000000000E+00,  
     1.8750000000E+00,  
     2.0000000000E+00,  
     2.1250000000E+00,  
     3.0000000000E+00,  
     4.0000000000E+00,  
     5.0000000000E+00,  
     6.0000000000E+00,  
     7.0000000000E+00,  
     10.0000000000E+00,  
     15.0000000000E+00,  
     20.0000000000E+00,  
     40.0000000000E+00 };

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

  *n_data = *n_data + 1;

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

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

void abram2_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    ABRAM2_VALUES returns some values of the Abramowitz2 function.
//
//  Discussion:
//
//    The function is defined by:
//
//      ABRAM2(x) = Integral ( 0 <= t < infinity ) t^2 * exp( -t^2 - x / t ) dt
// 
//    The data was reported by McLeod.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 August 2004
//
//  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.
//
//    Allan McLeod,
//    Algorithm 757:
//    MISCFUN: A software package to compute uncommon special functions,
//    ACM Transactions on Mathematical Software,
//    Volume 22, Number 3, September 1996, pages 288-301.
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 20

  double fx_vec[N_MAX] = { 
     0.44213858162107913430E+00,  
     0.43923379545684026308E+00,  
     0.42789857297092602234E+00,  
     0.38652825661854504406E+00,  
     0.26538204413231368110E+00,  
     0.16848734838334595000E+00,  
     0.13609200032513227112E+00,  
     0.11070330027727917352E+00,  
     0.82126019995530382267E-01,  
     0.74538781999594581763E-01,  
     0.67732034377612811390E-01,  
     0.35641808698811851022E-01,  
     0.17956589956618269083E-01,  
     0.94058737143575370625E-02,  
     0.50809356204299213556E-02,  
     0.28149565414209719359E-02,  
     0.53808696422559303431E-03,  
     0.44821756380146327259E-04,  
     0.46890678427324100410E-05,  
     0.20161544850996420504E-08 };

  double x_vec[N_MAX] = { 
     0.0019531250E+00,  
     0.0078125000E+00,  
     0.0312500000E+00,  
     0.1250000000E+00,  
     0.5000000000E+00,  
     1.0000000000E+00,  
     1.2500000000E+00,  
     1.5000000000E+00,  
     1.8750000000E+00,  
     2.0000000000E+00,  
     2.1250000000E+00,  
     3.0000000000E+00,  
     4.0000000000E+00,  
     5.0000000000E+00,  
     6.0000000000E+00,  
     7.0000000000E+00,  
     10.0000000000E+00,  
     15.0000000000E+00,  
     20.0000000000E+00,  
     40.0000000000E+00 };

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

  *n_data = *n_data + 1;

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

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

void agm_values ( int *n_data, double *a, double *b, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    AGM_VALUES returns some values of the AGM.
//
//  Discussion:
//
//    The AGM is defined for nonnegative A and B.
//
//    The AGM of numbers A and B is defined by setting
//
//      A(0) = A,
//      B(0) = B
//
//      A(N+1) = ( A(N) + B(N) ) / 2
//      B(N+1) = sqrt ( A(N) * B(N) )
//
//    The two sequences both converge to AGM(A,B).
//
//    In Mathematica, the AGM can be evaluated by
//
//      ArithmeticGeometricMean [ a, b ]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 February 2008
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    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 *A, *B, the argument ofs the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 15

  double a_vec[N_MAX] = {
        22.0, 
        83.0, 
        42.0, 
        26.0, 
         4.0, 
         6.0, 
        40.0,
        80.0,
        90.0,
         9.0,
        53.0,
         1.0,
         1.0,
         1.0,
         1.5 }; 
  double b_vec[N_MAX] = {
        96.0,
        56.0,
         7.0,
        11.0,
        63.0,
        45.0,
        75.0,
         0.0,
        35.0,
         1.0,
        53.0,
         2.0,
         4.0,
         8.0,
         8.0 };    
  double fx_vec[N_MAX] = {
        52.274641198704240049,
        68.836530059858524345,
        20.659301196734009322,
        17.696854873743648823,
        23.867049721753300163,
        20.717015982805991662,
        56.127842255616681863,
         0.000000000000000000,
        59.269565081229636528,
        3.9362355036495554780,
        53.000000000000000000,
        1.4567910310469068692,
        2.2430285802876025701,
        3.6157561775973627487,
        4.0816924080221632670 };
 
  if ( *n_data < 0 )
  {
    *n_data = 0;
  }

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *a = 0.0;
    *b = 0.0;
    *fx = 0.0;
  }
  else
  {
    *a = a_vec[*n_data-1];
    *b = b_vec[*n_data-1];
    *fx = fx_vec[*n_data-1];
  }

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

void airy_ai_values ( int *n_data, double *x, double *ai )

//****************************************************************************80
//
//  Purpose:
//
//    AIRY_AI_VALUES returns some values of the Airy Ai(x) function.
//
//  Discussion:
//
//    The Airy functions Ai(X) and Bi(X) are a pair of linearly independent
//    solutions of the differential equation:
//
//      W'' - X * W = 0;
//
//    In Mathematica, the function can be evaluated by:
//
//      AiryAi[x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    11 August 2004
//
//  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 *X, the argument of the function.
//
//    Output, double *AI, the value of the Airy AI function.
//
{
# define N_MAX 11

  double ai_vec[N_MAX] = { 
     0.3550280538878172E+00,  
     0.3292031299435381E+00,  
     0.3037031542863820E+00,  
     0.2788064819550049E+00,  
     0.2547423542956763E+00,  
     0.2316936064808335E+00,  
     0.2098000616663795E+00,  
     0.1891624003981501E+00,  
     0.1698463174443649E+00,  
     0.1518868036405444E+00,  
     0.1352924163128814E+00 };

  double x_vec[N_MAX] = { 
     0.0E+00,  
     0.1E+00,  
     0.2E+00,  
     0.3E+00,  
     0.4E+00,  
     0.5E+00,  
     0.6E+00,  
     0.7E+00,  
     0.8E+00,  
     0.9E+00,  
     1.0E+00 };

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

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *x = 0.0;
    *ai = 0.0;
  }
  else
  {
    *x  = x_vec[*n_data-1];
    *ai = ai_vec[*n_data-1];
  }

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

void airy_ai_int_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    AIRY_AI_INT_VALUES returns some values of the integral of the Airy function.
//
//  Discussion:
//
//    The function is defined by:
//
//      AIRY_AI_INT(x) = Integral ( 0 <= t <= x ) Ai(t) dt
//
//    The data was reported by McLeod.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 August 2004
//
//  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.
//
//    Allan McLeod,
//    Algorithm 757:
//    MISCFUN: A software package to compute uncommon special functions,
//    ACM Transactions on Mathematical Software,
//    Volume 22, Number 3, September 1996, pages 288-301.
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 20

  double fx_vec[N_MAX] = { 
     -0.75228838916610124300E+00,  
     -0.57348350185854889466E+00,  
     -0.76569840313421291743E+00,  
     -0.65181015505382467421E+00,  
     -0.55881974894471876922E+00,  
     -0.56902352870716815309E+00,  
     -0.47800749642926168100E+00,  
     -0.46567398346706861416E+00,  
     -0.96783140945618013679E-01,  
     -0.34683049857035607494E-03,  
      0.34658366917927930790E-03,  
      0.27657581846051227124E-02,  
      0.14595330491185717833E+00,  
      0.23631734191710977960E+00,  
      0.33289264538612212697E+00,  
      0.33318759129779422976E+00,  
      0.33332945170523851439E+00,  
      0.33333331724248357420E+00,  
      0.33333333329916901594E+00,  
      0.33333333333329380187E+00 };

  double x_vec[N_MAX] = { 
     -12.0000000000E+00,  
     -11.0000000000E+00,  
     -10.0000000000E+00,  
      -9.5000000000E+00,  
      -9.0000000000E+00,  
      -6.5000000000E+00,  
      -4.0000000000E+00,  
      -1.0000000000E+00,  
      -0.2500000000E+00,  
      -0.0009765625E+00,  
       0.0009765625E+00,  
       0.0078125000E+00,  
       0.5000000000E+00,  
       1.0000000000E+00,  
       4.0000000000E+00,  
       4.5000000000E+00,  
       6.0000000000E+00,  
       8.0000000000E+00,  
      10.0000000000E+00,  
      12.0000000000E+00 };         

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

  *n_data = *n_data + 1;

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

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

void airy_ai_prime_values ( int *n_data, double *x, double *aip )

//****************************************************************************80
//
//  Purpose:
//
//    AIRY_AI_PRIME_VALUES returns some values of the Airy function Ai'(x).
//
//  Discussion:
//
//    The Airy functions Ai(X) and Bi(X) are a pair of linearly independent
//    solutions of the differential equation:
//
//      W'' - X * W = 0;
//
//    In Mathematica, the function can be evaluated by:
//
//      AiryAiPrime[x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    11 August 2004
//
//  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 *X, the argument of the function.
//
//    Output, double *AIP, the derivative of the Airy AI function.
//
{
# define N_MAX 11

  double aip_vec[N_MAX] = { 
     -0.2588194037928068E+00,  
     -0.2571304219075862E+00,  
     -0.2524054702856195E+00,  
     -0.2451463642190548E+00,  
     -0.2358320344192082E+00,  
     -0.2249105326646839E+00,  
     -0.2127932593891585E+00,  
     -0.1998511915822805E+00,  
     -0.1864128638072717E+00,  
     -0.1727638434616347E+00,  
     -0.1591474412967932E+00 };

  double x_vec[N_MAX] = { 
     0.0E+00,  
     0.1E+00,  
     0.2E+00,  
     0.3E+00,  
     0.4E+00,  
     0.5E+00,  
     0.6E+00,  
     0.7E+00,  
     0.8E+00,  
     0.9E+00,  
     1.0E+00 };

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

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *x = 0.0;
    *aip = 0.0;
  }
  else
  {
    *x  = x_vec[*n_data-1];
    *aip = aip_vec[*n_data-1];
  }

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

void airy_bi_values ( int *n_data, double *x, double *bi )

//****************************************************************************80
//
//  Purpose:
//
//    AIRY_BI_VALUES returns some values of the Airy Bi(x) function.
//
//  Discussion:
//
//    The Airy functions Ai(X) and Bi(X) are a pair of linearly independent
//    solutions of the differential equation:
//
//      W'' - X * W = 0;
//
//    In Mathematica, the function can be evaluated by:
//
//      AiryBi[x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    11 August 2004
//
//  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 *X, the argument of the function.
//
//    Output, double *BI, the value of the Airy BI function.
//
{
# define N_MAX 11

  double bi_vec[N_MAX] = { 
     0.6149266274460007E+00,   
     0.6598616901941892E+00,   
     0.7054642029186612E+00,   
     0.7524855850873156E+00,   
     0.8017730000135972E+00,   
     0.8542770431031555E+00,   
     0.9110633416949405E+00,   
     0.9733286558781659E+00,   
     0.1042422171231561E+01,   
     0.1119872813134447E+01,   
     0.1207423594952871E+01 };

  double x_vec[N_MAX] = { 
     0.0E+00,  
     0.1E+00,  
     0.2E+00,  
     0.3E+00,  
     0.4E+00,  
     0.5E+00,  
     0.6E+00,  
     0.7E+00,  
     0.8E+00,  
     0.9E+00,  
     1.0E+00 };

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

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *x = 0.0;
    *bi = 0.0;
  }
  else
  {
    *x  = x_vec[*n_data-1];
    *bi = bi_vec[*n_data-1];
  }

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

void airy_bi_int_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    AIRY_BI_INT_VALUES returns some values of the integral of the Airy function.
//
//  Discussion:
//
//    The function is defined by:
//
//      AIRY_BI_INT(x) = Integral ( 0 <= t <= x ) Bi(t) dt
//
//    The data was reported by McLeod.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 August 2004
//
//  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.
//
//    Allan McLeod,
//    Algorithm 757:
//    MISCFUN: A software package to compute uncommon special functions,
//    ACM Transactions on Mathematical Software,
//    Volume 22, Number 3, September 1996, pages 288-301.
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 20

  double fx_vec[N_MAX] = { 
      0.17660819031554631869E-01,  
     -0.15040424806140020451E-01,  
      0.14756446293227661920E-01,  
     -0.11847304264848446271E+00,  
     -0.64916741266165856037E-01,  
      0.97260832464381044540E-01,  
      0.50760058495287539119E-01,  
     -0.37300500963429492179E+00,  
     -0.13962988442666578531E+00,  
     -0.12001735266723296160E-02,  
      0.12018836117890354598E-02,  
      0.36533846550952011043E+00,  
      0.87276911673800812196E+00,  
      0.48219475263803429675E+02,  
      0.44006525804904178439E+06,  
      0.17608153976228301458E+07,  
      0.73779211705220007228E+07,  
      0.14780980310740671617E+09,  
      0.97037614223613433849E+11,  
      0.11632737638809878460E+15 };

  double x_vec[N_MAX] = { 
     -12.0000000000E+00,  
     -10.0000000000E+00,  
      -8.0000000000E+00,  
      -7.5000000000E+00,  
      -7.0000000000E+00,  
      -6.5000000000E+00,  
      -4.0000000000E+00,  
      -1.0000000000E+00,  
      -0.2500000000E+00,  
      -0.0019531250E+00,  
       0.0019531250E+00,  
       0.5000000000E+00,  
       1.0000000000E+00,  
       4.0000000000E+00,  
       8.0000000000E+00,  
       8.5000000000E+00,  
       9.0000000000E+00,  
      10.0000000000E+00,  
      12.0000000000E+00,  
      14.0000000000E+00 };

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

  *n_data = *n_data + 1;

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

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

void airy_bi_prime_values ( int *n_data, double *x, double *bip )

//****************************************************************************80
//
//  Purpose:
//
//    AIRY_BI_PRIME_VALUES returns some values of the Airy function Bi'(x).
//
//  Discussion:
//
//    The Airy functions Ai(X) and Bi(X) are a pair of linearly independent
//    solutions of the differential equation:
//
//      W'' - X * W = 0;
//
//    In Mathematica, the function can be evaluated by:
//
//      AiryBiPrime[x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    11 August 2004
//
//  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 *X, the argument of the function.
//
//    Output, double *BIP, the derivative of the Airy BI function.
//
{
# define N_MAX 11

  double bip_vec[N_MAX] = { 
     0.4482883573538264E+00,   
     0.4515126311496465E+00,   
     0.4617892843621509E+00,   
     0.4800490287524480E+00,   
     0.5072816760506224E+00,   
     0.5445725641405923E+00,   
     0.5931444786342857E+00,   
     0.6544059191721400E+00,   
     0.7300069016152518E+00,   
     0.8219038903072090E+00,   
     0.9324359333927756E+00 };

  double x_vec[N_MAX] = { 
     0.0E+00,  
     0.1E+00,  
     0.2E+00,  
     0.3E+00,  
     0.4E+00,  
     0.5E+00,  
     0.6E+00,  
     0.7E+00,  
     0.8E+00,  
     0.9E+00,  
     1.0E+00 };

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

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *x = 0.0;
    *bip = 0.0;
  }
  else
  {
    *x  = x_vec[*n_data-1];
    *bip = bip_vec[*n_data-1];
  }

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

void airy_cai_values ( int *n_data, complex <double> *x, 
  complex <double> *cai )

//****************************************************************************80
//
//  Purpose:
//
//    AIRY_CAI_VALUES returns some values of the Airy Ai(x) for complex argument.
//
//  Discussion:
//
//    The Airy functions Ai(X) and Bi(X) are a pair of linearly independent
//    solutions of the differential equation:
//
//      W'' - X * W = 0;
//
//    In Mathematica, the function can be evaluated by:
//
//      AiryAi[x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    13 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, complex <double> *X, the argument of the function.
//
//    Output, complex <double> *CAI, the value of the Airy AI function.
//
{
# define N_MAX 10

  complex <double> cai_vec[N_MAX] = { 
   complex <double> ( 0.1352924163128814, + 0.0000000000000000  ), 
   complex <double> ( 0.1433824486882056, - 0.1092193342707378  ), 
   complex <double> ( 0.2215404472324631, - 0.2588711788891803  ), 
   complex <double> ( 0.4763929771766866, - 0.3036484220291284  ), 
   complex <double> ( 0.5983692170633874, - 0.08154602160771214 ),
   complex <double> ( 0.5355608832923521, + 0.00000000000000000 ),
   complex <double> ( 0.5983692170633874, + 0.08154602160771214 ),
   complex <double> ( 0.4763929771766866, + 0.3036484220291284  ), 
   complex <double> ( 0.2215404472324631, + 0.2588711788891803  ), 
   complex <double> ( 0.1433824486882056, + 0.1092193342707378  ) };

  complex <double> x_vec[N_MAX] = { 
   complex <double> (  1.0000000000000000, + 0.0000000000000000 ), 
   complex <double> (  0.8090169943749474, + 0.5877852522924731 ), 
   complex <double> (  0.3090169943749474, + 0.9510565162951536 ), 
   complex <double> ( -0.3090169943749474, + 0.9510565162951536 ), 
   complex <double> ( -0.8090169943749474, + 0.5877852522924731 ), 
   complex <double> ( -1.0000000000000000, + 0.0000000000000000 ), 
   complex <double> ( -0.8090169943749474, - 0.5877852522924731 ), 
   complex <double> ( -0.3090169943749474, - 0.9510565162951536 ), 
   complex <double> (  0.3090169943749474, - 0.9510565162951536 ), 
   complex <double> (  0.8090169943749474, - 0.5877852522924731 ) };

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

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *x = 0.0;
    *cai = 0.0;
  }
  else
  {
    *x  = x_vec[*n_data-1];
    *cai = cai_vec[*n_data-1];
  }

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

void airy_cbi_values ( int *n_data, complex <double> *x, 
  complex <double> *cbi )

//****************************************************************************80
//
//  Purpose:
//
//    AIRY_CBI_VALUES returns some values of the Airy Bi(x) for complex argument.
//
//  Discussion:
//
//    The Airy functions Ai(X) and Bi(X) are a pair of linearly independent
//    solutions of the differential equation:
//
//      W'' - X * W = 0;
//
//    In Mathematica, the function can be evaluated by:
//
//      AiryAi[x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    13 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, complex <double> *X, the argument of the function.
//
//    Output, complex <double> *CBI, the value of the Airy BI function.
//
{
# define N_MAX 10

  complex <double> cbi_vec[N_MAX] = { 
    complex <double> ( 1.207423594952871,  + 0.0000000000000000 ), 
    complex <double> ( 0.9127160108293936, + 0.3800456133135556 ), 
    complex <double> ( 0.6824453575635721, + 0.3343047153635002 ), 
    complex <double> ( 0.5726265660086474, + 0.3988641086982559 ), 
    complex <double> ( 0.2511841251049547, + 0.3401447690712719 ), 
    complex <double> ( 0.1039973894969446, + 0.0000000000000000 ), 
    complex <double> ( 0.2511841251049547, - 0.3401447690712719 ), 
    complex <double> ( 0.5726265660086474, - 0.3988641086982559 ), 
    complex <double> ( 0.6824453575635721, - 0.3343047153635002 ), 
    complex <double> ( 0.9127160108293936, - 0.3800456133135556 ) };

  complex <double> x_vec[N_MAX] = { 
    complex <double> (  1.0000000000000000, + 0.0000000000000000 ), 
    complex <double> (  0.8090169943749474, + 0.5877852522924731 ), 
    complex <double> (  0.3090169943749474, + 0.9510565162951536 ), 
    complex <double> ( -0.3090169943749474, + 0.9510565162951536 ), 
    complex <double> ( -0.8090169943749474, + 0.5877852522924731 ), 
    complex <double> ( -1.0000000000000000, + 0.0000000000000000 ), 
    complex <double> ( -0.8090169943749474, - 0.5877852522924731 ), 
    complex <double> ( -0.3090169943749474, - 0.9510565162951536 ), 
    complex <double> (  0.3090169943749474, - 0.9510565162951536 ), 
    complex <double> (  0.8090169943749474, - 0.5877852522924731 ) };

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

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *x = 0.0;
    *cbi = 0.0;
  }
  else
  {
    *x  = x_vec[*n_data-1];
    *cbi = cbi_vec[*n_data-1];
  }

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

void airy_gi_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    AIRY_GI_VALUES returns some values of the Airy Gi function.	
//
//  Discussion:
//
//    The function is defined by:
//
//      AIRY_GI(x) = Integral ( 0 <= t < infinity ) sin ( x*t+t^3/3) dt / pi
//
//    The data was reported by McLeod.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 August 2004
//
//  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.
//
//    Allan McLeod,
//    Algorithm 757:
//    MISCFUN: A software package to compute uncommon special functions,
//    ACM Transactions on Mathematical Software,
//    Volume 22, Number 3, September 1996, pages 288-301.
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 20

  double fx_vec[N_MAX] = { 
      0.20468308070040542435E+00,  
      0.18374662832557904078E+00,  
     -0.11667221729601528265E+00,  
      0.31466934902729557596E+00,  
     -0.37089040722426257729E+00,  
     -0.25293059772424019694E+00,  
      0.28967410658692701936E+00,  
     -0.34644836492634090590E+00,  
      0.28076035913873049496E+00,  
      0.21814994508094865815E+00,  
      0.20526679000810503329E+00,  
      0.22123695363784773258E+00,  
      0.23521843981043793760E+00,  
      0.82834303363768729338E-01,  
      0.45757385490989281893E-01,  
      0.44150012014605159922E-01,  
      0.39951133719508907541E-01,  
      0.35467706833949671483E-01,  
      0.31896005100679587981E-01,  
      0.26556892713512410405E-01 };

  double x_vec[N_MAX] = { 
      -0.0019531250E+00,  
      -0.1250000000E+00,  
      -1.0000000000E+00,  
      -4.0000000000E+00,  
      -8.0000000000E+00,  
      -8.2500000000E+00,  
      -9.0000000000E+00,  
     -10.0000000000E+00,  
     -11.0000000000E+00,  
     -13.0000000000E+00,  
       0.0019531250E+00,  
       0.1250000000E+00,  
       1.0000000000E+00,  
       4.0000000000E+00,  
       7.0000000000E+00,  
       7.2500000000E+00,  
       8.0000000000E+00,  
       9.0000000000E+00,  
      10.0000000000E+00,  
      12.0000000000E+00 };

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

  *n_data = *n_data + 1;

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

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

void airy_hi_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    AIRY_HI_VALUES returns some values of the Airy Hi function.
//
//  Discussion:
//
//    The function is defined by:
//
//      AIRY_HI(x) = Integral ( 0 <= t < infinity ) exp(x*t-t^3/3) dt / pi
//
//    The data was reported by McLeod.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 August 2004
//
//  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.
//
//    Allan McLeod,
//    Algorithm 757:
//    MISCFUN: A software package to compute uncommon special functions,
//    ACM Transactions on Mathematical Software,
//    Volume 22, Number 3, September 1996, pages 288-301.
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 20

  double fx_vec[N_MAX] = { 
     0.40936798278458884024E+00,  
     0.37495291608048868619E+00,  
     0.22066960679295989454E+00,  
     0.77565356679703713590E-01,  
     0.39638826473124717315E-01,  
     0.38450072575004151871E-01,  
     0.35273216868317898556E-01,  
     0.31768535282502272742E-01,  
     0.28894408288051391369E-01,  
     0.24463284011678541180E-01,  
     0.41053540139998941517E+00,  
     0.44993502381204990817E+00,  
     0.97220515514243332184E+00,  
     0.83764237105104371193E+02,  
     0.80327744952044756016E+05,  
     0.15514138847749108298E+06,  
     0.11995859641733262114E+07,  
     0.21472868855967642259E+08,  
     0.45564115351632913590E+09,  
     0.32980722582904761929E+12 };

  double x_vec[N_MAX] = { 
      -0.0019531250E+00,  
      -0.1250000000E+00,  
      -1.0000000000E+00,  
      -4.0000000000E+00,  
      -8.0000000000E+00,  
      -8.2500000000E+00,  
      -9.0000000000E+00,  
     -10.0000000000E+00,  
     -11.0000000000E+00,  
     -13.0000000000E+00,  
       0.0019531250E+00,  
       0.1250000000E+00,  
       1.0000000000E+00,  
       4.0000000000E+00,  
       7.0000000000E+00,  
       7.2500000000E+00,  
       8.0000000000E+00,  
       9.0000000000E+00,  
      10.0000000000E+00,  
      12.0000000000E+00 };

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

  *n_data = *n_data + 1;

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

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

void arccos_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    ARCCOS_VALUES returns some values of the arc cosine function.
//
//  Discussion:
//
//    In Mathematica, the function can be evaluated by:
//
//      ArcCos[x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 June 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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 12

  double fx_vec[N_MAX] = { 
    1.6709637479564564156, 
    1.5707963267948966192, 
    1.4706289056333368229, 
    1.3694384060045658278, 
    1.2661036727794991113, 
    1.1592794807274085998, 
    1.0471975511965977462, 
    0.92729521800161223243, 
    0.79539883018414355549, 
    0.64350110879328438680, 
    0.45102681179626243254, 
    0.00000000000000000000 };

  double x_vec[N_MAX] = { 
    -0.1, 
     0.0, 
     0.1, 
     0.2, 
     0.3, 
     0.4, 
     0.5, 
     0.6, 
     0.7, 
     0.8, 
     0.9, 
     1.0 };       

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

  *n_data = *n_data + 1;

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

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

void arccosh_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    ARCCOSH_VALUES returns some values of the hyperbolic arc cosine function.
//
//  Discussion:
//
//    In Mathematica, the function can be evaluated by:
//
//      ArcCosh[x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 June 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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 15

  double fx_vec[N_MAX] = { 
    0.0000000000000000000,
    0.14130376948564857735,
    0.44356825438511518913,
    0.62236250371477866781,
    0.75643291085695958624,
    0.86701472649056510395,
    0.96242365011920689500,
    1.3169578969248167086,
    1.7627471740390860505,
    1.8115262724608531070,
    2.0634370688955605467,
    2.2924316695611776878,
    2.9932228461263808979,
    5.2982923656104845907,
    7.6009022095419886114 };

  double x_vec[N_MAX] = { 
       1.0,
       1.01,
       1.1,
       1.2,
       1.3,
       1.4,
       1.5,
       2.0,
       3.0,
       3.1415926535897932385,
       4.0,
       5.0,
      10.0,
     100.0,
    1000.0 };       

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

  *n_data = *n_data + 1;

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

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

void arcsin_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    ARCSIN_VALUES returns some values of the arc sine function.
//
//  Discussion:
//
//    In Mathematica, the function can be evaluated by:
//
//      ArcSin[x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 June 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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 12

  double fx_vec[N_MAX] = { 
    -0.10016742116155979635, 
     0.00000000000000000000, 
     0.10016742116155979635, 
     0.20135792079033079146, 
     0.30469265401539750797, 
     0.41151684606748801938, 
     0.52359877559829887308, 
     0.64350110879328438680, 
     0.77539749661075306374, 
     0.92729521800161223243, 
     1.1197695149986341867, 
     1.5707963267948966192 };

  double x_vec[N_MAX] = { 
    -0.1, 
     0.0, 
     0.1, 
     0.2, 
     0.3, 
     0.4, 
     0.5, 
     0.6, 
     0.7, 
     0.8, 
     0.9, 
     1.0 };       

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

  *n_data = *n_data + 1;

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

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

void arcsinh_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    ARCSINH_VALUES returns some values of the hyperbolic arc sine function.
//
//  Discussion:
//
//    In Mathematica, the function can be evaluated by:
//
//      ArcSinh[x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 June 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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 20

  double fx_vec[N_MAX] = { 
    -2.3124383412727526203,
    -0.88137358701954302523,
     0.00000000000000000000,
     0.099834078899207563327,
     0.19869011034924140647,
     0.29567304756342243910,
     0.39003531977071527608,
     0.48121182505960344750,
     0.56882489873224753010,
     0.65266656608235578681,
     0.73266825604541086415,
     0.80886693565278246251,
     0.88137358701954302523,
     1.4436354751788103425,
     1.8184464592320668235,
     2.0947125472611012942,
     2.3124383412727526203,
     2.9982229502979697388,
     5.2983423656105887574,
     7.6009027095419886115 };

  double x_vec[N_MAX] = { 
       -5.0,
       -1.0,
        0.0,
        0.1,
        0.2,
        0.3,
        0.4,
        0.5,
        0.6,
        0.7,
        0.8,
        0.9,
        1.0,
        2.0,
        3.0,
        4.0,
        5.0,
       10.0,
      100.0,
     1000.0 };       

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

  *n_data = *n_data + 1;

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

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

void arctan_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    ARCTAN_VALUES returns some values of the arc tangent function.
//
//  Discussion:
//
//    In Mathematica, the function can be evaluated by:
//
//      ArcTan[x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 June 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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 11

  double fx_vec[N_MAX] = { 
    0.00000000000000000000, 
    0.24497866312686415417, 
    0.32175055439664219340, 
    0.46364760900080611621, 
    0.78539816339744830962, 
    1.1071487177940905030, 
    1.2490457723982544258, 
    1.3258176636680324651, 
    1.3734007669450158609, 
    1.4711276743037345919, 
    1.5208379310729538578 };

  double x_vec[N_MAX] = { 
    0.00000000000000000000, 
    0.25000000000000000000, 
    0.33333333333333333333, 
    0.50000000000000000000, 
    1.0000000000000000000, 
    2.0000000000000000000, 
    3.0000000000000000000, 
    4.0000000000000000000, 
    5.0000000000000000000, 
    10.000000000000000000, 
    20.000000000000000000 };       

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

  *n_data = *n_data + 1;

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

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

void arctan_int_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    ARCTAN_INT_VALUES returns some values of the inverse tangent integral.
//
//  Discussion:
//
//    The function is defined by:
//
//      ARCTAN_INT(x) = Integral ( 0 <= t <= x ) arctan ( t ) / t dt
//
//    The data was reported by McLeod.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Allan McLeod,
//    Algorithm 757:
//    MISCFUN: A software package to compute uncommon special functions,
//    ACM Transactions on Mathematical Software,
//    Volume 22, Number 3, September 1996, pages 288-301.
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 20

  double fx_vec[N_MAX] = { 
      0.19531241721588483191E-02,  
     -0.39062433772980711281E-02,  
      0.78124470192576499535E-02,  
      0.15624576181996527280E-01,  
     -0.31246610349485401551E-01,  
      0.62472911335014397321E-01,  
      0.12478419717389654039E+00,  
     -0.24830175098230686908E+00,  
      0.48722235829452235711E+00,  
      0.91596559417721901505E+00,  
      0.12749694484943800618E+01,  
     -0.15760154034463234224E+01,  
      0.24258878412859089996E+01,  
      0.33911633326292997361E+01,  
      0.44176450919422186583E+01,  
     -0.47556713749547247774E+01,  
      0.50961912150934111303E+01,  
      0.53759175735714876256E+01,  
     -0.61649904785027487422E+01,  
      0.72437843013083534973E+01 };

  double x_vec[N_MAX] = { 
       0.0019531250E+00,  
      -0.0039062500E+00,  
       0.0078125000E+00,  
       0.0156250000E+00,  
      -0.0312500000E+00,  
       0.0625000000E+00,  
       0.1250000000E+00,  
      -0.2500000000E+00,  
       0.5000000000E+00,  
       1.0000000000E+00,  
       1.5000000000E+00,  
      -2.0000000000E+00,  
       4.0000000000E+00,  
       8.0000000000E+00,  
      16.0000000000E+00,  
     -20.0000000000E+00,  
      25.0000000000E+00,  
      30.0000000000E+00,  
     -50.0000000000E+00,  
     100.0000000000E+00 };       

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

  *n_data = *n_data + 1;

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

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

void arctanh_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    ARCTANH_VALUES returns some values of the hyperbolic arc tangent function.
//
//  Discussion:
//
//    In Mathematica, the function can be evaluated by:
//
//      ArcTanh[x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 June 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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 15

  double fx_vec[N_MAX] = { 
    -0.54930614433405484570,
     0.00000000000000000000,
     0.0010000003333335333335,
     0.10033534773107558064,
     0.20273255405408219099,
     0.30951960420311171547,
     0.42364893019360180686,
     0.54930614433405484570,
     0.69314718055994530942,
     0.86730052769405319443,
     1.0986122886681096914,
     1.4722194895832202300,
     2.6466524123622461977,
     3.8002011672502000318,
     7.2543286192620472067 };

  double x_vec[N_MAX] = { 
    -0.5,
     0.0,
     0.001,
     0.1,
     0.2,
     0.3,
     0.4,
     0.5,
     0.6,
     0.7,
     0.8,
     0.9,
     0.99,
     0.999,
     0.999999 };       

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

  *n_data = *n_data + 1;

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

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

void bei0_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    BEI0_VALUES returns some values of the Kelvin BEI function of order NU = 0.
//
//  Discussion:
//
//    The function is defined by:
//
//      BER(NU,X) + i * BEI(NU,X) = exp(NU*Pi*I) * J(NU,X*exp(-PI*I/4))
//
//    where J(NU,X) is the J Bessel function.
//
//    In Mathematica, BEI(NU,X) can be defined by:
//
//      Im [ Exp [ NU * Pi * I ] * BesselJ [ NU, X * Exp[ -Pi * I / 4 ] ] ]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    28 June 2006
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 11

  double fx_vec[N_MAX] = {
    0.0000000000000000, 
    0.06249321838219946, 
    0.2495660400366597, 
    0.5575600623030867, 
    0.9722916273066612, 
    1.457182044159804, 
   1.937586785266043, 
   2.283249966853915, 
   2.292690322699300, 
   1.686017203632139, 
   0.1160343815502004 };
  double x_vec[N_MAX] = {
    0.0, 
    0.5, 
    1.0, 
    1.5, 
    2.0, 
    2.5, 
    3.0, 
    3.5, 
    4.0, 
    4.5, 
    5.0 };

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

  *n_data = *n_data + 1;

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

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

void bei1_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    BEI1_VALUES returns some values of the Kelvin BEI function of order NU = 1.
//
//  Discussion:
//
//    The function is defined by:
//
//      BER(NU,X) + i * BEI(NU,X) = exp(NU*Pi*I) * J(NU,X*exp(-PI*I/4))
//
//    where J(NU,X) is the J Bessel function.
//
//    In Mathematica, BEI(NU,X) can be defined by:
//
//      Im [ Exp [ NU * Pi * I ] * BesselJ [ NU, X * Exp[ -Pi * I / 4 ] ] ]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    28 June 2006
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 11

  double fx_vec[N_MAX] = {
    0.0000000000000000, 
    0.1711951797170153, 
    0.3075566313755366, 
    0.3678649890020899, 
    0.2997754370020335, 
    0.03866844396595048, 
   -0.4874541770160708, 
   -1.344042373111174, 
   -2.563821688561078, 
   -4.105685408400878, 
   -5.797907901792625 };
  double x_vec[N_MAX] = {
    0.0, 
    0.5, 
    1.0, 
    1.5, 
    2.0, 
    2.5, 
    3.0, 
    3.5, 
    4.0, 
    4.5, 
    5.0 };

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

  *n_data = *n_data + 1;

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

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

void bell_values ( int *n_data, int *n, int *c )

//****************************************************************************80
//
//  Purpose:
//
//    BELL_VALUES returns some values of the Bell numbers.
//
//  Discussion:
//
//    The Bell number B(N) is the number of restricted growth functions on N.
//
//    Note that the Stirling numbers of the second kind, S^m_n, count the
//    number of partitions of N objects into M classes, and so it is
//    true that
//
//      B(N) = S^1_N + S^2_N + ... + S^N_N.
//
//    The Bell numbers were named for Eric Temple Bell.
//
//    In Mathematica, the function can be evaluated by
//
//      Sum[StirlingS2[n,m],{m,1,n}]
//
//  Definition:
//
//    The Bell number B(N) is defined as the number of partitions (of
//    any size) of a set of N distinguishable objects.  
//
//    A partition of a set is a division of the objects of the set into 
//    subsets.
//
//  Examples:
//
//    There are 15 partitions of a set of 4 objects:
//
//      (1234), 
//      (123) (4), 
//      (124) (3), 
//      (12) (34), 
//      (12) (3) (4), 
//      (134) (2), 
//      (13) (24), 
//      (13) (2) (4), 
//      (14) (23), 
//      (1) (234),
//      (1) (23) (4), 
//      (14) (2) (3), 
//      (1) (24) (3), 
//      (1) (2) (34), 
//      (1) (2) (3) (4).
//
//    and so B(4) = 15.
//
//  First values:
//
//     N         B(N)
//     0           1
//     1           1
//     2           2
//     3           5
//     4          15
//     5          52
//     6         203
//     7         877
//     8        4140
//     9       21147
//    10      115975
//
//  Recursion:
//
//    B(I) = sum ( 1 <= J <=I ) Binomial ( I-1, J-1 ) * B(I-J)
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    06 February 2003
//
//  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, int *N, the order of the Bell number.
//
//    Output, int *C, the value of the Bell number.
//
{
# define N_MAX 11

  int c_vec[N_MAX] = { 
    1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975 };

  int n_vec[N_MAX] = { 
     0,  1,  2,  3,  4, 5,  6,  7,  8,  9,  10 };

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

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *n = 0;
    *c = 0;
  }
  else
  {
    *n = n_vec[*n_data-1];
    *c = c_vec[*n_data-1];
  }

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

void ber0_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    BER0_VALUES returns some values of the Kelvin BER function of order NU = 0.
//
//  Discussion:
//
//    The function is defined by:
//
//      BER(NU,X) + i * BEI(NU,X) = exp(NU*Pi*I) * J(NU,X*exp(-PI*I/4))
//
//    where J(NU,X) is the J Bessel function.
//
//    In Mathematica, BER(NU,X) can be defined by:
//
//      Re [ Exp [ NU * Pi * I ] * BesselJ [ NU, X * Exp[ -Pi * I / 4 ] ] ]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    27 June 2006
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 11

  double fx_vec[N_MAX] = {
    1.0000000000000000, 
    0.9990234639908383, 
    0.9843817812130869, 
    0.9210721835462558, 
    0.7517341827138082, 
    0.3999684171295313, 
   -0.2213802495986939, 
   -1.193598179589928, 
   -2.563416557258580, 
   -4.299086551599756, 
   -6.230082478666358 };
  double x_vec[N_MAX] = {
    0.0, 
    0.5, 
    1.0, 
    1.5, 
    2.0, 
    2.5, 
    3.0, 
    3.5, 
    4.0, 
    4.5, 
    5.0 };

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

  *n_data = *n_data + 1;

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

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

void ber1_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    BER1_VALUES returns some values of the Kelvin BER function of order NU = 1.
//
//  Discussion:
//
//    The function is defined by:
//
//      BER(NU,X) + i * BEI(NU,X) = exp(NU*Pi*I) * J(NU,X*exp(-PI*I/4))
//
//    where J(NU,X) is the J Bessel function.
//
//    In Mathematica, BER(NU,X) can be defined by:
//
//      Re [ Exp [ NU * Pi * I ] * BesselJ [ NU, X * Exp[ -Pi * I / 4 ] ] ]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    28 June 2006
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 11

  double fx_vec[N_MAX] = {
     0.0000000000000000, 
    -0.1822431237551121, 
    -0.3958682610197114, 
    -0.6648654179597691, 
    -0.9970776519264285, 
    -1.373096897645111, 
   -1.732644221128481, 
   -1.959644131289749, 
   -1.869248459031899, 
   -1.202821631480086, 
    0.3597766667766728 };
  double x_vec[N_MAX] = {
    0.0, 
    0.5, 
    1.0, 
    1.5, 
    2.0, 
    2.5, 
    3.0, 
    3.5, 
    4.0, 
    4.5, 
    5.0 };

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

  *n_data = *n_data + 1;

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

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

void bernoulli_number_values ( int *n_data, int *n, double *c )

//****************************************************************************80
//
//  Purpose:
//
//    BERNOULLI_NUMBER_VALUES returns some values of the Bernoulli numbers.
//
//  Discussion:
//
//    The Bernoulli numbers are rational.
//
//    If we define the sum of the M-th powers of the first N integers as:
//
//      SIGMA(M,N) = sum ( 0 <= I <= N ) I**M
//
//    and let C(I,J) be the combinatorial coefficient:
//
//      C(I,J) = I! / ( ( I - J )! * J! )
//
//    then the Bernoulli numbers B(J) satisfy:
//
//      SIGMA(M,N) = 1/(M+1) * sum ( 0 <= J <= M ) C(M+1,J) B(J) * (N+1)**(M+1-J)
//
//    In Mathematica, the function can be evaluated by:
//
//      BernoulliB[n]
//
//  First values:
//
//   B0  1                   =         1.00000000000
//   B1 -1/2                 =        -0.50000000000
//   B2  1/6                 =         1.66666666666
//   B3  0                   =         0
//   B4 -1/30                =        -0.03333333333
//   B5  0                   =         0
//   B6  1/42                =         0.02380952380
//   B7  0                   =         0
//   B8 -1/30                =        -0.03333333333
//   B9  0                   =         0
//  B10  5/66                =         0.07575757575
//  B11  0                   =         0
//  B12 -691/2730            =        -0.25311355311
//  B13  0                   =         0
//  B14  7/6                 =         1.16666666666
//  B15  0                   =         0
//  B16 -3617/510            =        -7.09215686274
//  B17  0                   =         0
//  B18  43867/798           =        54.97117794486
//  B19  0                   =         0
//  B20 -174611/330          =      -529.12424242424
//  B21  0                   =         0
//  B22  854,513/138         =      6192.123
//  B23  0                   =         0
//  B24 -236364091/2730      =    -86580.257
//  B25  0                   =         0
//  B26  8553103/6           =   1425517.16666
//  B27  0                   =         0
//  B28 -23749461029/870     = -27298231.0678
//  B29  0                   =         0
//  B30  8615841276005/14322 = 601580873.901
//
//  Recursion:
//
//    With C(N+1,K) denoting the standard binomial coefficient,
//
//    B(0) = 1.0
//    B(N) = - ( sum ( 0 <= K < N ) C(N+1,K) * B(K) ) / C(N+1,N)
//
//  Special Values:
//
//    Except for B(1), all Bernoulli numbers of odd index are 0.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 August 2004
//
//  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, int *N, the order of the Bernoulli number.
//
//    Output, double *C, the value of the Bernoulli number.
//
{
# define N_MAX 10

  double c_vec[N_MAX] = { 
      0.1000000000000000E+01,  
     -0.5000000000000000E+00,  
      0.1666666666666667E+00,  
      0.0000000000000000E+00,  
     -0.3333333333333333E-01,  
     -0.2380952380952380E-01,  
     -0.3333333333333333E-01,  
      0.7575757575757575E-01,  
     -0.5291242424242424E+03,  
      0.6015808739006424E+09 };

  int n_vec[N_MAX] = { 
     0,  1,  2,  3,  4, 6,  8, 10, 20, 30 };

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

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *n = 0;
    *c = 0.0E+00;
  }
  else
  {
    *n = n_vec[*n_data-1];
    *c = c_vec[*n_data-1];
  }

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

void bernoulli_poly_values ( int *n_data, int *n, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    BERNOULLI_POLY_VALUES returns some values of the Bernoulli polynomials.
//
//  Discussion:
//
//    In Mathematica, the function can be evaluated by:
//
//      BernoulliB[n,x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 August 2004
//
//  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, int *N, the order of the Bernoulli polynomial.
//
//    Output, double *X, the argument of the Bernoulli polynomial.
//
//    Output, double *FX, the value of the Bernoulli polynomial.
//
{
# define N_MAX 27

  double fx_vec[N_MAX] = { 
      0.1000000000000000E+01,  
     -0.3000000000000000E+00,  
      0.6666666666666667E-02,  
      0.4800000000000000E-01,  
     -0.7733333333333333E-02,  
     -0.2368000000000000E-01,  
      0.6913523809523810E-02,  
      0.2490880000000000E-01,  
     -0.1014997333333333E-01,  
     -0.4527820800000000E-01,  
      0.2332631815757576E-01,  
     -0.3125000000000000E+00,  
     -0.1142400000000000E+00,  
     -0.0176800000000000E+00,  
      0.0156800000000000E+00,  
      0.0147400000000000E+00,  
      0.0000000000000000E+00,  
     -0.1524000000000000E-01,  
     -0.2368000000000000E-01,  
     -0.2282000000000000E-01,  
     -0.1376000000000000E-01,  
      0.0000000000000000E+01,  
      0.1376000000000000E-01,  
      0.2282000000000000E-01,  
      0.2368000000000000E-01,  
      0.1524000000000000E-01,  
      0.0000000000000000E+01 };

  int n_vec[N_MAX] = { 
     0, 
     1, 
     2, 
     3, 
     4, 
     5, 
     6, 
     7, 
     8, 
     9, 
    10, 
     5, 
     5, 
     5, 
     5, 
     5, 
     5, 
     5, 
     5, 
     5, 
     5, 
     5, 
     5, 
     5, 
     5, 
     5, 
     5 };

  double x_vec[N_MAX] = { 
      0.2E+00,  
      0.2E+00,  
      0.2E+00,  
      0.2E+00,  
      0.2E+00,  
      0.2E+00,  
      0.2E+00,  
      0.2E+00,  
      0.2E+00,  
      0.2E+00,  
      0.2E+00,  
     -0.5E+00,  
     -0.4E+00,  
     -0.3E+00,  
     -0.2E+00,  
     -0.1E+00,  
      0.0E+00,  
      0.1E+00,  
      0.2E+00,  
      0.3E+00,  
      0.4E+00,  
      0.5E+00,  
      0.6E+00,  
      0.7E+00,  
      0.8E+00,  
      0.9E+00,  
      1.0E+00 };

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

  *n_data = *n_data + 1;

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

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

void bernstein_poly_values ( int *n_data, int *n, int *k, double *x, double *b )

//****************************************************************************80
//
//  Purpose:
//
//    BERNSTEIN_POLY_VALUES returns some values of the Bernstein polynomials.
//
//  Discussion:
//
//    The Bernstein polynomials are assumed to be based on [0,1].
//
//    The formula for the Bernstein polynomials is
//
//      B(N,I)(X) = [N!/(I//(N-I)!)] * (1-X)**(N-I) * X**I
//
//    In Mathematica, the function can be evaluated by:
//
//      Binomial[n,i] * (1-x)^(n-i) * x^i
//
//  First values:
//
//    B(0,0)(X) = 1
//
//    B(1,0)(X) =      1-X
//    B(1,1)(X) =                X
//
//    B(2,0)(X) =     (1-X)**2
//    B(2,1)(X) = 2 * (1-X)    * X
//    B(2,2)(X) =                X**2
//
//    B(3,0)(X) =     (1-X)**3
//    B(3,1)(X) = 3 * (1-X)**2 * X
//    B(3,2)(X) = 3 * (1-X)    * X**2
//    B(3,3)(X) =                X**3
//
//    B(4,0)(X) =     (1-X)**4
//    B(4,1)(X) = 4 * (1-X)**3 * X
//    B(4,2)(X) = 6 * (1-X)**2 * X**2
//    B(4,3)(X) = 4 * (1-X)    * X**3
//    B(4,4)(X) =                X**4
//
//  Special values:
//
//    B(N,I)(X) has a unique maximum value at X = I/N.
//
//    B(N,I)(X) has an I-fold zero at 0 and and N-I fold zero at 1.
//
//    B(N,I)(1/2) = C(N,K) / 2**N
//
//    For a fixed X and N, the polynomials add up to 1:
//
//      Sum ( 0 <= I <= N ) B(N,I)(X) = 1
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    19 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    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, int *N, the degree of the polynomial.
//
//    Output, int *K, the index of the polynomial.
//
//    Output, double *X, the argument of the polynomial.
//
//    Output, double *B, the value of the polynomial B(N,K)(X).
//
{
# define N_MAX 15

  double b_vec[N_MAX] = { 
     0.1000000000000000E+01,  
     0.7500000000000000E+00,  
     0.2500000000000000E+00,  
     0.5625000000000000E+00,  
     0.3750000000000000E+00,  
     0.6250000000000000E-01,  
     0.4218750000000000E+00,  
     0.4218750000000000E+00,  
     0.1406250000000000E+00,  
     0.1562500000000000E-01,  
     0.3164062500000000E+00,  
     0.4218750000000000E+00,  
     0.2109375000000000E+00,  
     0.4687500000000000E-01,  
     0.3906250000000000E-02 };

  int k_vec[N_MAX] = { 
    0, 
    0, 1, 
    0, 1, 2, 
    0, 1, 2, 3, 
    0, 1, 2, 3, 4 };

  int n_vec[N_MAX] = { 
    0, 
    1, 1, 
    2, 2, 2, 
    3, 3, 3, 3, 
    4, 4, 4, 4, 4 };

  double x_vec[N_MAX] = { 
     0.25E+00,  
     0.25E+00,  
     0.25E+00,  
     0.25E+00,  
     0.25E+00,  
     0.25E+00,  
     0.25E+00,  
     0.25E+00,  
     0.25E+00,  
     0.25E+00,  
     0.25E+00,  
     0.25E+00,  
     0.25E+00,  
     0.25E+00,  
     0.25E+00 };

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

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *n = 0;
    *k = 0;
    *x = 0.0;
    *b = 0.0;
  }
  else
  {
    *n = n_vec[*n_data-1];
    *k = k_vec[*n_data-1];
    *x = x_vec[*n_data-1];
    *b = b_vec[*n_data-1];
  }

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

void bessel_i0_int_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    BESSEL_I0_INT_VALUES returns some values of the Bessel I0 integral.
//
//  Discussion:
//
//    The function is defined by:
//
//      I0_INT(x) = Integral ( 0 <= t <= x ) I0(t) dt
//
//    The data was reported by McLeod.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    29 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Allan McLeod,
//    Algorithm 757:
//    MISCFUN: A software package to compute uncommon special functions,
//    ACM Transactions on Mathematical Software,
//    Volume 22, Number 3, September 1996, pages 288-301.
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 20

  double fx_vec[N_MAX] = { 
      0.19531256208818052282E-02,  
     -0.39062549670565734544E-02,  
      0.62520348032546565850E-01,  
      0.12516285581366971819E+00,  
     -0.51051480879740303760E+00,  
      0.10865210970235898158E+01,  
      0.27750019054282535299E+01,  
     -0.13775208868039716639E+02,  
      0.46424372058106108576E+03,  
      0.64111867658021584522E+07,  
     -0.10414860803175857953E+08,  
      0.44758598913855743089E+08,  
     -0.11852985311558287888E+09,  
      0.31430078220715992752E+09,  
     -0.83440212900794309620E+09,  
      0.22175367579074298261E+10,  
      0.58991731842803636487E+10,  
     -0.41857073244691522147E+11,  
      0.79553885818472357663E+12,  
      0.15089715082719201025E+17 };

  double x_vec[N_MAX] = { 
       0.0019531250E+00,  
      -0.0039062500E+00,  
       0.0625000000E+00,  
       0.1250000000E+00,  
      -0.5000000000E+00,  
       1.0000000000E+00,  
       2.0000000000E+00,  
      -4.0000000000E+00,  
       8.0000000000E+00,  
      18.0000000000E+00,  
     -18.5000000000E+00,  
      20.0000000000E+00,  
     -21.0000000000E+00,  
      22.0000000000E+00,  
     -23.0000000000E+00,  
      24.0000000000E+00,  
      25.0000000000E+00,  
     -27.0000000000E+00,  
      30.0000000000E+00,  
      40.0000000000E+00 };       

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

  *n_data = *n_data + 1;

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

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

void bessel_i0_spherical_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    BESSEL_I0_SPHERICAL_VALUES returns some values of the Spherical Bessel function i0.
//
//  Discussion:
//
//    In Mathematica, the function can be evaluated by:
//
//      Sqrt[Pi/(2*x)] * BesselI[1/2,x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    06 January 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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 21

  double fx_vec[N_MAX] = { 
    1.001667500198440E+00,
    1.006680012705470E+00,
    1.026880814507039E+00,
    1.061089303580402E+00,
    1.110132477734529E+00,
    1.175201193643801E+00,
    1.257884462843477E+00,
    1.360215358179667E+00,
    1.484729970750144E+00,
    1.634541271164267E+00,
    1.813430203923509E+00,
    2.025956895698133E+00,
    2.277595505698373E+00,
    2.574897010920645E+00,
    2.925685126512827E+00,
    3.339291642469967E+00,
    3.826838748926716E+00,
    4.401577467270101E+00,
    5.079293155726485E+00,
    5.878791279137455E+00,
    6.822479299281938E+00 };

  double x_vec[N_MAX] = { 
     0.1E+00,  
     0.2E+00,  
     0.4E+00,  
     0.6E+00,  
     0.8E+00,  
     1.0E+00,  
     1.2E+00,  
     1.4E+00,  
     1.6E+00,  
     1.8E+00,  
     2.0E+00,  
     2.2E+00,  
     2.4E+00,  
     2.6E+00,  
     2.8E+00,  
     3.0E+00,  
     3.2E+00,  
     3.4E+00,  
     3.6E+00,  
     3.8E+00,  
     4.0E+00 };

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

  *n_data = *n_data + 1;

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

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

void bessel_i0_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    BESSEL_I0_VALUES returns some values of the I0 Bessel function.
//
//  Discussion:
//
//    The modified Bessel functions In(Z) and Kn(Z) are solutions of
//    the differential equation
//
//      Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0.
//
//    The modified Bessel function I0(Z) corresponds to N = 0.
//
//    In Mathematica, the function can be evaluated by:
//
//      BesselI[0,x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    20 August 2004
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 20

  double fx_vec[N_MAX] = { 
     0.1000000000000000E+01,  
     0.1010025027795146E+01,  
     0.1040401782229341E+01,  
     0.1092045364317340E+01,  
     0.1166514922869803E+01,  
     0.1266065877752008E+01,  
     0.1393725584134064E+01,  
     0.1553395099731217E+01,  
     0.1749980639738909E+01,  
     0.1989559356618051E+01,  
     0.2279585302336067E+01,  
     0.3289839144050123E+01,  
     0.4880792585865024E+01,  
     0.7378203432225480E+01,  
     0.1130192195213633E+02,  
     0.1748117185560928E+02,  
     0.2723987182360445E+02,  
     0.6723440697647798E+02,  
     0.4275641157218048E+03,  
     0.2815716628466254E+04 };

  double x_vec[N_MAX] = { 
     0.00E+00,  
     0.20E+00,  
     0.40E+00,  
     0.60E+00,  
     0.80E+00,  
     0.10E+01,  
     0.12E+01,  
     0.14E+01,  
     0.16E+01,  
     0.18E+01,  
     0.20E+01,  
     0.25E+01,  
     0.30E+01,  
     0.35E+01,  
     0.40E+01,  
     0.45E+01,  
     0.50E+01,  
     0.60E+01,  
     0.80E+01,  
     0.10E+02 };

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

  *n_data = *n_data + 1;

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

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

void bessel_i1_spherical_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    BESSEL_I1_SPHERICAL_VALUES returns some values of the Spherical Bessel function i1.
//
//  Discussion:
//
//    In Mathematica, the function can be evaluated by:
//
//      Sqrt[Pi/(2*x)] * BesselI[3/2,x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    06 January 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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 21

  double fx_vec[N_MAX] = { 
    0.03336667857363341E+00,
    0.06693371456802954E+00,
    0.1354788933285401E+00,
    0.2072931911031093E+00,
    0.2841280857128948E+00,
    0.3678794411714423E+00,
    0.4606425870674146E+00,
    0.5647736480096238E+00,
    0.6829590627779635E+00,
    0.8182955028627777E+00,
    0.9743827435800610E+00,
    1.155432469636406E+00,
    1.366396525527973E+00,
    1.613118767572064E+00,
    1.902515460838681E+00,
    2.242790117769266E+00,
    2.643689828630357E+00,
    3.116811526884873E+00,
    3.675968313148932E+00,
    4.337627987747642E+00,
    5.121438384183637E+00 };

  double x_vec[N_MAX] = { 
     0.1E+00,  
     0.2E+00,  
     0.4E+00,  
     0.6E+00,  
     0.8E+00,  
     1.0E+00,  
     1.2E+00,  
     1.4E+00,  
     1.6E+00,  
     1.8E+00,  
     2.0E+00,  
     2.2E+00,  
     2.4E+00,  
     2.6E+00,  
     2.8E+00,  
     3.0E+00,  
     3.2E+00,  
     3.4E+00,  
     3.6E+00,  
     3.8E+00,  
     4.0E+00 };

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

  *n_data = *n_data + 1;

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

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

void bessel_i1_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    BESSEL_I1_VALUES returns some values of the I1 Bessel function.
//
//  Discussion:
//
//    The modified Bessel functions In(Z) and Kn(Z) are solutions of
//    the differential equation
//
//      Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0.
//
//    In Mathematica, the function can be evaluated by:
//
//      BesselI[1,x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    20 August 2004
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 20

  double fx_vec[N_MAX] = { 
     0.0000000000000000E+00,  
     0.1005008340281251E+00,  
     0.2040267557335706E+00,  
     0.3137040256049221E+00,  
     0.4328648026206398E+00,  
     0.5651591039924850E+00,  
     0.7146779415526431E+00,  
     0.8860919814143274E+00,  
     0.1084810635129880E+01,   
     0.1317167230391899E+01,  
     0.1590636854637329E+01,  
     0.2516716245288698E+01,  
     0.3953370217402609E+01,  
     0.6205834922258365E+01,  
     0.9759465153704450E+01,  
     0.1538922275373592E+02,  
     0.2433564214245053E+02,  
     0.6134193677764024E+02,  
     0.3998731367825601E+03,  
     0.2670988303701255E+04 };

  double x_vec[N_MAX] = { 
     0.00E+00,  
     0.20E+00,  
     0.40E+00,  
     0.60E+00,  
     0.80E+00,  
     0.10E+01,  
     0.12E+01,  
     0.14E+01,  
     0.16E+01,  
     0.18E+01,  
     0.20E+01,  
     0.25E+01,  
     0.30E+01,  
     0.35E+01,  
     0.40E+01,  
     0.45E+01,  
     0.50E+01,  
     0.60E+01,  
     0.80E+01,  
     0.10E+02 };

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

  *n_data = *n_data + 1;

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

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

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

//****************************************************************************80
//
//  Purpose:
//
//    BESSEL_IN_VALUES returns some values of the In Bessel function.
//
//  Discussion:
//
//    The modified Bessel functions In(Z) and Kn(Z) are solutions of
//    the differential equation
//
//      Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0.
//
//    In Mathematica, the function can be evaluated by:
//
//      BesselI[n,x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 August 2004
//
//  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, int *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

  double fx_vec[N_MAX] = { 
     0.5016687513894678E-02,  
     0.1357476697670383E+00,  
     0.6889484476987382E+00,  
     0.1276466147819164E+01,  
     0.2245212440929951E+01,  
     0.1750561496662424E+02,  
     0.2281518967726004E+04,  
     0.3931278522104076E+08,  
     0.2216842492433190E-01,  
     0.2127399592398527E+00,  
     0.1033115016915114E+02,  
     0.1758380716610853E+04,  
     0.2677764138883941E+21,  
     0.2714631559569719E-03,  
     0.9825679323131702E-02,  
     0.2157974547322546E+01,  
     0.7771882864032600E+03,  
     0.2278548307911282E+21,  
     0.2752948039836874E-09,  
     0.3016963879350684E-06,  
     0.4580044419176051E-02,  
     0.2189170616372337E+02,  
     0.1071597159477637E+21,  
     0.3966835985819020E-24,  
     0.4310560576109548E-18,  
     0.5024239357971806E-10,  
     0.1250799735644948E-03,  
     0.5442008402752998E+19 };

  int nu_vec[N_MAX] = { 
     2,  2,  2,  2, 
     2,  2,  2,  2, 
     3,  3,  3,  3, 
     3,  5,  5,  5, 
     5,  5, 10, 10, 
    10, 10, 10, 20, 
    20, 20, 20, 20 };

  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;
    *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
}
//****************************************************************************80

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

//****************************************************************************80
//
//  Purpose:
//
//    BESSEL_IX_VALUES returns some values of the Ix Bessel function.
//
//  Discussion:
//
//    This set of data considers the less common case in which the
//    index of the Bessel function In is actually not an integer.
//    We may suggest this case by occasionally replacing the symbol
//    "In" by "Ix".
//
//    The modified Bessel functions In(Z) and Kn(Z) are solutions of
//    the differential equation
//
//      Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0.
//
//    In Mathematica, the function can be evaluated by:
//
//      BesselI[n,x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    02 March 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

  double fx_vec[N_MAX] = {
    0.3592084175833614E+00,  
    0.9376748882454876E+00,  
    2.046236863089055E+00,   
    3.053093538196718E+00,   
    4.614822903407601E+00,   
    26.47754749755907E+00,   
    2778.784603874571E+00,   
    4.327974627242893E+07,   
    0.2935253263474798E+00,  
    1.099473188633110E+00,   
    21.18444226479414E+00,   
    2500.906154942118E+00,   
    2.866653715931464E+20,   
    0.05709890920304825E+00, 
    0.3970270801393905E+00,  
    13.76688213868258E+00,   
    2028.512757391936E+00,   
    2.753157630035402E+20,   
    0.4139416015642352E+00,  
    1.340196758982897E+00,   
    22.85715510364670E+00,   
    2593.006763432002E+00,   
    2.886630075077766E+20,   
    0.03590910483251082E+00, 
    0.2931108636266483E+00,  
    11.99397010023068E+00,   
    1894.575731562383E+00,   
    2.716911375760483E+20 };

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

  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
}
//****************************************************************************80

void bessel_j0_int_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    BESSEL_J0_INT_VALUES returns some values of the Bessel J0 integral.
//
//  Discussion:
//
//    The function is defined by:
//
//      J0_INT(x) = Integral ( 0 <= t <= x ) J0(t) dt
//
//    The data was reported by McLeod.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    29 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Allan McLeod,
//    Algorithm 757:
//    MISCFUN: A software package to compute uncommon special functions,
//    ACM Transactions on Mathematical Software,
//    Volume 22, Number 3, September 1996, pages 288-301.
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 20

  double fx_vec[N_MAX] = { 
      0.97656242238978822427E-03,  
      0.39062450329491108875E-02,  
     -0.62479657927917933620E-01,  
      0.12483733492120479139E+00,  
     -0.48968050664604505505E+00,  
      0.91973041008976023931E+00,  
     -0.14257702931970265690E+01,  
      0.10247341594606064818E+01,  
     -0.12107468348304501655E+01,  
      0.11008652032736190799E+01,  
     -0.10060334829904124192E+01,  
      0.81330572662485953519E+00,  
     -0.10583788214211277585E+01,  
      0.87101492116545875169E+00,  
     -0.88424908882547488420E+00,  
      0.11257761503599914603E+01,  
     -0.90141212258183461184E+00,  
      0.91441344369647797803E+00,  
     -0.94482281938334394886E+00,  
      0.92266255696016607257E+00 };

  double x_vec[N_MAX] = { 
       0.0009765625E+00,  
       0.0039062500E+00,  
      -0.0625000000E+00,  
       0.1250000000E+00,  
      -0.5000000000E+00,  
       1.0000000000E+00,  
      -2.0000000000E+00,  
       4.0000000000E+00,  
      -8.0000000000E+00,  
      16.0000000000E+00,  
     -16.5000000000E+00,  
      18.0000000000E+00,  
     -20.0000000000E+00,  
      25.0000000000E+00,  
     -30.0000000000E+00,  
      40.0000000000E+00,  
     -50.0000000000E+00,  
      75.0000000000E+00,  
     -80.0000000000E+00,  
     100.0000000000E+00 };       

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

  *n_data = *n_data + 1;

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

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

void bessel_j0_spherical_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    BESSEL_J0_SPHERICAL_VALUES returns some values of the Spherical Bessel function j0.
//
//  Discussion:
//
//    In Mathematica, the function can be evaluated by:
//
//      Sqrt[Pi/(2*x)] * BesselJ[1/2,x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 August 2004
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 21

  double fx_vec[N_MAX] = { 
      0.9983341664682815E+00,  
      0.9933466539753061E+00,  
      0.9735458557716262E+00,  
      0.9410707889917256E+00,  
      0.8966951136244035E+00,  
      0.8414709848078965E+00,  
      0.7766992383060220E+00,  
      0.7038926642774716E+00,  
      0.6247335019009407E+00,  
      0.5410264615989973E+00,  
      0.4546487134128408E+00,  
      0.3674983653725410E+00,  
      0.2814429918963129E+00,  
      0.1982697583928709E+00,  
      0.1196386250556803E+00,  
      0.4704000268662241E-01,  
     -0.1824191982111872E-01,  
     -0.7515914765495039E-01,  
     -0.1229223453596812E+00,  
     -0.1610152344586103E+00,  
     -0.1892006238269821E+00 };

  double x_vec[N_MAX] = { 
     0.1E+00,  
     0.2E+00,  
     0.4E+00,  
     0.6E+00,  
     0.8E+00,  
     1.0E+00,  
     1.2E+00,  
     1.4E+00,  
     1.6E+00,  
     1.8E+00,  
     2.0E+00,  
     2.2E+00,  
     2.4E+00,  
     2.6E+00,  
     2.8E+00,  
     3.0E+00,  
     3.2E+00,  
     3.4E+00,  
     3.6E+00,  
     3.8E+00,  
     4.0E+00 };

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

  *n_data = *n_data + 1;

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

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

void bessel_j0_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    BESSEL_J0_VALUES returns some values of the J0 Bessel function.
//
//  Discussion:
//
//    In Mathematica, the function can be evaluated by:
//
//      BesselJ[0,x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    10 August 2004
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 21

  double fx_vec[N_MAX] = { 
     -0.1775967713143383E+00,  
     -0.3971498098638474E+00,  
     -0.2600519549019334E+00,  
      0.2238907791412357E+00,  
      0.7651976865579666E+00,  
      0.1000000000000000E+01,  
      0.7651976865579666E+00,  
      0.2238907791412357E+00,  
     -0.2600519549019334E+00,  
     -0.3971498098638474E+00,  
     -0.1775967713143383E+00,  
      0.1506452572509969E+00,  
      0.3000792705195556E+00,  
      0.1716508071375539E+00,  
     -0.9033361118287613E-01,  
     -0.2459357644513483E+00,  
     -0.1711903004071961E+00,  
      0.4768931079683354E-01,  
      0.2069261023770678E+00,  
      0.1710734761104587E+00,  
     -0.1422447282678077E-01 };

  double x_vec[N_MAX] = { 
     -5.0E+00,  
     -4.0E+00,  
     -3.0E+00,  
     -2.0E+00,  
     -1.0E+00,  
      0.0E+00,  
      1.0E+00,  
      2.0E+00,  
      3.0E+00,  
      4.0E+00,  
      5.0E+00,  
      6.0E+00,  
      7.0E+00,  
      8.0E+00,  
      9.0E+00,  
     10.0E+00,  
     11.0E+00,  
     12.0E+00,  
     13.0E+00,  
     14.0E+00,  
     15.0E+00 };

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

  *n_data = *n_data + 1;

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

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

void bessel_j1_spherical_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    BESSEL_J1_SPHERICAL_VALUES returns some values of the Spherical Bessel function j1.
//
//  Discussion:
//
//    In Mathematica, the function can be evaluated by:
//
//      Sqrt[Pi/(2*x)] * BesselJ[3/2,x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 August 2004
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 21

  double fx_vec[N_MAX] = { 
     0.3330001190255757E-01,  
     0.6640038067032223E-01,  
     0.1312121544218529E+00,  
     0.1928919568034122E+00,  
     0.2499855053465475E+00,  
     0.3011686789397568E+00,  
     0.3452845698577903E+00,  
     0.3813753724123076E+00,  
     0.4087081401263934E+00,  
     0.4267936423844913E+00,  
     0.4353977749799916E+00,  
     0.4345452193763121E+00,  
     0.4245152947656493E+00,  
     0.4058301968314685E+00,  
     0.3792360591872637E+00,  
     0.3456774997623560E+00,  
     0.3062665174917607E+00,  
     0.2622467779189737E+00,  
     0.2149544641595738E+00,  
     0.1657769677515280E+00,  
     0.1161107492591575E+00  };

  double x_vec[N_MAX] = { 
     0.1E+00,  
     0.2E+00,  
     0.4E+00,  
     0.6E+00,  
     0.8E+00,  
     1.0E+00,  
     1.2E+00,  
     1.4E+00,  
     1.6E+00,  
     1.8E+00,  
     2.0E+00,  
     2.2E+00,  
     2.4E+00,  
     2.6E+00,  
     2.8E+00,  
     3.0E+00,  
     3.2E+00,  
     3.4E+00,  
     3.6E+00,  
     3.8E+00,  
     4.0E+00 };

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

  *n_data = *n_data + 1;

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

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

void bessel_j1_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    BESSEL_J1_VALUES returns some values of the J1 Bessel function.
//
//  Discussion:
//
//    In Mathematica, the function can be evaluated by:
//
//      BesselJ[1,x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 August 2004
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 21

  double fx_vec[N_MAX] = { 
      0.3275791375914652E+00,   
      0.6604332802354914E-01,   
     -0.3390589585259365E+00,   
     -0.5767248077568734E+00,   
     -0.4400505857449335E+00,   
      0.0000000000000000E+00,   
      0.4400505857449335E+00,   
      0.5767248077568734E+00,   
      0.3390589585259365E+00,   
     -0.6604332802354914E-01,   
     -0.3275791375914652E+00,   
     -0.2766838581275656E+00,   
     -0.4682823482345833E-02,   
      0.2346363468539146E+00,   
      0.2453117865733253E+00,   
      0.4347274616886144E-01,   
     -0.1767852989567215E+00,   
     -0.2234471044906276E+00,   
     -0.7031805212177837E-01,   
      0.1333751546987933E+00,   
      0.2051040386135228E+00 };

  double x_vec[N_MAX] = { 
     -5.0E+00,  
     -4.0E+00,  
     -3.0E+00,  
     -2.0E+00,  
     -1.0E+00,  
      0.0E+00,  
      1.0E+00,  
      2.0E+00,  
      3.0E+00,  
      4.0E+00,  
      5.0E+00,  
      6.0E+00,  
      7.0E+00,  
      8.0E+00,  
      9.0E+00,  
     10.0E+00,  
     11.0E+00,  
     12.0E+00,  
     13.0E+00,  
     14.0E+00,  
     15.0E+00 };

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

  *n_data = *n_data + 1;

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

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

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

//****************************************************************************80
//
//  Purpose:
//
//    BESSEL_JN_VALUES returns some values of the Jn Bessel function.
//
//  Discussion:
//
//    In Mathematica, the function can be evaluated by:
//
//      BesselJ[n,x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    29 April 2001
//
//  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, int *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 20

  double fx_vec[N_MAX] = { 
      0.1149034849319005E+00,  
      0.3528340286156377E+00,  
      0.4656511627775222E-01,  
      0.2546303136851206E+00,  
     -0.5971280079425882E-01,  
      0.2497577302112344E-03,  
      0.7039629755871685E-02,  
      0.2611405461201701E+00,  
     -0.2340615281867936E+00,  
     -0.8140024769656964E-01,  
      0.2630615123687453E-09,  
      0.2515386282716737E-06,  
      0.1467802647310474E-02,  
      0.2074861066333589E+00,  
     -0.1138478491494694E+00,  
      0.3873503008524658E-24,  
      0.3918972805090754E-18,  
      0.2770330052128942E-10,  
      0.1151336924781340E-04,  
     -0.1167043527595797E+00 };

  int nu_vec[N_MAX] = { 
     2,  2,  2,  2, 
     2,  5,  5,  5, 
     5,  5, 10, 10, 
    10, 10, 10, 20, 
    20, 20, 20, 20 };

  double x_vec[N_MAX] = { 
      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;
    *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
}
//****************************************************************************80

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

//****************************************************************************80
//
//  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 GNU LGPL 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

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

  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
}
//****************************************************************************80

void bessel_k0_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    BESSEL_K0_VALUES returns some values of the K0 Bessel function.
//
//  Discussion:
//
//    The modified Bessel functions In(Z) and Kn(Z) are solutions of
//    the differential equation
//
//      Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0.
//
//    The modified Bessel function K0(Z) corresponds to N = 0.
//
//    In Mathematica, the function can be evaluated by:
//
//      BesselK[0,x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    15 August 2004
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 20

  double fx_vec[N_MAX] = { 
     0.2427069024702017E+01,   
     0.1752703855528146E+01,   
     0.1114529134524434E+01,   
     0.7775220919047293E+00,   
     0.5653471052658957E+00,   
     0.4210244382407083E+00,   
     0.3185082202865936E+00,   
     0.2436550611815419E+00,   
     0.1879547519693323E+00,   
     0.1459314004898280E+00,   
     0.1138938727495334E+00,   
     0.6234755320036619E-01,   
     0.3473950438627925E-01,   
     0.1959889717036849E-01,   
     0.1115967608585302E-01,   
     0.6399857243233975E-02,   
     0.3691098334042594E-02,   
     0.1243994328013123E-02,   
     0.1464707052228154E-03,   
     0.1778006231616765E-04 };

  double x_vec[N_MAX] = { 
      0.1E+00,   
      0.2E+00,   
      0.4E+00,   
      0.6E+00,   
      0.8E+00,   
      1.0E+00,   
      1.2E+00,   
      1.4E+00,   
      1.6E+00,   
      1.8E+00,   
      2.0E+00,   
      2.5E+00,   
      3.0E+00,   
      3.5E+00,   
      4.0E+00,   
      4.5E+00,   
      5.0E+00,   
      6.0E+00,   
      8.0E+00,   
     10.0E+00  };

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

  *n_data = *n_data + 1;

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

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

void bessel_k0_int_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    BESSEL_K0_INT_VALUES returns some values of the Bessel K0 integral.
//
//  Discussion:
//
//    The function is defined by:
//
//      K0_INT(x) = Integral ( 0 <= t <= x ) K0(t) dt
//
//    The data was reported by McLeod.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    29 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Allan McLeod,
//    Algorithm 757:
//    MISCFUN: A software package to compute uncommon special functions,
//    ACM Transactions on Mathematical Software,
//    Volume 22, Number 3, September 1996, pages 288-301.
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 20

  double fx_vec[N_MAX] = { 
     0.78587929563466784589E-02,  
     0.26019991617330578111E-01,  
     0.24311842237541167904E+00,  
     0.39999633750480508861E+00,  
     0.92710252093114907345E+00,  
     0.12425098486237782662E+01,  
     0.14736757343168286825E+01,  
     0.15606495706051741364E+01,  
     0.15673873907283660493E+01,  
     0.15696345532693743714E+01,  
     0.15701153443250786355E+01,  
     0.15706574852894436220E+01,  
     0.15707793116159788598E+01,  
     0.15707942066465767196E+01,  
     0.15707962315469192247E+01,  
     0.15707963262340149876E+01,  
     0.15707963267948756308E+01,  
     0.15707963267948966192E+01,  
     0.15707963267948966192E+01,  
     0.15707963267948966192E+01  };

  double x_vec[N_MAX] = { 
       0.0009765625E+00,  
       0.0039062500E+00,  
       0.0625000000E+00,  
       0.1250000000E+00,  
       0.5000000000E+00,  
       1.0000000000E+00,  
       2.0000000000E+00,  
       4.0000000000E+00,  
       5.0000000000E+00,  
       6.0000000000E+00,  
       6.5000000000E+00,  
       8.0000000000E+00,  
      10.0000000000E+00,  
      12.0000000000E+00,  
      15.0000000000E+00,  
      20.0000000000E+00,  
      30.0000000000E+00,  
      50.0000000000E+00,  
      80.0000000000E+00,  
     100.0000000000E+00 };       

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

  *n_data = *n_data + 1;

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

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

void bessel_k1_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    BESSEL_K1_VALUES returns some values of the K1 Bessel function.
//
//  Discussion:
//
//    The modified Bessel functions In(Z) and Kn(Z) are solutions of
//    the differential equation
//
//      Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0.
//
//    The modified Bessel function K1(Z) corresponds to N = 1.
//
//    In Mathematica, the function can be evaluated by:
//
//      BesselK[1,x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    15 August 2004
//
//  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 *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 20

  double fx_vec[N_MAX] = { 
     0.9853844780870606E+01,   
     0.4775972543220472E+01,   
     0.2184354424732687E+01,   
     0.1302834939763502E+01,   
     0.8617816344721803E+00,   
     0.6019072301972346E+00,   
     0.4345923910607150E+00,   
     0.3208359022298758E+00,   
     0.2406339113576119E+00,   
     0.1826230998017470E+00,   
     0.1398658818165224E+00,   
     0.7389081634774706E-01,   
     0.4015643112819418E-01,   
     0.2223939292592383E-01,   
     0.1248349888726843E-01,   
     0.7078094908968090E-02,   
     0.4044613445452164E-02,   
     0.1343919717735509E-02,   
     0.1553692118050011E-03,   
     0.1864877345382558E-04 };

  double x_vec[N_MAX] = { 
      0.1E+00,   
      0.2E+00,   
      0.4E+00,   
      0.6E+00,   
      0.8E+00,   
      1.0E+00,   
      1.2E+00,   
      1.4E+00,   
      1.6E+00,   
      1.8E+00,   
      2.0E+00,   
      2.5E+00,   
      3.0E+00,   
      3.5E+00,   
      4.0E+00,   
      4.5E+00,   
      5.0E+00,   
      6.0E+00,   
      8.0E+00,   
     10.0E+00  };

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

  *n_data = *n_data + 1;

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

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

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

//****************************************************************************80
//
//  Purpose:
//
//    BESSEL_KN_VALUES returns some values of the Kn Bessel function.
//
//  Discussion:
//
//    The modified Bessel functions In(Z) and Kn(Z) are solutions of
//    the differential equation
//
//      Z^2 * W'' + Z * W' - ( Z^2 + N^2 ) * W = 0.
//
//    In Mathematica, the function can be evaluated by:
//
//      BesselK[n,x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 August 2004
//
//  Author:
//
//    John Burkard