
% Iterating the Henon map near fixed points to check 
% for contraction/expansion of an e-ball: fig 2.9

clear

henon_map = inline('[a-x.^2+b*y x]','a','b','x','y');

a=0*0.27
b=0.4

%fixed points:
xs1 = (b-1+sqrt((1-b)^2+4*a))/2
ys1 = xs1;
xs2 = (b-1-sqrt((1-b)^2+4*a))/2
ys2 = xs2;

figure(1)
clf
lw=3;
lw=2;
plot(xs1,ys1,'r+','LineWidth',lw)
%axis([-1 0.25 -1 0.25]); 
%axis([0.15 0.45 0.15 0.45]); 
%axis([0.299 0.301 0.299 0.301]); 
e1=0.1;
e2=0.1;
ax1=1.1*[-e1+xs1 e1+xs1 -e1+ys1 e1+ys1]; 
ax2=[-3*e2+xs2 3*e2+xs2 -3*e2+ys2 3*e2+ys2]; 
ax=1.1*[min(ax1(1),ax2(1)) max(ax1(2),ax2(2)) min(ax1(3),ax2(3)) max(ax1(4),ax2(4))];

axis(ax);
%axis(ax1);
%axis(ax2);

hold on
plot(xs2,ys2,'b+','LineWidth',lw)
set(gca,'FontSize',[18]);
xlabel('x_n'); ylabel('y_n');
title(['Henon map: a=',num2str(a),', b=',num2str(b)])

% build circles around f.p.
Ncircles=5;
col=jet(2*Ncircles);
Npts=200;
ang=0:2*pi/Npts:2*pi;
ang=ang';
x1= xs1+e1*cos(ang);
y1= ys1+e1*sin(ang);
x2= xs2+e2*cos(ang);
y2= ys2+e2*sin(ang);
 
plot(x1,y1,'r-','LineWidth',lw)
plot(x2,y2,'b-','LineWidth',lw)


for i=1:Ncircles
 fprintf('Press any key to continue...\n');
 pause
 H1=henon_map(a,b,x1,y1);
 H2=henon_map(a,b,x2,y2);
 x1=H1(:,1);
 y1=H1(:,2);
 x2=H2(:,1);
 y2=H2(:,2);
 if lw>1  
%  lw=lw-1;
  lw=lw+0.25;
 end
 plot(x1,y1,'r-','LineWidth',lw)
 plot(x2,y2,'b-','LineWidth',lw)
 plot(x1,y1,'-','LineWidth',lw,'Color',col(2*i,:))
 plot(x2,y2,'-','LineWidth',lw,'Color',col(2*i-1,:))
end

