function fresnel_test ( )

%*****************************************************************************80
%
%% fresnel_test() tests fresnel().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    07 July 2025
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../fresnel' );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'fresnel_test():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Test fresnel().\n' );

  fresnel_cos_values_test ( );
  fresnel_sin_values_test ( );
  fresnel_cos_plot ( );
  fresnel_sin_plot ( );
  fresnel_phase_plot ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'fresnel_test():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../fresnel' );

  return
end
function fresnel_cos_plot ( )

%*****************************************************************************80
%
%% fresnel_cos_plot() tests fresnel_cos().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    07 July 2025
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'fresnel_cos_plot():\n' );
  fprintf ( 1, '  Plot (X,C(X)), where C(X) is the Fresnel cosine integral.\n' );

  a = -5.0;
  b = +5.0;
  n = 501;
  x = linspace ( a, b, n );
  y = zeros ( n, 1 );

  for i = 1 : n
    y(i) = fresnel_cos ( x(i) );
  end

  clf ( );
  plot ( x, y, 'b-', 'linewidth', 2 );
  grid ( 'On' );
  xlabel ( '<-- X -->' );
  ylabel ( '<-- C(X) -->' );
  title ( 'C(X) = Fresnel cosine integral' );
  filename = 'fresnel_cos_plot.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
  close ( );

  return
end
function [ n_data, x, fx ] = fresnel_cos_values ( n_data )

%*****************************************************************************80
%
%% fresnel_cos_values() returns values of the Fresnel cosine integral function.
%
%  Discussion:
%
%    The Fresnel cosine integral is defined by:
%
%      C(X) = integral ( 0 <= T <= X ) cos ( PI * T^2 / 2 ) dT
%
%    In Mathematica, the function can be evaluated by:
%
%      FresnelC[x]
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    16 November 2015
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Input:
%
%    integer N_DATA.  The user sets N_DATA to 0 before the first call.  
%    Thereafter, it should simply be the value returned by the previous call.
%
%  Output:
%
%    integer N_DATA.  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.
%
%    real X, the argument of the function.
%
%    real FX, the value of the function.
%
  n_max = 16;

  fx_vec = [ ...
     0.0000000000000000E+00, ...
     0.1999210575944531E+00, ...
     0.3974807591723594E+00, ...
     0.5810954469916523E+00, ...
     0.7228441718963561E+00, ...
     0.7798934003768228E+00, ...
     0.7154377229230734E+00, ...
     0.5430957835462564E+00, ...
     0.3654616834404877E+00, ...
     0.3336329272215571E+00, ...
     0.4882534060753408E+00, ...
     0.6362860449033195E+00, ...
     0.5549614058564281E+00, ...
     0.3889374961919690E+00, ...
     0.4674916516989059E+00, ...
     0.6057207892976856E+00 ];

  x_vec = [ ...
     0.0E+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 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  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);
    fx = fx_vec(n_data);
  end

  return
end
function fresnel_cos_values_test ( )

