function [ p, np ] = pendulum_ode_period_ode45 ( )

%*****************************************************************************80
%
%% pendulum_ode_period_ode45() uses ode45() to detect the period of pendulum_ode().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    21 December 2021
%
%  Author:
%
%    John Burkardt
%
%  Output:
%
%    real p: an estimate of the period.
%
%    integer np: the number of full periods that were observed.
%

%
%  Get parameter values.
%
  [ g, l, m, t0, y0, tstop ] = pendulum_parameters ( );
%
%  Solve the ODE.
%
  f = @ pendulum_deriv;
  tspan = [ t0, tstop ];
  options = odeset ( 'Events', @ pendulum_event );
  [ t, y, te, ye, ie ] = ode45 ( f, tspan, y0, options );
%
%  Count the number of events.
%
  if ( length ( te ) < 3 )
    np = NaN;
    p = NaN;
  else
    np = length ( te );
    p = ( te(np) - te(np-2) );
    np = floor ( np / 2 );
  end

  return
end
function [ value, isterminal, direction ] = pendulum_event ( t, y )

%*****************************************************************************80
%
%% pendulum_event() watches for repetition of the initial value for pendulum_ode().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    21 December 2021
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real t: the current time.
%
%    real y(2): the current solution.
%
%  Output:
%
%    real value: defines an event when y(1) = 0.0.
%
%    logical isterminal: is true if the event should cause
%    the calculation to halt.
%
%    logical direction: is 0 if all zeroes are to be located,
%    +1 or -1 if events only occur when value is increasing or decreasing.
%
  [ g, l, m, t0, y0, tstop ] = pendulum_parameters ( );

  value = y(1);
  isterminal = false;
  direction = 0;

  return
end

