function y0 = henon_heiles_initial ( e, lam, t0, u, v, vp )

%*****************************************************************************80
%
%% henon_heiles_initial(): consistent initial condition for henon_heiles_ode().
%
%  Discussion:
%
%    It must be the case that 0 <= 2e-vp^2-u^2-v^2-2u^2v+2v^3/3.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    23 April 2021
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Michel Henon, Carl Heiles,
%    The applicability of the third integral of motion: 
%    Some numerical experiments,
%    The Astronomical Journal,
%    Volume 69, pages 73-79, 1964.
%
%  Input:
%
%    real E: proposed initial energy.
%
%    real LAM: a parameter.
%
%    real T0: the initial time.
%
%    real U, V, VP: proposed initial values for three coordinates.
%
%  Output:
%
%    real Y0[4]: a consistent initial condition.
%
  upsq = 2.0 * e - vp.^2 - u.^2 - v.^2  ...
    - 2.0 * lam * ( u.^2 .* v - v.^3 / 3.0 );

  if ( upsq < 0.0 )
    fprintf ( 1, 'henon_heiles_choose_xp - Fatal error!\n' );
    fprintf ( 1, '  2e-vp^2-u^2-v^2-2u^2v+2v^3/3 < 0\n' );
    error ( 'henon_heiles_choose_xp - Fatal error!' );
  end

  up = sqrt ( upsq );

  y0 = [ u, up, v, vp ];

  return
end

