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:
%
%    21 March 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_data ( );

  t1 = 0.5 * m1 * l1^2 * z(:,2).^2;
  t2 = 0.5 * 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) );
  v1 = m1 * g * l1 * cos ( z(:,1) );
  v2 = m2 * g * l1 * cos ( z(:,3) ) ...
     + m2 * g * l2 * cos ( z(:,3) );

  e = t1 + t2 + v1 + v2;

  return
end

