function double_pendulum_test ( )

%*****************************************************************************80
%
%% double_pendulum_test tests the double pendulum problem.
%
%  Discussion:
%
%    Options to explore include:
%    * varying the ODE solver;
%    * varying n, the number of steps for fixed step methods;
%    * specifying the small, medium, or large initial condition;
%    * varying parameters in double_precision_data(), particularly m2 and l2;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 March 2020
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../scripts' );

  f = @ double_pendulum_deriv;
  tspan = [ 0.0, 50.0 ];
  z0 = [ 0.25; 0.0; 0.0; 0.0 ];
% z0 = [ pi/3.0; 0.0; 0.0; 2.0 ];
% z0 = [ pi; 0.0; pi; 5.0 ];
  n = 100000;

% [ t, z ] = euler ( f, tspan, z0, n );
% [ t, z ] = rk4 ( f, tspan, z0, n );
% [ t, z ] = midpoint_fixed ( f, tspan, z0, n );
  [ t, z ] = ode45 ( f, tspan, z0 );

  [ g, m1, m2, l1, l2 ] = double_pendulum_data ( );

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

  x2 = x1 + l2 * sin ( z(:,3) );
  y2 = y1 - l2 * cos ( z(:,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', 'y1' );
  title ( 'double pendulum (x1,y1)' )
  filename = 'double_pendulum_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', 'y2' );
  title ( 'double pendulum (x2,y2)' )
  filename = 'double_pendulum_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 = 'double_pendulum_phase.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Energy plot.
%
  figure ( 4 );
  e = double_pendulum_energy ( z );
  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 = 'double_pendulum_energy.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end

