% iteration of a 2d map acting on the unit disc. barnett 9/18/07

henon = 1;

if henon
  %a = -0.09;
  a=0;
  b = 0.4;
  f = @(x) [a-x(1,:).^2+b*x(2,:);  x(1,:)];        % Henon map w/ params a,b
else
  %A = [.8 .5; 0 1.3];
  A = [.7 2; 0 .7];
  %A = [4.5 8;  -2 -3.5];
  f = @(x) A*x;
end
% notice that f takes each column 2-vector in x to its 2-vector f(x)

N = 6;

[rr, tt] = meshgrid(0:0.1:1, (0:0.01:1)*2*pi); si = size(rr);
%p = [-0.6 -0.6]; eps = 0.5;
p = [0 0]; eps = 0.5;         % choose disc center and radius

disc = [p(1) + eps*rr(:)'.*cos(tt(:)'); p(2) + eps*rr(:)'.*sin(tt(:)')];  % define a list of points

figure;
for n=1:N
  surf(reshape(disc(1,:),si), reshape(disc(2,:),si), n*ones(si));   % show it
  shading interp; axis equal;
  axis([-5 5 -5 5 0 N]); axis([-3 3 -3 3 0 N]);
  view(0,90);
  disc = f(disc);                            % apply f to all points in disc
  caxis([1 N]); grid off; drawnow; hold on; pause(0.5);          % wait
end

% done
