function rigid_body_midpoint ( n )

%*****************************************************************************80
%
%% rigid_body_midpoint() tests rigid_body_ode with the midpoint method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    26 May 2021
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N: the number of steps to take.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'rigid_body_midpoint():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Solve rigid_body_ode using the midpoint method.\n' );
  fprintf ( 1, '  rigid_body_ode() models motion on the surface of a sphere.\n' );

  [ i1, i2, i3, t0, y0, tstop ] = rigid_body_parameters ( );

  tspan = [ t0, tstop ];
  f = @rigid_body_deriv;

  [ t, xyz ] = midpoint ( f, tspan, y0, n );

  n = size ( t, 1 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of equal steps n = %d\n', n );

  if ( 500 < n )
    ilo = n - 500;
  else
    ilo = 1;
  end

  figure ( 1 );
  plot ( t(ilo:n), xyz(ilo:n,:), 'linewidth', 3 );
  grid ( 'on' );
  xlabel ( '<-- X -->', 'Fontsize', 16 );
  ylabel ( '<-- Y -->', 'Fontsize', 16 );
  title ( 'rigid body, midpoint: (X(t),Y(t),Z(t))' );
  legend ( 'X(t)', 'Y(t)', 'Z(t)' );
  filename = 'rigid_body_midpoint_plot.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  figure ( 2 );
  plot3 ( xyz(ilo:n,1), xyz(ilo:n,2), xyz(ilo:n,3), 'linewidth', 3 );
  grid ( 'on' );
  xlabel ( '<-- X -->', 'Fontsize', 16 );
  ylabel ( '<-- Y -->', 'Fontsize', 16 );
  zlabel ( '<-- Z -->', 'Fontsize', 16 );
  title ( 'rigid body, midpoint: (x,y,z)(t)', 'fontsize', 24 );
  view ( 3 );
  filename = 'rigid_body_midpoint_plot3d.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  [ h1, h2 ] = rigid_body_conserved ( xyz );

% h1 = xyz(:,1).^2 + xyz(:,2).^2 + xyz(:,3).^2;

  figure ( 3 );
  clf ( );
  hold ( 'on' );
    plot ( t, h1, 'r-', 'linewidth', 3 );
    plot ( [ t(1), t(n) ], [h1(1),h1(1) ], 'b--', 'linewidth', 3 );
    plot ( [ t(1), t(n) ], [0.0, 0.0 ], 'k-', 'linewidth', 3 );
  hold ( 'off' );
  grid ( 'on' );
  xlabel ( '<-- T -->', 'Fontsize', 16 );
  ylabel ( '<-- H1(T) -->', 'Fontsize', 16 );
  title ( 'rigid_body midpoint: H1(T)' );
  filename = 'rigid_body_midpoint_h1.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

%  h2 = (xyz(:,1)/i1).^2 + (xyz(:,2)/i2).^2 + (xyz(:,3)/i3).^2;

  figure ( 4 );
  clf ( );
  hold ( 'on' );
    plot ( t, h2, 'r-', 'linewidth', 3 );
    plot ( [ t(1), t(n) ], [h2(1),h2(1) ], 'b--', 'linewidth', 3 );
    plot ( [ t(1), t(n) ], [0.0, 0.0 ], 'k-', 'linewidth', 3 );
  hold ( 'off' );
  grid ( 'on' );
  xlabel ( '<-- T -->', 'Fontsize', 16 );
  ylabel ( '<-- H2(T) -->', 'Fontsize', 16 );
  title ( 'rigid_body midpoint(): H2(T)' );
  filename = 'rigid_body_midpoint_h2.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end
