function y1 = kepler( t, y )
%
%  y1 = kepler( t, y )
%
%  Problema di Keplero. 
%  y = [q1 q2 p1 p2]';
%  Se nargin==1, calcola il valore
%  	dell'Hamiltoniana e del momento angolare in y ;
%  se nargin==2, calcola f .
%
if nargin==0
   epsil = 0.5; 
   y1    = zeros(1,4);
   y1(4) = sqrt((1+epsil)/(1-epsil));
   y1(1) = 1-epsil;
elseif nargin==1
   y    = t;
   roq  = y(:,1).^2+y(:,2).^2;
   rop  = y(:,3).^2+y(:,4).^2;
   sroq = sqrt(roq);
   y1   = [.5*rop-1./sroq y(:,1).*y(:,4)-y(:,2).*y(:,3)];
else % nargin==2
   y1(:,1:2) = y(:,3:4);
   roq       = (y(:,1).^2+y(:,2).^2);
   roq1      = roq.^(3/2);
   sroq      = sqrt(roq);
   y1(:,3)   = -y(:,1)./roq1;
   y1(:,4)   = -y(:,2)./roq1;
end
return      
   
