
% Iterating the Henon map
% Stability of a period-2 orbit Ex 2.13

%clear
fs=18;

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

a=input('a [0.43] = ');
b=input('b [0.4] = ');

if(isempty(a))
 a=0.43
end

if(isempty(b))
 b=0.4
end

x21 = (1-b+sqrt((1-b)^2-4*(-a+(1-b)^2)))/2
x22 = (1-b-sqrt((1-b)^2-4*(-a+(1-b)^2)))/2

figure(1);clf
figure(2);clf

figure(1)
hold on
plot(x21,x22,'k+','LineWidth',4,'MarkerSize',14)
plot(x22,x21,'k+','LineWidth',4,'MarkerSize',14)
axis([-0.45 1 -0.45 1]); 
set(gca,'FontSize',[fs]);
xlabel('x_n'); ylabel('y_n');
fprintf('Position windows and press any key to continue...\n');
pause;

Npts=10*0+2;
NICs=25;
col=jet(2*NICs);
fs=15;

for j=1:NICs
 
% ICs=input('ICs [x0 y0] = ');
% x0=ICs(1);
% y0=ICs(2);

 if j>1
  fprintf('Press any key to continue...\n');
  pause;
 end
 x0=x21+0.10*cos(j/NICs*2*pi+0.1);
 y0=x22+0.10*sin(j/NICs*2*pi+0.1);

 X=x0;
 Y=y0;
 xo=x0;
 yo=y0;
 
 index=[0];
 
 for i=1:Npts
  index=[index i];

%  H=henon_map(a,b,xo,yo);
%  xo=H(1);
%  yo=H(2);

  H=henon_map(a,b,xo,yo);
  xn=H(1);
  yn=H(2);
  X=[X xn];
  Y=[Y yn];
  xo=xn;
  yo=yn;
 end
 
 figure(1)
 hold on
 plot(X,Y,'o','LineWidth',1,'MarkerFaceColor',col(2*j,:),'Color',col(2*j,:))
 plot(x21,x22,'k+','LineWidth',4,'MarkerSize',14)
 plot(x22,x21,'k+','LineWidth',4,'MarkerSize',14)
 plot(X(1),Y(1),'o','LineWidth',3,'MarkerFaceColor','w','Color',col(2*j,:))
 %axis tight;
 axis image
 set(gca,'FontSize',[fs]);
 xlabel('x_n'); ylabel('y_n');
 title(['Henon map: a=',num2str(a),', b=',num2str(b),', (x_0,y_0)=(',num2str(x0),',',num2str(y0),')']);
 
% fprintf('Press any key to continue...');
% pause
 
 figure(2)
 hold on
 
 disp=150;
 if disp>Npts-1
  disp=Npts-1;
 end
% subplot(2,1,2)
 plot(index(1:disp),X(1:disp),'-o','LineWidth',1,'MarkerSize',10,'MarkerFaceColor',col(2*j,:))
 axis tight
 V=axis;
 if max([abs(V(3)) abs(V(4))]) > 5
  axis([0 NICs -1 0.0]); 
 end
 set(gca,'FontSize',[fs]);
 xlabel('n');ylabel('x_n');
 title(['Henon map: a=',num2str(a),', b=',num2str(b),', (x_0,y_0)=(',num2str(x0),',',num2str(y0),')']);
% subplot(2,1,2)
% plot(index(1:disp),X(1:disp),'-o','LineWidth',1,'MarkerSize',10,'MarkerFaceColor',col(j,:))
% plot(index(Npts-disp:Npts),X(Npts-disp:Npts),'r-')
% axis tight
% set(gca,'FontSize',[fs]);
% xlabel('n');ylabel('x_n');
 
end
