function double_pendulum_test ( )

%*****************************************************************************80
%
%% double_pendulum_test tests the double pendulum problem.
%
%  Discussion:
%
%    Options to explore include:
%    * vary the ODE solver;
%    * vary n, the number of steps for fixed step methods;
%    * specify the small, medium, or large initial condition;
%    * keep m1 = 1, but set m2 = 5, or m2 = 1/5.  
%    * keep l1 = 1, but set l2 = 5, or l2 = 1/5.
%
%    What happens if, at the initial condition:
%    * pendulum 1 is not moving:
%      z0 = [ 0, 0, 0, 5 ]
%    * pendulum 1 and pendulum 2 are at the same angle:
%      z0 = [ pi/3, 0, pi/3, 0 ]
%    * pendulum 1 and pendulum 2 are at opposite angles:
%      z0 = [ pi/3, 0, -pi/3, 0 ]
%    * the angular velocity is just large enough that the pendulum will
%      reach the top, and keep going.  (Then theta1 and theta2 will not
%      necessarily stay between -pi and pi.)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 November 2020
%
%  Author:
%
%    John Burkardt
%
  [ g, m1, m2, l1, l2 ] = double_pendulum_parameters ( );

  f = @ double_pendulum_deriv;
  t0 = 0.0;
  tstop = 50.0;
  tspan = [ t0, tstop ];
%
%  Choose the initial condition.
%
  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 ( f, tspan, z0, n );
  [ t, z ] = ode45    ( f, tspan, z0 );
%
%  Convert theta1 to (x1,y1) coordinates.
%
  x1 =   l1 * sin ( z(:,1) );
  y1 = - l1 * cos ( z(:,1) );
%
%  Convert theta2 to (x2,y2) coordinates.
%
  x2 = x1 + l2 * sin ( z(:,3) );
  y2 = y1 - l2 * cos ( z(:,3) );
%
%  Plot the (x,y) coordinates of m1 and m2 over time.
%
  figure ( 1 );
  clf ( );
  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 );
  clf ( );
  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 );
  clf ( );
  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 );
  clf ( );
  e = double_pendulum_energy ( z );
  clf ( );
  hold ( 'on' );
  plot ( t, e, 'r-', 'linewidth', 3 );
  plot ( [t0,tstop], [0.0, 0.0], 'k--', 'linewidth', 2 );
  hold ( 'off' );
  grid ( 'on' );
  xlabel ( '<--- t --->' );
  ylabel ( '<--- energy --->' );
  title ( 'double pendulum energy plot' )
  legend ( 'Energy', 'X-axis' );
  filename = 'double_pendulum_energy.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end

