% interactive explore basin of attraction for 2d iterated map. barnett 9/18/07
% fast vectorized version...

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

nmax = 100;                                          % max # iterations
p = [-2 -2]; win = [4 4];                  % initial view (bot L, wid, hei)
m = 300;                                  % resolution: slow if make large

while (1) % ............                   don't stop zooming
  xs = p(1) + win(1)*(0:m)/m; ys = p(2) + win(2)*(0:m)/m;
  [xx, yy] = meshgrid(xs, ys); si = size(xx);
  x = [xx(:)'; yy(:)'];                     % pack into 2-by-(# points) list
  
  figure;
  for n=1:nmax
    x = f(x);                                    % iterate all points at once
    c = reshape(sum(x.^2, 1)>1e6 | isnan(x(1,:)), si);
    imagesc(xs, ys, c); set(gca, 'ydir', 'normal'); caxis([0 1]);
    axis equal tight; colormap(1-gray(256)); xlabel('x'); ylabel('y');
    title(sprintf('Henon map basin with a=%g, b=%g', a, b)); drawnow;
  end
  disp('done');

  % code for interactive grabbing a rectangle... from 'help rbbox'
  k = waitforbuttonpress;
  pt1 = get(gca,'CurrentPoint');    % button down detected
  finalRect = rbbox;                   % return figure units
  pt2 = get(gca,'CurrentPoint');    % button up detected
  pt1 = pt1(1,1:2);              % extract x and y
  pt2 = pt2(1,1:2);
  p = min(pt1,pt2);             % calculate location of bottom left corner
  win = abs(pt1-pt2);         % and window dimensions (wid, hei)
  x = [p(1) p(1)+win(1) p(1)+win(1) p(1) p(1)];
  y = [p(2) p(2) p(2)+win(2) p(2)+win(2) p(2)];
  hold on
  axis manual
  plot(x,y);                            % redraw in dataspace units
end  % ..............


% timing test - v slow, about 2e4 per sec! also if vectorize with 100
%tic; x = [0; 0]; for n=1:1e5, x = f(x); end; toc       % 4.5 sec
%tic; x = zeros(2,100); for n=1:1e5, x = f(x); end; toc    % 5 sec!

% done
