The Julia set is the set of complex numbers z that do not diverge under the following iteration:
$$ z = z^2 + c $$
Where c is a constant complex number.
Different values of z reach infinity at different rates. A colourful fractal is produced by colouring the complex plane based on the number of iterations required to reach infinity.
function colour = julia(c, total_iterations, image_size, limits)
% Calculates julia set by iterating z = z^2 + c, where z and c are complex,
% and recording when z reaches infinity.
%
% Inputs:
%
% c Fixed complex number, of form a + bi
% total_iterations Number of iterations of z = z^2 + c
% image_size 2D vector with number of complex coordinates in
% x and y directions
% limits Vector with 4 elements: min x, max x, min y, max y
%
% Outputs:
%
% colour Matrix of doubles, with size equal to image_size.
% Plotting this matrix will produce a julia set
im_step = (limits(4) - limits(3)) / (image_size(2) - 1);
re_step = (limits(2) - limits(1)) / (image_size(1) - 1);
reals = limits(1) : re_step : limits(2); % Real numbers
imags = limits(3) : im_step : limits(4); % Imaginary numbers
z = bsxfun(@plus, reals(:), (imags(:) * 1i)'); % Complex coordinates
colour = inf(size(z)); % Colour of Julia set
for iteration = 1:total_iterations
index = isinf(z);
% Only perform calculation on the z values that are less than infinity
z(~index) = z(~index).^2 + c;
% Colour depends on number of iterations to reach infinity
colour(index & isinf(colour)) = iteration;
end
colour = colour'; % Transpose so that plot will have reals on the x axis
end
Example:
image_size = [1000 1000]; % Image size
limits = [-1.5 1.5 -1 1]; % Real and imaginary limits
c = -0.4 + 0.6i; % Fixed complex number
total_iterations = 300;
colour = julia(c, total_iterations, image_size, limits);
I = imagesc(colour);
colormap(hot)
myaa('publish')
Note: myaa is a 3rd party anti-aliasing function.
