% HW6: estimating the correlation dimension of the Henon attractor

a = 1.4; b = 0.3;
f = @(x) [a-x(1,:).^2+b*x(2,:);  x(1,:)]; % Henon map (vec func) w/ params a,b

N = 10000;

% your task: add code here which fills x with a 2-by-N array of points on the
% Henon attractor.
% ...

x = zeros(2, N+1);
x(:,1) = [0;0];
for n=1:N
  x(:,n+1) = f(x(:,n));
end
N = N+1;      % note the rest only works if x has size N

% estimate cor dim (takes a few sec, since an O(N^2) process)

r = 0.2;
c = 0;     % counts
tic;
for n=1:N
  c = c + numel(find(sum((kron(x(:,n), ones(1,N)) - x).^2, 1) < r^2));
end
toc
disp(sprintf('C(r) = %g', c/N^2));

% r= 0.01: C(r) = 0.00131879
% r = 0.03: C(r) = 0.00440619
% r = 0.1: C(r) = 0.0187927
% r = 0.2: C(r) = 0.0436107

d = log(0.00440619/0.00131879) / log(3)

d = log(0.0187927/0.00440619) / log(10/3) 
% the one I gave them: gives 1.2047266523

d = log(0.0436107/0.00440619) / log(20/3)
