function yprime = pendulum_elastic_deriv ( t, y )

%*****************************************************************************80
%
%% pendulum_elastic_deriv(): right hand side of pendulum_elastic_ode().
%
%  Discussion:
%
%    Y1 is the angular displacement;
%    Y2 is the angular velocity;
%    Y3 is the strain
%    Y4 is the time derivative of strain
%
%    G is the gravitational coefficient.
%    K is the spring constant;
%    L is the length of the pendulum;
%    M is the mass;
%
%    y1' = y2
%    y2' = ?
%    y3' = y4;
%    y4' = ?
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    23 November 2020
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul D Williams,
%    The RAW Filter: An Improvement to the Robert–Asselin Filter 
%    in Semi-Implicit Integrations,
%    Monthly Weather Review,
%    Volume 139, pages 1996-2007, June 2011.
%
%  Input:
%
%    real T, the current time.
%
%    real Y(4,1), the current state values.
%
%  Output:
%
%    real YPRIME(4,1), the time derivatives of the current state values.
%
  [ g, k, l, m, t0, y0, tstop ] = pendulum_elastic_parameters ( );

  wlo = sqrt ( g / l );
  whi = sqrt ( k / m );

  th = y(1);
  thdot = y(2);
  nu = y(3);
  nudot = y(4);

  yprime(1,1) = thdot;
  yprime(2,1) = ( - wlo^2 * sin ( th ) - 2.0 * thdot * nudot ) / ( 1.0 + nu );
  yprime(3,1) = nudot;
  yprime(4,1) = - wlo^2 * ( 1.0 - cos ( th ) ) ...
                - whi^2 * nu + ( 1.0 + nu ) * thdot^2;

  return
end

