function p = vanderpol_period ( mu )

%*****************************************************************************80
%
%% vanderpol_period() estimates the period of vanderpol_ode().
%
%  Discussion:
%
%    Five approximations are available for the period.  These estimates
%    seen to disagree for low values of mu.  The formula of Urabe
%    is evaluated here.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 December 2021
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Mary Cartwright,
%    Van der Pol's equation for relaxation oscillations,
%    Contributions to the theory of nonlinear oscillations, II,
%    Annals of Mathematical Studies 29,
%    Princeton, 1952, pages 3-18.
%
%    John D Cook,
%    Calculating the period of Van der Pol Oscillators,
%    https://www.johndcook.com/blog/2019/12/26/van-der-pol-period/
%    26 December 2019.
%
%    Anatolii Dorodnitsyn,
%    Asymptotic solution of Van der Pol's equation,
%    Prikladnaia Matematika i Mekhanika, 
%    (Journal of Applied Mathematics and Mechanics),
%    American Mathematical Society Translation 88.
%
%    Roger Grimshaw,
%    Nonlinear Ordinary Differential Equations,
%    CRC Press, 1991, pages 160-163.
%
%    Minoru Urabe,
%    Journal of Science of Hiroshima University,
%    Series A, Volume 24, Number 2, October 1960, pages 197-199.
%
%  Input:
%
%    real MU: the parameter.
%
%  Output:
%
%    real P: the estimated period.
%
  if ( mu == 0.0 )

    p1 = 2.0 * pi;
    p2 = 2.0 * pi;
    p3 = 2.0 * pi;
    p4 = 2.0 * pi;
    p5 = 2.0 * pi;

  else
%
%  Cook
%
    p1 = ( 3.0 - 2.0 * log ( 2.0 ) ) * mu;
%
%  Cartwright
%
    p2 = ( 3.0 - 2.0 * log ( 2.0 ) ) * mu + 2.0 * pi / mu^(1.0/3.0);
%
%  Grimshaw.
%
    alpha = 2.338;
    p3 = ( 3.0 - 2.0 * log ( 2.0 ) ) * mu ...
      + 2.0 * alpha / mu^(1.0/3.0);
%
%  Dorodnicyn.
%
    alpha = 2.338107;
    b0 = 0.1723;
    d = 0.4889;

    p4 = ( 3.0 - 2.0 * log ( 2.0 ) ) * mu ...
      + 3.0 * alpha / mu^(1.0/3.0) ...
      - 22.0 / 9.0 *  log ( mu ) / mu ...
      + ( 3.0 * log ( 2.0 ) - log ( 3.0 ) - 1.0 / 6.0 + b0 - 2.0 * d ) / mu;
%
%  Urabe
%
    alpha = 2.338107;
    b0 = 0.1723;
    d = 0.4889;

    p5 = ( 3.0 - 2.0 * log ( 2.0 ) ) * mu ...
      + 3.0 * alpha / mu^(1.0/3.0) ...
      - 1.0 / 3.0 *  log ( mu ) / mu ...
      + ( 3.0 * log ( 2.0 ) - log ( 3.0 ) - 1.5 + b0 - 2.0 * d ) / mu;

  end

  p = p5;

  return
end

