function [y1,T] = cassini(t,y)
%    
%   Usage:
%   [y0,T] = cassini()    or    y' = cassini(t,y)   or  H = cassini(y)
%   
%   y = [q p]. Cassini problem with polynomial Hamiltonian
%   H(q,p) = (q^2+p^2)^2-alfa^2(q^2-p^2),   alfa^2 = 10.
%   T is the period of the orbit starting at y0.
%
%   See pages 134-136 of
%   L.Brugnano, F.Iavernaro.
%   Line Integral Methods for Conservative Problems.
%   Chapman and Hall/CRC - 2015
%   Series: Monographs and Research Notes in Mathematics.
%  
alfa2 = 10;
if nargin==1 % Hamiltonian
    y = t;
    q = y(:,1); 
    p = y(:,2); 
    q2 = q.^2;
    p2 = p.^2;
    y1 = (q2+p2).^2-alfa2*(q2-p2);
elseif nargin==2 % f(y)
    q  = y(:,1); 
    p  = y(:,2);
    c1 = 4*(q.^2+p.^2);
    alfa2 = 2*alfa2;
    y1 = [alfa2*p+p.*c1 alfa2*q-q.*c1];
elseif nargin==0 % y0 and period
    y1=[0 1e-6]; T = 3.131990057003955;  
else
    error
end
return