% ### EXfractal1.m ###     10.08.14
% Creates a (color) fractal based upon the Mandelbrot set
% [modified from source code of Steve Lord/Birkhauser]

% Uses the difference map (i.e., discrete ODE) 
%   z_{n+1} = (z_n)^2 + c  (z,c are both complex)
% where c = x+i*y (i.e., some point in the complex plane)

% Basic idea is to choose some initial z (e.g., z_0= 0) and 'iterate' the 
% map (i.e., increment n some # of times). For that c, see if the value
% 'diverges' (i.e., blows up towards infinity). If not, we say that that
% point is part of the Mandelbrot set. If so, we can color code the
% divergence rate (so to get a prettier color image)

% NOTE: At best, we only have an approximation of the Mandelbrot set, since
% we can only iterate the map a finite # of times (and a limited # at
% that!). Basically, if |z|>2, the point will likely diverge. In general,
% this can be seen after ~30 iterations.

% Some values to try (successively zooming in):
% - xB=[-2.4 1.2]; yB=[-1.5 1.5]; scale=0.005;  % {BROAD VIEW}
% - xB=[-1.4 -1.05]; yB=[0.15 0.5]; scale=0.0005;
% - xB=[-1.31 -1.22]; yB=[0.3 0.45]; scale=0.0001;

clear; figure(1); clf;
% -------------------------------------
Nmax = 30;      % # of iterations to run to check for divergence {30}
maxZ= 2;        % threshold |z| to determine divergence {2}

% spatial scale of complex plane (x+i*y) to consider
xB= [-2.4 1.2]; % x-range 
yB= [-1.5 1.5]; % y-range 
scale = 0.005;  % spacing between adjacent points

fileN= './mandelbrot1.jpg';    % filename to save image to

% -------------------------------------
% generate a grid of points in the specified complex plane, and determine
% the associated c value (i.e., simply x+i*y)
[x,y]=meshgrid(xB(1):scale:xB(2), yB(1):scale:yB(2));
c = x+ i*y;

% Initialize matricies for iterated value (z) and number of iterations
% before blowing up (k); thus everything will be done at once in a single
% while loop (rather than nested for loops for each x,y pair)
z = zeros(size(c));
k = zeros(size(c));
% +++++
% loop to iterate the grid and check for divergence (flagging the array k);
% loop stops once N reaches Nmax or all points have been determined to have diverged
N = 1;  % initialize index 
while N<Nmax && ~all(k(:))
    z = z.^2 +c; % Square w, add z [can also try variations: z.^3 +z +c, z.^5 +c, ...]
    % Any k locations for which abs(w)>maxZ at this iteration and no
    % previous iteration get assigned the value of N (indicates rate of 
    % divergence and therefore means to color code)
    k(~k & abs(z)>maxZ) = N;    % clever way to update values conditionally
    N = N+1;    % Increment iteration count
end

% If any k's are equal to 0 (i.e. the corresponding w's never blew up) set
% them to the final iteration number 
k(k==0) = Nmax;

% Display the matrix as a surface (color by default)
s=pcolor(x,y,k); 
tb = colorbar('peer',gca);  set(get(tb,'ylabel'),'String', 'k'); % create colorbar and label
% plot in B&W??
if (1==0), colormap(bone);   end
% Turn off the edges of the surface (because the cells are so small, the
% edges would drown out any useful information if we left them black)
set(s,'edgecolor','none')
% Adjust axis limits, ticks, and tick labels
axis([xB(1) xB(2) yB(1) yB(2)])
fs=15; xlabel('Re z','FontSize',fs); ylabel('Im z','FontSize',fs);

% +++++
% save picture to file as a jpg w/ a user=specified resolution
REZ= '-r180';     % resolution for exporting colormaps to jpg
print('-djpeg',REZ,[fileN]);