% ### EXfractal2.m ###     10.08.14
% Creates a (color) fractal based upon the Julia set [http://mathworld.wolfram.com/JuliaSet.html]

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

% Some values to try:
% --- (successively zooming in)
% c= -0.297491-0.641051*i, xB=[-1.35 1.35], yB=[-1.05 1.05], scale=0.005 {BROAD VIEW}
% c= -0.297491-0.641051*i, xB=[-0.7 -0.2], yB=[-0.05 0.4], scale=0.0005
% c= -0.297491-0.641051*i, xB=[-0.6 -0.48], yB=[0.24 0.36], scale=0.00005
% ---
% c= -0.726895347709+ 0.1888871290438*i, xB=[-1.5 1.5], yB=[-1 1], scale=0.005 {BROAD VIEW}
% ---
% c= -0.123 + 0.745*i, xB=[-1.5 1.5], yB=[-1 1], scale=0.005 {Douady's Rabbit Fractal}

clear; figure(1); clf; 
% -------------------------------------
Nmax = 100;      % # of iterations to run to check for divergence {30}
maxZ= 2;        % threshold |z| to determine divergence {2}
c= -0.297491-0.641051*i;   % complex const. value (changing leads to different fractals)
% spatial scale of complex plane (x+i*y) to consider
xB= [-1.5 1.5]; % x-range 
yB= [-1 1]; % y-range 
scale = 0.005;  % spacing between adjacent points
fileN= './fractal1.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));
z = x+ i*y;
zI= z;  % Initialize matricies for iterative loop
k = zeros(size(z));
% +++++
% loop to iterate the grid and check for divergence (flagging the array k)
N = 1;  % initialize index 
while N<Nmax && ~all(k(:))
    zI = zI.^2+c;
    % Any k locations for which abs(zI)>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(zI)>maxZ) = N;    % clever way to update values conditionally
    N = N+1;    % Increment iteration count
end
k(k==0) = Nmax; % For k's still =0, set to N (i.e., these are points in the Julia set)
% +++++
% Visualize
s=pcolor(x,y,k);  % Display the matrix as a surface (color by default)
tb = colorbar('peer',gca);  set(get(tb,'ylabel'),'String', 'k'); % create colorbar and label
if (1==0), colormap(bone);   end    % plot on B&W?? 
set(s,'edgecolor','none')   % Turn off the edges
axis([xB(1) xB(2) yB(1) yB(2)]); % Adjust axis limits, ticks, and tick labels
fontsize=15; xlabel('Re z','FontSize',fontsize); ylabel('Im z','FontSize',fontsize);
% +++++
% save picture to file as a jpg w/ a user=specified resolution
REZ= '-r180';     % resolution for exporting colormaps to jpg
print('-djpeg',REZ,[fileN]);