%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% plot up the Mandelbro set! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
clear all;
close all;

% set parameters and initial conditions
x_min = -3.0;
x_max = 2.0;
y_min = -2.5;
y_max = 2.5;
% set up xy_mag so that if zero, 
% two mouse button clicks are needed and a box is defined
xy_mag = 4.0;

% create some altermative color scales
c_max = 128;
loc_cmap1 = [1+zeros(1,c_max+1); 1:-1/c_max:0; 1:-1/c_max:0]';
loc_cmap2 = [1:-1/c_max:0; zeros(1,c_max+1); 0:1/c_max:1]';
brotmap12 = [ 0 0 0;
    loc_cmap1;
    loc_cmap2];
c_max = 256;
loc_cmap0 = [1:-1/c_max:0; 1:-1/c_max:0; 1+zeros(1,c_max+1)]';
brotmap0 = [ 0 0 0;
    loc_cmap0];
c_max = 128;
loc_cmap3 = [1:-1/c_max:0; 1:-1/c_max:0; 1+zeros(1,c_max+1)]';
loc_cmap4 = [zeros(1,c_max+1); zeros(1,c_max+1); 1:-(1-0.25)/c_max:0.25]';
brotmap34 = [ 0 0 0;
    loc_cmap3;
    loc_cmap4];
brotmapjet2 = [ 0 0 0;
    jet;
    flipud(jet)];

% continually loop:
% find the mouse location (left click), re-calculate, and re-plot
% exit loop is right-hand mouse button is clicked
go = true;
while go
    % call function
    zz = fun_brot(x_min,x_max,y_min,y_max);
    [ylen, xlen] = size(zz);
    % plot
    figure;
    cmap = colormap(brotmapjet2);
    imagesc(zz);
    h=gca;
    set(h,'XTick',[1 xlen],'XTickLabel',{num2str(x_min), num2str(x_max)});
    set(h,'YTick',[1 ylen],'YTickLabel',{num2str(y_max), num2str(y_min)});
    xlabel('Real part');
    ylabel('Imaginary part');
    %
    if (xy_mag == 0.0),
        [xlims,ylims,button] = ginput(2);
    else
        [xlims,ylims,button] = ginput(1);
    end
    %rect = getrect;
    if (button(1)==3),
        go = false;
    else
        %
        if (xy_mag == 0.0),
            xlims = sort(xlims);
            x_min_new = x_min + (x_max-x_min)*xlims(1)/xlen;
            x_max_new = x_min + (x_max-x_min)*xlims(2)/xlen;
            ylims = sort(ylims);
            y_min_new = y_min + (y_max-y_min)*ylims(1)/ylen;
            y_max_new = y_min + (y_max-y_min)*ylims(2)/ylen;
        else
            x_min_new = x_min + (x_max-x_min)*xlims/xlen-(x_max-x_min)/xy_mag/2;
            x_max_new = x_min + (x_max-x_min)*xlims/xlen+(x_max-x_min)/xy_mag/2;
            y_min_new = y_min + (y_max-y_min)*ylims/ylen-(y_max-y_min)/xy_mag/2;
            y_max_new = y_min + (y_max-y_min)*ylims/ylen+(y_max-y_min)/xy_mag/2;            
        end
        %
        x_min = x_min_new;
        x_max = x_max_new;
        y_min = y_min_new;
        y_max = y_max_new;
        %
        close;
    end
end
