
% Iterating the Henon map

clear

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

a=1.4;b=0.3;
x0=0.1;
y0=0.2;
Ntransient=1000;
Npts=10000;

%Transient iterations to be droped 
for i=1:Ntransient
 H=henon_map(a,b,x0,y0);
 x0=H(1);
 y0=H(2);
end

%Iterations to keep:
X=zeros(Npts+1,1);
Y=zeros(Npts+1,1);
X(1)=x0;
Y(1)=y0;
for i=1:Npts
 H=henon_map(a,b,X(i),Y(i));
 X(i+1)=H(1);
 Y(i+1)=H(2);
end

figure(1);clf
plot(X,Y,'.','MarkerSize',0.1)
axis tight
title(['Henon attractor'])

%construct bins for histogram
myeps = 0.05; % Size of bins
xmin=min(X); xmax=max(X);
ymin=min(Y); ymax=max(Y);
xbins=xmin:myeps:xmax+myeps;
ybins=ymin:myeps:ymax+myeps;

%Compute histogram
H=histcounts2(X,Y,xbins,ybins);
Z=1*(H>0); % non-empty bins
Ne=sum(sum(Z>0)); % just sum the non-empty bins
[NBx,NBy]=size(Z);NB=NBx*NBy; % number of total bins

figure(2);clf
surf(ybins(1:end-1),xbins(1:end-1),H)
view([90 -90])
axis tight
colorbar
colormap(jet(1024))
title(['Histogram with epsilon = ',num2str(myeps)])

figure(3);clf
surf(ybins(1:end-1),xbins(1:end-1),Z)
view([90 -90])
axis tight
colorbar
colormap(jet(1024))
title(['Non-empty bins: ',num2str(Ne),' out of ',num2str(NB),' bins'])





