clear xi yi xs ys

mu=h^2;
per=2;

% matrice di linearizzazione in [x0; y0]
A=[1-mu*cos(x0) 1; -mu*cos(x0) 1]
[V,D]=eig(A^per);
% assicurarsi che la componente y sia positiva
     if V(2,1) < 0
       V(:,1)=-V(:,1);
     end
     if V(2,2) < 0
       V(:,2)=-V(:,2);
     end
 % autovalori ed autovettori
     lam(1)=D(1,1)
     V(:,1)
     lam(2)=D(2,2)
     V(:,2)
% scegliere un punto lungo la separatrice instabile del linearizzato 
eps=1e-3
if lam(1) > lam(2)
     x1=x0 + eps*V(1,1);
     y1=y0 + eps*V(2,1);
else
     x1=x0 + eps*V(1,2);
     y1=y0 + eps*V(2,2);
end

% trovare la sua immagine per la mappa iterata per volte
% riportarsi all rettangolo 0< x < 2*pi ; 0< y < 2*pi
xb=x1;yb=y1;
for it=1:per
   yb=yb-mu*sin(xb);yb=yb-(2*pi)*floor(yb/(2*pi)); 
   xb=xb+yb;xb=xb-(2*pi)*floor(xb/(2*pi));
end
x2=xb;
y2=yb; 
% oppure a -pi <y<pi
%if y2 > pi
%y2=y2-2*pi;
%end
%segmento che unisce (x1,y1) con la sua immagine
np=100
yi=y1:(y2-y1)/np:y2;
xi=x1:(x2-x1)/np:x2;
nni=size(xi);ni=nni(2);
% scegliere un punto lungo la separatrice stabile del linearizzato 
eps=1e-3
if lam(1) < lam(2)
     x3=x0 + eps*V(1,1);
     y3=y0 + eps*V(2,1);
else
     x3=x0 + eps*V(1,2);
     y3=y0 + eps*V(2,2);
end

% trovare la sua immagine inversa per la mappa iterata per volte 
% riportarsi al solito rettangolo
xb=x3;yb=y3;
for it=1:per
   xb=xb-yb;xb=xb-(2*pi)*floor(xb/(2*pi));
   yb=yb+mu*sin(xb);yb=yb-(2*pi)*floor(yb/(2*pi));
end
x4=xb;
y4=yb;
%if y4 > pi
%y4=y4-2*pi;
%end

%segmento che unisce (x3,y3) con la sua immagine
ys=y3:(y4-y3)/np:y4;
xs=x3:(x4-x3)/np:x4;
nns=size(xs);ns=nns(2);

%inizializzazione della grafica
hold off
plot(xi,yi,'r')
hold on
plot(xs,ys,'g')
w=0.05;
pause
axis([x0-w x0+w y0-w y0+w])

% iterazione della mappa e della mappa inversa
nit=4;
for j=1:6
for ii=1:nit
% mappa standard in forma vettoriale, iterata per volte
for it=1:per
   yi=yi-mu*sin(xi);yi=yi-(2*pi)*floor(yi./(2*pi));
   xi=xi+yi;xi=xi-(2*pi)*floor(xi./(2*pi));
%  yi(yi>pi)=yi(yi>pi)-2*pi;
end
% tracciare la separatrice instabile
%plot(xi,yi,'.r')
for i=2:ni
   if abs(xi(i)-xi(i-1)) < .1 & abs(yi(i)-yi(i-1)) <.1
      plot([xi(i-1) xi(i)],[yi(i-1) yi(i)],'r')
   else
%        plot([xi(i-1) xi(i)],[yi(i-1) yi(i)],':r')
        plot([xi(i-1) xi(i)],[yi(i-1) yi(i)],'.r')
   end
end
% inversa della mappa standard^per in forma vettoriale
for it=1:per
   xs=xs-ys;xs=xs-(2*pi)*floor(xs./(2*pi));
   ys=ys+mu*sin(xs);ys=ys-(2*pi)*floor(ys./(2*pi));
%   ys(ys>pi)=ys(ys>pi)-2*pi;
end
%tracciare la separatrice stabile
%plot(xs,ys,'.g')
for i=2:ns
   if abs(xs(i)-xs(i-1)) < .1 & abs(ys(i)-ys(i-1)) <.1
      plot([xs(i-1) xs(i)],[ys(i-1) ys(i)],'g')
   else
%      plot([xs(i-1) xs(i)],[ys(i-1) ys(i)],':g')
      plot([xs(i-1) xs(i)],[ys(i-1) ys(i)],'.g')
   end
end
ii
pause
end
% scelta finestra
w=input('finestra?')
if w <3
axis([x0-w x0+w y0-w y0+w])
else
axis([0 2*pi 0 2*pi])
end
end




