clear lp lm da d x y
x(1)=x0
y(1)=y0
per=1
h=1
mu=h^2;
nit=input('assegna numero iterazioni ')

An=eye(2);
d0=[1 ; 1]*1/sqrt(2)
xb=x0;yb=y0;

for n=1:nit
 
% iterazioni della mappa, e jacobiana
    for p=1:per
        A=[1-mu*cos(xb) 1; -mu*cos(xb) 1];
        An=A*An;
        yb=yb-mu*sin(xb);
        xb=xb+yb;
        xb=xb-2*pi*floor(xb/(2*pi));
        yb=yb-2*pi*floor(yb/(2*pi));
%        if yb > pi 
%          yb=yb-2*pi;
%        end
    end
    ll=eig(An);da(n)=det(An);
    if abs(ll(1)) > abs(ll(2))
       lp(n)=ll(1);lm(n)=ll(2);
    else
       lp(n)=ll(2);lm(n)=ll(1);
    end
    vv=An*d0;
    d(n)=sqrt(vv(1)^2+vv(2)^2);

end
hold off
plot(log10(abs(lp)),'.r')   
hold on
lm(abs(lm)<1e-20)=1e-20;
plot(log10(abs(lm)),'.g')
da(abs(da)<1e-20)=1e-20;
plot(log10(abs(da)),'.k')
plot(log10(d),'b')
nn=1:1:nit;
pp=polyfit(nn,log(d),1);lce=pp(1)
pp=polyfit(nn,log10(d),1);
plot(nn,polyval(pp,nn),'k')
ylabel('log10')
xlabel('iterazione numero')
