% NEWTON_FRACTAL  An implemenation of Newton fractal for z^3 - 1
% It demonstrates how to work with vectorization, logical matrices,
% complex numbers and bitmap images construction.
%
% Compatible with MATLAB and GNU Octave 4.0

% (C) 2015 Alexey Voskov
% License: public domain
function newton_fractal
% Prepare initial approximations (mesh grid)
dx = 0.005;
xv = [-3:dx:(-dx), dx:dx:3];
[A,B] = meshgrid(xv, xv);
X = A + B*1i;
X0 = X - 1e-5;
% Newton method iterations
IterNum = zeros(size(X)); % Number of iterations
while max(max(abs(X0 - X))) > eps
    ii = abs(X0 - X) > eps;
    IterNum(ii) = IterNum(ii) + 1;
    Xii = X(ii);
    X0(ii) = Xii;
    X(ii) = Xii - (Xii.^3 - 1)./(3*Xii.^2);    
end
% Find all roots of the cubic equation and find convergence region
% of the complex plane for each region
x0 = roots([1 0 0 -1]);
R = nan(size(X));
for i = 1:numel(x0)
    R(abs(X - x0(i)) < 1e-10) = i;
end
% Show the fractal
picVariant = menu('Select the color scheme', 'Simple', 'Shaded');
if picVariant == 1
    % Variant 1: simple three-colored image
    image(R);
    colormap([0 0 1; 0 1 0; 1 0 0]);
else
    % Variant 2: take into account the number of iterations
    IterNumMin = min(min(IterNum));
    IterNumMax = max(max(IterNum));
    col = 1 + (20 - round(20 * (IterNum - IterNumMin) / (IterNumMax - IterNumMin)));
    image(21*(R-1) +col)
    cmz = zeros(1,21); cmc = (1:21)/21;
    colormap_custom = [cmc cmz cmz; cmz cmc cmz; cmz cmz cmc]';
    colormap(colormap_custom);
end
end
