Filled Julia Set Generator

% This script generates numerically a Filled Julia/Mandelbrot set of a
% function f(z)=z^2+c, where 'c' is a complex constant (or c=z in
% Mandelbrot).
% TODO: Parallel version (with parfor).

%----------------- Parameters:

c=complex(-0.8,0.156); % Complex constant that goes on f(z)=z^2+c (for
% Julia only). Try: -0.8,0.156 or -0.74543,0.11301

reaxis=[-1.5  1.5]; % Real axis range.
imaxis=[-1  1]; % Imaginary axis range.

re_res=1500; % Horizontal resolution, over the real axis.
im_res=1000; % Vertical resolution, over the imaginary axis.
iter=100; % Number of iterations to perform.

mandelbrot=false; % true (1) if Mandelbrot, false (0) if Julia.

%----------------- Code:

h=waitbar(0,'Starting...'); % Progress bar

[Re,Im]=meshgrid(linspace(reaxis(1),reaxis(2),re_res), ...
                 linspace(imaxis(2),imaxis(1),im_res));
z=complex(Re,Im); % Creates a complex rectangular plane/grid matrix where
clear Re Im;      % the iterated calculations will be performed

if mandelbrot
    zcons=z;
else
    zcons=c;
end

itcnt=zeros(im_res,re_res); % Iterations count after becoming NaN/Inf

waitbar(0,h,'Calculating...');

try
    iter_fun=sym('z^2'); % Function to be performed (z or c is added later)

    for k=1:iter
        z=double(subs(iter_fun,z))+zcons; % Performs the function
        jset=logical(isnan(z)); % Checks for NaN values (isnan or isinf?)
        itcnt=itcnt+jset; % Counts iterations after becoming NaN
        waitbar(k/(iter+2),h); % Refresh progress bar
    end
catch
    err_msg = lasterror.message;
    if strfind (err_msg,"sym")
        no_sym_flag=1;
    end
end

if no_sym_flag
    for k=1:iter
        z=z.^2+zcons; % Performs the function
        jset=logical(isnan(z)); % Checks for NaN values (isnan or isinf?)
        itcnt=itcnt+jset; % Counts iterations after becoming NaN
        waitbar(k/(iter+2),h); % Refresh progress bar
    end
end

%----------------- Post-processing:

waitbar(k/(iter+1),h,'Post-processing...');

itper=(1-itcnt/iter)/(1-min(min(itcnt))/iter); % Normalized inetation
% counts (range: [0,1]).
zfix=abs(z); % Normalized z data (without NaN and Inf).
zfix(isnan(zfix))=0;
itper(itper==1)=0;
prettyfig=zfix+itper; % Combination of zfix and itper, displaying purpose.
prettyfig(prettyfig>1)=1;

waitbar(1,h,'Done...');

%----------------- Displaying the images:
imtool(jset);
imtool(prettyfig);

%----------------- Saving the images:
numb=dir('prettyfig*.png'); % Detecting largest file number
if size(numb,1) ~= 0
    numb=struct2cell(numb); numb=numb(1,:);
    y=regexp(numb,'[0-9]{3}');
    numb=cell2mat(numb');
    numb=num2str(nanmax(str2num(numb(:,y{1}:y{1}+2))+1),'%03u');
else
    numb='000';
end
imwrite(uint8(-1+prettyfig*64),pink,['prettyfig' numb '.png']) % Strangely this
imwrite(uint8(-1+jset*64),gray,['juliaset' numb '.png']) % opens an empty figure

%----------------- Clearing used variables:
% clear c im_res imaxis itcnt iter itper jset k prettyfig re_res
% clear reaxis z zfix

close(h);