%author: ankur kumar
%Julia - Mandelbrot set fractal
%source: Wikipedia: http://en.wikipedia.org/wiki/Julia_set
%quadratic polynomial: F(z) = z^2 + c

golden = (1+sqrt(5))/2; %golden ratio
clear j %make sure j=sqrt(-1)
%c  values
% c = (golden-2)+j*(golden-1);
% c= 0.45+0.1428j;
%c = j;
c=0.285+0.01*j;

%use backwards iteration z^2(n+1) = z(n) - c
n = 10000; %number of iterations
z = ones(1,n);
%initial z value
% z(1,1) = 1 - golden*j;
z(1,1) = 1;
for i=2:n
    pick = rand(1)>0.5; %randomly pick one of the two square roots
    roots = [-sqrt(z(1,i-1)) sqrt(z(1,i-1))]; %the two square roots to choose from
    temp = roots(pick+1);
    z(1,i) = temp - c;
end
data = [real(z)' imag(z)']; %set up the data for graphing
map = jet(n);
figure(1)
hold on
%plot in parallel coordinates
for i=1:length(data(:,1))
    plot(data(i,:),'color',map(i,:));
end
%the axes
axisx=[1;2];
axisy=[0;0];
for i=1:2
    axisx=[axisx,[i;i]];
    axisy=[axisy,[-5;5]];
end
line(axisx, axisy, 'color', 'b');
set(gca, 'xtick',[1 2])
set(gca, 'xticklabel', {'Re', 'Im'})
hold off
figure(2)
map = jet(length(data(:,1)));
hold on
for i=1:length(data(:,1))
    plot(data(i,1),data(i,2),'color',map(i,:));
end
hold off

