function e = double_pendulum_energy ( z )

%*****************************************************************************80
%
%% double_pendulum_energy: enegy for the double pendulum.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 November 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real Z(:,4): theta1, dtheta1/dt, theta2, dtheta2/dt.
%
%  Output:
%
%    real E(:): the energy.
%
  [ g, m1, m2, l1, l2 ] = double_pendulum_parameters ( );

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

  return
end

