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).^3+i; 
    % 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