function dydt = pendulum_double_deriv ( t, y )

%*****************************************************************************80
%
%% pendulum_double_deriv(): right hand side of pendulum_double_ode().
%
%  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 Y:
%
%      y(1) =   theta(1)
%      y(2) = d theta(1) dt
%      y(3) =   theta(2)
%      y(4) = d theta(2) dt
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    24 April 2021
%
%  Author:
%
%    Alexander Erlich (alexander.erlich@gmail.com)
%    Modifications by John Burkardt
%
%  Input:
%
%    real T: the current time.
%
%    real Y(4): the current solution.
%
%  Output:
%
%    real DYDT(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, t0, y0, tstop ] = pendulum_double_parameters ( );

  dydt = zeros(4,1);

  dydt(1) = y(2);

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

  dydt(3) = y(4);

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

  return
end

