m=mandelbrot(width, iterations)

function m=mandelbrot(width, iterations)
% MANDELBROT Generate Mandelbrot set with WIDTH pixels on the X axis.
%
% Algorithms based on the Wikipedia article:
% http://en.wikipedia.org/wiki/Mandelbrot_set

    % Height is the even number >= two-thirds of the width
    height = floor(width * 2.5 / 3);
    if mod(height, 2) ~= 0, height = height + 1; end
    
    m = zeros(height, width, 3);
    
    % Display a progress bar
    msg = sprintf('Computing %d x %d Mandelbrot set...', width, height);
    h = waitbar(0, msg);
    col = 1;
    
    % Compute the values of the Mandelbrot set using the popular
    % "Escape Time" algorithm. Each color in the image represents the
    % number of steps required for the iterative series to diverge to
    % infinity. The color black means the series never diverged. The black
    % pixels represent the points in the Mandelbrot set. See the Wikipedia
    % article for much more detail:
    %
    % http://en.wikipedia.org/wiki/Mandelbrot_set
    %
    % The Mandelbrot set has real (x) coordinate from -2 to +1, and
    % imaginary coordinate from -1.5 to +1.5. Since the set is perfectly
    % symmetrical about the X axis, we compute half the data (the
    % negative imaginary plane, Y from -1.5 to 0) and then mirror it to
    % produce the final image.
    for xp = -2.2:3.2/width:1
        row = 1;
        % Compute escape time for the negative imaginary values.
        for yp = -1.25:2.5/height:0 % Mirror around X-axis
            
            k = 0;
            x = 0;
            y = 0;
            
            % Optimization: check to see if the point can be
            % deterministically shown to be in the cardioid or the bulb.
            p = sqrt((xp-.25)^2 + yp^2);
            if ((xp < (p - (2*p^2) + 0.25)) || ...
                (xp+1)^2 + yp^2 < (1/16))
                k = iterations;
            end
            
            % Iterate: determine if the series converges to C (here 2) or
            % diverges to infinity.
            while ( k < iterations && x*x + y*y <= (2*2))
                [k, x, y] = iterate(k, x, y, xp, yp);
            end
            % If the series diverged, color the pixel according to the
            % number of steps required for divergence. Otherwise, it
            % converged, so color the pixel black.
            if k < iterations
                m(row, col, :) = iterationColor(k, x, y, xp, yp, iterations);
            end
            row = row + 1;
        end
        col = col + 1;
        waitbar(col / width);
    end
    % Close the colorbar
    close(h);
    % Mirror the negative imaginary values into the positive half-plane.
    m(row-1:height, :, :) = m(row-2:-1:1, :, :);
    % Convert color data from HSV to RGB. I used HSV because it's easier to
    % navigate, but MATLAB's image command expects RGB colors.
    m = hsv2rgb(m);
    m = uint8(floor(m * 255));
end

function [k, x, y] = iterate(k, x, y, xp, yp)
% ITERATE Solve one iterative step of the equation (Z(n) = Z(n-1)^2 + C).
    xt = x*x - y*y + xp;
    y = 2*x*y + yp;
    x = xt;
    k = k + 1;    
end

function hsv = iterationColor(k, x, y, xp, yp, iterations)
% ITERATIONCOLOR Smoothly interpolate pixel color based on escape count.

    % Iterate three more times to reduce error term.
    for l = 1:3
        [k, x, y] = iterate(k, x, y, xp, yp);
    end
    % Normalize count by taking the log log of norm of Z
    c = normalizeCount(k, x, y);
    % Interpolate the color around the top of the HSV cone. Start at clr
    % degrees. Wrap around using MOD. Set S (saturation) and V (value) to 
    % 1 for bright colors. CLR = 2/3 to start at blue, 1/2 to start at
    % red.
    clr = (1/2);
    hsv = [mod((c / iterations) + clr, 1) 1 1];
end

function nC = normalizeCount(k, x, y)
% NORMALIZECOUNT Fractional iteration count for smooth color interpolation.

    % Take the log log of the norm of Z: log(log(|Z(n)|))
    nC = k + 1 - log(log(sqrt(x*x+y*y))) / log(2);
end