function pendulum_comparison_ode_test ( )

%*****************************************************************************80
%
%% pendulum_comparison_ode_test() compares linear and nonlinear pendulum models.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    16 April 2022
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../pendulum_comparison_ode' );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'pendulum_comparison_ode_test():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  pendulum_comparison_ode() compares linear and nonlinear pendulum models.\n' );

  [ g, l, m, t0, y0 ] = pendulum_parameters ( );
%
%  We set the initial deflect to 9/10 of the "upright" deflection of pi.
%
  y0(1) = 0.9 * pi;
  y0(2) = 0.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  parameters:\n' );
  fprintf ( 1, '    gravity (m/s^2) = %g\n', g );
  fprintf ( 1, '    length (m)      = %g\n', l );
  fprintf ( 1, '    mass (g)        = %g\n', m );
  fprintf ( 1, '    t0 (s)          = %g\n', t0 );
  fprintf ( 1, '    u0 (rad)        = %g\n', y0(1) );
  fprintf ( 1, '    v0 (rad/s)      = %g\n', y0(2) );

  tstop = 20.0;
  tspan = [ t0, tstop ];

  f1 = @ pendulum_deriv;
  [ t1, y1 ] = ode45 ( f1, tspan, y0 );

  f2 = @ pendulum_nonlinear_deriv;
  [ t2, y2 ] = ode45 ( f2, tspan, y0 );
%
%  Plot u = theta.
%
  figure ( 1 );
  clf ( );
  plot ( t1, y1(:,1), 'r-', ...
         t2, y2(:,1), 'b-', 'linewidth', 2 );
  grid ( 'on' );
  xlabel ( '<-- Time -->' );
  ylabel ( '<-- Angular deflection (radians) -->' );
  title ( 'Linear/nonlinear pendulum, theta' );
  legend ( 'linear', 'nonlinear' );
  filename = 'pendulum_comparison_theta.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved in "%s"\n', filename );
  close ( );
%
%  Plot v = thetadot.
%
  figure ( 2 );
  clf ( );
  plot ( t1, y1(:,2), 'r-', ...
         t2, y2(:,2), 'b-', 'linewidth', 2 );
  grid ( 'on' );
  xlabel ( '<-- Time -->' );
  ylabel ( '<-- Angular velocity (radians/sec) -->' );
  title ( 'Linear/nonlinear pendulum, thetadot' );
  legend ( 'linear', 'nonlinear' );
  filename = 'pendulum_comparison_thetadot.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved in "%s"\n', filename );
  close ( );
%
%  Phase plots.
%
  figure ( 3 );
  clf ( );
  plot ( y1(:,1), y1(:,2), 'r-', ...
         y2(:,1), y2(:,2), 'b-', 'linewidth', 2 );
  grid ( 'on' );
  xlabel ( '<-- Theta -->' );
  ylabel ( '<-- Thetadot -->' );
  title ( 'Linear/nonlinear pendulum, phase plot' );
  legend ( 'linear', 'nonlinear' );
  axis ( 'equal' );
  filename = 'pendulum_comparison_phase.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved in "%s"\n', filename );
  close ( );
%
%  Plot energy.
%
  h1 = pendulum_conserved ( y1 );
  h2 = pendulum_nonlinear_conserved ( y2 );

  figure ( 4 );
  clf ( );
  hold ( 'on' );
    plot ( t1, h1, 'r-', 'linewidth', 2 );
    plot ( t2, h2, 'b-', 'linewidth', 2 );
    plot ( [t0, tstop ], [ 0.0, 0.0 ], 'k-', 'linewidth', 2 );
  hold ( 'off' );
  grid ( 'on' );
  xlabel ( '<-- Time -->' );
  ylabel ( '<-- Energy -->' );
  title ( 'Linear/nonlinear pendulum, energy' );
  legend ( 'linear', 'nonlinear' );
  filename = 'pendulum_comparison_energy.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved in "%s"\n', filename );
  close ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'pendulum_comparison_ode_test():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../pendulum_comparison_ode' );

  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

