function pendulum_double_ode_test ( )

%*****************************************************************************80
%
%% pendulum_double_ode_test tests pendulum_double_ode().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    11 November 2020
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../pendulum_double_ode' )

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'pendulum_double_ode_test:\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Test pendulum_double_ode using ode45.\n' );

  [ g, m1, m2, l1, l2, t0, y0 ] = pendulum_double_parameters ( );

  f = @ pendulum_double_deriv;
  tspan = [ t0, 50.0 ];

  [ t, theta ] = ode45 ( f, tspan, y0 );

  x1 =   l1 * sin ( theta(:,1) );
  y1 = - l1 * cos ( theta(:,1) );

  x2 = x1 + l2 * sin ( theta(:,3) );
  y2 = y1 - l2 * cos ( theta(:,3) );
%
%  Solution plot.
%
  figure ( 1 );
  plot ( t, x1, 'g', ...
         t, y1, 'r', ...
         'LineWidth', 2 );
  grid ( 'on' );
  xlabel ( 'Time' );
  ylabel ( 'x1(t), y1(t)' );
  legend ( 'x1(t)', 'y1(t)' );
  title ( 'double pendulum (x1,y1)' )
  filename = 'pendulum_double_x1y1.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  figure ( 2 );
  plot ( t, x2, 'g', ...
         t, y2, 'r', ...
         'LineWidth', 2 );
  grid ( 'on' );
  xlabel ( 'Time' );
  ylabel ( 'x2(t), y2(t)' );
  legend ( 'x2(t)', 'y2(t)' );
  title ( 'double pendulum (x2,y2)' )
  filename = 'pendulum_double_x2y2.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Phase plane plot.
%
  figure ( 3 );
  hold ( 'on' );
  plot ( x1, y1, 'r-', 'linewidth', 3 );
  plot ( x2, y2, 'b-', 'linewidth', 1 );
  hold ( 'off' );
  grid ( 'on' );
  xlabel ( '<--- x(t) --->' );
  ylabel ( '<--- y(t) --->' );
  legend ( '(x1,y1)', '(x2,y2)' );
  axis ( 'equal' );
  title ( 'double pendulum phase' )
  filename = 'pendulum_double_phase.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Energy plot.
%
  figure ( 4 );
  e = pendulum_double_conserved ( theta );
  clf ( )
  hold ( 'on' );
  plot ( t, e );
  plot ( [ t(1),t(end)], [0.0,0.0], 'k-', 'linewidth', 2 );
  hold ( 'off' );
  grid ( 'on' );
  xlabel ( '<--- t --->' );
  ylabel ( '<--- energy --->' );
  title ( 'double pendulum energy plot' )
  filename = 'pendulum_double_energy.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'pendulum_double_ode_test:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../pendulum_double_ode' )

  return
end
