function h = arenstorf_conserved ( y )

%*****************************************************************************80
%
%% arenstorf_conserved() evaluates conservation for arenstorf_ode().
%
%  Discussion:
%
%    The angular momentum should also be conserved, but I haven't been able
%    to evaluate and verify this.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    10 June 2021
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real Y(4,N): the current solution.
%
%  Output:
%
%    real H(N): the value of energy.
%
  q1 = y(1,:);
  q2 = y(2,:);
  p1 = y(3,:);
  p2 = y(4,:);

  [ mu1, mu2, t0, y0, tstop ] = arenstorf_parameters ( );

  r1 = sqrt ( ( q1 + mu1 ).^2 + q2.^2 );
  r2 = sqrt ( ( q1 - mu2 ).^2 + q2.^2 );

  h = 0.5 * ( p1.^2 + p2.^2 ) ...
    - 0.5 * ( q1.^2 + q2.^2 ) ...
    - mu1 ./ r2 ...
    - mu2 ./ r1;

  return
end