%*****************************************************************************80
%
%% fresnel_cos_values_test() tests fresnel_cos().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'fresnel_cos_values_test():\n' );
  fprintf ( 1, '  fresnel_cos_values() stores values of\n' );
  fprintf ( 1, '  the Fresnel cosine integral C(X).\n' );
  fprintf ( 1, '  fresnel_cos(x) computes the value directly.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '        X               C(X)                      C(X)\n' );
  fprintf ( 1, '                        Tabulated                 Computed\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( true )

    [ n_data, x, fx1 ] = fresnel_cos_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = fresnel_cos ( x );

    fprintf ( 1, '  %12f  %24.16f  %24.16f\n', x, fx1, fx2 );

  end

  return
end
function fresnel_phase_plot ( )

%*****************************************************************************80
%
%% fresnel_phase_plot() tests fresnel_cos() and fresnel_sin().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    07 July 2025
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'fresnel_phase_plot():\n' );
  fprintf ( 1, '  Plot (C(X),S(X)), where C(X) and S(X) are \n' );
  fprintf ( 1, '  the Fresnel cosine and sine integrals.\n' );

  a = -5.0;
  b = +5.0;
  n = 501;
  x = linspace ( a, b, n );
  c = zeros ( n, 1 );
  s = zeros ( n, 1 );

  for i = 1 : n
    c(i) = fresnel_cos ( x(i) );
    s(i) = fresnel_sin ( x(i) );
  end

  clf ( );
  plot ( c, s, 'b-', 'linewidth', 2 );
  grid ( 'On' );
  xlabel ( '<-- C(X) -->' );
  ylabel ( '<-- S(X) -->' );
  title ( 'Fresnel Phase Plot' );
  filename = 'fresnel_phase_plot.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
  close ( );

  return
end
function fresnel_sin_plot ( )

%*****************************************************************************80
%
%% fresnel_sin_plot() tests fresnel_sin().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    07 July 2025
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'fresnel_sin_plot():\n' );
  fprintf ( 1, '  Plot (X,S(X)), where S(X) is the Fresnel sine integral.\n' );

  a = -5.0;
  b = +5.0;
  n = 501;
  x = linspace ( a, b, n );
  y = zeros ( n, 1 );

  for i = 1 : n
    y(i) = fresnel_sin ( x(i) );
  end

  clf ( );
  plot ( x, y, 'b-', 'linewidth', 2 );
  grid ( 'On' );
  xlabel ( '<-- X -->' );
  ylabel ( '<-- S(X) -->' );
  title ( 'S(X) = Fresnel sine integral' );
  filename = 'fresnel_sin_plot.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
  close ( );

  return
end
function [ n_data, x, fx ] = fresnel_sin_values ( n_data )

%*****************************************************************************80
%
%% fresnel_sin_values() returns some values of the Fresnel sine integral function.
%
%  Discussion:
%
%    The Fresnel sine integral is defined by
%
%      S(X) = integral ( 0 <= T <= X ) sin ( pi * T^2 / 2 ) dT
%
%    In Mathematica, the function can be evaluated by:
%
%      FresnelS[x]
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    16 November 2015
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Input:
%
%    integer N_DATA.  The user sets N_DATA to 0 before the first call.  
%    Thereafter, it should simply be the value returned by the previous call.
%
%  Output:
%
%    integer N_DATA.  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.
%
%    real X, the argument of the function.
%
%    real FX, the value of the function.
%
  n_max = 16;

  fx_vec = [ ...
     0.0000000000000000E+00, ...
     0.4187609161656762E-02, ...
     0.3335943266061318E-01, ...
     0.1105402073593870E+00, ...
     0.2493413930539178E+00, ...
     0.4382591473903548E+00, ...
     0.6234009185462497E+00, ...
     0.7135250773634121E+00, ...
     0.6388876835093809E+00, ...
     0.4509387692675831E+00, ...
     0.3434156783636982E+00, ...
     0.4557046121246569E+00, ...
     0.6196899649456836E+00, ...
     0.5499893231527195E+00, ...
     0.3915284435431718E+00, ...
     0.4963129989673750E+00 ];

  x_vec = [ ...
     0.0E+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 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  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);
    fx = fx_vec(n_data);
  end

  return
end
function fresnel_sin_values_test ( )

%*****************************************************************************80
%
%% fresnel_sin_values_test() tests fresnel_sin().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    07 July 2025
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'fresnel_sin_values_test():\n' );
  fprintf ( 1, '  fresnel_sin_values() stores values of\n' );
  fprintf ( 1, '  the Fresnel sine integral S(X).\n' );
  fprintf ( 1, '  fresnel_sin(x) computes the value directly.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '        X               S(X)                      S(X)\n' );
  fprintf ( 1, '                        Tabulated                 Computed\n' );

  n_data = 0;

  while ( true )

    [ n_data, x, fx1 ] = fresnel_sin_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = fresnel_sin ( x );

    fprintf ( 1, '  %12f  %24.16f  %24.16f\n', x, fx1, fx2 );

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% timestamp() prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

