function dzdt = double_pendulum_deriv ( t, z )

%*****************************************************************************80
%
%% double_pendulum_deriv: double pendulum ordinary differential equation.
%
%  Discussion:
%
%    In the double pendulum problem, a rod of length l1 is fixed at
%    one end (0,0), and forms an angle theta1 with the downward vertical,
%    so that its endpoint is at (x1,y1) = (l1*cos(theta1),l1*sin(theta1)).
%    A weight of mass m1 is attached to this end of the first rod.
%
%    A second rod, of length l2 is also attached to this end of the
%    first rod.  It forms an angle theta2 with the downward vertical.
%    A weight of mass m2 is attached to the free end of the second rod.
%    The position of this weight is 
%    (x2,y2) = (x1,y1) + (l2*cos(theta2),l2*sin(theta2)).
%
%    Gravity has a force coefficient of G.
%
%    The unknowns are the angles and their time derivatives.
%    We store the unknowns in a vector Z:
%
%      z(1) =   theta(1)
%      z(2) = d theta(1) dt
%      z(3) =   theta(2)
%      z(4) = d theta(2) dt
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2020
%
%  Author:
%
%    Alexander Erlich (alexander.erlich@gmail.com)
%    Modifications by John Burkardt
%
%  Input:
%
%    real T: the current time.
%
%    real Z(4): the current solution.
%
%  Output:
%
%    real DZDT(4): the right hand side of the ODE.
%
%  Local:
%
%    real G: the gravitational force coefficient.
%
%    real M1, M2: the masses of pendulums 1 and 2.
%
%    real L1, L2: the lengths of pendulums 1 and 2.
%
  [ g, m1, m2, l1, l2 ] = double_pendulum_data ( );

  dzdt = zeros(4,1);

  dzdt(1) = z(2);

  dzdt(2) = - ...
  ( ...
    (...
      g * ( 2.0 * m1 + m2 ) * sin ( z(1) ) + ...
      m2 * ...
      ( ...
        g * sin ( z(1) - 2.0 * z(3) ) + 2.0 * ...
        ( ...
          l2 * z(4)^2 + l1 * z(2)^2 * cos ( z(1) - z(3) ) ...
        ) ...
        * sin ( z(1) - z(3) ) ...
      ) ...
    ) ...
    / ...
    ( ...
     2.0 * l1 * ( m1 + m2 - m2 * cos ( z(1) - z(3) )^2 ) ...
    ) ...
  );

  dzdt(3) = z(4);

  dzdt(4) = ...
  ( ...
    ( ...
      ( m1 + m2 ) * ( l1 * z(2)^2 + g * cos ( z(1) ) ) ...
      + l2 * m2 * z(4)^2 * cos ( z(1) - z(3) ) ...
    ) * ...
    sin ( z(1) - z(3) ) ...
  ) ...
  / ...
  ( ...
    l2 * ( m1 + m2 - m2 * cos ( z(1) - z(3) )^2 ) ...
  );

  return
end

