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

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:

• Calculate the step-size and use it to create a vector with the correct number of steps using colon

• 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.

StackExchange Code Review Q#145752, answer score: 2