% Basin of attraction for the period 2 0rbit of the Henon map
% RCG Sep'22

a=1.28;
b=-0.3;
Nt=300;           % # of iterates
Nx=700;Lx=2.5;Ny=Nx;Ly=Lx;
distp=1e-2;       % close to P2
distM=1e+5;       % cap before going to Inf

plottransient=20; % plot during evolution?
mycolmap=[[1 1 1]*0.8;1 0 0;0 0 0];
figure(1);colormap(mycolmap);
if(plottransient) figure(2);clf;colormap(jet(2048));end

%period-1 points:
p1x1=((b-1)+sqrt((1-b)^2+4*a))/2;
p1x2=((b-1)-sqrt((1-b)^2+4*a))/2;
p1y1=p1x1;
p1y2=p1x2;

%period-2 points:
p2x1=((1-b)+sqrt((1-b)^2-4*(-a+(1-b)^2)))/2;
p2x2=((1-b)-sqrt((1-b)^2-4*(-a+(1-b)^2)))/2;
p2y1=(a-p2x1^2)/(1-b);
p2y2=(a-p2x2^2)/(1-b);

%Henon map (stores result in a 2-list)
H = @(X,Y) {a-X.^2+b*Y,X};

%Space of ICs
x=linspace(-Lx,Lx,Nx);
y=linspace(-Ly,Ly,Ny);
[X0,Y0]=meshgrid(x,y);
X=X0;Y=Y0;

fprintf('it=')
for ii=1:Nt  % iterate
 fprintf('%d,',ii)
 XY=H(X,Y);X=XY{1};Y=XY{2};
 DD=X.^2+Y.^2;X(DD>distM)=distM;Y(DD>distM)=-distM; %avoid Inf and NaN
 if(ii<=plottransient)  % plot during evolution
  figure(2)
  D1=(X-p2x1).^2+(Y-p2y1).^2;  % distance to 1st point
  D2=(X-p2x2).^2+(Y-p2y2).^2;  % distance to 2nd point
  D=sqrt(min(D1,D2));
  surf(X0,Y0,log(D)/log(10))
  hold on
  plot([p1x1,p1x2],[p1y1,p1y2],'mx','MarkerSize',10,'LineWidth',2)  
  plot([p2x1,p2x2],[p2y1,p2y2],'b+','MarkerSize',10,'LineWidth',2)
  hold off
  shading interp;view([0 90]);axis image;
  colorbar
  title(['Log_{10} of distance to P2 after ',num2str(ii),' iterates'])
  drawnow
 end
end
fprintf('\n')

%build distance matrix
figure(1)
D1=(X-p2x1).^2+(Y-p2y1).^2;  % distance to 1st point
D2=(X-p2x2).^2+(Y-p2y2).^2;  % distance to 2nd point
D=min(D1,D2);
B=0*D;         % default: stay around
B(D>distM)= 1; % went to Inf
B(D<distp)=-1; % converged to P2
surf(X0,Y0,B)
hold on
plot([p1x1,p1x2],[p1y1,p1y2],'mx','MarkerSize',10,'LineWidth',2)
plot([p2x1,p2x2],[p2y1,p2y2],'b+','MarkerSize',10,'LineWidth',2)
hold off
shading interp;view([0 90]);axis image;
title('Basins of attraction for P2 for Henon Map')
