function y = ochs ( t )

%*****************************************************************************80
%
%% pendulum_nonlinear_exact(): exact solution of nonlinear pendulum ODE.
%
%  Discussion:
%
%    This solution satisfies the pendulum ODE:
%
%    u' = v
%    v' = - g / l * sin ( u )
%
%    u(t0) = u0
%    v(t0) = v0
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    21 November 2020
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Karlheinz Ochs,
%    A comprehensive analytical solution of the nonlinear pendulum,
%    European Journal of Physics,
%    Volume 32, pages 479-490, 2011.
%
%    Paul D. Williams, 
%    Achieving seventh-order amplitude accuracy in leapfrog integrations, 
%    Monthly Weather Review,
%    Volume 141, Number 9, pages 3037-3051, September 2013.
%
%  Input:
%
%    real T(:): the current time.
%
%  Output:
%
%    real Y(2): the exact solution at time T.
%
  nt = length ( t );

  [ g, l, m, t0, y0, tstop ] = pendulum_nonlinear_parameters ( );

  theta0 = y0(1);
  thetadot0 = y0(2);
  
  k0 = sin ( theta0 / 2.0 );

  omega = sqrt ( g / l );
  ep = 4.0 * g / l;
  e0 = thetadot0^2 + ep * k0^2;
  k = sqrt ( e0 / ep );
  chi = 1.0 / k;
%
%  Equation 21 of the reference.
%
%  bob(t) is a linear function of t.
%
  bob = r8_sign ( thetadot0 ) * k * omega * ( t - t0 ) ...
    + jacobi_snk_inverse ( k0, chi );

  theta = 2.0 * r8_sign ( jacobi_cnk ( bob, chi ) ) .* ...
    asin ( jacobi_snk ( bob, chi ) );

  thetadot = r8_sign ( thetadot0 ) * sqrt ( e0 ) .* jacobi_dnk ( bob, chi );
%
%  Pack the results into Y.
%
  y = zeros ( nt, 2 );
  y(:,1) = theta;
  y(:,2) = thetadot;

  return
end

