function e = pendulum_double_conserved ( y )

%*****************************************************************************80
%
%% pendulum_double_conserved(): conserved quantity for pendulum_double_ode().
%
%  Discussion:
%
%    This quantity is the total energy.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    24 April 2021
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real Y(4): theta1, dtheta1/dt, theta2, dtheta2/dt.
%
%  Output:
%
%    real E: the energy.
%
  [ g, m1, m2, l1, l2, t0, y0, tstop ] = pendulum_double_parameters ( );

  e = 0.5 * ( m1 + m2 ) * l1^2 * y(:,2).^2 ...
    + 0.5 * m2 * l2^2 * y(:,4).^2 ...
    + m2 * l1 * l2 * y(:,2) .* y(:,4) .* cos ( y(:,1) - y(:,3) ) ...
    - ( m1 + m2 ) * g * l1 * cos ( y(:,1) ) ...
    - m2 * g * l2 * cos ( y(:,3) );

  return
end

