% interactive explore Mandelbrot set. barnett 10/8/07
% fast vectorized... but using matlab's C-#s.
% 10/28/07: R-button gives iteration plot

nmax = 200;                                          % max # iterations
p = [-2.2 -1.5]; win = [3 3];                  % initial view (bot L, wid, hei)
m = 300;                                  % resolution: slow if make large
disp('Left button: drag rectangle to zoom, Right button: show iterations');

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);
  c = xx + 1i*yy;                     % pack into complex c list
  x = zeros(si);                        % iteration values
  its = zeros(si);                              % number of its array
  a = ones(si);
  
  figure;
  for n=1:nmax           % calculation (iteration) loop
    x(a) = x(a).^2 + c(a);                      % iterate relevant points
    a = ~(abs(x)>4 |isnan(x));
    d = find(its==0 & ~a);
    its(d) = n;
    if mod(n,10)==0
      imagesc(xs, ys, mod(its, 64)); set(gca, 'ydir', 'normal'); caxis([0 63]);
      axis equal tight;colormap([0 0 0;jet(63)]);xlabel('Re z');ylabel('Im z');
      title(sprintf('Mandelbrot iter n=%d', n)); drawnow
    end
  end

  % code for interactive grabbing a rectangle... from 'help rbbox'

  f = gcf;
  [x, y, b] = ginput(1);
  while (b~=1)
    c = x + 1i*y          % make into c
    figure;
    z = 0;
    for n=1:1e3
      z = z^2 + c; plot(z, '.'); hold on;
    end
    drawnow; title(sprintf('iterations for c=%g', c));
    figure(f); [x, y, b] = ginput(1);
  end
  %pt1 = get(gca,'CurrentPoint');    % button down detected
  pt1 = [x y];
  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
