function y1 = kepler2( t, q )
%
%   Usage:
%   q'' = M*kepler2(t,q)  or  H = kepler2([q,p])  or  M = kepler2() == 1 
%
%   Kepler problem in second-order form. 
%
%   See pages 131-134  of
%   L.Brugnano, F.Iavernaro.
%   Line Integral Methods for Conservative Problems.
%   Chapman and Hall/CRC - 2015
%   Series: Monographs and Research Notes in Mathematics.
%
if nargin==0 % M
   y1 = 1;   % is a scalar matrix, the scalar suffices
elseif nargin==1   % Hamiltonian
   q = t(:,1:2); 
   p = t(:,3:4);
   sroq = sqrt( q(:,1).^2 + q(:,2).^2 );
   rop  = p(:,1).^2+p(:,2).^2;
   y1   = .5*rop - 1./sroq;
else               % M * nabla U(q)
   roq1 = ( q(:,1).^2+q(:,2).^2 ).^(3/2);
   y1   = [-q(:,1)./roq1 -q(:,2)./roq1];
end
return      
   
