Function for plotting Julia Set

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. This is a very nice, well commented and well structured code!

Your comments in the start are exactly how they should be. You already know this, but I'll say it anyway. Writing a function where you have comments like this in the start aids the user because help functionName will provide the information needed in order to use the function.

There's is only one minor thing I would do differently. You know the end points, and how many steps there are. You have two options:

1. Calculate the step-size and use it to create a vector with the correct number of steps using colon
2. Just use what you know directly in linspace

Instead of calculating re_step and do reals = limits(1) : re_step : limits(2); you can use linspace:

reals = linspace(limits(1), limits(2), image_size(1));
imags = linspace(limits(4), limits(3), image_size(2));

The following comment is incorrect: "% Transpose so that plot will have reals on the x axis"

You are working with complex numbers, and in those cases, ' is the complex conjugate transpose, not the regular transpose. I don't know if you want the transpose or the complex conjugate transpose, but it's nice to keep in mind.

• Thanks Stewie, I tried to include your suggestions from last time. I didn't know about the transpose acting differently for complex numbers. Is there a reason you'd prefer linspace over using colons? – Vermillion Nov 1 '16 at 16:53
• In this case, simplicity. You know the end points, and how many steps there are. You have two options: 1. Calculate the step-size and use it to create a vector with the correct number of steps using colon, or 2. Just use what you know directly in linspace. – Stewie Griffin Nov 1 '16 at 19:50
• That's a good point, it's more direct – Vermillion Nov 1 '16 at 22:47
• r.e. the transpose: colour is a matrix of reals anyway (number of iterations) so don't think this note really applies, although it's worth being aware of. Also, you can avoid the colour = colour'; line all together by using z = bsxfun(@plus, reals', imags * 1i); I think. – Wolfie Apr 25 '17 at 13:06