% Julia set
close all
clear all

c=0;
%c=-0.17+i*0.78;
%c=0.38+i*0.32;
%c=0.32+i*0.043;
%c=-0.5+0.5*i;
%c=-1;
%c=i;
a=real(c);
b=imag(c);
%The square
xmin=-2;
xmax=2;
ymin=-2;
ymax=2;

Nrpointsx=200;
Nrpointsy=200;

x=linspace(xmin,xmax,Nrpointsx);
y=linspace(xmin,xmax,Nrpointsy);
Nx=length(x);
Ny=length(y);

Nriter=100;
Infi=2;
z=zeros(Nx,Ny);

for i=1:Nx
    for j=1:Ny
        xn=x(i);
        yn=y(j);
        flag=0; % checks if the new values are higher or lower than Infi
        k=1;
        while (k<=Nriter) && (flag==0)
            k=k+1;
            xnew=xn^2-yn^2+a;
            ynew=2*xn*yn+b;
            xn=xnew;
            yn=ynew;
            r=sqrt(xn^2+yn^2);
            flag=(r>Infi);
        end
        %if flag==1
        %    z(i,j)=1;
        %end
        z(i,j)=k;
    end
end

z=z';
%cmap = rot90(hsv(Nriter),2);
cmap = rot90(hot(Nriter),2);
colormap(cmap);
image(x,y,z)
axis tight off

%spy(z)